%單頻信號頻率估計主函數
%
%*********************************************************
clear all;
close all;
clc;
%*******************************************************
 fs=100e6;  %采樣頻率100MHz
 deltt=1/fs; %采樣間隔
 numMC=100; %Monte Carlo仿真次數
 %****************************************
 for ii=1:numMC;   %%開始蒙特卡洛仿真
 T=4e-6;     %脈寬設置1us
 SNR=0;     %設定信噪比
 
 %產生單頻信號+噪聲
 N=fix(T/deltt);  %信號點數，在本算法中要求N為偶數
 if mod(N,2)==1   %如果N不是偶數，則減1變成偶數
     N=N-1   
 end 
 
%這是N可能已經變化，相應改變信號的脈寬
T=N*deltt; %脈寬
%***********************************************
n=0:N-1;

randnoise=randn(1,N)+i*randn(1,N); %產生I,Q路的噪聲，功率各為1，總功率2
%fc=(2*rand(1)-1)*35e6;      %隨機產生頻率為-35~35MHz的信號
fc=10.567e6;
%phase=(2*rand(1)-1)*pi;     %隨機產生初始相位-pi~pi的信號
phase=pi/6;
signal_A=10^(SNR/20)*sqrt(2);  %噪聲功率為2，電壓幅度值sqrt(2),所以信號幅度*sqrt(2)
%signal_A=1;
signalN=signal_A*exp(i*(2*pi*fc*n*deltt+phase))+randnoise;  %信號+噪聲
%************************************************************************
%調用單頻信號的頻率與相位的估計子函數SFM_estimation
[f_e(ii) phase_e(ii)]=SFM_estimation(signalN,N,deltt);
%**************計算誤差
f_error(ii)=(f_e(ii)-fc)/fc;   %估計值-真值；
phase_error(ii)=(phase_e(ii)-phase)/phase; %估計值-真值；
%將phase_error限制在-pi~pi之間
if phase_error(ii)>pi
    phase_error(ii)=phase_error(ii)-2*pi;
end;
if phase_error(ii)<-pi
    phase_error(ii)=phase_error(ii)+2*pi;
end;
end;  %end for;
%**********************************************************
%對估計的f_e進行平均
f_e_mean=mean(f_e);
phase_e_mean=mean(phase_e);
%均方根誤差
f_e_rms=sqrt(mean((f_e-f_e_mean).^2));  %方均根
phase_e_rms=sqrt(mean((phase_e-phase_e_mean).^2));  %方均根
% %平均絕對誤差
% f_e_mae=mean(abs(f_e-f_e_mean));
% phase_e_mae=mean(phase_e-phase_e_mean);


%畫圖
figure(1)  %產生的信號圖
plot(n*deltt,real(signalN),'r',n*deltt,imag(signalN),'b');
xlabel('時間');
ylabel('幅值');
title('仿真產生的信號');
legend('實部I，虛部Q');
%
figure(2)   %信號頻譜
signalN_fft=abs(fft(signalN));
f=0:fs/N:(N-1)*fs/N;
plot(f,signalN_fft);
xlabel('頻率/Hz');
ylabel('幅值');
title('仿真信號的頻譜');
%
figure(3)   %頻率估計誤差
plot(f_error);
xlabel('次數');
ylabel('誤差值/Hz');
title('相位法頻率估計誤差');
%
figure(4)   %相位估計誤差
plot(phase_error);
xlabel('次數');
ylabel('誤差值/弧度');
title('相位法相位估計誤差');
%求均方根誤差