GPS信号跟踪

GPS课程设计实验报告（3）学院姓名LSC班级学号指导教员一、试验名称：GPS信号跟踪二、试验目的：1. 熟悉GPS信号跟踪的基本概念；2. 掌握载波跟踪环和码跟踪环的基本思想及处理流程；3. 训练利用课堂所学知识解决实际问题的能力。

三、试验内容1. 编写GPS信号跟踪子程序，其中载波环及码环阶数自定，鉴别器算法自选，通道数自选（简单起见，可为1）。

2. 在实验二的基础上，调用以上信号跟踪子程序对信号进行跟踪处理，解调出导航电文并作图。

3. 画出每一个通道码环鉴别器和载波环鉴别器的输出随时间变化的曲线。

四、试验原理4.1载波和码元跟踪GPS信号的跟踪示意图如下，传统锁相环的输入通常是连续波形或频率调制信号，压控振荡器VCO的频率用来跟随输入信号的频率。

在一台GPS接收机中，输入是GPS信号，接收机的锁相环必须跟随（或跟踪）这个信号。

然而，GPS信号是双向编码信号。

如以前提到的，GPS接收机和GPS卫星的运动产生多普勒效应，使得载波和码元频率会发生改变。

为了跟踪GPS信号，必须首先剥离C/A 码。

结果，就需要两个环路来跟踪GPS信号，一个用来跟踪C/A码，另一个用来跟踪载波频率。

这两个环路必须成对一起使用。

捕获程序找到了C/A码的起始点，码元环路产生相位超前C/A码近似1/2个基波，同时产生相位滞后近似1/2个基波的滞后码。

超前码与滞后码同输入C/A码进行卷积运算，产生两个结果，这两个卷积结果经过动态均值滤波器之后再经过平方运算，比较得到的这两个卷积结果的平方项，用以产生一个本地C/A 码速率的控制信号，使之与输入信号的C/A码匹配。

当C/A码从输入信号中剥离后，载波频率环路接收到的是一个相位只被导航电文调制的连续信号。

捕获过程得到载波频率的初始值，压控振荡器VCO根据捕获过程得到的初始值产生一个载波频率，这个载波频率信号分成两路：一路同相信号和一路正交相信号(相位偏移90°)。

这两路信号与输入信号相关运算，相关结果经过滤波器后，通过反正切比较器，比较它们之间的相位，用以调节本地振荡器的载波频率，来跟随输入的连续信号。

产生的载波频率同时用来从输入信号中剥离载波。

4.2利用锁相环跟踪GPS信号4.2.1二阶锁相环：二阶锁相环意味着传递函数H(s)的分母是s的二阶函数。

构成这种二阶锁相环的传递函数：2s s 101)t (t 44t 81n n s n k k C ωζωζω++= 2s s 2n 102)t (t 44t 41n n s k k C ωζωω++=）（ 其中，ωn 表示固有频率，ζ是抑制因子 (1) 设置码元环路和载波环路的带宽与增益。

环路增益包括相位探测器增益和VCO 增益。

(2)抑制因子ζ=0.707， 此值常被认为是最理想值。

(3) 由上式求出固有频率。

(4) 选择码元环路增益(k0k1)为50，载波环路增益为4π×100，这个设置也是若干种设置中的一种，常量C1、C2可由上式求。

这四步为两个环路提供了必要参数，通过调整码元环路的相位和载波频率相位，使其跟踪输入信号。

因为C/A 码是1 ms 长，因此，环路要在每毫秒更新一次，C/A 码也必须每毫秒重新产生一次， 且新产生的C/A 码相位必须与前一时刻相位连续。

典型输出数据设置如图所示， 输出幅值随时间变化，由于跟踪环的传输效应，最终使幅值达到稳定。

五、Matlab 代码及仿真结果5.1追踪子程序function[diata_dll,diata_fll]=tracksignal(iniphcode,inifd,iniph,snr,Code_Meth od_flag,Carrier_Method_flag,codew,codeb,carrierw,carrierfllb,carrierp llb)svnum = 10; %卫星号iniphcode = 120; %生成信号源的码相位inifd = 4600; %生成信号源的载波多普勒频率iniph = 0; %生成信号源的载波初相位snr = 0; %生成信号源的信噪比global time_unit; % 数据跳变时间单位global time; % 数据发送时间global time_cyc;% 一个完整扩频码周期global fs; % 采样率110def τωk k n 22211012102110121022def )(n n n n s s s k k s k k s k k s k k s H ωξωωξωττττττ++++++=global nn; % 一个完整扩频周期采样点数global kk; % 数据总采样点global F_if; % 载波中频global CA_freq; % PN码速率global tc;global CA ; % 扩频码基玛global F_Carrier; % L1波段载波频率%参数设置time = 100*(10^(-3));time_unit = 20*(10^(-3));time_cyc = 1*(10^(-3));fs = 5*(10^6);nn = time_cyc*fs;kk = (time/time_cyc)*nn;F_if = 1.25*(10^6);F_Carrier = 1575.42*(10^6);CA_freq = 1.023*(10^6);%生成C/A以供使用（未用试验一的C/A码）PN = codegen(svnum);CA = [];k = 5;for n = 1:length(PN)CA((1+k*(n-1)):k*n) = PN(n)*ones(1,k);endtc = 1/(k*CA_freq);loop_time = time/time_cyc;%生成信号源 %Signal_Source是生成的信号源, buffer_bit_data是随机生成的数据位，用于与最后解调的数据进行比对[Signal_Source,Phase_signal,buffer_bit_data]=CreateSource(iniphco de,inifd,iniph,snr);%捕获[fd_ac,f_ac_code] = acqu(Signal_Source);%定义跟踪中用到的参数IPSum_old = 0.01;QPSum_old = 0.01;%码跟踪环滤波器参数设置Code_Method_flag=2; %码跟踪环鉴相法选择标志diffoffside = 0.5; %码鉴相时，正负偏移半个码片k0 =10^(-6); %码跟踪环鉴相器增益k1=10^(-3); %码跟踪环NCO增益k1 = 50/k0; %码跟踪环NCO增益codew = 20 %自然圆频率codeb = 2 %阻尼系数offside = f_ac_code;theta_code_old = 0;offside_old = f_ac_code;CodeErr_old = 0;Bk_DLL = [];Track_Code_Buffer =[];%载波跟踪环滤波器参数设置Last_Phase = 0;Control_Buffer = [];ts = 1/fs; %采样时间间隔Carrier_Method_flag =3; %fll,pll,fll->pll的方法选择标志dem_flag = 0; %fll->pll的切换标志add = 0; %fll->pll的切换过程中用到的变量carrierw = 20 %自然圆频率carrierfllb = 0.707 %阻尼系数carrierpllb = 0.706Track_Freq_Buffer = [];track_dopplar_old = 0;%FLL环参数FLLinput_old=0;FLLoutput_old=0;track_freq_fll = 0;Sita_fll = 0;Bk_FLL = [];%PLL环参数PLLinput_old=0;PLLoutput_old=0;track_freq_pll = 0;Sita_pll = 0;Bk_PLL = [];Buffer_Data =[];adj_buffer = [];ALL_Buffer_Data = [];count_buffer = [];Demodulate_I = [];Local_Ph_Buffer = [];%跟踪循环for i = 1:1:loop_timeSignal = Signal_Source((i-1)*nn+1:i*nn);Signalph= Phase_signal((i-1)*nn+1:i*nn);t = [0:nn-1]*ts;track_dopplar = fd_ac+track_freq_fll+track_freq_pll;Track_Freq_Buffer = [Track_Freq_Buffer track_dopplar]; Local_I = cos(2*pi*(F_if+track_dopplar)*t + Last_Phase); Local_Q = sin(2*pi*(F_if+track_dopplar)*t + Last_Phase); Iph = 2*pi*(F_if+track_dopplar)*t + Last_Phase;Local_Ph_Buffer = [Local_Ph_Buffer Iph];Last_Phase=Last_Phase2*pi*(F_if+track_dopplar)*time_cyc; %%上一次积分结束点的相位Carrier_I = Local_I;Carrier_Q = Local_Q;ph_code_p = offside;fd_code_p = track_dopplar;CA_Code_p = CAcode(ph_code_p,fd_code_p,i);lc_p = CA_Code_p.*Signal;ph_code_e = offside+diffoffside;fd_code_e = track_dopplar;CA_Code_e = CAcode(ph_code_e,fd_code_e,i);lc_e = CA_Code_e.*Signal;ph_code_l = offside-diffoffside;fd_code_l = track_dopplar;CA_Code_l = CAcode(ph_code_l,fd_code_l,i);lc_l = CA_Code_l.*Signal;Local_P_I = lc_p.*Carrier_I;Local_P_Q = lc_p.*Carrier_Q;Local_E_I = lc_e.*Carrier_I;Local_E_Q = lc_e.*Carrier_Q;Local_L_I = lc_l.*Carrier_I;Local_L_Q = lc_l.*Carrier_Q;IPSum = sum(Local_P_I);QPSum = sum(Local_P_Q);IESum = sum(Local_E_I);QESum = sum(Local_E_Q);ILSum = sum(Local_L_I);QLSum = sum(Local_L_Q);if Code_Method_flag==1%鉴相器theta_code = k0*((IESum.^2+QESum.^2)-(ILSum.^2+QLSum.^2));%码环路滤波器[CodeErr] = CodeloopFilter(codew,codeb,theta_code,theta_code_old,CodeErr_old); % CodeErr是经过滤波器输出的码相位误差的估计值%码环NCOoffside=offside_old+60*k1*CodeErr; %码NCO的输出theta_code_old = theta_code; %将当前结果保存，用于下一个循环的码跟踪CodeErr_old = CodeErr; %将当前结果保存，用于下一个循环的码跟踪offside_old = offside; %将当前结果保存，用于下一个循环的码跟踪Bk_DLL = [Bk_DLL theta_code]; %记录跟踪过程中的码环鉴想器的输出Track_Code_Buffer = [Track_Code_Buffer offside]; %记录跟踪过程中的码环NCO的数出elseif Code_Method_flag==2%鉴想器theta_code = ((IESum.^2+QESum.^2)-(ILSum.^2+QLSum.^2))/((IESum.^2+QESum.^2)+(ILSum .^2+QLSum.^2));theta_code=10.^(-3)*(1-sqrt(1-theta_code.^2))/(2*theta_code);%码环路滤波器[CodeErr] = CodeloopFilter(codew,codeb,theta_code,theta_code_old,CodeErr_old); % CodeErr是经过滤波器输出的码相位误差的估计值%码环NCOoffside=offside_old+500*CodeErr; %码NCO的输出theta_code_old = theta_code; %将当前结果保存，用于下一个循环的码跟踪CodeErr_old = CodeErr; %将当前结果保存，用于下一个循环的码跟踪offside_old = offside; %将当前结果保存，用于下一个循环的码跟踪Bk_DLL = [Bk_DLL theta_code]; %记录跟踪过程中的码环鉴想器的输出Track_Code_Buffer = [Track_Code_Buffer offside]; %记录跟踪过程中的码环NCO的数出endif Carrier_Method_flag == 1 % fll跟踪环路%鉴想器real_Q = IPSum_old*QPSum-QPSum_old*IPSum;real_I = IPSum_old*IPSum+QPSum_old*QPSum;a=real_Q/real_I;theta_fll = atan(real_Q/real_I);Bk_FLL = [Bk_FLL theta_fll];FLLinput = theta_fll/(2*pi*time_cyc);%环路滤波器FLLoutput = CarrierLoopFilter(carrierw,carrierfllb/2,FLLinput,FLLinput_old,FLLout put_old);%FLL环NCOSita_fll = Sita_fll+FLLoutput;track_freq_fll = -Sita_fll; %FLL环跟踪到的多普勒频率 (由于反正切主值区间造成的)FLLinput_old=FLLinput; %将当前结果保存，用于下一个循环的载波跟踪FLLoutput_old=FLLoutput; %将当前结果保存，用于下一个循环的载波跟踪IPSum_old = IPSum; %将当前结果保存，用于下一个循环的载波跟踪QPSum_old = QPSum; %将当前结果保存，用于下一个循环的载波跟踪elseif Carrier_Method_flag == 2 %costa跟踪环路theta_pll = atan(QPSum/IPSum);PLLinput = theta_pll/(2*pi*time_cyc);Bk_PLL = [Bk_PLL theta_pll];%LoopFilterPLLoutput=CarrierLoopFilter(carrierw,carrierpllb/2,PLLinput,PLLinput_old,PLLout put_old);track_freq_pll = -PLLoutput;PLLinput_old=PLLinput; %将当前结果保存，用于下一个循环的载波跟踪PLLoutput_old=PLLoutput; %将当前结果保存，用于下一个循环的载波跟踪elseif Carrier_Method_flag == 3 %fll跟踪环路->costa if dem_flag == 0real_Q = IPSum_old*QPSum-QPSum_old*IPSum;real_I = IPSum_old*IPSum+QPSum_old*QPSum;theta_fll = atan(real_Q/real_I);FLLinput = theta_fll/(2*pi*time_cyc);Bk_FLL = [Bk_FLL theta_fll];%LoopFilterFLLoutput = CarrierLoopFilter(carrierw,carrierfllb/2,FLLinput,FLLinput_old,FLLout put_old);;%NCOSita_fll = Sita_fll+FLLoutput;track_freq_fll = -Sita_fll;FLLinput_old=FLLinput; %将当前结果保存，用于下一个循环的载波跟踪FLLoutput_old=FLLoutput; %将当前结果保存，用于下一个循环的载波跟踪IPSum_old = IPSum; %将当前结果保存，用于下一个循环的载波跟踪QPSum_old = QPSum; %将当前结果保存，用于下一个循环的载波跟踪elseif dem_flag == 1theta_pll = atan(QPSum/IPSum);PLLinput = theta_pll/(2*pi*time_cyc);Bk_PLL = [Bk_PLL theta_pll];%LoopFilterPLLoutput = CarrierLoopFilter(carrierw,carrierpllb/2,PLLinput,PLLinput_old,PLLout put_old);track_freq_pll = -PLLoutput;PLLinput_old=PLLinput; %将当前结果保存，用于下一个循环的载波跟踪PLLoutput_old=PLLoutput; %将当前结果保存，用于下一个循环的载波跟踪endendadj_flag = track_dopplar - track_dopplar_old; %相邻两次跟踪到的多普勒频率值之差，用以判断是否FLL跟踪的频率已经足够精确,从而转入PLLtrack_dopplar_old = track_dopplar;adj_buffer = [adj_buffer adj_flag];outdata = sign(real(IPSum));ALL_Buffer_Data = [ALL_Buffer_Data outdata];if adj_flag < 1 %看相邻两次跟踪到的多普勒频率之差是否小于1Hz add = add+1;elseadd = 0;endif add >= 2 %看是否有连续两次跟踪到的多普勒频率之差小于1Hz，若有，则认为频率跟踪已很稳定而精确，可以转入PLLdem_flag = 1;endif dem_flag == 1count_time = i;count_buffer = [count_buffer count_time];Buffer_Data = [Buffer_Data outdata]endend%位同步与数据解调Buffer_Data_out = framecheck2(Buffer_Data,count_buffer); %位同步l_i_d = time/time_unit;l_o_d = length(Buffer_Data_out);l_zeros = l_i_d - l_o_d;Buffer_Data_out = [zeros(1,l_zeros) Buffer_Data_out];%计算跟踪精度Track_Code_Buffer;Track_Freq_Buffer;l_dll = length(Track_Code_Buffer);l_fll = length(Track_Freq_Buffer);diata_dll =sqrt(sum((Track_Code_Buffer(40:l_dll)-iniphcode).^2)/length(Track_Cod e_Buffer(40:l_dll))); %码跟踪环跟踪精度diata_fll = sqrt(sum((Track_Freq_Buffer(40:l_fll)-inifd).^2)/length(Track_Freq_Bu ffer(40:l_fll))); %载波跟踪环的跟踪精度%画图显示跟踪结果subplot(1,2,1);plot(Track_Code_Buffer);xlabel('时间(ms)');ylabel('码相位跟踪值');title('码跟踪结果');grid onsubplot(1,2,2);plot(Track_Freq_Buffer);xlabel('时间(ms)');ylabel('多普勒频率跟踪结果(Hz)')title('载波跟踪结果');grid onfigure;plot([1:length(buffer_bit_data)],buffer_bit_data,'b*',[1:length(B uffer_Data_out)],Buffer_Data_out,'ro');set(gca,'xtick',[1:1:5]);set(gca,'xticklabel',{'1','2','3','4','5'});set(gca,'ytick',[-1:1:1]);set(gca,'yticklabel',{'-1','0','1'});xlabel('数据位');ylabel('解调结果')title('数据解调输出结果');grid on5.2试验截图。