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