p7_1.asv

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

%p7 1.m **performs acquisition on collected data对数据进行捕获***
clear
%*******初始化数据*******
svnum=input('enter satellite number=');  %接收一个卫星号存入svnum
intodat=10001;%输入读取数据起始点
fs=5e6;%***采样频率
ts=1/fs;%***采样时间
n=fs/1000;%***数据块长度1ms
nn=[0:n-1];%***total no.of pts
fc=1.25e6;%***无多普勒频移的中心频率
nsat=length(svnum);%***total number of satellites to be processed

%*******读取数据文件*********
fid=fopen('G:\GNSS\GNSS_signal_records\GPS_and_GIOVE_A-NN-fs16_3676-if4_1304.bin','r');
fseek(fid,intodat-1,'bof');
x2=fread(fid,12*n,'schar');%从文件的第inodat个字节处读取6*n个'signed char'类型数据
yy=zeros(21,n);

%*******开始捕获*******
%     2.25 @ 2009       %
%**********************
code=digitizg(n,fs,0,svnum);%量化 C/A 码
xf=fft(x2(1:n)');
for i=[1:21];%****find coarse freq 1kHz 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')));%找出最高峰值
[amp ccn]=max(max(abs(yy)));%ccn为yy中最大值所在列数
pt_init=ccn;%初始点
cfrq=fc+1000*(crw-11);%捕获的粗略频率

%*******产生5ms剥离C/A码的数据*******
z5=x2(pt_init:pt_init+5*n-1);%从输入数据的C/A码起始点读取5ms
za5=z5'.*[code code code code code];%create cw from 5 seconds of data

%*******计算分辨率为400HZ的中频频率*******
for i=[1:3];
fr=cfrq-400+(i-1)*400;%在粗捕获点上移动步长400HZ
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;%中频频率

%*******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;
%*******调整相位并获取可能的相移*******
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;%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)])
figure

plot(codegen(svnum));
axis([0 1200 0 1.5])
format
pt_init
format long e
frr