acquisition.m

MATLAB编写的GPS 信号的产生捕获程序

clear
clc
svnum=input('enter satellite numner=');
intodat=10001;   %input('enter initial pt into data(multiple of n)=');
fs=12e6;   %Sampling freq
ts=1/fs;  %sampling time
n=fs/1000;   %data pu in 1ms
nn=[0:n-1];  % total no. of pts
fc=3.563e6;   %center freq without doppler
nsat=length(svnum);

[fid,message]=fopen('d:/compactdata_20050407_142600.bin','r','b');
%fseek(fid,6*n,'bof');
x2=fread(fid,n*6,'bit1');
yy=zeros(21,n);

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*n); % generate local code
    lcf=fft(lc);
    yy(i,:)=ifft(xf .* conj(lcf)); % circular correlation
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         

% ***** gerenate 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 KHz apart *****
for i=[1:3];
    fr=cfrq-400+(i-1)*400;
    guance=exp(j*2*pi*fr*ts*nn);
    mfrq0(i)=sum(za5(1:12000) .* 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(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


plot(abs(yy(crw,1:n)))
title(['GPS=' num2str(svnum) 'max at' num2str(pt_init)])
figure 
plot(abs(yy(:,ccn)),'*')
title(['GPS=' num2str(svnum) 'Freq='  num2str(frr)])

format
pt_init
format long e
frr