第七章功率谱估计(2009II)

中国石油大学（北京）电子信息工程系第七章平稳随机信号的自相关和功率谱估计7.1概述7.2估计准则7.3均值与自相关的估计7.4经典谱估计经典谱估计实际求出的自相关函数观测到有限个样本数据：傅里叶变换功率谱估计中国石油大学（北京）电子信息工程系经典谱估计的质量估计值的均值中国石油大学（北京）电子信息工程系主瓣宽度

渐近无偏中国石油大学（北京）电子信息工程系能分辨出两个分立的线间的最小频率间隔

（估计谱＝理想谱*三角窗谱）

石油大学（北京）电子信息工程系频率分辨率：periodgram_PSD.m,f1=50,f2=51产生两个问题：（1）主瓣宽度，功率在主瓣内扩散，分辨率降低。中国石油大学（北京）电子信息工程系（2）存在旁瓣，功率“泄漏”，产生假峰，淹没弱信号。改进2.方差假定是高斯零均值过程经典功率谱估计不是一致估计中国石油大学（北京）电子信息工程系，如何改善？选择好的窗（1）主瓣窄；（2）旁瓣小6.4经典谱估计的改进任务：改进周期图的估计的性能目标：主要是改进方差的性能方法：平滑与平均1.平滑（Smoothing)中国石油大学（北京）电子信息工程系用对的加窗来实现窗函数中国石油大学（北京）电子信息工程系中国石油大学（北京）电子信息工程系中国石油大学（北京）电子信息工程系2.平均（Average）

理论依据:k个独立同分布随即变量:

将一个较长的信号分成若干段，对每一段求功率谱，然后平均之。中国石油大学（北京）电子信息工程系Welch平均法调用格式：[Pxx,f]=pwelch(x,window,overlap,nfft,Fs)式中x是随机信号，Fs是抽样频率，Nfft是对x作FFT时的长度，window是选用的窗函数,Noverlap是估计x的功率谱时每一段叠合的长度。缺省时，Nfft=256,noverlap=0,window=Hanning(Nfft),Fs=2。输出的Px是估计出的功率谱，按上述调用格式给出的是幅平方值，f是频率轴坐标。

中国石油大学（北京）电子信息工程系例3N=1000;n=0:N-1;Fs=1000;x=sin(2*pi*40/Fs*n)+sin(2*pi*100/Fs*n)+2*randn(1,N);nfft=1024;window=boxcar(100);overlap=20;[Pxx,f]=pwelch(x,window,overlap,nfft,Fs)中国石油大学（北京）电子信息工程系中国石油大学（北京）电子信息工程系window=boxcar中国石油大学（北京）电子信息工程系window=hamming中国石油大学（北京）电子信息工程系f1=80f2=100f1=90f2=100（1）分辨率低。由于不可避免的加窗影响，使得真实谱的功率在窗口主瓣内扩散，降低了分辨率。分辨率正比于，为所观测到的数据长度；优点：物理概念明确，可用FFT实现，计算快速、简单。缺点：中国石油大学（北京）电子信息工程系经典功率谱估计的特点：（2）功率“泄漏”到旁瓣，较大的旁瓣有可能“淹没”较弱的成分，或是产生的假的峰值。当分析的数据较短时，这些影响更为突出；（3）方差性能不好，不是的一致估计，越大，谱曲线起伏越大；（4）周期图的“平滑”与“平均”主要是用来改善方差性能，但这降低了分辨率，偏差大，没有一个窗函数能使估计的谱在方差、偏差和分辨率各方面都得到改善。中国石油大学（北京）电子信息工程系

为克服传统谱估计方法的这些缺陷，人们开始寻求另外的谱估计方法，三十多年来，以参数模型为基础逐步形成了一套新的谱估计理论：现代谱估计。

中国石油大学（北京）电子信息工程系现代谱估计的提出主要是针对经典谱估计的分辨率低和方差性能不好的问题。现代谱估计的内容极其丰富，涉及的学科及应用领域也相当广泛，至今，每年都有大量的论文出现。目前尚难对现代谱估计的方法做出准确的分类。从现代谱估计的方法上，大致可分为参数模型谱估计（AR模型、MA模型、ARMA模型、PRONY指数模型等）和非参数模型谱估计（熵谱估计法、最大似然法、特征分解法、MUSIC法等）。中国石油大学（北京）电子信息工程系第七章功率谱估计（续）7.5功率谱估计的参数模型7.6AR模型的正则方程与参数计算7.7AR模型谱估计的性质与阶次选择7.8AR模型的Burg算法7.9AR模型的修正协方差法7.10Pisarenko谐波分解法中国石油大学（北京）电子信息工程系基本思想：经典谱估计将观测到有限个样本数据：以外的数据令其为0，与实际不符，中国石油大学（北京）电子信息工程系参数模型谱估计是对随机信号建立参数模型，认为是由白噪声通过此模型产生的，再由参数模型计算功率谱。7.5功率谱估计的参数模型即是对建立的数学模型。2.由观测到有限个样本数据,估计的参数1.假定所研究的平稳过程是由一白噪声序列激励一线性系统所产生的输出；中国石油大学（北京）电子信息工程系3.功率谱估计随机信号通过LSI系统的输入输出关系中国石油大学（北京）电子信息工程系三种参数模型：差分方程系统函数中国石油大学（北京）电子信息工程系

AR(Auto—Regressive，自回归)模型全极点模型若：并假定：则：中国石油大学（北京）电子信息工程系

MA(Moving—Average,移动平均)模型若：则：全零点模型中国石油大学（北京）电子信息工程系

ARMA(Auto-RegressiveMoving-Average,自回归－移动平均）模型极—零模型ARMA模型如果：不全为零则：中国石油大学（北京）电子信息工程系AR模型:全极模型，线性，用的最多，被研究的也最多。MA模型：全零模型，看起来简单， 但是非线性；ARMA模型：极－零模型，二者的综合。中国石油大学（北京）电子信息工程系关于模型：

任何一个有限方差的平稳ARMA或MA过程都可以表示成唯一的、阶数可能无穷大的AR过程。

上述定理是重要的，因为如果在三种模型中选择了一个错误的模型，我们仍然可以通过一个很高的阶数获得一个合理的逼近。因此，一个ARMA模型可以用一个足够高阶的AR模型来近似。而AR建模比ARMA建模和MA建模在计算上要简单得多。正是由于这个缘故，AR模型得到了广泛应用。中国石油大学（北京）电子信息工程系KaySM,MarpleSL.SpectrumAnalysis:amodern Perspective.Proc.IEEE,69(Nov):1380-1419,1981MakhoulJ.LinearPrediction:atutorialreview.Proc.IEEE,62(April):561-580,1975KaySM.ModernSpectrumEstimation:TheoryandApplication.19884MarpleSL.DigitalSpectrumAnalysiswithApplication.1987

推荐如下参考文献：7.6AR模型的正则方程与参数计算未知参数中国石油大学（北京）电子信息工程系估计自相关：求未知参数：求解方法：由AR模型的差分方程入手：两边同乘，求均值，得根据观测到有限个样本：7.6.1正则方程中国石油大学（北京）电子信息工程系

和的互相关， 是白噪声中国石油大学（北京）电子信息工程系正则方程（NormalEq.）中国石油大学（北京）电子信息工程系正则方程（丁）p148(4.3.9)中国石油大学（北京）电子信息工程系又称Yule-Walker方程中国石油大学（北京）电子信息工程系利用Yule-Walker方程，可求解出AR模型的参数：有快速算法中国石油大学（北京）电子信息工程系7.6.2Levinson-Durbin快速算法Yule-Walker方程中国石油大学（北京）电子信息工程系思路：利用Toeplitz矩阵特点，由低阶高阶系数：反射系数m正在计算的阶次在阶次为m时系数的序号中国石油大学（北京）电子信息工程系p=2递推公式递推过程中，要始终保持：中国石油大学（北京）电子信息工程系k阶Levinson算法的计算步骤：Step1：由估计Step2：解Yule-walker方程 得Step3：中国石油大学（北京）电子信息工程系用FFT计算中国石油大学（北京）电子信息工程系AR2005.m7.7AR模型谱估计稳定性及阶次选择7.7.1

AR模型的稳定性能否保证始终是稳定的？取决于R的性质中国石油大学（北京）电子信息工程系可以证明（胡第10章）若正定，则保证求解的模型的稳定性。（自相关函数是估计出的）中国石油大学（北京）电子信息工程系7.7.2AR模型阶次p的选择Levinson递推给出：（1）最终预测误差准则（2）信息论准则递减、恒正p中国石油大学（北京）电子信息工程系MATLAB中，PowerSpectralDensity(PSD)estimateviaYule-Walker'smethod.函数Pyulear,levinson和aryule中国石油大学（北京）电子信息工程系fb1=100.;fb2=500.;fb3=1000.;fb4=1500.例1中国石油大学（北京）电子信息工程系7.8AR模型的Burg算法先由观测到的有限长序列，估计自相关序列，Levinson递推算法：再由递推公式计算出参数较短时，误差较大。中国石油大学（北京）电子信息工程系求出其它参数。直接从观测到的有限长序列，令前后向预测误差功率之和为最小，估计出反射系数，再利用Levinson递推公式：Burg递推算法：中国石油大学（北京）电子信息工程系（1）前向与后向预测前向预测：即利用x(n)之前的p个值对x(n)作预测。误差功率前向预测误差序列ForwardPrediction中国石油大学（北京）电子信息工程系对同样一组数据，我们可以实现后向预测：BackwardPrediction后向预测误差功率后向预测误差序列前后向预测误差功率＝min

求出其它参数。Burg递推算法：令Levinson递推公式：中国石油大学（北京）电子信息工程系可以推导出前、后向预测误差序列有如下的递推关系：初始条件推导见板书中国石油大学（北京）电子信息工程系（2）Burg递推公式前向预测误差功率

中国石油大学（北京）电子信息工程系后向预测误差功率

中国石油大学（北京）电子信息工程系前后向预测误差功率

平稳过程，时间平均代替集平均？中国石油大学（北京）电子信息工程系有限长数据：P阶中国石油大学（北京）电子信息工程系递推步骤:1、令：求出2、求时的参数3、求出，再求4、用Levinson算法，求时的5、重复上述过程，直到Burg算法：一个公认的较好的算法。Burg2005.mfb1=500.;fb2=800.;fb3=1500.;fb4=1600.;Burg算法的特点：1.同时使用前向后后向预测，即使最小2.在每一级，仅对最小，然后套用自相关法的Levinson递推算法，影响分辨率；3.直接用数据递推，方法简单。4.对于白噪声加正弦信号，有时可能会出现谱峰分裂及谱线偏移现象。

中国石油大学（北京）电子信息工程系同Burg算法注意：这是Marple算法和Burg算法的最大区别。Burg算法仅：7.9AR模型的修正协方差法——Marple法中国石油大学（北京）电子信息工程系上述最小化的结果是得到一个协方差方程：注意：该矩阵不是Toeplitz矩阵，因此不能用Levinson算法求解。中国石油大学（北京）电子信息工程系Marple于1983年给出的求解上式的快速递归算法。所以，该算法称作“改进的协方差法，或Marple算法。该算法的估计性能最好，但计算复杂。参考文献：MarpleS.L.,Anewautogressivespectrumanalysisalgorithm.IEEETrans.OnASSP,1980,12:441~454中国石油大学（北京）电子信息工程系中国石油大学（北京）电子信息工程系7.10Pisarenko谐波分解法信号：p个正弦+白噪声下面推导x(n)满足的方程先设x(n)为单个正弦波：中国石油大学（北京）电子信息工程系Z变换特征多项式单个正弦波三角恒等式中国石油大学（北京）电子信息工程系一对共轭复根P个正弦波，P对共轭复根由此决定正弦频率中国石油大学（北京）电子信息工程系对应的差分方程考虑加上白噪声将代入ARMA中国石油大学（北京）电子信息工程系用矢量表示中国石油大学（北京）电子信息工程系是自相关阵的特征值，而是对应特征值的特征向量，求解转化成自相关阵的特征值分解—Pisarenko谐波分解中国石油大学（北京）电子信息工程系当y(n)由p个频率为wi，功率为Pi，i=1,…,p

的正弦+白噪声组成，中国石油大学（北京）电子信息工程系上述两个矩阵，各列线性相关，故秩=1，两矩阵之和的秩=2。若p个频率互不相同，则p个这样的矩阵之和的秩=2p，R为非负定矩阵，其2p个特征值中国石油大学（北京）电子信息工程系Ry的维数为2p+1，其中2p个特征值都是信号和噪声的贡献，而另一个则是噪声的贡献，其值为（最小特征值）是最小特征值所对应的特征向量。求出，再解下列方程，求出2p个根，从而得到具体计算步骤见p169(丁）PISAR_1_SIN_REAL.M与本章内容有关的MATLAB文件：1.pyulear.m用AR模型的自相关法估计信号的功率谱，其基本调用格式是：

[Px,F]=pyulear(x,order,Nfft,Fs)2.pburg.m用AR模型的Burg算法估计信号的功率谱，其基本调用格式是：

[Px,F]=pburg(x,order,Nfft,Fs)（一）、有关功率谱估计的MATLAB文件中国石油大学（北京）电子信息工程系3.pcov.m用AR模型方差方法估计信号的功率谱，其基本调用格式是：

[Px,F]=pcov(x,order,Nfft,Fs)4.pmcov.m用AR模型的改进的方差方法估计信号的功率谱，其基本调用格式是：

[Px,F]=pmcov(x,order,Nfft,Fs)中国石油大学（北京）电子信息工程系

5.pmem.m最大熵功率谱估计，其估计性能类似pyulear,其基本调用格式是：

[Px,F]=pmem(x,order,Nfft,Fs)6.pmusic.m用自相关矩阵分解的MUSIC算法估计信号的功率谱，其基本调用格式是：

[Px,F]=pmusic(x,order,Nfft,Fs)中国石油大学（北京）电子信息工程系7.peig.m用自相关矩阵分解的特征向量法估计信号的功率谱，其基本调用格式是：

[Px,F]=peig(x,order,Nfft,Fs),[Px,F，V,E]=peig(x,order,Nfft,Fs),x:信号向量，order:模型的阶次，Fs:抽样频率，Nfft:对x作FFT时的长度。Px:估计出的功率谱，F是频率轴坐标。对peig,输出的E是由自相关矩阵的特征值所组成的向量，V是由特征向量组成的矩阵。V的列向量张成了噪声子空间，V的行数减去列数即是信号子空间的维数。

中国石油大学（北京）电子信息工程系（二）有关AR模型参数估计的文件：8.aryule.m用自相关法(即Yule-Walker法)估计AR模型的参数，其基本调用格式是：[a,E]=aryule(x,order),[a,E，k]=aryule(x,order)中国石油大学（北京）电子信息工程系9.arburg.m用Burg算法估计AR模型的参数，其基本调用格式是：[a,E]=arburg(x,order)[a,E，k]=arburg(x,order)10.arcov.m用方差方法估计AR模型的参数，其基本调用格式是：[a,E]=arcov(x,order)中国石油大学（北京）电子信息工程系11.armcov.m用改进的方差方法估计AR模型的参数，其基本调用格式是：[a,E]=armcov(x,order)x：信号向量；ord