clf;
%fm=100;f0=2.4e10;B=2.4e8;FR=320R;FD=160v;
Fs = 130000;            % Sampling frequency                    
Ts = 1/Fs;             % Sampling period       
N = 512;             % Length of signal
n=0:N-1;
%t1 = (0:1:N-1)*Ts;        % Time vector
t1=n/Fs;
%y = wgn(m, n, p);
% Form a signal containing a 50 Hz sinusoid of amplitude 0.7 and a 120 Hz sinusoid of amplitude 1.
S = sqrt(2)*sqrt(5)*sin(2*pi*50000*t1) + sqrt(2)*sqrt(3)*sin(2*pi*14400*t1);
% Corrupt the signal with zero-mean white noise with a variance of 4.
NOISE=wgn(1, 512, 3);
X = S+NOISE;
sigPower = sum(abs(S).^2)/length(S);           %求出信号功率
noisePower=sum(abs(X-S).^2)/length(X-S); %求出噪声功率
SNR=10*log10(sigPower/noisePower) ;
% Plot the noisy signal in the time domain. It is difficult to identify the frequency components by looking at the signal X(t).
 %对信号进行快速Fourier变换
Y = fft(X); %求得Fourier变换后的振幅
mag1=abs(2*Y/N);
P1 = mag1(1:N/2);
P1=(abs(P1));
%P1(2:end-1) = 2*P1(2:end-1);
f1=(1:N/2)*Fs/N;    %频率序列
figure(1);
plot(1000*t1,X);title('');xlabel('t (milliseconds)');ylabel('X(t)');
figure(2);
plot(f1,P1,'r');  xlabel('频率/Hz');ylabel('dbv');title('： FFT','color','r');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CA-CFAR %%%%%%%%%%%%%%%%%%%%%%%%% 
T_std=10;
M=N/2;P_fa=10.^(-5); 
R=16;%参考单元数 
n=R/2; 
L_slipper=n;%滑窗长度 
L_move=1;%滑窗间隔 
%L_num=floor((M-L_slipper)/L_move)-1;%滑窗次数 
L_num=251;
Z=zeros(1,L_num); 
Z1=zeros(1,n); 
Z2=zeros(1,n);
S1=zeros(1,256-n-3);
%%%%%%%%%
for i=n+3:256-n-3
for j=0:n-1
%Z(i)=P1(i+3+j).^2+P1(i-3-j).^2;
Z1(j+1)=P1(i-3-j).^2;%前窗
Z2(j+1)=P1(i+3+j).^2;
Z11=mean(Z1);%前窗均值
Z21=mean(Z2);
Z12=std(Z1,1);%前窗标准差
Z22=std(Z2,1);
end
if Z12<T_std&&Z22<T_std;
        Z(i)=sum(Z1)+sum(Z2);
        T=P_fa.^(-1/R)-1 ;
         S1(i)=Z(i).*T; 
    else if Z12<T_std&&Z22>T_std;
            Z(i)=sum(Z1);
            T=P_fa.^(-1/(R/2))-1 ;
            S1(i)=Z(i).*T; 
        else if Z12>T_std&&Z22<T_std;
                Z(i)=sum(Z2);
                T=P_fa.^(-1/(R/2))-1 ;
                S1(i)=Z(i).*T; 
            else if Z12>T_std&&Z22>T_std;
                        [a_m1,b_m1]=max(Z1);
                        Z1(b_m1)=[];
                        [a_mi1,b_mi1]=min(Z1);
                        Z1(b_mi1)=[];
                    [a_m2,b_m2]=max(Z2);
                        Z2(b_m2)=[];
                        [a_mi2,b_mi2]=min(Z2);
                        Z2(b_mi2)=[];
                        Z(i)=sum(Z1)+sum(Z2);
                        T=P_fa.^(-1/(R-4))-1 ;
                        S1(i)=Z(i).*T; 
                end
            end
        end
        
 
 
end 
%Z(i)=Z(i)-P1((i-1)*L_move+R/2+1);
    
end 
%%%%%%%%%%%%%%

figure(3);plot(f1,20.*log10(P1.^2)); 
hold on;
plot(f1(n+3:256-n-3),20.*log10(S1(n+3:256-n-3)),'r');