tic

clear
format short
c=3*10.^8;
L=1;
N=256;
M=80;

fm1=0.35;%信號1起始頻率
Km1=0.05;%信號2調頻斜率
fm2=0.25;%信號2起始頻率
Km2=0.12;%信號2調頻斜率
f0=fm1+Km1/2;
lamta=c/f0;
len=lamta/2;
doature=10;
P1=doature*pi/180;


t=1:N;

%%%%%%%%%%%%%%%%%%
%CS的方法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%構造原子庫
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
atom=zeros(51*201,N);
for l=1:51
    fa=(l-1)*0.01;%0~0.5
    for k=1:201
        Ka=k*0.001;    %-0.05~0.05
         atom((l-1)*201+k,:)=1.*exp(j.*2*pi*(Ka/2*t.^2/N+fa*t));
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%設置信噪比個數及每個信噪比試驗次數
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

time=1;
A=(1/sqrt(M))*randn(M,N);       %構造觀測矩陣

sig1=exp(j*2*pi*(Km1/2*(t.^2)/N+fm1*t));
sig2=exp(j*2*pi*(Km2/2*(t.^2)/N+fm2*t)); 

for ii=1:6;
  snr=(ii-1)*4-10;

 Amp=sqrt(2*10^(snr/10));


%生成信號

  

%控制循環次數 
 for times=1:time
     
 noise=randn(1,N)+j*randn(1,N);
z=Amp*sig1+Amp*sig2+noise; 

%求出觀測值y
y1=A*z';


%解優化問題

cvx_begin
 variable x(51*201)
 minimize(sum_square(real(y1)-A*real(atom')*x)+0.01*norm(x,1))
cvx_end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求解起始頻率和調頻斜率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[m,n]=max(abs(x));
e_fm1=floor(n/201)/100;
e_Km1=(n-floor(n/201)*201)/1000;
x(n)=0;
[m,n]=max(abs(x));
e_fm2=floor(n/201)/100;
e_Km2=(n-floor(n/201)*201)/1000;
if e_fm1>e_fm2
ef1(times)=e_fm1;
eK1(times)=e_Km1;
ef2(times)=e_fm2;
eK2(times)=e_Km2;
else
ef1(times)=e_fm2;
eK1(times)=e_Km2;
ef2(times)=e_fm1;
eK2(times)=e_Km1;
end
end
  
f1_mean(ii)=mean(ef1);
K1_mean(ii)=mean(eK1);
f2_mean(ii)=mean(ef2);
K2_mean(ii)=mean(eK2);

f1_resm(ii)=rmse(ef1,fm1);
K1_resm(ii)=rmse(eK1,Km1);
f2_resm(ii)=rmse(ef2,fm2);
K2_resm(ii)=rmse(eK2,Km2);

end







%%%%%%%%%%%%%%%%%顯示結果
 i=1:6;
 snr=(i-1)*4-10;

figure(1)
plot(snr,f1_resm,'r-o',snr,f2_resm,'b-*');
xlabel('SNR/dB','fontsize',12,'fontweight','b');
ylabel('RMSE','fontsize',12,'fontweight','b');
title('起始頻率');
legend('信號1-起始頻率','信號2-起始頻率')

figure(2)
plot(snr,K1_resm,'r-o',snr,K2_resm,'b-*');
xlabel('SNR/dB','fontsize',12,'fontweight','b');
ylabel('RMSE','fontsize',12,'fontweight','b');
title('調頻斜率');

legend('信號1-調頻斜率','信號2-調頻斜率')

toc