%-----------函数ML_Estimation-------------clc;clearall;Nt=1;%仿真次数M=8;%MPSKNsym=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值fortime=1:Nt%仿真循环Nt次d=randint(1,Nsym,[0M-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;fori=1:length(hrcos)sum3=sum3+hrcos(i)^2;endh_k=hrcos/sqrt(sum3);m_k=conv(b_k,h_k);%序列成型L=length(m_k);%成型后的数据总长度K=Nss*Nsym;%s1=0;s2=0;fork=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;forl=1:K%公式中两个求和因子计算s1=s1+real(conj(r_k(l))*m_k(l));s2=s2+abs(r_k(l))^2;endrho_ml(k,time)=(Nss^2)*(s1^2)/(K^2)/(s2/(K-1.5)-(Nss*(s1^2)/K/(K-1.5)));%代入公式计算信噪比估计值%Pe=[Pemean((rho_ml-rho(k)).^2)/(rho(k).^2)];endendNMSE=2./rho/Nsym+1/Nss/Nsym;%克拉美罗界计算CRBfigure(3);semilogy(snr,NMSE,'-b*');holdon,gridon;fori=1:length(snr)sum1=0;forj=1:Ntsum1=sum1+(rho_ml(i,j)-rho(i)).^2;endpe(i)=sum1./(rho(i).^2)/Nt/K;%归一化均方误差NMSEendsemilogy(snr,pe,'-r*');gridon;