tracksignal.asv

GPS的CA码捕获跟踪源代码

function [diata_dll,diata_fll]=tracksignal(iniphcode,inifd,iniph,snr,Code_Method_flag,Carrier_Method_flag,codew,codeb,carrierw,carrierfllb,carrierpllb)

svnum = 10;    %卫星号
iniphcode  = 120;  %生成信号源的码相位
inifd = 4600;      %生成信号源的载波多普勒频率
iniph = 0;         %生成信号源的载波初相位
snr = 0;           %生成信号源的信噪比


global time_unit; % 数据跳变时间单位
global time; % 数据发送时间
global time_cyc;% 一个完整扩频码周期
global fs; % 采样率
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以供使用%%%%%%%%%%
PN = codegen(svnum);
CA = [];
k = 5;
for n = 1:length(PN)
CA((1+k*(n-1)):k*n) = PN(n)*ones(1,k);
end

tc = 1/(k*CA_freq); 
loop_time = time/time_cyc;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%生成信号源 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Signal_Source是生成的信号源, buffer_bit_data是随机生成的数据位，用于与最后解调的数据进行比对
[Signal_Source,Phase_signal,buffer_bit_data]=CreateSource(iniphcode,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;          %自然圆频率

Track_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_time   

     % =================part signal source========= %

     Signal = Signal_Source((i-1)*nn+1:i*nn);
%    Signalph= Phase_signal((i-1)*nn+1:i*nn);

     %==============Local carrier===================%
     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_Phase + 2*pi*(F_if+track_dopplar)*time_cyc;   %%上一次积分结束点的相位
     Carrier_I = Local_I;
     Carrier_Q = Local_Q;
     

     % =================Create Localcode============== %
     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;

    % ==========================multiply================================ %
     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;


     % ========================integration============================= %
     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);
 
     % =======================Code Control==================================
     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是经过滤波器输出的码相位误差的估计值
    
        %码环NCO
        offside=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是经过滤波器输出的码相位误差的估计值
    
        %码环NCO
        offside=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的数出
     end
     % ====================Carrier control=============================%
     if 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,FLLoutput_old);
     
        %FLL环NCO
        Sita_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];
        
        %LoopFilter       
        PLLoutput = CarrierLoopFilter(carrierw,carrierpllb/2,PLLinput,PLLinput_old,PLLoutput_old); 
        
        track_freq_pll = -PLLoutput;
        
        PLLinput_old=PLLinput;      %将当前结果保存，用于下一个循环的载波跟踪
        PLLoutput_old=PLLoutput;    %将当前结果保存，用于下一个循环的载波跟踪 

      elseif Carrier_Method_flag == 3        %%%%%%%%%%%%%%%%%%%%fll跟踪环路->costa
          if dem_flag == 0 
             real_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];
             
             %LoopFilter
             FLLoutput = CarrierLoopFilter(carrierw,carrierfllb/2,FLLinput,FLLinput_old,FLLoutput_old);;
             
             %NCO
             Sita_fll = Sita_fll+FLLoutput;
             track_freq_fll = -Sita_fll;
             
             FLLinput_old=FLLinput;      %将当前结果保存，用于下一个循环的载波跟踪
             FLLoutput_old=FLLoutput;    %将当前结果保存，用于下一个循环的载波跟踪 
             IPSum_old = IPSum;          %将当前结果保存，用于下一个循环的载波跟踪
             QPSum_old = QPSum;          %将当前结果保存，用于下一个循环的载波跟踪
             
          elseif dem_flag == 1
             theta_pll = atan(QPSum/IPSum);
             PLLinput = theta_pll/(2*pi*time_cyc);
             Bk_PLL = [Bk_PLL theta_pll];
             
             %LoopFilter
             PLLoutput = CarrierLoopFilter(carrierw,carrierpllb/2,PLLinput,PLLinput_old,PLLoutput_old); 
             
             track_freq_pll = -PLLoutput;
             
             PLLinput_old=PLLinput;      %将当前结果保存，用于下一个循环的载波跟踪
             PLLoutput_old=PLLoutput;    %将当前结果保存，用于下一个循环的载波跟踪 
          end
     end

     adj_flag = track_dopplar - track_dopplar_old;  %相邻两次跟踪到的多普勒频率值之差，用以判断是否FLL跟踪的频率已经足够精确,从而转入PLL
     track_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;
     else
         add = 0;
     end
     
     if add >= 2          %看是否有连续两次跟踪到的多普勒频率之差小于1Hz，若有，则认为频率跟踪已很稳定而精确，可以转入PLL
         dem_flag = 1;
     end
    
     if dem_flag == 1
         count_time = i;
         count_buffer = [count_buffer count_time];
         Buffer_Data = [Buffer_Data outdata]
     end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%位同步与数据解调%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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_Code_Buffer(40:l_dll)));  %码跟踪环跟踪精度
diata_fll = sqrt(sum((Track_Freq_Buffer(40:l_fll)-inifd).^2)/length(Track_Freq_Buffer(40:l_fll)));      %载波跟踪环的跟踪精度

%Track_Ph_Buffer = Local_Ph_Buffer-Phase_signal;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画图显示跟踪结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(1,2,1);
plot(Track_Code_Buffer);
xlabel('时间(ms)');
ylabel('码相位跟踪值');
title('码跟踪结果');
grid on

subplot(1,2,2);
plot(Track_Freq_Buffer);
xlabel('时间(ms)');
ylabel('多普勒频率跟踪结果(Hz)')
title('载波跟踪结果');
grid on

figure;
plot([1:length(buffer_bit_data)],buffer_bit_data,'b*',[1:length(Buffer_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('数据解调输出结果');
legend('数据','解调输出数据');
grid on


% figure;
% plot(Bk_DLL);
% title('Track Code 输入控制')
% grid on

% figure;
% plot(Track_Ph_Buffer);
% title('Track differ Phase')
% grid on