统计信号处理 实验二

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=