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