matlab信噪比估计程序

%----------- 函数ML_Estimation -------------

clc;

clear all;

Nt=1; %仿真次数

M=8; %MPSK

Nsym=64; %信源仿真符号数

Nss=16; %上采样点数

Es=1; %每符号能量

snr=1:30; %仿真信噪比范围

rho=10.^(snr/10); %实际信噪比大小

%sigma=sqrt(1/2)*(10.^(-snr/20)); %噪声根方差

No=Es./rho;

rho_ml=zeros(length(snr),Nt); %估计信噪比存储器,一列对应一次仿真，一行对应一个SNR值

for time=1:Nt %仿真循环Nt次

d=randint(1,Nsym,[0 M-1]); %产生Nsym个随机数作为信源符号

a_n=exp(j*(2*pi/M*d+pi/M)); %构成MPSK调制符号

figure(1);

plot_astrology(real(a_n),imag(a_n));

b_k=zeros(1,Nss*Nsym);

n=1:Nss:length(b_k);

b_k(n)=a_n;

%b_k=upsample(a_n,Nss); %上采样，每符号取Nss个采样点

%hrcos=firrcos(256,1,1,16,'sqrt'); %成型滤波器，采用RRC方式，阶数127，滚将系数0.5

%gk=conv(hrcos,hrcos);

hrcos=rcosflt(1,1,16,'sqrt',0.5,3);

sum3=0;

for i=1:length(hrcos)

sum3=sum3+hrcos(i)^2;

end

h_k=hrcos/sqrt(sum3);

m_k=conv(b_k,h_k); %序列成型

L=length(m_k); %成型后的数据总长度

K=Nss*Nsym;

%s1=0; s2=0;

for k=1:length(snr) %计算每个信噪比对应的信噪比估计值

%r_k=m_k+No(k)*(randn(1,L)+j*randn(1,L));

%r_k=m_k+nb(k)/sqrt(2)*(randn(1,L)+j*randn(1,L));

%r_k=sqrt(Es)*m_k+sqrt(No(k))/sqrt(2)*(randn(1,L)+j*randn(1,L)); %往信号中加入高斯白噪声

r_k=awgn(m_k,snr(k),'measured');

y_k=conv(r_k,h_k);

y_n=y_k(112+1:Nss:(length(y_k)-112));

figure(2);

plot_astrology(real(y_n),imag(y_n));

s1=0; s2=0;

for l=1:K %公式中两个求和因子计算

s1=s1+real(conj(r_k(l))*m_k(l));

s2=s2+abs(r_k(l))^2;

end

rho_ml(k,time)=(Nss^2)*(s1^2)/(K^2)/(s2/(K-1.5)-(Nss*(s1^2)/K/(K-1.5))); %代入公式计算信噪比估计值

% Pe=[Pe mean((rho_ml-rho(k)).^2)/(rho(k).^2)];

end

end

NMSE=2./rho/Nsym+1/Nss/Nsym; %克拉美罗界计算 CRB

figure(3);

semilogy(snr,NMSE,'-b*');

hold on,grid on;

for i=1:length(snr)

sum1=0;

for j=1:Nt

sum1=sum1+(rho_ml(i,j)-rho(i)).^2;

end

pe(i)=sum1./(rho(i).^2)/Nt/K; %归一化均方误差 NMSE

end

semilogy(snr,pe,'-r*');

grid on;