function [dm, iap]=msq( sam, fc, cyc_dig, ncs, Wn, sw_iap, vsearch, law, hrm)
% function [dm, iap]=msq( sam, fc, cyc_dig, ncs, Wn, sw_iap, vsearch, law, hrm)
%
% fc=75;                           % carrier frequency
% cyc_dig=2;                       % cycles per digit ( or Q)
% ncs=10;                          % number of periods (cirsums) to process
% Wn=40;                           % demod filter cutoff
% sw_iap=1;                        % iap align circums
% vsearch=-5:0.2:5;                % doppler search (m/s)
% law = [ 1 0 0 0 0 0 1 1 0 1 1];  % 2033 oct
% hrm =3;                          % For 3rd harmonic AST 84 Hz hrm=3,
% hrm =1;                          % all others hrm=1 

sam=sam(:);

if nargin==1
  fc=75; cyc_dig=2;
  ncs=10; Wn=40; sw_iap=1;
  vsearch=-5:0.2:5;
  law = [ 1 0 0 0 0 0 1 1 0 1 1];
  hrm=1;
elseif nargin==2
  cyc_dig=2;
  ncs=10; Wn=40; sw_iap=1;
  vsearch=-5:0.2:5;
  law = [ 1 0 0 0 0 0 1 1 0 1 1];
  hrm=1;
elseif nargin==3
  ncs=10; Wn=40; sw_iap=1;
  vsearch=-5:0.2:5;
  law = [ 1 0 0 0 0 0 1 1 0 1 1];
  hrm=1;
elseif nargin==4
  Wn=40; sw_iap=1;
  vsearch=-5:0.2:5;
  law = [ 1 0 0 0 0 0 1 1 0 1 1];
  hrm=1;
elseif nargin==5
  sw_iap=1;
  vsearch=-5:0.2:5;
  law = [ 1 0 0 0 0 0 1 1 0 1 1];
  hrm=1;
elseif nargin==6
  vsearch=-5:0.2:5;
  law = [ 1 0 0 0 0 0 1 1 0 1 1];
  hrm=1;
elseif nargin==7
  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
  hrm=1;
elseif nargin==8
  hrm=1;
end

nvs=length(vsearch);
nr=length(law)-1;             % length of sequence register
L = 2.^nr - 1;                % digits per period
sam_cyc=4;                    % samples per cycle
sam_dig=sam_cyc*cyc_dig;      % samples per digit
fs=sam_cyc*fc;                % sampling frequency
N=sam_cyc*cyc_dig*L;          % samples per period
T=N/fs;                       % period length
dt=1/fs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[firb, fira]=butter(9, 2*Wn/fs);                 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The original filtering
%Wn=40;
%uwarp = 2*fs*tan(pi*Wn/fs);
%[firb, fira]=besself(9, uwarp); 
%[firb, fira]=impinvar( firb, fira, fs);
%clear uwarp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% more filter designs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Wnx=2*Wn/fs;
%uwarp = 4*tan(pi*Wnx/2);
%Wnx = uwarp;
%%[zero,pole,gain]   = buttap(9);
%%[zero,pole,gain]   = besselap(9);
%%[zero,pole,gain]   = besselap(21);
%[zero,pole,gain]   = besselap(9);
%[ssa,ssb,ssc,ssd]  = zp2ss(zero,pole,gain);
%[ssa,ssb,ssc,ssd]  = lp2lp(ssa,ssb,ssc,ssd,Wnx);
%[ssa,ssb,ssc,ssd]  = bilinear(ssa,ssb,ssc,ssd,2,2*Wn/fs);
%fira = poly(ssa);
%firb = poly(ssa-ssb*ssc)+(ssd-1)*fira;
%clear Wnx uwarp zero pole gain ssa 
%clear ssb ssc ssd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

uwarp = 2*fs*tan(pi*Wn/fs);
[zero,pole,gain]   = besselap(9);
[ssa,ssb,ssc,ssd]  = zp2ss(zero,pole,gain);
[ssa,ssb,ssc,ssd]  = lp2lp(ssa,ssb,ssc,ssd,uwarp);
[ssa,ssb,ssc,ssd]  = bilinear(ssa,ssb,ssc,ssd,fs);
fira = poly(ssa);
firb = poly(ssa-ssb*ssc)+(ssd-1)*fira;
clear Wnx uwarp zero pole gain ssa 
clear ssb ssc ssd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% estimate filter response
%%[hrsp, frind]=freqz( firb, fira, [0:2*Wn], fs);
%%plot( frind, 20*log10(abs(hrsp)));
%%grid
%%pause
%%[grpd, frind]=grpdelay( firb, fira, [0:2*Wn], fs);
%%plot( frind, grpd/fs);
%%grid
%%pause
%
%%taxis=0:dt:1;
%%offset=20;
%%sigin=zeros(size(taxis));
%%sigin(offset+[1:4*cyc_dig])=1;
%%sigout=filter(firb, fira, sigin);
%
%%[dmpk, pki]=max(abs(sigout));
%%y1=abs(sigout(pki-1));
%%y2=abs(sigout(pki));
%%y3=abs(sigout(pki+1));
%%tte=taxis(pki) + dt*(y1-y3)/(y1-2*y2+y3)/2;
%
%%clf
%%plot( taxis, sigin)
%%hold
%%plot( taxis, sigout)
%%plot( tte, sigout(pki), 'ro');
%
%%dly=tte-dt*offset; %delay estimate for use in travel time measurements
%
%%clear grpd frind
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sremset( law, hrm); 
cspeed=1484;
dm=zeros(nvs,N); iap=zeros(1,ncs);
for vi=1:nvs
  vopen=vsearch(vi);
  alpha=1-vopen/cspeed;

  fd=alpha*fc;
  n=1:ncs*N; n=n-1;
  arg=rem(fd*n/fs, 1);
  ds=exp(-j*2*pi*arg).*sam.';
  %ds=filter( firb, fira, ds);
  nh=n/alpha;
  ind=find(nh>=max(n)); nh(ind)=[];
  frac=nh-floor( nh);
  dsi=(1-frac).*ds( 1+floor( nh)) + frac.*ds( 2+floor( nh));
  %dsi=ds;
  ndsi=length(dsi);
  if ndsi< N*ncs
    dsi( ndsi+1:N*ncs)=zeros(1, N*ncs-ndsi);
  end

  %dsi=filter( firb, fira, dsi);

  circ=reshape(dsi, N, ncs);
  cor=zeros(size(circ'));
  for ip=1:ncs
    cirsum=circ(:,ip);
    cirsum=reshape( cirsum, sam_dig, L);
    corr=zeros(size(cirsum));
    for ci=1:sam_dig
      corr(ci,:)=srem( cirsum(ci,:)).'; 
    end
    cor(ip,:)=reshape( corr, 1, N);
  end

  cor=cor.';
  cor=cor(:);
  cor=filter( firb, fira, cor);
  %cor=filtfilt( firb, fira, cor);
  %cor=cor/fs;
  %cor=cor.';
  cor=reshape( cor, N, ncs);
  cor=cor.';

  %cor=cor.';
  %cor=[cor(:,1); cor(:); cor(:,ncs)];
  %cor=filter( firb, fira, cor);
  %%    cor=filtfilt( firb, fira, cor);
  %cor=cor(N+1:(ncs+1)*N)/fs;
  %cor=cor.';
  %cor=reshape( cor, N, ncs);
  %cor=cor.';

  if (ncs>1)
    if (sw_iap==1) 
      iapt=0.0; iap(1)=iapt;
      dm(vi,:)=cor(1,:);
      for ip=2:ncs
        iapt=iapt+angle(sum(cor(ip,:).*conj(cor(ip-1,:))));
        dm(vi,:)=dm(vi,:)+cor(ip,:)*exp( -j*iapt);
        iap(ip)=iapt;
      end
    else
      dm(vi,:)=sum(cor);
    end
  else
    dm(vi,:)=cor;
  end

end

return;



