catch.m

很好的GPS信号的捕获和跟踪程序

%%初始状态设置%%%%%
svnum=input('enter satellite number = ');
intodat=10001;      
fs=5e6;             %采样频率
ts=1/fs;            %采样时间
n=fs/1000;          %1ms的数据点 
nn=[0:n-1];         %所有的数据点
fc=1.25e6;          %未考虑多普勒效应的中心频率
nsat=length(svnum); %要捕获的卫星数
%%%%%%%%%%输入数据文件%%%%%%%%%%%%%%
fid=fopen('D:\Program Files\MATLAB\R2008a\work\GPSRawData.dat','r');
fseek(fid,intodat-1,'bof');
x2=fread(fid,6*n,'schar');
yy=zeros(21,n);
%%%%%%%%%%%%%%执行捕获%%%%%%%%%%%%%%%
code=digitizg(n,fs,0,svnum);
xf=fft(x2(1:n)');
for i=[1:21];
    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')));
[amp ccn]=max(max(abs(yy)));
pt_init=ccn;
cfrq=fc+1000*(crw-11);
%%%%%%%%%%%%%%%%%%%%% 产生5ms剥离了C/A码的数据%%%%%%%%
z5=x2(pt init:pt_init+5*n-1);
za5=z5'.*[code code code code code];

%%%%%%%%%%%%%%%找出分辨率400khz的中频%%%%%%%%%%%%%%
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);
mfrq=cfrq+200*(mrw-2);
fr=mfrq;
%%%%%%%%%%%%%%%%%找出精频%%%%%%%%%
zb5=za5.*exp(j*2*pi*fr*ts*[0:5*n-1]);
zc5=diff(-angle(sum(reshape(zb5,n,5))));
zc5fix=zc5;
%%%%%%%%%%%%%%%%%%%调整相位，去除可能的相位偏移%%%%%%%%%%%55
threshold=2.3*pi\5;
for i=1:4;
    if abs(zc5(i))>threshold;
        zc5(i)=zc5fix(i)-2*pi;
        if abs(zc5(i))>threshold;
            zc5(i)=zc5fix(i)+2*pi;
            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
dfrq=mean(zc5)*1000/(2*pi);
frr=fr+dfrq;
plot(abs(yy(crw,1:n)));
title(['GPS = 'num2str(svnum) 'max at' num2str(pt init)]);
figure
plot(abs(yy(:,ccn)),'*');
format
pt init
format long e
frr

%%%%%%%%%%%%%%%%%%%%%%%%digitizg.m生成C/A码并数字化%%%%%%%
function code2=digitizg(n,fs,offset,svnum);
%code - 戈等码
%n - 采样点
%fs - 采样频率
%offset - 延迟时间，必须小于1/fs,才不会产生偏移
%svnum - 卫星编号
gold rate =1.023e6;
ts=1/fs;
tc=1/gold rate;
cmd1=codegen(svnum);
code in =cmd1;
%%%%%%%%%%%%%%%%形成数字化所用的16个C/A码%%%%%%%%
code a =[code in code in code in code in];
code a =[code a code a];
code a =[code a code a];

%%%%%%%%%%%%%%%%%%%%%%%数字化%%%%%%%%%%%%%%
b=[1:n];
c= ceil((ts*b+offset)/tc);
code=code a(c);

%%%%%%%%%%%%%%调整初始点%%%%%%%%%%%%%%%%%%%

if offset>=0;
    code2=[code(1) code(1:n-1)];
else
    code2=[code(n) code(1:n-1)];
end

%%%%%%%%%%%%%%%%%%codegen.m产生32个C/A码中的一个%%%%%%%%%%%%%%
function[ca used]=codegen(svnum);
%caused ; 包含期望输入序列的向量
%g2s:向量存放g2码正确的向量
%C/A码：（if sv#19，使用G2为g2s(19)=471的偏离）               
%svnum：卫星编号
            
g2s=[5;6;7;8;17;18;139;140;141;251;252;254;255;256;257;258;269;470;471;472;473;474;509;512;513;514;515;516;859;860;861;862];
g2shift=g2s(svnum,1);
%%%%%%%%%%%%%产生G1码中的一个%%%%%%%%%%%%%%%%
reg = -1*ones(1,10);
for i=1:1023,
    gl(i)=reg(10);
    save1=reg(3)*reg(10);
    reg(1,2:10)=reg(1:1:9);
    reg(1)=save1;
end
%%%%%%%%%%%%%%%产生G2码中的一个%%%%%%%%%%%%%%%%%%%%

reg=-1*ones(1,10);
for i=1:1023,
    g2(i)=reg(10);
    save2=reg(2)*reg(3)*reg(6)*reg(8)*reg(9)*reg(10);
    reg(1,2:10)=reg(1:1:9);
    reg(1)=save2;
end

%%%%%%%%%%%%%%%%5移动 G2码%%%%%%%%%%%

g2tmp(1,1:g2shift)=g2(1,1023-g2shift+1:1023);
g2tmp(1,g2shift+1:1023)=g2(1,1:1023-g2shift);
g2=g2tmp;

%%%%%%%%%%%%%%%%%%通过G1与G2相乘形成单个采样C/A码%%%%%%%%%%%%%%%
ss_ca = g1.*g2;
caused=-ss_ca;