acq3.m

GPS中频信号频率和ca码搜索程序

%**** initial condition*****
svnum=input('enter satellite number=');
%function [pt_init,frr]=acq3(svnum);
tic
intodat=16368*0;
fs=16.368e6;  %sampling freq    down frequency 20.4916MHz
ts=1/fs;   %samplint time
n=fs/1000;  %data pt in 1 ms
nn=[0:n-1];  %total numbers of pts
fc=4.124e6;   %center freq without doppler
nsat=length(svnum); %total number of satellites to be precessed

%*****input data file****
%fid=fopen('d:\gps\mendraw.dat','r');
fid=fopen('d:\gps\mendraw1558.dat','r');
fseek(fid,0,'bof');
x2=fread(fid,6*n,'schar');

yy=zeros(21,n);

%********start acquisition******
code=digitizg(n,fs,0,svnum); %digitize c/a code
xf=fft(x2(1:n)');
for i=[1:21];   %find coarse freq 1 khz resolution
    fr=fc-10000+(i-1)*1000;
    lc=code.*exp(j*2*pi*fr*ts*nn);
    lcf=fft(lc);
    yy(i,:)=ifft(xf.*conj(lcf));
end
[amp crw]=max(max(abs(yy')));  %find highest peak
[amp ccn]=max(max(abs(yy)));
pt_init=ccn; %initial point
cfrq=fc+1000*(crw-11);  %coarse freq
                                                %******************test*********
%code2=codeshift(code,pt_init);
%wave2=cos(2*pi*cfrq*ts*nn);
%x3=x2(1:n)';
%x4=code2.*x3;
%x5=x4.*wave2;
%figure(1);
%plot(x4(1:100))                                        %continue wave
%figure(2);
%plot(x5(1:100))
%sum_out=sum(x5);
%*****generate 5 ms data by stripping c/a code***
z5=x2(pt_init:pt_init+5*n-1);   %take 5 ms data starting with c/a code
za5=z5'.*[code code code code code]; % create cw from 5 sets of data

%******find medium freq resolution 400 hz aprat***
for i=[1:3];
    fr=cfrq-100+(i-1)*400;
    mfrq0(i)=sum(za5(1:n).*exp(j*2*pi*fr*ts*nn));
    mfrq1(i)=abs(mfrq0(i));
end
[mamp mrw]=max(mfrq1); % find highest peak
mfrq=cfrq+200*(mrw-2);
fr=mfrq;   %medium freq

%*****find fine freq******
zb5=za5.*exp(j*2*pi*fr*ts*[0:5*n-1]);  %one DFT component
zc5=diff(-angle(sum(reshape(zb5,n,5)))); %find difference angle
zc5fix=zc5;

%   adjust phase and take out possible phase shift
threshold=2.3*pi/5;
for i=1:4;
    if abs(zc5fix(i))>threshold;
        zc5(i)=zc5fix(i)-2*pi;
        if abs(zc5(i))>threshold;
            zc5(i)=zc5fix(i)+2*pi;   %end
            if abs(zc5(i))>2.2*pi/5;
                zc5(i)=zc5fix(i)-pi;
                if abs(zc5(i))>threshold;
                    zc5(i)=zc5fix(i)-3*pi;
                    if abs(zc5(i))>threshold;
                        zc5(i)=zc5fix(i)+pi; %end
                    end
                end
            end
        end
    end
end
dfrq=mean(zc5)*1000/(2*pi);
frr=fr+dfrq; %fine freq
toc     %compute the time
figure(3);
plot(abs(yy(crw,1:n)))
xlabel('GPS= svnum  max at pt_init')
title(num2str(pt_init))
figure(4);
plot(abs(yy(:,ccn)),'*')
%plot(abs(yy(:,ccn)))
xlabel('GPS =  num2str(svnum)  Freq = num2str(frr)')
title(num2str(frr))
format
pt_init
format long e 
frr
cfrq
mfrq
%figure;
%mesh(abs(yy));