多种功率谱估计郑州大学随机信号处理大作业

随机信号处理大作业多种功率谱估计的算法实现及性能比较一、引言二、原理及过程Xx(ejw)即ep(n)=ep-1(n-1)+k,eg-1(n)=1,2,…,p(3.4)①由初始条件ed(n)=eg(n)=x(n)和式(3.4)求出k₁③由k₁和式(3.3)求出ef(n),e{(n),再由式(3.4)估计k₂即e=[1eiw.…ei(N-1)w]T若令a;=1(j=M+1,M+2,…,N),则计算功率谱人为的设置信噪比范围为-15db~15db,通过分别调用上面的求功率谱的函数，求特定信噪比下1000次功率谱估计的均方误差，分别绘制采用上述功率谱估计时，均方误差①求古典法在-15db下的1000次功率谱估计的均方误差②同上求-14db~15db下的1000次功率谱估计的均方误差③绘制均方误差与信噪比的关系④同理，求出采用上文中现代谱估计时，均方误差与信噪比的关系结果分析原始信号为信噪比和均方误差的关系，并分析各种功率谱估计方法的性能优劣。1,古典法N不小于130时可以识别出两个波峰n*nnn=128f1=0.20002,直接解yole-walker方程n=130f1=0.2000f2=2.2126P=10时，无论N取多少，均不能识别两个波峰P=40N=128f1=0.1998fP=20N=256f1=0.1995p=40n=256f1=0.1995f2=0.2135P=20N=128f1=0.1992P=40N=128f1=0.1998f2=0.2140P=404Burg算法N=256f1=0.1995f2=0.2135P=10时，N取任何值都不能识别出两个波峰P=10N=128f1=0.2031P=10N=256f1=0.2027P=20时，N不小于55时可识别出两个波峰P=20N=128f1=0.2003f2=0.2138P=20N=256f1=0.2005f2=0.2125P=40时，N大于50时可识别出两个波峰P=40N=128f1=0.2005f2=0.2138P=40N=256f1=0.2008f2=0.2130在相同的条件下，f1=0.2,f2=0.213,N=128时，观察图形发现，用古典法只能估计一个频率点，用AR模型估计时，当阶次P达到一定时，仍能分辨出两频率点，所以只要选择合适的阶次，就能用AR模型估计两接近的频率点，而N越大，图像越精细，越易于分辨两频率点，对于古典法，当N合适时，仍能分辨相近的频率点。由图易知，分辨率的关系为：古典法<直接解Yole-walker方程法<快速递推法<Burg算法。基于AR模型的现代功率谱估计质量明显优于经典谱估计，基于AR模型的现代谱估计是逐步改进，性能递增，精度提高的，直接解Yole-Walker方程运算量大，且不一定有解，采用Levinson递推可以大大减小运算量，但由于自相关函数的估计默认取样点以外的值为零，引入了一定的误差，而Burg算法引入前后向预测误差，来估计功率谱，避免了自相关函数带来的误差，因而精度较高，因此，Burg算法求得的AR模型最稳定，而且Burg算法不需要自相关函数，所以性能优于自相关法，由图可知，在P=20,N=55时即可分辨出两个峰值，所以在短数据时，Burg算法优势明显，具有较高但是Burg算法进行谱估计会出现谱线分裂、谱峰偏移等问题，而且对于低信噪比情况，用AR模型很难准确估计出淹没在噪声中的正弦波频率，而改进的MUSIC算法采用特征分解技术，在估计正弦信号和噪声叠加的信号中明显优于AR模型功率谱估计。6信噪比与均方误差的关系RR起初在刚接触这门课时，甚至有些反感老师的严格要求，大概是习惯了其他老师的“温柔乡”吧，现在回过头，在开始写这篇报告时，对老师有一种说不出的感觉，我想那就是感谢的羞涩表达吧。虽然老师的表达有时不靠谱，但老师比那些看似靠谱的老师要强很多，他牢记灌输知识的使命，又深谙人心之道，不时来点心灵鸡汤，努力的鞭挞我们前进，将知识灌输给我们，这让我不由得想起了：“小时候，不爱吃饭的我被妈妈训斥着吃饭的事”,想到这，对老师有一种莫名的亲切，简直深得我心。天啊，我也不知道要写这些，也许，这tm是发自肺腑，不由而衷吧。首先，让我对Matlab和Matlab的语法有了更深的了解，提高了我的编程能力，虽然以前开设过Matlab的实验课，却从没有真正入门Matlab,可以说，陈老师是我Matlab的启蒙老师啊!然后，让我真正理解了一部分信号处理信号分析的知识，加深了对概念的理解，以前学信号分析信号处理时，只是记忆公式和概念，很少去分析为什么,甚至也不去管为什么,陈老师是一位真正有学识的人，他不照本宣科，总能说出自己的理解，而这些很多是其他老师未曾提到的(或许是我没听到其他老师提),这对于融会贯通整个学科至关重要。最后，让我对一些问题有了更深层次的看法。程序附录1经典法(周期图法)N=130;w=n/N;wn=randn(1,N)xn=10*sin(2*pi*0.2*n+(pi/3))+5*sin(2*pi*0.213*n+(pi/4))+wnXk=fft(xn);'两个正弦信号与白噪声叠加的时域波形);f=(locs-1)/(2*N)2直接解Yole-Walker方程N=256;p1=40;x=10*sin(2*pi*f1*n+pi/3)+5*sin(2*pi*f2*n+pi/4)+randn(1,N);sum=sum+x(n)*x(n);sum=0;sum=x(n)*x(n+m)+sum;sum=0;%求自相关函数form=1:p1a(i,k)=R(i-k);end%a为p行p列的矩阵，对角线及以上部分为0B(i,i)=RO,G=(RO+b)^(1/2);f=(locs-1)/(2*2000)3Levinson-Durbin快速递推法N=256;n=1:N;p1=40;x=10*sin(2*pi*f1*n+pi/3)+5*sin(2*pi*f2*n+pi/4)+randn(1,N);'两个正弦信号与白噪声叠加的时域波形);%产生两个正弦信号与白噪声的sum=sum+x(n)*x(n);sum=0;sum=x(n)*x(n+m)+sum;sum=0;%求自相关函数form=2:p1forj=1:m-1sum=sum+aa(m-1,i)*R(m-i);aa(m,m)=-(R(m)+sum)/p(m-1);%定义levinson-durbin递推式中的km%由H(z)用freqz()求解功率谱G=p(p1)^(1/2);A=[1,aa(p1;:)];[H,w]=freqz(B,A,2000);%将pi分成2000份f=(locs-1)/(2*2000)4Brug算法xn=10*sin(2*pi*n*f1+pi/3)+5*sin(2*pi*n*f2+pi/4)wn=randn(1,N);%产生高斯白噪声中数组索引从1开始forj=1:Nsum=sum+xn(i)*xn(i);r(1)=sum/N;%因matlab中数组索引从1开始suma=0;sumb=0;suma=suma+ef(m-1,n)*eb(m-1,n-1);%初步构造KM中的分子sumb=sumb+ef(m-1,n)^2+eb(m-1,n-1)^2;%构造KM中的分母k(m-1)=-2*suma/sumb;%完全构造KM,根据前后向预测均方误差之和最小来求反射系suma=0;sumb=0;%每次循环前赋零值forj=1:m-1aa(m,m)=k(m);p(m)=p(m-1)*(1-k(m)G=p(P)^(1/2);A=[1,aa(P;)];title('Burg算法功率谱)f=(locs-1)/(2*2000)x=10*sin(2*pi*0.2*n+pi/3)+5*sum=sum+x(n)*x(n);end%求自相关函数Rx(m,n)=R(m-n);%将值赋值给矩阵Rx(m,m)=R0;%将主对角线的值变为R(0)列是对应的右特征向量，使得A*V=V*Dp=0;%特侦知的个数pev=e*V(:1:N-p);title('MUSIC算法求功率谱)xlabel('频率)axis([00.5010])%e矩阵%功率谱信噪比与均方误差的关系将求功率谱的各种算法改写成函数以便直接调用，由与采用不同功率谱估计方法时，仅需替换代码中的一部分，关系图谱已附，以下仅提供采用Burg时的代码和子程序代码：MSE1=0;MSE2=0;V=26;s=[1:26];c=-15;10;M1=zeros(1,v);%定义1行V列的全零矩阵forc=1:v%求30次信噪比和均方误差fori=1:500%循环500次求均方误差x=awgn(x1,c,measured');%把信噪比为c的高斯白噪声加到信号x1中pw1=Q(x);%我将burg算法改写成了Q函数，此处直接调用。对应的位置f1=(locs(1,1)-1)/50f2=(locs(1,2)-1)/500;M1(1,s)=MSE1;M2(1,s)=MSE2S(1,s)=c;xlabel('SNR/dB');xlabel('SNR/dB');Q函数(burg子程序)的代码：N=512;%定义采样点数P=40;%滤波器阶数的最大取值x1=10*sin(2*pi*n*f1+pi/3)+5*sin(2*pi*n*f2+pi/4);x=awgn(x1,10,'measured');%把信噪比为c的高斯白噪声加到信号x1中ef(1,i)=x(i);中数组索引从1开始fori=1:Nr(1)=sum/N;%因matlab中数组索引从1开始%Burg递推suma=0;sumb=0form=2:(P+1);%循环的阶次，因matlab中数组索引从1开始suma=suma+ef(m-1,n)*eb(m-1,n-1);%初步构造KM中的分子sumb=sumb+ef(m-1,n)^2+eb(m-1,n-1)^2;%构造KM中的分母k(m-1)=-2*suma/sumb;%完全构造KM,根据前后向预测均方误差之和最小来求反射系suma=0;sumb=0;%每次循环前赋零值forn