统计信号处理-实验二

《统计信号处理》实验二目的：1．掌握参数估计方法；2．掌握用计算机分析数据的方法。内容：假设一个某信号为两个正弦波的和：1122()sin(2)sin(2)stAftAft其中1A、2A为未知参数，1f、2f分别为两个正弦波的频率，已知12fHz，23fHz。在对这个信号进行观测的时候，接收到的信号为：()()()xtstnt其中()xt为观测到的信号，()(0,9)ntN，为白色观测噪声。假设我们总共接收了1秒的信号，然后用0.005ts的采样间隔测量到了200个观测值(0)x,()xt，…,(199)xt。1）分别用最大似然，最小二乘方法，推导出根据观测结果估计参数1A、2A的公式，并编程实现算法；假设A1=0.5，A2=3，得到1000组观测值作为样本(X1,...XN)最大似然估计：NiNiNiiiNNxxxxxP12112212121)21ln()|,...,,(ln))(t6sin())(t4sin(21iAiA令0ln1AP，0ln2AP得到0))}(t4sin(]))(t6sin())(t4sin())(t4({[sinN1212ixAiiAiii0))}(t6sin(]))(t6sin())(t4sin())(t6({[sinN1122ixAiiAiii则可写作矩阵乘法NiiNiiiiiiixixAAitititititit1121N12N1N1N12))(t6sin())(t4sin())(6(sin))(6sin())(4sin())(6sin())(4sin())(4(sin即bM利用matlab可以求得已知样本(X1,...XN)时A1，A2的最大似然估计。最小二乘估计：假设观测结果和待估计参量间为线性关系：x=H.*p+n𝐱=[𝑿1…𝑿𝑁]𝑚×𝑁𝑇，𝐩=[𝐴1…𝐴1𝐴2…𝐴2]2×𝑚，𝐇=[sin⁡(4πt(1))𝑠𝑖𝑛⁡(6𝜋𝑡(1))……𝑠𝑖𝑛⁡(4𝜋𝑡(𝑁))𝑠𝑖𝑛⁡(6𝜋𝑡(𝑁))]𝑁×2，𝐧=[𝒏1…𝒏𝑁]𝑚×𝑁𝑇其中N=200，m=1000最小二乘解𝒑̂lms=(𝐇𝑇𝐇)−1𝐇𝑇𝐱加权最小二乘解𝒑̂lms=(𝐇𝑇𝐖𝐇)−1𝐇𝑇𝐖𝐱运行程序，得到结果Etheta=(0.4998,3.0235)Ep=(0.4998,3.0235)在A1=0.5，A2=3，得到1000组观测值作为样本(X1,...XN)的前提下，求得最大似然估计与最小二乘估计平均结果完全相同，且都与假设值接近，可以认为是无偏估计。2）用Monte_Carlo法，计算出上面两种方法求出的参数的偏差和方差；平均偏差标准偏差方差ML估计A15.3211e-040.01682.8314e-04A22.4522e-040.00786.0134e-05LS估计A15.3211e-040.01682.8314e-04A22.4522e-040.00786.0134e-05两种方法方差相同，因此都是有效的。3）自由发挥：你认为还能通过Monte-Carlo仿真，对算法的那些特性进行研究？可以研究算法估计是否具有无偏性、有效性、一致性。无偏性、有效性已讨论。当m继续增大m=5000时得到的结果5000个点更加集中于实际值附近Etheta=(0.4956,2.9926)Ep=(0.4956,2.9926)m=10000时得到的结果Etheta=(0.4997,3.0027)Ep=(0.4997,3.0027)可见m增大到无穷大时，估计值会与实际值相同，故具有一致性。4）利用估计出的参数，得到()st的估计()st，并用Monte_Carlo法计算在0,,2,...,199tttt等各个时间点上对目标位置估计的方差和偏差；平均偏差0.009，方差0.08，估计效果理想。5）将噪声的分布改为在（-1，+1）区间分布，应用上面推导出的最大似然，最小二乘公式对参数进行估计，并计算估计的偏差和方差。对A1、A2估计：平均偏差标准偏差方差ML估计A16.8823e-060.01684.7366e-08A21.3544e-050.00781.8345e-07LS估计A16.8823e-060.01684.7366e-08A21.3544e-050.00781.8345e-07对目标位置x估计：最大似然估计与最小二乘估计结果仍然非常相近。6）假设s(t)是两个幅度未知、频率与上面正弦波一样的三角波的和，试推出这时的最大似然，最小二乘方法估计公式，并验证其性能。f1=2Hz，f2=3Hz，A1，A2未知假设A1=0.5，A2=3，用傅里叶级数表示三角波：𝑠1=0.5𝐴1+4𝜋2𝐴1∑1𝑘2𝐾𝑘=2𝑛−1cos(4𝜋kt),(n=1,2,…)𝑠2=0.5𝐴2+4𝜋2𝐴2∑1𝑘2𝐾𝑘=2𝑛−1cos(6𝜋kt),(n=1,2,…)NiNiNiiiNNxxxxxP12112212121)21ln()|,...,,(ln21ss令0ln1AP，0ln2AP得到：∑𝐴1[0.5+4𝜋2∑1𝑘2𝐾𝑘=2𝑛−1cos(4𝜋kt)]2+𝐴2[0.5+4𝜋2∑1𝑘2𝐾𝑘=2𝑛−1cos(6𝜋kt)][0.5𝑁=1+4𝜋2∑1𝑘2𝐾𝑘=2𝑛−1cos(4𝜋kt)]=∑[0.5+4𝜋2∑1𝑘2𝐾𝑘=2𝑛−1cos(4𝜋kt)]𝑁=1∑𝐴2[0.5+4𝜋2∑1𝑘2𝐾𝑘=2𝑛−1cos(6𝜋kt)]2+𝐴1[0.5+4𝜋2∑1𝑘2𝐾𝑘=2𝑛−1cos(6𝜋kt)][0.5𝑁=1+4𝜋2∑1𝑘2𝐾𝑘=2𝑛−1cos(4𝜋kt)]=∑[0.5+4𝜋2∑1𝑘2𝐾𝑘=2𝑛−1cos(6𝜋kt)]𝑁=1则可写作矩阵乘法[∑[.+∑=𝒏−()]=∑[.+∑=𝒏−𝐜𝐨𝐬(𝟔𝐤𝐭)][.+∑=𝒏−𝐜𝐨𝐬(𝐤𝐭)]=∑[.+∑=𝒏−𝐜𝐨𝐬(𝟔𝐤𝐭)][.+∑=𝒏−𝐜𝐨𝐬(𝐤𝐭)]=∑[.+∑=𝒏−(𝟔)]=][]=[∑[.+∑=𝒏−𝐜𝐨𝐬(𝐤𝐭)]=∑[.+∑=𝒏−𝐜𝐨𝐬(𝟔𝐤𝐭)]=]即bM利用matlab可以求得已知样本(X1,...XN)时A1，A2的最大似然估计。假设观测结果和待估计参量间为线性关系：x=H.*p+n𝐱=[𝑿1…𝑿𝑁]𝑚×𝑁𝑇，𝐩=[𝐴1…𝐴1𝐴2…𝐴2]2×𝑚，𝐇=[0.5+4𝜋2∑1𝑘𝐾𝑘=2𝑛−1𝑠(4𝜋𝑘𝑡(1))0.5+4𝜋2∑1𝑘𝐾𝑘=2𝑛−1𝑠(6𝜋𝑘𝑡(1))……0.5+4𝜋2∑1𝑘𝐾𝑘=2𝑛−1𝑠(4𝜋𝑘𝑡())0.5+4𝜋2∑1𝑘𝐾𝑘=2𝑛−1𝑠(6𝜋𝑘𝑡())]𝑁×2𝐧=[𝒏1…𝒏𝑁]𝑚×𝑁𝑇其中N=200，m=1000最小二乘解𝒑̂lms=(𝐇𝑇𝐇)−1𝐇𝑇𝐱加权最小二乘解𝒑̂lms=(𝐇𝑇𝐖𝐇)−1𝐇𝑇𝐖𝐱代码：%（1）clearall;clc;A1=0.5;A2=3;N=200;m=1000;s=zeros(m,N);n=zeros(m,N);theta=zeros(2,m);forj=1:mi=1:N;t=i/N;s(j,:)=A1*sin(4*pi*t)+A2*sin(6*pi*t);n(j,:)=3.*randn(1,N);x=s+n;A=sum((sin(4*pi*t)).^2);B=sum(sin(4*pi*t).*sin(6*pi*t));C=sum((sin(6*pi*t)).^2);U=sum(x(j,:).*sin(4*pi*t));%矩阵行求和V=sum(x(j,:).*sin(6*pi*t));M=[AB;BC];b=[U;V];theta(:,j)=M\b;%最大似然估计endsubplot(2,2,1);plot(theta(1,:),'.');title('ML估计A1值');gridon;subplot(2,2,2);plot(theta(2,:),'.');title('ML估计A2值');gridon;Etheta=mean(theta,2);p=zeros(2,m);%最小二乘估计H=zeros(N,2);H=[sin(4*pi*t);sin(6*pi*t)]';X=x';p=inv(H'*H)*H'*X;subplot(2,2,3);plot(p(1,:),'.');title('LS估计A1值');gridon;subplot(2,2,4);plot(p(2,:),'.');title('LS估计A2值');gridon;Ep=mean(p,2);%（2）forj=1:mEBmla1=sum(abs(theta(1,j)-Etheta(1)))/m;%平均偏差NBmla1=sqrt((sum(theta(1,j)-Etheta(1)))^2/(m-1));%标准偏差EBmla2=sum(abs(theta(2,j)-Etheta(2)))/m;NBmla2=sqrt((sum(theta(2,j)-Etheta(2)))^2/(m-1));EBlsa1=sum(abs(p(1,j)-Ep(1)))/m;NBlsa1=sqrt((sum(p(1,j)-Ep(1)))^2/(m-1));EBlsa2=sum(abs(p(2,j)-Ep(2)))/m;NBlsa2=sqrt((sum(p(2,j)-Ep(2)))^2/(m-1));Dmla1=NBmla1^2*(m-1)/m;%方差Dmla2=NBmla2^2*(m-1)/m;Dlsa1=NBlsa1^2*(m-1)/m;Dlsa2=NBlsa2^2*(m-1)/m;end%（4）a1=Etheta(1);a2=Etheta(2);a3=Ep(1);a4=Ep(2);s_hat1=zeros(m,N);s_hat2=zeros(m,N);Ex=zeros(N,1);Ex=mean(x,1);EBs1=zeros(1,N);DBs1=zeros(1,N);EBs2=zeros(1,N);DBs2=zeros(1,N);forj=1:mi=1:N;t=i/N;s_hat1(j,:)=a1*sin(4*pi*t)+a2*sin(6*pi*t);s_hat2(j,:)=a1*sin(4*pi*t)+a2*sin(6*pi*t);endx_hat1=s_hat1+n;x_hat2=s_hat1+n;fori=1:Nforj=1:mEBs1(i)=sum(abs(x_hat1(j,i)-Ex(i)))/m;%平均偏差DBs1(i)=(sum(x_hat1(j,i)-Ex(i)))^2/m;%方差EBs2(i)=sum(abs(x_hat1(j,i)-Ex(i)))/m;DBs2(i)=(sum(x_hat1(j,i)-Ex(i)))^2/m;endend%xlswrite('C:\Users\LENOVO\Desktop\EBs',EBs');%xlswrite('C:\Users\LENOVO\Desktop\DBs',DBs');subplot(2,2,1);plot(EBs1);title('ML估计平均偏差');gridon;subplot(2,2,2);plot(DBs1);title('ML估计方差');gridon;subplot(2,2,3);plot(EBs2);title('LS估计平均偏