carrier_trace.m

包括载波跟踪和伪码相位跟踪

% 本程序用于实现扩频系统中，信号捕获之后的跟踪过程。

% 输入信号为，捕获阶段得到的多普勒频率和码相位
% 程序将进一步找到更精确的频率和相位，实现载波的跟踪
% 目前实际多普勒频率1000，相位86
function Carrier_Trace(fd_estimate, Correl_Peak_location_real)

% 获取系统参数
[code_cycle,SNR,Ts,T,T_interp,fi,M,N,fd,f0,delay_time,fd_estimate_init,sigma_n,fs,Bandwidth,dot_insert,coef,coef_mod,B,A1]=set_parameter;

%%%%产生接受中频信号伪码，抽头为1，8
code_phase=[1,8];
[inputcode]=code_gen(code_phase);   % 1023大小的GOLD伪码
%%%%本地伪码
native_code=inputcode;
inputcode_2=[native_code,native_code];
EmulateIndex = 0;
k = 0;

half_dot_insert = dot_insert/2;

while (1)

    % 模拟输入的中频信号
    EmulateIndex = EmulateIndex + 1;
    [signal,signal_noise]=SignalGenerator(EmulateIndex,inputcode,fs,Ts,fi,delay_time,fd,SNR,Bandwidth,code_cycle,sigma_n,f0,N);

    IF_signal(1:N)=signal(1:N);
    IF_signal_noise(1:N)=signal_noise(1:N);
%     % 绘制信号及噪声
%     figure(1);
%     plot(IF_signal, 'b');  % 绘制生成的信号(含噪声)
%     hold on;
%     plot(IF_signal_noise, 'r');  % 信号的噪声
%     hold off;

    % 当前的载波NCO ××× fd_estimate在不断调整
    [NCO_I_Quanti,NCO_Q_Quanti]=Nco_gen(Ts,fi,fd_estimate,N);
    
    % 下变频
	Rece_Signal_Quan_down_I = IF_signal(1:N).* NCO_I_Quanti;
	Rece_Signal_Quan_down_Q = IF_signal(1:N).* NCO_Q_Quanti;
	Rece_Signal_Quan_down_I_noise=IF_signal_noise.*NCO_I_Quanti;
	Rece_Signal_Quan_down_Q_noise=IF_signal_noise.*NCO_Q_Quanti;

    %%%%经过低通滤波器，截止频率为1.023MHz
    [signal_I,signal_Q,signal_I_noise,signal_Q_noise]=LP_filter(N,N+1,code_cycle,fs,Rece_Signal_Quan_down_I,Rece_Signal_Quan_down_Q,Rece_Signal_Quan_down_I_noise,Rece_Signal_Quan_down_Q_noise,EmulateIndex);       %进过低通滤波器
    
%     % 绘制比较，低通滤波器的作用
%     figure(2);
%     plot(Rece_Signal_Quan_down_I, 'b'); % 原始
%     hold on
%     plot(signal_I, 'r');    % 
%     hold off
    
    %%%%采用线性内插
	for index_new=1:dot_insert
        index_previous=floor((index_new-1)*T_interp/Ts)+1;                                                 %产生内插抽取所需的原序列索引值
        Rece_Signal_down_I_1(index_new)=signal_I(index_previous)+round(((signal_I(index_previous+1)-signal_I(index_previous)))*mod(((index_new-1)*coef),coef_mod)/coef_mod);
        Rece_Signal_down_Q_1(index_new)=signal_Q(index_previous)+round(((signal_Q(index_previous+1)-signal_Q(index_previous)))*mod(((index_new-1)*coef),coef_mod)/coef_mod);
	end%%内插完成     
    % 下变频之后的两路信号，保存到Rece_Signal_down_I_1，和Rece_Signal_down_Q_1中，4096点
    
    % 获取本地即时码，此处依然是连续的两对GOLD码，
    Code_native = code_native(inputcode_2,fd_estimate,T_interp,T,f0,dot_insert);

    k = k+1;
    Ips(k) = sum(Code_native(Correl_Peak_location_real:(Correl_Peak_location_real+half_dot_insert-1)).*Rece_Signal_down_I_1(1:half_dot_insert));
    Qps(k) = sum(Code_native(Correl_Peak_location_real:(Correl_Peak_location_real+half_dot_insert-1)).*Rece_Signal_down_Q_1(1:half_dot_insert));
    
	FdQueue(k) = fd_estimate;
	PhaseQueue(k) = Correl_Peak_location_real;
    
    if(k==1)
        continue;
    end

    % 载波环锁定判决
    E(k) = (Ips(k)*Ips(k-1) + Qps(k)*Qps(k-1))/(Ips(k)*Ips(k) + Qps(k)*Qps(k));
    % 低通滤波-平滑E(k)
    y = 0.9*E(k-1) + 0.1*E(k);
    E(k) = y;
    
    if(E(k) > 0.95)
        % 已经锁定载波
        % 进入码相位跟踪
        if(PN_Phase_Trace(fd_estimate, Correl_Peak_location_real) ~= 0)
            break;
        end;
    end;
    
    Tid = T;   % 积分清除的周期为1023个GOLD码占的时间
    % 频率判决
    % 生成结果如：0   52.2740  -42.6017  -52.2780  144.9547   -7.2020
    Delta(k) = (Ips(k-1)*Qps(k)-Ips(k)*Qps(k-1))/(2*pi*Tid*(Ips(k)*Ips(k)+Qps(k)*Qps(k)));    % 公式1
    %Delta(k) = (atan(Qps(k)/Ips(k))-atan(Qps(k-1)/Ips(k-1)))/(2*pi*Tid);    % 公式2
    
    if(abs(Delta(k)) > 10)
        % 频率牵引---四相鉴频器矫正
        if(abs(Ips(k)) >= abs(Qps(k)))
            Beta(k) = sign0(Ips(k))*(Qps(k)-Qps(k-1))/sqrt(Ips(k)*Ips(k)+Qps(k)*Qps(k));
        else
            Beta(k) = (-1)*sign0(Qps(k))*(Ips(k)-Ips(k-1))/sqrt(Ips(k)*Ips(k)+Qps(k)*Qps(k));
        end;
        fd_estimate = fd_estimate + Beta(k);
    else
        % 相位判决
        if(abs(Qps(k)/Ips(k)) > 0.176)
            % FLL 鉴频器，频率跟踪
            Efk = atan((Ips(k-1)*Qps(k)-Ips(k)*Qps(k-1))/(2*pi*Tid*(Ips(k)*Ips(k)+Qps(k)*Qps(k))));
            fd_estimate = fd_estimate + Efk;
        else
            % PLL
            Epk = atan2(Qps(k), Ips(k));
            fd_estimate = fd_estimate + Epk;
            %Correl_Peak_location_real = Correl_Peak_location_real + int32(Epk);
        end;
    end;

    if(mod(length(Delta),1000)==1)
        % 绘制多普勒频率跟踪的结果
        figure(81);
        plot(FdQueue);
        % 绘制跟踪过程中的频率判决
        figure(82);
        plot(E);
    end;
    
end