costas.m

该程序包含一般数字PLL仿真maltab程序和数字COSTAS环程序

function [out,Fre,Alt_phase,Alter_frequency,phase_discri] = Costas(fs,T,fc,Signal); 
% [out,Alter_frequency,phase_discri] 
% fs 信号采样率 
% T 处理数据的时间（单位:毫秒） 
% Signal 去除P码后的纯载波信号 
% 使用方法：[out,Fre,Alt_phase,Alter_frequency,phase_discri] = Costas(81.84e+6,50,15.345e+6+100,data); 
% Costas(81.84e+6,50,15.345e+6+100,data); 
f = 10.23e+6; % 基准频率 
ts = 1/fs;    % time delay 
Crate = fs/f; % Crate倍的采样率 
num = T * 10230 * Crate;   % 总共要处理的点数 
Point_Num = floor(fs/1000);% 每次积累处理的点数  floor  
LocalUsed_n = [0:Point_Num-1]; 
MaxCount = floor(length(Signal)/Point_Num); % 循环处理的次数 
Pifc = 2*pi*fc;  
temp = 0; 
phase = 0; 
[c1,c2] = GetC1C2(0.707,30); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
load PCodes_U PCodes_U;                       %%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for i=1:MaxCount  
    exp_Local = exp(j*(Pifc*ts*LocalUsed_n+phase)); % 本地信号的复数形式 
    sin_Signal = imag(exp_Local); % 本地信号 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    sin_Signal = PCodes_U((1:Point_Num)+((i-1)*Point_Num)) .*sin_Signal; %%%%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    cos_Signal = real(exp_Local); % 本地信号的pi/2相移 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    cos_Signal = PCodes_U((1:Point_Num)+((i-1)*Point_Num)) .*cos_Signal;%%%%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    Used_Signal = Signal((1:Point_Num)+((i-1)*Point_Num));% 本次处理使用的信号 
    Q_Sin = Used_Signal .* sin_Signal; 
    I_Cos = Used_Signal .* cos_Signal; 
    
    Q = sum(Q_Sin); 
    I = sum(I_Cos); 
    if(i==1 && Q<0 && I<0) 
        phase = phase + pi; 
    end; 
    out(i) = I; 
    phase_discri(i) = atan(Q/I); 
    % phase_discri(i) = atan(I/Q); 
    % out((1:Point_Num)+((i-1)*Point_Num)) = Used_Signal.*sin(Pifc*ts*LocalUsed_n +phase_discri(i)); 
    Frequency_dis = c1*phase_discri(i)+temp; 
    temp = temp+c2*phase_discri(i); 
    Pifc = Pifc-Frequency_dis*2*pi;  
    Alter_frequency(i) = Pifc;  
    phase = Pifc*ts*Point_Num+phase; 
    Alt_phase(i) = mean(phase); 
end 
figure; 
plot(phase_discri); 
figure; 
plot(Alter_frequency/(2*pi)); 
Fre = Alter_frequency/(2*pi); 
% b_c1=fir1(3300,1/(2.0460e+007));  
% out = fftfilt(b_c1, out); 
figure ; 
plot(out);