costas_simulation_bak.m

本程序是用COstas环跟踪信号的处理过程

%--------------------------------------------------------------------------
% The goal of this simulation is to simulate carrier tracking loop with the
% way of costas loop. 
% Author:   Meng, xianwei, copyright of sun-create
% Data:     2006.9.12
%--------------------------------------------------------------------------
clear all;
clc;
%------------simulate costas loop-------------

% disp(' ');

fdel = 30;          %input('Enter the size of frequency step in Hertz>');%

fn = 50.7143;       %input('Enter  loop filter natural frequency in Hertz>');%10
lambda = 0.01;      %input('Enter lambda, the relative pole offset>');%0.1
zeta = 1/sqrt(2);

%------采样频率------------------
fs = 4000*1.023;
tstop = 18;          %----结束时间----

npts = fs * tstop + 1;
t = (0:(npts-1))/fs;

nsettle = fix(npts/10);
tsettle = nsettle/fs;

freq = fs;
fin  = [zeros(1, nsettle), fdel * ones(1,npts-nsettle)];    %--输入的频率

phin = [zeros(1, nsettle), 2*pi*fdel*t(1:npts - nsettle)];  %--输入的相位

len = length(phin);
Acostx = zeros(1,len);
carrier_freq = 2000;
% nsettle = 0;
nsettle = 64;
randn('state',0);
for k=1:len
    if ( k<=nsettle )
        Acostx(k) = cos(2*pi*carrier_freq*t(k) + 0);
    else
%         if(k / 100 =0)
            
            Acostx(k) = cos(2*pi*carrier_freq*t(k) + 0.008*randn(1)*fdel);%*(-1);
%         end;
    end
end

% Generate original signal xref and a noisy version x adding 
% a standard Gaussian white noise. 
% sqrt_snr = 8; init = 2055615866;
%-----信号加入噪声---------------------------------------------

Acost = Acostx + 0.2*randn(size(t));;%
dumptime =90;          %积分清零时间频率计数门限
Carrier_Tracking_loop(Acost,carrier_freq,fs,t,dumptime);

% % Acost(11001:19000) = Acost(11001:19000) .* ones(1,8000) * (-1);
% % Acost(71001:79000) = Acost(71001:79000) .* ones(1,8000) * (-1);
% % Acost(91001:99000+10000) = Acost(91001:99000+10000) .* ones(1,8000+10000) * (-1);
% %Acost((91001+10000):99000+10000) = Acost((91001+10000):99000+10000) .* ones(1,8000+10000) * (-1);
% %---------initialize---------
% w2b = 0; w2c = 0; s5 =0; phivco =0 ;
% twopi = 2 * pi;
% twofs = 2 * fs;
% 
% %--define 2*fs;
% 
% %---set loop gain-----------zeta = 1/sqrt(2);
% G = pi*fn*(zeta + sqrt(zeta*zeta - lambda));
% 
% %---set filter parameter----Wn = 2*pi*fn
% a = 2*pi*fn/(zeta + sqrt(zeta*zeta - lambda));
% 
% %---define constant---
% a1 = a*(1-lambda); a2 = a*lambda;
% 
% %---初始化向量----
% % phierror = zeros(1,npts);
% % fvco = zeros(1, npts);
% % vco_sin = zeros(1, npts);
% % vco_cos = zeros(1, npts);
% 
% phivco = 0;lpw2c = 0; Qpw2c=0;
% 
% dumptime = 90;
% dumptime_counter = 0;
% vertor_len = 0;
% 
% %---------costas loop simulation---------
% for i =1:npts
%     dumptime_counter = dumptime_counter + 1;
%     if (dumptime_counter>dumptime)
%         dumptime_counter = dumptime_counter - dumptime;
%     end;
%     vco_sin(i) = sin(2*pi*carrier_freq.*t(i) + phivco);
%     vco_cos(i) = cos(2*pi*carrier_freq.*t(i) + phivco);
%     
%     I_cos = vco_cos(i) * Acost(i);
%     Q_sin = vco_sin(i) * Acost(i);
%     s11 = Q_sin;
%     Qs11= I_cos;
%     %mul_tmp_Q(i) = Q_sin;
% 
%     %------低通滤波器-------Integrate and dump-------
%     Qps6  = Qs11;             %--V输入
%     Qpw1c = Qps6 + Qpw2c;   %-积分步骤1
%     Qpw2c = Qps6 + Qpw1c;   %--积分步骤2
%     Qs11 = Qpw1c;%/twofs;      %--V输出
% 
%     %------低通滤波器-------------------------------
%     lps6  = s11;             %--V输入
%     lpw1c = lps6 + lpw2c;   %-积分步骤1
%     lpw2c = lps6 + lpw1c;   %--积分步骤2
%     s11 = lpw1c;%/twofs;      %--V输出
%     
%     %------积分清零，信号抽取----------------------
%     if (dumptime_counter == dumptime)
%         vertor_len = vertor_len + 1; 
%         s11 = s11/dumptime;
%         Qs11= Qs11/dumptime;
%         Qft_tmp(vertor_len) = Qs11;
%         ft_tmp(vertor_len) = s11;    
%         %------------------------------------------
%    %---------------鉴相器---------------------------------
%         s2 = ft_tmp(vertor_len) * Qft_tmp(vertor_len);
%         s3 = G*s2;              %--滤波器增益乘信号
%         s4 = a1*s3;
%     
%         s4a = s4 - a2*s5;       %--环路滤波器积分输出
%         w1b = s4a + w2b;        %--滤波器积分步骤1
%         w2b = s4a + w1b;        %--滤波器积分步骤2
%         s5  = w1b/twofs;        %--产生滤波器输出
%    %----------VCO---------------------------------------- 
%         s6  = s3 + s5;          %--VCO积分输入
%         w1c = s6 + w2c;         %--VCO积分步骤1
%         w2c = s6 + w1c;         %--VCO积分步骤2
%         phivco = w1c/twofs;     %--VCO积分输出
%     
%         phierror(i) = s11;       %--相位差矩阵
%         fvco(vertor_len) = s6/twofs;     %--VCO输入矩阵
%         
%         %-------积分清零-----------------------
%         Qs11 = 0;
%         s11  = 0;
%         Qpw2c =0;
%         lpw2c = 0;
%     end;
%  
%     
% end;
% 
% % freqerror = fin - fvco;     %--频差矩阵/fs
% 
% plot([1:vertor_len]*dumptime./fs, Qft_tmp);
% hold on;
% plot([1:vertor_len]*dumptime./fs, ft_tmp,'r')
% grid on;
% xlabel('时间（t）');
% ylabel('I和Q路幅值');
% title('载波跟踪环（costas loop）模拟');
% set(gcf, 'Color', [1 1 1])




% %-----------------------Global variable defination-------------------------
% 
%        PRN_Fs = 1023000;
%        NavData_Fs = 50;
%        sampleTimes = 18;
%        Carrier_Frequency = 3000000;
%        Carrier_doppleFre_Shift = 0;             %4000;
%        Phase_Diff = pi/9;
%        Pha_pos = 1023;                          % C/A Code Phase
%        len = 120;                               %Signal plot length
% %--------------------------- Create signal.--------------------------------
% 
%        PRN_num = 1023;
%        
%        PRN_temp = double(0: (1/(sampleTimes*PRN_Fs)) :PRN_num*(1/(PRN_Fs)));
%        PRN_t = PRN_temp(1:(length(PRN_temp)-1));
%        %--------------------------GPS PRN signal---------------------------
%  
%        
%        %------------------------carrier signal-----------------------------
%        %Fs = 4000;   t = 0:1/Fs:1;
%        Carrier_sample_time = PRN_t;
%        Carrier = cos(2 * pi * (Carrier_Frequency + Carrier_doppleFre_Shift)...
%                         * Carrier_sample_time + Phase_Diff); 
%        %-----plot--
%        figure;
%        an=[1:21];
%        stem(an.*(1/(sampleTimes*PRN_Fs)),Carrier(1:21));
%        hold on; 
%        plot(an.*(1/(sampleTimes*PRN_Fs)),Carrier(1:21),'r');
%        grid on;
%        if (sampleTimes*1000000/Carrier_Frequency>2)
%             title('Carrier Signal Over Sampling Rate');
%        else
%            title('Carrier Signal UnderSampling Rate');
%        end;