% lms.m
%
% demonstrate linear maximal sequences
%

% states associated with generator polynomial
law = [ 1 0 1 1];                % 13 oct
law = [ 1 0 0 1 1 ];             % 23 oct
law = [ 1 0 0 1 0 1];            % 45 oct
law = [ 1 0 0 0 0 1 1];          % 103 oct
law = [ 1 0 0 0 0 0 1 1];        % 203 oct
law = [ 1 0 0 0 1 1 1 0 1];      % 435 oct
law = [ 1 0 0 0 0 1 0 0 0 1];    % 1021 oct
law = [ 1 0 0 0 0 0 1 1 0 1 1];  % 2033 oct
%law = [ 1 0 0 0 0 1 0 0 0 1];    % 1021 oct
law = [ 1 0 0 0 0 1 1];          % 103 oct
law=law(:);

nr=length(law)-1;             % length of sequence register
L = 2.^nr - 1;                % digits per period
thetamod=-atan(sqrt(L));      % modulation angle

% ssrg (simple shift register generator)
G=[ law(2:nr+1)'; [eye(nr-1) zeros(nr-1,1)] ];
Q=zeros(nr,L); Q(1,1)=1;
for il=1:L-1
   Q(:,il+1) = rem(G*Q(:,il),2);
end
Q=[Q(:,L-nr+2:L), Q(:,1:L-nr+1)]; %delay to match KM
Qp=2.^[0:nr-1]*Q;

% msrg (modular shift register generator)
G=[ [zeros(1,nr-1); eye(nr-1)] law(1:nr) ];
P=zeros(nr,L); P(1,1) = 1;
for il=1:L-1
   P(:,il+1)=rem(G*P(:,il),2);
end
Pp=2.^[0:nr-1]*P;

% modulation

j=sqrt(-1);
dg=(1-2*Q(1,:));
u=1/tan(thetamod); v=j;
dm=u+v*dg;

clf 
plot( angle( dm));
hold on
plot( angle( dm), 'x');
set( gca, 'xlim', [1 L]);
set( gca, 'fontname', 'times')
title('Baseband M-sequence')
hold off
pause

% FFT processing

R=[ 0; L];
M=[ u-v/L -u/L-v/L; u-v/L -u/L+v];
M=conj(M); U=M\R;

ref = U(1) + U(2)*dg;
Ref = conj(fft(ref));

t0=cputime;
Rm=ifft(fft(dm).*Ref)/L;
te=cputime-t0;

[Rmax imax]=max(abs(Rm));
subplot(211)
plot( 20*log10(abs(Rm)))
hold on
plot( imax, 20*log10(abs(Rm(imax))), 'go')
hold off
set( gca, 'xlim', [1 L]);
set( gca, 'fontname', 'times')
title('FFT Processing')
subplot(212)
plot( angle(Rm)/2/pi)
hold on
plot( imax, angle(Rm(imax))/2/pi, 'go')
hold off
set( gca, 'xlim', [1 L]);
set( gca, 'fontname', 'times')
fprintf(1,'Rmax: %.2f Angle: %.2f MaxIndex: %d etime: %f\n',...
        Rmax, angle(Rm(imax)), imax, te);
pause

% Hadamard processing

R=[ 0; L];
M=[ u-v/L -u/L-v/L; u-v/L -u/L+v];
M=conj(M); U=M\R;

H = hadamard(L+1);

t0=cputime;
Sm=sum(dm)/L;
A=Sm/(u-v/L);
  
sp(Qp)=dm;
saug=[ 0; sp.'];
hqr=H*saug;
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;
te=cputime-t0;

[Rmax imax]=max(abs(Rm));
subplot(211)
plot( 20*log10(abs(Rm)))
hold on
plot( imax, 20*log10(abs(Rm(imax))), 'go')
hold off
set( gca, 'xlim', [1 L]);
set( gca, 'fontname', 'times')
title('Hadamard Processing')
subplot(212)
plot( angle(Rm)/2/pi)
hold on
plot( imax, angle(Rm(imax))/2/pi, 'go')
hold off
set( gca, 'xlim', [1 L]);
set( gca, 'fontname', 'times')
fprintf(1,'Rmax: %.2f Angle: %.2f MaxIndex: %d etime: %f\n',...
        Rmax, angle(Rm(imax)), imax, te);
pause

% integer Hadamard ( faster)

R=[ 0; L];
M=[ u-v/L -u/L-v/L; u-v/L -u/L+v];
M=conj(M); U=M\R;

t0=cputime;
Sm=sum(dm)/L;
A=Sm/(u-v/L);

sp(Qp)=dm;
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;
te=cputime-t0;

[Rmax imax]=max(abs(Rm));
subplot(211)
plot( 20*log10(abs(Rm)))
hold on
plot( imax, 20*log10(abs(Rm(imax))), 'go')
hold off
set( gca, 'xlim', [1 L]);
set( gca, 'fontname', 'times')
title('Integer Hadamard Processing')
subplot(212)
plot( angle(Rm)/2/pi)
hold on
plot( imax, angle(Rm(imax))/2/pi, 'go')
hold off
set( gca, 'xlim', [1 L]);
set( gca, 'fontname', 'times')
fprintf(1,'Rmax: %.2f Angle: %.2f MaxIndex: %d etime: %f\n',...
        Rmax, angle(Rm(imax)), imax, te);









