白噪声中正弦组合的功率谱分析-郑大原创论文

通信工程孙大江郑州大学白噪声中正弦组合的功率谱分析专业：通信工程(1)学生姓名：孙大江学号:20092420128指导教师：陈恩庆完成时间：2020年4月7日通信工程孙大江摘要本文分别选取古典谱估计,现代参数方程谱估计,特征分解法谱估计中的各一种对淹没在白噪声中的正弦信号进行功率谱分析,并进行了相应的matlab仿真,仿真程序和结果见附录.关键词:功率谱,自相关,估计AbstractThispaperselectedtheclassicalspectrumestimation,parameterequationofmodernspectralestimation,featuredecompositionspectralestimationtoanalysisthesinesignalconcealedinthewhitenoise,andconductedtheimplementationofthecorrespondingmatlabsimulation,theSimulationprogramandresultsareSaveontheappendix.Keywords:powerspectrum,autocorrelation,estimation一.课题简介由于白噪声中的正弦组合是最常见的随机过程,而估计淹没在噪声中的正弦波频率即是信号处理中最有实际应用价值的技术之一,也是测试所有谱估计性能的基础.所以在本课题中以信号()20sin(20.2)2sin(20.3)()xnnnn为例进行研究,并分别以经典谱估计中的自相关法(间接法),和现代谱估计AR模型的Yule-Walker法(自相关法),以及分离信号子空间和噪声子空间的特征分解法(MUSIC)分别进行功率谱估计.最后再根据不同方法得到的不同结果进行分析.二.信号分析的演变在进行课题研究之前,我们首先明确一下为什么要进行信号谱分析.最初的信号研究知识局限在时域,但是随着科技的进步,时域分析已经完全满足不了科研需要,于是人们开始把目光转向了频域.最初,根据任何信号f(t)都可以在特定区间内展开成由完备正交集构成的无穷级数的理论,我们可以吧一个周期信号分解为有不同频率谐波组成的形式,并且分别在三角集和指数集中得到三角型傅里叶级数和指数型傅里叶级数njjntnn1f(t)Ae2,由于指数型运算的简便性,人们选择了指数型运算进行了深通信工程孙大江入发展,即以njnn1Ae2为振幅,以n为相角分别得到周期信号的频域幅度谱和相位谱对于非周期信号,由于相当于周期T区域无穷,所以基波角频率区域无穷小,个分量的幅度也相应的趋于无穷小,不能在进行幅度谱分析了,故引出新的概念频谱密度即单位频率的振幅,也就是我们通常所说的频谱,这里引出了傅里叶变换,即jtnTF(j)limF.Tf(t)edt根据傅里叶变换,我们可以进行频谱分析了下一步,引入能量谱的概念,信号能量22Ef(t)dt=F(j)df,其中2F(j)为单位频率的信号能量,也称为能量密度函数,即我们所说的能量谱,第四步,与能量谱类似引入功率谱概念T22TT21Plimf(t)dtT,由于随机信号不能用频谱表示(能量无限),但是由维纳辛钦定理,功率有限信号的功率谱函数语气自相关函数是一对傅里叶变换,所以我们可以由信号的自相关函数求得其功率谱,然后利用功率谱来描述随机信号的频域特性.三.本文中用到的三种谱分析方法Ⅰ.古典谱估计之相关法相关法谱估计是以相关函数为媒介来计算功率谱,又叫做间接法它的理论基础是维纳-辛钦定理,其具体步骤如下:第一步,由获得的N点数据构成的有限长序列xn(n)来估计自相关函数,即:N1xNNn01R(m)x(n)x(nm)N第二步,由自相关函数的傅里叶变换求功率谱,即通信工程孙大江1(1)ˆˆS()()MjjmxxmMeRme以上两步经历了两次截断,一次是估计xR(m)时仅利用了x(n)的N个观测值()Nxn,这相当于对x(n)加矩形窗截断.该窗是加在数据上的,一般称为加数据窗.另一次是估计ˆS()jxe时仅利用了从-(M-1)到(M-1)的xR(m),这相当于对xR()m加矩形窗截断,将xR()m截成(2M-1)长,这称为加延时窗.式中xR(m)和ˆS()jxe分别表示对xR()m和S()jxe的估值.由于MN,使得上式的运算量不是很大,因此在FFT问世之前,相关法是最常用的谱估计方法.FFT问世以后,情况有所变化,这里只进行基本思想介绍,不再进行matlab程序仿真.因为阶段后的()Nxn可视作能量信号,由相关卷积定理x1R(m)[()*()]NNNxmxm对上式两边取(2N-1)点DFT,并将时域卷积变为频域乘积,有221212111ˆS()[()()]()NNjxNNNeXkXkXk于是可得用FFT求xR(m)方案如下:(1)对N长()Nxn充(N-1)个零,成为(2N-1)长的21()Nxn.(2)求(2N-1)点的FFT,得:222121210()()NnkNNNnXkxnW(3)求2211()NNXk(4)求(2N-1)点的IFFT,得:12x2121(1)11R(m)()2N-1NNmkNNkNXkW上面的相关运算中,充零的目的是为了能用圆周卷积代替线性卷积,以便采用快速卷计算法.Ⅱ,参数模型法谱估计之levinson-durbin快速递推法Levinson-durbin递推算法是解尤勒-沃克方程的快速有效的算法,这种方法利用自相关矩阵具有的一些好的性质,是运算量大大减少,在介绍Levinson-durbin递推算法之前,先大概了解一下尤勒-沃克方程:如前所述,P阶ＡＲ模型的系统函数为通信工程孙大江1G()1miiiHzaz可以看出，Ｐ阶ＡＲ模型有Ｐ＋１个待定系数，由上式，可得白噪声激励得到的系统输出1()()()piixnaxniGn可以理解为，用ｎ时刻之前的ｐ个值的线性组合1()piiaxni来预测ｎ时刻的值()xn预测误差为()Gn．在均方误差最小准则下，组合系数12,,,paaa的选择应使预测误差()Gn的均方值最小．经过一系列的运算，最终可以得到ＡＲ模型的正则方程121(),1,2,()(),0pixixpixiaRmimpRmaRiGm也就是尤勒－沃克（ｙｕｌｅ－ｗａｌｋｅｒ）方程．在上面方程中令ｐ＝１，得一阶ＡＲ模型的尤勒-沃克方程为111(0)(1)(1)(1)(1)(0)0xxxxRaRRaR可以解得1(1)(1)(0)xxRaR进而求得2211(1)(0)(0)[1(1)](0)xxxxRRRaR在一阶的基础上进行递推，将阶次为ｍ时的第ｍ个系数()mam定义为反射系数，用mk表示，于是可以将计算ｍ阶ＡＲ模型参数的Levinson-durbin递推算法表示如下：通信工程孙大江11111121[()()()]()()(),1,2,3,,1(1)mxmximmmmmmmmmRmaiRmikaiaikamiimk式中1,2,mp，20[()](0)xExnR．Ⅲ.特征分解法谱估计之MUSIC方法前面分别讨论了古典法和参数模型法的一个例子，下面讨论用特征分解法对白噪声中的正弦波频率进行谱估计．特征分解技术的主要思想是，把自相关矩阵中的信息空间分成两个子空间，即信号子空间和噪声子空间，这两个子空间中的矢量的函数在正弦波频率上有尖锐的峰，据此可以估计正弦波的频率．输入信号为()1()(),0,1,2,,1iiMjiixnAennN其向量形式为1iMjiiiAexeω由向量x可以求出自相关矩阵为21[]MHHiiiNNiEPxsωR=R+RxxeeI其中，1MHiiiiPsRee，2NNωRI，分别为信号自相关矩阵和噪声自相关矩阵．也就是说，数据自相关矩阵可以分解为信号自相关矩阵和噪声自相关矩阵之和．sR和ωR都是N行N列方阵,秩分别为M和N,即()rankMsR()rankNωR.如果已知xR或者其估计,则可以通过分解由xR得到sR,再进一步分解由sR得到ie,从而求出正弦波频率的估计.通信工程孙大江省略分解sR的过程,定义N个特征向量中非零的M12v,v,,v为主特征向量M+1M+2nv,v,,v为噪声特征向量,因为信号子空间与噪声子空间相互正交,所以信号向量ie与噪声子空间的向量M+1M+2nv,v,,v都正交,与它们的线性组合也正交,有1()0,1,2,NHijjjMaiMev式中(1)[1],(1,2,,)iijjNiieeiMe为M个正弦信号的频率.令(1)()[1]iijjNTeee当(1)()()[1]iijjNTiieeeee时,应有1()[]()0NHHjjjjMaevve即21()0NHjjjMaev因此可以定义一中类似于功率谱的函数211ˆ()()xNHjjjMPaev若令1ja,则有21111ˆ()()()()()xNNHHHjjjjMjMPevveev上式取峰值的M个就是M个正弦信号频率的估计.理论上这M个应使上式所示的函数值为无限大,但由于存在估计误差,所以ˆ()xP为有限值,但呈现尖锐的峰值.也就是说,ˆ()xP的M个最大值所对应的频率就是正弦信号频率的估计.通信工程孙大江四.程序和结果时域图如下:通信工程孙大江(1)间接法程序clear;n=1:1000;x=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.3*n)+randn(size(n));N=1000;forh=1:NR(h)=0;endform=1:Ns=0;fori=1:N-(m-1)s=s+x(i)*x(i+m-1);endR(m)=s/N;%估计自相关函数,P89页4.1.1endw=0:0.001:0.5;Sx=0;M=1000;fork=1:(M-1)Sx=Sx+R(k)*exp(-j*w*k)+R(k)*exp(j*w*k);%由自相关函数的傅里叶变换求功率谱,P89页4.1.2endSx=log10(abs(Sx));plot(w,abs(Sx));title('间接法谱估计');gridon;通信工程孙大江通信工程孙大江(2)Levinson法程序clc;clearallfs=100;%设采样频率为100，N=150;%数据长度,改变数据长度会导致分辨率的变化；n=0:0.01:2;%变化范围x=sqrt(20)*sin(2*pi*0.2*n)+sqrt(10)*sin(2*pi*0.3*n)+randn(size(n));form=1:N%接下来的循环计算自相关函数的无偏估计,仍使用P89也方法R(m)=0;endform=1:Ns=0;forn=1:N-(m-1)s=s+x(n)*x(n+m-1);endR(m)=s/N;endforM=1:N-5%接下来的循环定义初值----阶数M要小于数据长度Na(M,M)=0;FPE(M)=0;P(M)=0;end%接下来初始化：由Levinson公式a(1,1)=-R(2)/R(1);%计算一阶模型参数a(1,1),见P66页公式3.3.15；P(1)=R(1)+a(1,1)*R(2);%计算一阶模型的误差功率P(1),p66页底部,误差功率P=G^2；FPE(1)=P(1)*(N+1+1)/(N-1-1)%计算一阶模型的最终预测误差准则函数FPE(1)；sum=0;%接下来进行递推forM=2:N-5fork=1:M-1sum=sum+a(M-1,