clc
clear all
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L = 64;
N = L/2;
n = 2;
T = 1/10^n;
fs = 1/T;
f_bias = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%添加AWGN%%%%%%%%%%%%%%%%%%%%%%%%%
SNR = -5:1:15;
times = 400; % 设定循环次数
sigma = zeros(1,length(SNR));
var_Fitz = zeros(1,length(SNR));
var_Kay = zeros(1,length(SNR));
var_LR = zeros(1,length(SNR));
f=[];
for t = 1:length(SNR)
sigma(t) = sqrt(1/(2*10^(SNR(t)/10)));
for k = 1:times
pilot_input = ones(1,2*L);
[pilot] = qpsk(pilot_input);
pilot_frequency = [pilot];
for n = 1:length(pilot_frequency)
s(n) = pilot_frequency(n)*exp(j*2*pi*f_bias*n/fs);%上变频
end %大概是因为没有成形的缘故，处理的信号是离散时间的
r = s+sigma(t)*(randn(1,length(s))+j*randn(1,length(s)));%叠加噪声后的接收信号
%r_conj = zeros(1,length(r));%一般赋值时先给出空间在赋值,也可有可无
Ip = 1:L;%前置，IP是什么？
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Fitz%%%%%%%%%%%%%%%%%%%%%%%%%
sum_arg_cor = 0;
for m = 1:N
sum_z = 0;
for n = m:L
sum_z = sum_z+r(Ip(n))*conj(r(Ip(n)-m+1));
end;
Correlation(m) = sum_z/(L-m);
sum_arg_cor = sum_arg_cor+angle(Correlation(m));
end
f_est = sum_arg_cor/(pi*N*(N+1)*T);
f_est = sum_arg_cor/(pi*N*(N-1)*T);
var_Fitz(t) = var_Fitz(t)+((f_est-f_bias)*T)^2/times;
%——————————————————Kay算法—————————————————
f_est_Kay = 0;
for k = 1:L-1
R = 3*L/2/(L^2-1)*(1-((2*k-L)/L)^2);
f_est_Kay = f_est_Kay+R*angle(r(k+1)*conj(r(k)))/2/pi/T;
end
var_Kay(t) = var_Kay(t)+((f_est_Kay-f_bias)*T)^2/times;
%——————————————————L&R算法——————————————————
sum_cor = 0;
for m = 1:N
sum_z = 0;
for n = m:L
sum_z = sum_z+r(Ip(n))*conj(r(Ip(n)-m+1));
end;
Correlation(m) = sum_z/(L-m);
sum_cor = sum_cor+Correlation(m);
end
sum_arg =angle(sum_cor);
f_estimation = sum_arg/(pi*N*T);
var_LR(t) = var_LR(t)+((f_estimation-f_bias)*T)^2/times;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%MCRB%%%%%%%%%%%%%%%%%%%%%%%
MCRB=zeros(1,length(SNR)) ;
for k=1:length(SNR)
EsN0(k)=10^(SNR(k)/10);%snr和ESN之间的转化。。
MCRB(k)=3/(2*pi^2*L^3*EsN0(k));
end
figure(1)
semilogy(SNR,MCRB,'--','LineWidth',3);hold on
semilogy(SNR,var_Fitz,'-square','LineWidth',1.5);hold on
semilogy(SNR,var_Kay,'->r','LineWidth',1.5);hold on
semilogy(SNR,var_LR,'-dg','LineWidth',1.5);hold on
xlabel('Es/N0, dB');
ylabel('归一化载波频率估计方差');%最小均方差准则
title('三种时域频率估计算法性能比较');
legend('CRB','Fitz算法','Kay算法','L&R算法');
