function Rm=srem( circ)
% Rm=srem( circ)
%
% Hadamard process linear maximal sequences
%

global sr_L sr_u sr_v sr_Qp sr_Pp sr_U
L=sr_L; u=sr_u; v=sr_v; Qp=sr_Qp; Pp=sr_Pp; U=sr_U;

circ=circ(:);
 
Sm=sum(circ)/L;
A=Sm/(u-v/L);

sp=[];
sp(Qp)=circ;
saug=[ 0; sp.'];

[l,m] = size(saug);
if l <= 1, hqr = saug;
else 
   n = 1;
   while n < l,
      ln = l/n;
      indx = zeros(n,ln); indx(:) = 1:l;
      indx1 = indx(:,1:2:ln); indx1 = indx1(:);
      indx2 = indx1 + n;
      x1 = saug(indx1,:); x2 = saug(indx2,:);
      saug(indx1,:) = x1 + x2;
      saug(indx2,:) = x1 - x2;
      n = 2*n;
   end
   hqr = saug;
end
hqr=hqr(Pp+1);

%%bias correction
%Rb=(hqr/A + u)/v;
%Rm=u*conj(U(1)) - u*conj(U(2))/L - v*conj(U(1))/L + v*conj(U(2))*Rb/L;
%Rm=A*Rm;
%bias correction (same as above, avoid divide by zero problems);
Rb=(hqr + A*u)/v;
Rm=A*(u*conj(U(1)) - u*conj(U(2))/L - v*conj(U(1))/L) + v*conj(U(2))*Rb/L;

return



