x=read_complex_binary('radio192kHz.bin');
dcf=real(x);
gps=imag(x);
fe=192e3;
time=[0:length(x)-1]'/fe;
dcf=dcf.*exp(j*2*pi*(77500+1.08)*time);
lpf=firls(250,[0 720 790 fe/2]*2/fe,[1 1 0 0]);
dcf=filter(lpf,1,dcf);
x=dcf(1:59:end);
time=time(1:59:end);

fe=192000/59;
N=50;

xa=conv(abs(x),ones(N,1)/N)(N:end-N);
subplot(311);plot([0:length(xa)-1]/fe,xa);
xlabel('time (s)');ylabel('amplitude (a.u.)')

xp=angle(x); 
subplot(312);plot(time,xp)
xlabel('time (s)');ylabel('phase (rad)')

load lfsr.dat
np=192000/59*(120/77500)
oldP=0;
for k=1:length(lfsr)
  P=round(k*np);
  if (lfsr(k)==1) longlfsr(oldP+1:P)=ones(P-oldP,1);
  else longlfsr(oldP+1:P)=zeros(P-oldP,1);
  endif
  oldP=P;
end
yc=xcorr(longlfsr-mean(longlfsr),xp-mean(xp));
yc=yc(1:floor(length(yc)/2));
subplot(313)
plot([0:length(yc)-1]/fe,abs(yc))
xlabel('time (s)');ylabel('xcorr(phase,code)')

hold on
plot([0:length(gps)-1]/192000,(gps-min(gps))*500,'r');

figure
plot([0:length(yc)-1]/fe,abs(yc))
xlabel('time (s)');ylabel('xcorr(phase,code)')
