统计信号处理 实验二[学习分享]

1、统计信号处理实验二目的：1 掌握参数估计方法；2 掌握用计算机分析数据的方法。内容：假设一个某信号为两个正弦波的和：其中、为未知参数，、分别为两个正弦波的频率，已知，。在对这个信号进行观测的时候，接收到的信号为： 其中为观测到的信号，为白色观测噪声。假设我们总共接收了1秒的信号，然后用的采样间隔测量到了200个观测值,，,。 1） 分别用最大似然，最小二乘方法，推导出根据观测结果估计参数、的公式，并编程实现算法；假设A1=0.5，A2=3，得到1000组观测值作为样本(X1,.XN)最大似然估计： 令，得到则可写作矩阵乘法即 利用matlab可以求得已知样本(X1,.XN)时A1，A2的最大似

2、然估计。最小二乘估计：假设观测结果和待估计参量间为线性关系：x=H.*p+nx=X1XNmNT，p=A1A1A2A22m，H=sin(4t(1)sin(6t(1)sin(4t(N)sin(6t(N)N2，n=n1nNmNT其中N=200，m=1000最小二乘解plms=(HTH)-1HTx加权最小二乘解plms=(HTWH)-1HTWx运行程序，得到结果Etheta=(0.4998, 3.0235)Ep=(0.4998, 3.0235)在A1=0.5，A2=3，得到1000组观测值作为样本(X1,.XN)的前提下，求得最大似然估计与最小二乘估计平均结果完全相同，且都与假设值接近，可以认为是无偏

3、估计。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个点更加集中于实际值附近E

4、theta=(0.4956, 2.9926)Ep=(0.4956, 2.9926)m=10000时得到的结果Etheta=(0.4997, 3.0027)Ep=( 0.4997, 3.0027)可见m增大到无穷大时，估计值会与实际值相同，故具有一致性。4） 利用估计出的参数，得到的估计，并用Monte_Carlo法计算在等各个时间点上对目标位置估计的方差和偏差；平均偏差0.009，方差0.08，估计效果理想。5） 将噪声的分布改为在（-1，+1）区间分布，应用上面推导出的最大似然，最小二乘公式对参数进行估计，并计算估计的偏差和方差。对A1、A2估计：平均偏差标准偏差方差ML估计A16.8823

5、e-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，用傅里叶级数表示三角波：s1=0.5A1+42A1k=2n-1K1k2cos4kt, (n=1,2,)s2=0.5A2+42A

6、2k=2n-1K1k2cos6kt, (n=1,2,)令，得到：i=1NA10.5+4/2k=2n-1K1k2cos4kt2+A20.5+4/2k=2n-1K1k2cos6kt*0.5+4/2k=2n-1K1k2cos4kt=i=1Nxi0.5+4/2k=2n-1K1k2cos4kti=1NA20.5+4/2k=2n-1K1k2cos6kt2+A10.5+4/2k=2n-1K1k2cos6kt*0.5+4/2k=2n-1K1k2cos4kt=i=1Nxi0.5+4/2k=2n-1K1k2cos6kt则可写作矩阵乘法i=1N0.5+4/2k=2n-1K1k2cos4kt2i=1N0.5+4/2k

7、=2n-1K1k2cos6kt*0.5+4/2k=2n-1K1k2cos4kti=1N0.5+4/2k=2n-1K1k2cos6kt*0.5+4/2k=2n-1K1k2cos4kti=1N0.5+4/2k=2n-1K1k2cos6kt2A1A2=i=1Nxi0.5+4/2k=2n-1K1k2cos4kti=1Nxi0.5+4/2k=2n-1K1k2cos6kt即 利用matlab可以求得已知样本(X1,.XN)时A1，A2的最大似然估计。假设观测结果和待估计参量间为线性关系：x=H.*p+nx=X1XNmNT，p=A1A1A2A22m，H=0.5+4/2k=2n-1K1k2cos4kt(1)0

8、.5+4/2k=2n-1K1k2cos6kt(1)0.5+4/2k=2n-1K1k2cos4kt(N)0.5+4/2k=2n-1K1k2cos6kt(N)N2n=n1nNmNT其中N=200，m=1000最小二乘解plms=(HTH)-1HTx加权最小二乘解plms=(HTWH)-1HTWx代码：%（1）clear all;clc;A1=0.5;A2=3;N=200;m=1000;s=zeros(m,N);n=zeros(m,N);theta=zeros(2,m);for j=1:m i=1:N; t=i/N; s(j,:)=A1*sin(4*pi*t)+A2*sin(6*pi*t); n(j

9、,:)=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=A B;B C; b=U;V; theta(:,j)=Mb; %最大似然估计endsubplot(2,2,1);plot(theta(1,:),.);title(ML估计A1值);grid on;subplot(2,2,2);plot(theta(2,:),.);tit

10、le(ML估计A2值);grid on;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值);grid on;subplot(2,2,4);plot(p(2,:),.);title(LS估计A2值);grid on;Ep=mean(p,2);%（2）for j=1:m EBmla1=sum(abs(theta(1,j)-Etheta(1)/m; %平均偏差 NB

11、mla1=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=NBmla12*(m-1)/m; %方差 Dmla

12、2=NBmla22*(m-1)/m; Dlsa1=NBlsa12*(m-1)/m; Dlsa2=NBlsa22*(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);for j=1:m i=1:N; t=i/N; s_hat1(j,:)=a1*sin(4*pi*t)+a2*sin(

13、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;for i=1:N for j=1:m EBs1(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:UsersLENOVODesktopEBs

14、,EBs);% xlswrite(C:UsersLENOVODesktopDBs,DBs);subplot(2,2,1);plot(EBs1);title(ML估计平均偏差);grid on;subplot(2,2,2);plot(DBs1);title(ML估计方差);grid on; subplot(2,2,3);plot(EBs2);title(LS估计平均偏差);grid on;subplot(2,2,4);plot(DBs2);title(LS估计方差);grid on;%（5）clear all;clc;A1=0.5;A2=3;N=200;m=1000;s=zeros(m,N);n

15、=unifrnd(-1,1,m,N); %产生由(-1,1)上均匀分布的随机数组成的m*N矩阵theta=zeros(2,m);for j=1:m i=1:N; t=i/N; s(j,:)=A1*sin(4*pi*t)+A2*sin(6*pi*t); 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=A B;B C; b=U;V; theta(:,j)=M

16、b; %最大似然估计endsubplot(2,2,1);plot(theta(1,:),.);title(ML估计A1值);grid on;subplot(2,2,2);plot(theta(2,:),.);title(ML估计A2值);grid on;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;figure(1);subplot(2,2,3);plot(p(1,:),.);title(LS估计A1值);grid on;subplot(2,

17、2,4);plot(p(2,:),.);title(LS估计A2值);grid on;Ep=mean(p,2);for j=1:m EBmla1=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)-

18、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=NBmla12*(m-1)/m; %方差 Dmla2=NBmla22*(m-1)/m; Dlsa1=NBlsa12*(m-1)/m; Dlsa2=NBlsa22*(m-1)/m;enda1=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,

19、N);DBs1=zeros(1,N);EBs2=zeros(1,N);DBs2=zeros(1,N);for j=1:m i=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;for i=1:N for j=1:m EBs1(i)=sum(abs(x_hat1(j,i)-Ex(i)/m; %平均偏差 DBs1(i)=(sum(x_hat1(j,i)-Ex(i)2/m; %方差 EBs

20、2(i)=sum(abs(x_hat1(j,i)-Ex(i)/m; DBs2(i)=(sum(x_hat1(j,i)-Ex(i)2/m; endendfigure(2);subplot(2,2,1);plot(EBs1);title(ML估计平均偏差);grid on;subplot(2,2,2);plot(DBs1);title(ML估计方差);grid on; subplot(2,2,3);plot(EBs2);title(LS估计平均偏差);grid on;subplot(2,2,4);plot(DBs2);title(LS估计方差);grid on;%（6）clear all;clc;

21、A1=0.5;A2=3;N=200;m=1000;K=99;s=zeros(m,N);s1=zeros(1,N);s2=zeros(1,N);n=zeros(m,N);theta=zeros(2,m);H=zeros(N,2);for j=1:m for i=1:N; t=i/N; k=1:2:K; s1(i)=A1/2+4*A1/(pi*pi)*sum(k.*k).-1.*cos(4*pi*k*t); %傅里叶级数表示三角波 s2(i)=A2/2+4*A2/(pi*pi)*sum(k.*k).-1.*cos(6*pi*k*t); s(j,:)=s1(i)+s2(i); n(j,:)=3.*r

22、andn(1,N); x=s+n; A=sum(0.5+4/(pi*pi)*sum(k.*k).-1.*cos(4*pi*k*t).2); %A、B、C为M矩阵元素 B=sum(0.5+4/(pi*pi)*sum(k.*k).-1.*cos(4*pi*k*t)*(0.5+4/(pi*pi)*sum(k.*k).-1.*cos(6*pi*k*t); C=sum(0.5+4/(pi*pi)*sum(k.*k).-1.*cos(6*pi*k*t).2); U=sum(x(j,:).*(0.5+4/(pi*pi)*sum(k.*k).-1.*cos(4*pi*k*t); %矩阵行求和 V=sum(x(

23、j,:).*(0.5+4/(pi*pi)*sum(k.*k).-1.*cos(6*pi*k*t); H=(0.5+4/(pi*pi)*sum(k.*k).-1.*cos(4*pi*k*t);(0.5+4/(pi*pi)*sum(k.*k).-1.*cos(6*pi*k*t); %H为已知的参数矩阵 end M=A B;B C; b=U;V; theta(:,j)=Mb; %最大似然估计endsubplot(2,2,1);plot(theta(1,:),.);title(ML估计A1值);grid on;subplot(2,2,2);plot(theta(2,:),.);title(ML估计A2

24、值);grid on;Etheta=mean(theta,2);p=zeros(2,m); %最小二乘估计X=x;p=inv(H*H)*H*X;subplot(2,2,3);plot(p(1,:),.);title(LS估计A1值);grid on;subplot(2,2,4);plot(p(2,:),.);title(LS估计A2值);grid on;Ep=mean(p,2);for j=1:m EBmla1=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