costas.m

一个关于科斯塔斯环的MATLAB仿真程序

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);