%testmsq.m

%ATOC parameters
fc=75;                           % carrier frequency
cyc_dig=2;                       % cycles per digit ( or Q)
Pn=40;                           % periods per transmission
nper=1;                          % periods per circulating sum
nper=8;                          % periods per circulating sum (AVATOC popup)
ncs=Pn/nper;                     % number of circulating sums
law = [ 1 0 0 0 0 0 1 1 0 1 1];  % 2033 oct

Wn=40;                           % dmod filter bandwidth (Hz)

harm=1;                          % 75 Hz
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;
taxis=0:dt:T-dt;
clear nr sam_cyc 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

snrdb=-10;                        % input signal-to-noise ratio
%snrdb=5;                         % input signal-to-noise ratio
%path(path,'/net/wgasa/mad/msq/msq');
sg( 0.0, 0.0, fc, cyc_dig, Pn, snrdb, law, harm);

sname=sprintf( 'proto.sam');
fid=fopen( sname, 'r', 'ieee-le');
line=fgetl( fid);
while strcmp( line, '*DATA*') == 0
  line=fgetl( fid);
end

start=ftell(fid); fseek(fid,0,'eof'); stop=ftell(fid);
fseek(fid,start,'bof'); sizs=(stop-start)/2;
data=fread( fid, sizs, 'short'); data=data(:);
fclose( fid);
clear ans fid start stop sizs line

rmsin=std(data);
%fprintf(1, 'rms input power: %6.2f (bits)\n', rmsin);

circ=zeros(ncs*N,1);
for ip=1:ncs
  skip=(ip-1)*N;
  blah=data(1+nper*skip:nper*N+nper*skip);
  blah=reshape(blah,N,nper);
  if ( nper>1)
    blah=sum(blah.');
  else
    blah=blah.';
  end
  circ(1+skip:N+skip)=blah.';
end
clear ip skip blah

sw_iap=0;
vsearch=0.0;
%path( path, '/net/wgasa/mad/matlab/msq');

dmod=zeros(1,N); iap=zeros(ncs,1);
eltim=cputime;
[dmod iap]=msq( circ, fc, cyc_dig, ncs, Wn, sw_iap, vsearch, law, harm);
eltmmsq=cputime-eltim;

%complete matched filter
%mtfla=1; mtflb=ones(sam_dig,1);
%dmod=filter(mtflb, mtfla, dmod);

[dmpk, pki]=max(abs(dmod));
dmmd=median(abs(dmod));
dmang=180*angle(dmod(pki))/pi;

%For complex demodulates
%x=normrnd(0,5,100000,1);
%y=normrnd(0,5,100000,1);
%z=x+i*y;
%median(abs(z))/sqrt(log(2))
%std(z)
%sqrt(log(2))

tmax=taxis(pki);

md2rms=sqrt(log(2));
dmrms=dmmd/md2rms;
snr=20*log10(dmpk/dmrms);

y1=abs(dmod(pki-1));
y2=abs(dmod(pki));
y3=abs(dmod(pki+1));

tte=taxis(pki) + dt*(y1-y3)/(y1-2*y2+y3)/2;
clear y1 y2 y3

fprintf(1, 'input snr:  %6.2f (dB re 1 quanta)\n', snrdb);
fprintf(1, 'rms input power: %6.2f (quanta) %6.1f (bits) %6.2f (dB)\n', ...
            rmsin, log10(rmsin)/log10(2), 20*log10(rmsin));

exsg=20*log10(ncs*nper*L/sqrt(2));
exng=10*log10(ncs*nper*L/2)+10*log10(2*Wn/fs);
exsnrg=exsg-exng;
fprintf(1,'expected signal gain: %5.1f\n', exsg);
fprintf(1,'expected noise gain: %5.1f\n', exng);
fprintf(1,'expected snr gain: %5.1f\n', exsnrg);
fprintf(1,'tte: %8.5f tmax: %6.3f ang: %7.2f\n', tte, tmax, dmang);
fprintf(1,'dmpk: %5.1f (dB) dmrms: %5.1f (dB) snr: %5.1f (dB)\n',...
        20*log10(dmpk), 20*log10(dmrms), snr);
fprintf(1,'eltmmsq: %8.5f\n', eltmmsq);

clf
subplot( 211)
plot( taxis, 20*log10(abs(dmod)), 'k');
hold on
plot( taxis(pki), 20*log10(abs(dmod(pki))), 'ro');
hold off
grid on; zoom on
set( gca, 'xlim', [0 T])
ylabel( 'Amplitude (dB)');

subplot( 212)
plot( taxis, 180*angle(dmod)/pi, 'k');
hold on
plot( taxis(pki), 180*angle(dmod(pki))/pi, 'ro');
hold off
grid on; zoom on
set( gca, 'xlim', [0 T])
xlabel( 'Delay time (s)');
ylabel( 'Phase (degrees)');

pause
subplot( 211)
set( gca, 'xlim', [0 1])
subplot( 212)
set( gca, 'xlim', [0 1])
drawnow

orient landscape
%print -dpsc proto.ps

return



