myaciquition.m

利用频域并行的方法对GPS的C码捕获的方法进行了仿真。直接解压就可以了。

clear
%**** initial condition ********
%svnum = input('enter satellite number = ');
svnum  =5;
fs = 5e6; %sample rate
ts = 1/fs;
n= fs/1000;% signal sample numer 
nn= [0:n-1];
t=nn*ts;
fc=1.25e6;% carrier frequency 
df=8000;% Doppler frequency
nsat = length(svnum);

%********* simulated signal ***********
%%********* produce simulated signal *******
code2 = digitizg(5000,5000000,0,svnum);
xx = digitizg(30000,5000000,0,svnum);
%code1 = digitizg(5000,5000000,0,9)
data= code2.*cos(2*pi*(fc+df)*t);% input signal
%x2=x2.*code1;
%plot(x2);
%yy = zeros(21,n);
%axis([0 100 -1.4 1.5]);
%*********** start acquisition *****
code = digitizg(n,fs,0,5);  %generate local  C/A code
xf = fft(data(1:n));%fft of input signal
for i=[1:21];  %** find coarse freq 1KHz resolution
    fr = fc-10000+(i-1)*1000;
    lc = code.*exp(j*2*pi*fr*ts*nn);%generate local code
    lcf = fft(lc);
    yy(i,:) = ifft(xf.*conj(lcf));
end
[amp crw]=max(max(abs(yy')));
[amp ccn]=max(max(abs(yy)));
figure(2)
mesh(abs(yy));

pt_init = ccn;
cfrq = fc+1000*(crw-11);%coarse freq
%*** generate 5ms data by stripping C/A code **********

z5= xx(pt_init:pt_init+5*n-1);% take 5ms 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 KHz apart *****
for i = [1:3];
    fr=cfrq-400+(i-1)*400;
    mfrq0(i)=sum(za5(1:5000) .* 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
temp=reshanpe(zb5,n,5);
temp1=sum(temp);
%zc5=diff(-angle(sum(reshape(zb5,n,5)))); % find difference angle
temp2=-angle(temp1);
zc5=diff(temp2);
zc5fix=zc5;
% ***** Adjust phase and take out possible phase shift *****
threshold=2.3*pi/5;
for i=1:4;
    if abs(zc5(i))>threshold;% for angle adjustment
        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; % for pi phase shift correction
                        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
figure(1)
plot(abs(yy(crw,1:n)))
figure(3)

plot(abs(yy(:,ccn)))
title(['Gps = ' num2str(svnum)  'max at'  num2str(pt_init)]);grid;