第五章随机信号的功率谱估计

1、第五章第五章 随机信号的功率谱估计随机信号的功率谱估计u 功率谱估计的经典和现代方法功率谱估计的经典和现代方法u AR模型法的功率谱估计u AR模型法的主要性质u Yule-Walker方程的Levinson Durbin求解算法u 格型滤波器u AR模型参数提取算法u 噪声对AR谱估计的影响u ARMA和MA模型法简介u 白噪声中正弦波频率的估计1u随机信号不能直接进行傅立叶变换，但平稳随机过程的自相关函数与功率谱密度之间存在傅立叶变换关系；u随机信号的功率谱与确定信号的频谱作用类似，是随机信号的频域特征分析；u随机信号的功率谱估计在通信系统分析、噪声监测、信号检测、 模式识别、机械故障诊断

2、等领域应用广泛； u依据有限的N个样本观测数据对平稳随机过程的功率谱密度进行估计。u随机信号的谱估计分为：1.1. 经典谱估计经典谱估计2.2. 现代谱估计现代谱估计关于信号的频域分析关于信号的频域分析2经典谱估计经典谱估计 u基本思想：基本思想：以傅立叶变换为基础，附以平均、加窗、平滑等预处理或后处理；u主要方法主要方法：周期图法（BT法）、间接周期图法； u优点：优点：简单易行、计算效率高；u缺点：缺点：存在分辨率低、旁瓣效应、估计精度不高，短数据时更为突出；u适用范围：适用范围：长数据序列；u相互关系：相互关系： 二者均等效为进行矩形加窗处理，均可通过加窗函数改善性能； 二者均适合采用F

3、FT算法； 二者存在频率分辨率低的致命缺点，且窗函数改善无效。32|( )|( )limTTFSTFTIFT( ),( )( )()220, TTTTTfttftftF其 他确定信号 的功率谱密度：( )f t平稳离散随机信号x(n)的功率谱密度：1() () ( )()2 ()()jmxxjmxxmRmE x nm x nSedSRm e21( )(0)= ( )()2xxx nRE xnSd的功率：x(n)的自相关函数与功率谱密度存在傅立叶变换关系( )xS表示功率在频域分布的谱密度！4如果随机信号x(n)是各态遍历的，集平均自相关函数可以由一个取样时间序列的时间平均自相关函数替代： 1(

4、)() ( )21limNxNnNRmx nm x nN101( )() ( ), 1NmxnRmx mn x nmNN (有偏、渐进无偏)101( )() ( ), 1NmxnRmx mn x nmNNm (无偏估计) 实际中信号自相关函数需要有限样值 估计：(0), (1), (1)xxx N 称为取样自相关函数取样自相关函数。( )xRm5 自相关法自相关法(B-T(B-T法）法）直接周期图法直接周期图法 1958年，Blackman 和 Tukey提出。先估计信号的自相关函数，再求出信号的功率谱密度估计 ： (0), (1), (1)( )( )xNxxxxx NR mS 间接周期图间

5、接周期图(Periodogram)(Periodogram)法法: :取样自相关函数实际上是下x(n)与x(-n)的卷和,即1( )( )()xRmx nxnN*211( )( )( )|( )|xSXXXNN间接周期图法有两种理解：(1)对信号进行加窗处理得xN(n)，再进行离散傅立叶变换得X()，再求模的平方得功率谱密度；(2)对信号xN(n)进行周期延托，再计算功率谱密度。6( )( )( ) (0), (1), (1)Nx nxnx nxxx N2101()( )NjnxnSx n eN直接相关法和间接周期图法得到的谱估计相同。1(1)( )( )Nj mxxmNE SE R m e2

6、( )( )j mNxmgm R m e1(1)( )Nj mxmNR m e21,|-1()0, NmNgm其中其它21( )( )( )( )2xNE SGSS(有偏，渐进无偏)对x(n)加窗处理：( )X21|( )|XNFFT( )Nxn7主要缺点：主要缺点： u相关图法主观认为未观测数据都等于0，造成频谱能量的泄漏；u假设数据是以N为周期的周期性延拓，把不真实的信息加于时间序列之上，频率分辨率低。常用改进方法：常用改进方法： 改进数据加窗方法，降低谱的旁瓣泄露：将矩形窗改为其他窗函数，如：汉宁(Hanning)窗 、汉明(Harmming)窗、布拉克曼(Blackman)窗、三角窗(

7、Bartlett)、凯塞窗(Kaiser)等。 对B-T法的相关函数加窗：()()()jmxxmSw m Rm e窗函数8u 将长度为N的序列分K段，每段长为M，分别对每段进行谱估计，再进行总平均，得平均周期图(如下图)。u 如各段数据相互独立，则所得估计的方差为原来不分段时的 。缺点是点数减少，分辨率进一步恶化。 修正周期图法(分段平均)(Welch，1967)( )11()()KkxxkSSK01N 2k M1k MkKM(1)( )xS(2)( )xS()( )KxS为改善点数减少恶化频率分辨率的缺陷，常对数据进行交叠分段，并进行功率谱平均，如下图： 01N 3k M1k MkKM2k

8、21 K9现代谱估计现代谱估计u数学基础：以随机过程或信号的的参数模型为基础，又称为参数模型法或参数法。u主要发展：从非工程领域(如实验数据和观测数据的处理、统计学)的时间序列分析到工程领域的现代谱估计；现代谱估计始于上世纪60年代，主要经历了以下阶段： 线性预测滤波 最大熵谱估计(Burg,1967) 自回归(AR)谱估计方法(1968，Parzen) Pisarenko谐波分解 多重信号分类算法(MUSIC,1981，Schimit) HOS方法10功率谱估计的参数法（现代谱估计）功率谱估计的参数法（现代谱估计）参数法谱估计有理谱模型线谱模型（谐波）AR模型ARMA模型MA模型uAR模型

9、(全极点模型)uMA模型(全零点模型)uARMA模型00( )( )( )qkkkpkkkb zB zH zA za z011( )( )pkkkH zA za zqkkkzbzBzH0)()(11任何平稳随机信号x(n)都可以看成由白噪声序列 激励一个因果和稳定的线性时不变系统H(z)产生的输出。 ( )u n( )H z( )u n( )x n平稳随机信号模型 任何有限方差的广义平稳过程可以分为完全随机的部分和确定的部分，完全随机部分对应的功率谱为连续； 确定的随机过程完全可以通过过去无限个取样值加以预测； 任何ARMA过程或AR过程可以用无限阶的MA过程来表示；任何ARMA过程或MA 过

10、程可以用无限阶的AR过程来表示。谱分解定理及其推论是谱分解定理及其推论是参数法谱估计的理论基础参数法谱估计的理论基础12u 功率谱估计的经典和现代方法u AR模型法的功率谱估计模型法的功率谱估计u AR模型法的主要性质u Yule-Walker方程的Levinson Durbin求解算法u 格型滤波器u AR模型参数提取算法u 噪声对AR谱估计的影响u ARMA和MA模型法简介u 白噪声中正弦波频率的估计第五章第五章 随机信号的功率谱估计随机信号的功率谱估计13AR模型参数估计法的功率谱估计：模型参数估计法的功率谱估计：21( )( -( )(0,)(pkku nx na x n kNu n

11、式中，222-2111( )( )( )( )1|1|uppkjkkkkkH zSH zSa za eu基本原理：根据随机采样样本x(0), x(1) x(N1) 估计随机时间序列的功率谱密度：u求解方法(分三步)： 模型阶数p不确定时数学上很难处理，因此先假定p，求模型参数。 阶数p已知时，对模型两边同求某种统计特征以将随机变量转化为确定性的量。14对各种阶数下的模型进行比较，应用某种准则估计选择最好的模型(得阶数p、ak及 )。2u AR(p)模型的Yule - Walker方程组：2121(0)(1)(2)( )(1)(0)(1)(-1)0(2)(1)(0)(-2)0( )(-1)(-2

12、)(0)0pRRRR paRRRR paRRRR paR pR pR pR 15u功率谱估计步骤：N个样值x(0),x(1)x(N-1)模型参数和误差功率212 ( ,)pa aa估计自相关函数 (0), (1), ( )RRR p功率谱密度2-21( )/ |1|kpjkkSa e16u AR模型阶数的选择：阶数对估计性能的影响n阶数选得太低，功率谱被平滑太厉害，无法分辨真实峰(P130图4.3)；n阶数选得太高，谱分辨率提高，但产生虚假峰(P130图4.4)。n也可用以下matlab源程序验证：N=512;Nfft=1024;Fs=2*pi;n=0:N-1;xn=cos(0.3*pi*n)

13、+ cos(0.32*pi*n)+randn(size(n);order=10;figure(1)pburg(xn,order,Nfft,Fs)title(Burg Algorithm,p=10)order=100;figure(2)pburg(xn,order,Nfft,Fs)title(Burg Algorithm,p=100)选择方法：1.实验方法：观察拟合误差法。2.分析方法： 最终预测误差（FPE）准则。 Akaike信息准则（AIC）。 判别自回归传输函数（CAT）准则。17221( )(, )1kkNkFPR kfkNk定义：最终预测误差N为观测数据长度。使上式最小化的阶数k即为

14、最优阶数最优阶数。FPE准则得到模型阶数一般偏低。(1) 最终预测误差(Final Prediction Error, FPE)准则(2) Akaike信息准则(Akaike Information Criteria, AIC）2( )ln2kAIC kNk使上式最小化的阶数k即为最优阶数。AIC准则得到模型阶数一般偏高。(3) 判别自回归传输函数准则(Criteria Autoregressive Transfer Function, CAT)18估计性能：估计性能：精确分析很困难，只能给出大样本理论的近似关系。u估值的均值（N，p ）：u估值的方差（N，p ）： ( )( )E SS224

15、( ) , 0, ( )2( ) , pSNVar SpSothersN22221111CAT( ) (kjjjjkNkNNj其中）最小化上式得最优阶数。19第五章第五章 随机信号的功率谱估计随机信号的功率谱估计u 功率谱估计的经典和现代方法u AR模型法的功率谱估计u AR模型法的主要性质模型法的主要性质u Yule-Walker方程的Levinson Durbin求解算法u 格型滤波器u AR模型参数提取算法u 噪声对AR谱估计的影响u ARMA和MA模型法简介u 白噪声中正弦波频率的估计20性质性质1：隐含着自相关函数的外推：隐含着自相关函数的外推估计功率谱与自相关函数的关系：令：利用谱

16、分解定理，并考虑功率谱与自相关函数的关系：-0101 ( )( )( )1, 1ppkkkkkkH zA zA za za za 为因果系统21-( )( )( ) ()mxxmSzR m zA z A z 211- ( )( )(1/ )mxmA zR m zAz ( )(0) lim ( ) 1, 0(- )( ) zH zhH zmh mm因为滤波器是因果的，且故当时,212*2()( ) , 0(1/ )hmmfor mAz 左边：211-0 ( )( )* ( )()pmmkmkA zR m zaR ma R m k 右边： 20 ( ) (), 0pkkma R m kfor m

17、故有： 201( )( )(), 0 (1)pkkR mma R mkfor ma 设2111-(- ) , 0-(- ) , 1,2,-(- ) , ( ) ()kkkpkpkpka R m ka R m kma R m kmpR mRpmm自相关函数偶对称22结论： AR谱估计是将有限个自相关函数值按照Yule - Walker方程进行外推后进行傅立叶变换得到的结果。 由于AR谱估计将自相关函数进行了外推，克服了经典谱估计方法加窗导致分辨率低和旁瓣“泄漏”的问题，因此，AR谱估计有高的分辨率。ker( )( )1( )Yule WalR mR mmNR m此式与方程相同，只是用代替了。即：

18、对于的不是简单假设为零，而是按表达式外推，更合理！23性质性质2：与最大熵谱估计等效：与最大熵谱估计等效最大熵谱估计的提出(Burg)：u经典谱估计方法具有分辨率低和旁瓣“泄漏”的问题。其根本原因是自相关函数加窗，这样克服这些问题必须对自相关函数进行外推。uBurg提出以最大熵作为自相关函数外推的准则，其合理性在于这样处理后，时间序列的随机性最大、对自相关函数的约束最少、功率谱最平坦。信息熵 (Entropy)是信息论中的一个非常重要的概念，是信号含信息量的度量，表征了信息发生的随机程度。信息量：事件X，事件 发生时(概率 )，包含的信息 以e为底的信息量单位为nat (奈特)；以2为底的单位

19、bit (比特)kXxkP1()loglogkkkI xPP 24信息量的以下几个性质是显然的：01, ()0 (kkPI x1 .若 则肯定发生的事件发生时，不带来任何信息)；02 . ()001 ,kkI xP任何事件发生时，都不会造成信息的丢失；0,()()kjkjPPI xI x3 .若则。熵：熵：平均信息量称为熵，记为H(X)()()logkkkkkkxxH XP I xPP 设X为连续型随机变量，则：()( ) ( )( )log( ) ( )H Xp x I x dxp xp x dxE I x 25121121 ( )exp()(2 )det ( )ln( ) ln(exp()

20、ln(2 )ln(det)2 TNTxNxTxNxp xx RxRH xEp xNEx R xxRxx 设 维零均值正态随机矢量，则1ln(2 )ln( )lndet)2 ln(2 )ln(det)2( )ln( )det( ) (d1( )ln()2et)xxTxNxxxxxxNRNNRH xSRSHExSSdHRxxR也就是说，定义：信号的功率谱熵为：26已知 共2p+1个样本相关函数，使( ),0, 1,xR kkp( )maxxH S问题：求 ( )xS()()1max( )maxln( )21( )( )( ), 0, 1,2xxxxSSj kxxxH SSdR kSedR kkp约

21、束优化问题：约束条件：与相关函数相匹配与相关函数相匹配11( )ln( )( )( )22pj kxxkxxkpJ SSdR kSed拉格朗日约束多项式：27( )110( ) ( )kxxkppj kj kxkkkpkpJ SSSee与与ARAR功率谱等价功率谱等价实际上，由Fejer-Riesz定理(证明需用泛函知识)：( )( )ppkj kkkkpkpjz eW zzWe 220( )( )piiiW zA za z若 ，则一定可以找到一个 满足( )0W0( )piiiA za z而且若 的根全部在单位园内，则A(z)是唯一确定的。故有：( )0A z 282011( ) xppj

22、kj kkkkpkSea e1( )( ), 0, 1,21ln( )( ), 0, 1,2j kj kSedR kkpSedC kkp相关函数匹配：约束条件倒谱匹配：1ln( )max2Sd11( )ln( )( )( )221 ( )ln( )2pj kkkpqj kkkqJSdR kSedC kSed更深一步，功率谱(spectrum) 和 倒谱(cepstrum)( )Sln( )S( )qpj kj kkkkqkpSee与与ARMA功率谱等价功率谱等价29结论： AR谱估计相当于对自相关函数以序列信息熵最大为原则进行外推后进行傅立叶变换的结果； AR谱估计相当于在p+1个自相关函数值

23、确定的情况下，以功率谱密度最平坦为准则得到的估计结果； AR谱估计相当于在p+1个自相关函数值约定情况下，以功率谱熵最大为准则得到的估计结果。性质性质3：与最佳线性预测谱估计等效：与最佳线性预测谱估计等效1( )()pkkx na x nk 2, 1, 2 Min( )( ( )- ( ) kakpnEx nx n最佳问题的定义：最佳线性预测谱估计：用随机时间序列x(n)前p个时刻的值的线性组合来预测当前值，即：30min11() , 0( )() , 0pkkpkka R mkmR ma R mkmu 最佳线性预测滤波器权系数等于AR模型参数，最小预测误差功率min等于AR模型中激励噪声方差

24、2。u 上式还与AR(p)模型的Yule -Walker方程组相同。若二者自相关函数相同，则解必然相同。线性预测结果：预测误差滤波器：预测误差滤波器：1-1 ( )( )()( )( ) )1)(pkkpkkkAe nx na x nkE zX zA za zz 式中，于是22min-01|()|( )( )2jxnA eSdE e n31结论： AR谱估计相当于用时间序列前p个时刻的值在最小均方误差准则下来预测当前值并作为外推后值进行谱估计的结果。 预测误差滤波器传输函数A(z)是模型参数确定的最小相位多项式。(1)x(0)x( )x n(1)x n()x np( )x nN已知 个数据外推

25、( +1)x n性质性质4：等效于输入序列的最佳白化处理：等效于输入序列的最佳白化处理u AR模型的输入信号为白噪声，线性预测相当于AR模型的逆系统模型；u AR谱估计与最佳线性预测谱估计等效，若最佳线性预测输出的误差信号接近于白噪声，则AR谱估计与输入序列的白化处理等效；u 通过“谱平坦度” 衡量误差预测滤波器输出的误差信号的白化程度。32谱平坦度：谱平坦度：随机时间序列x(n)的谱平坦度定义为：-1expln ( )2011( )2xxxxSdSd( )1xxxS是的几何均值与算数均值之比。越大，谱越平坦。特别地，白噪声的。-, 1, 2 -1expln ( )2 1( )2keeakpe

26、SdMaxSd问题的描述：问题的描述：预测误差滤波器A(z)为最小相位滤波器，输入时间序列x(n)经线性预测后的输出误差序列为e(n)，求解以下优化问题：33问题的求解：问题的求解：2-1ln |A()| d02-211ln ( )ln |( )|( )221 ln ( )2exxSdASdSd-11expln ( )expln ( )2211( )( )221( )(0)21(0)( )2exeeexxxxeeSdSdSdSdSdRRSdu 若要使e最大，则Re(0)要最小；u 预测误差谱平坦度最大(白化度最大)等效于预测误差功率最小；u AR谱估计器等效于对输入时间序列的白化处理。34结论

27、：AR谱估计等效于预测误差最佳白化处理；根据预测误差的白化程度可以判断时间随机序列对AR(p)模型的符合程度。AR模型法谱估计主要性质小结：1. 相当于根据Yule - Walker方程组对自相关函数值外推后进行傅立叶变换的结果；2. 相当于对随机时间序列以最大熵准则外推后估计信号的功率谱密度；3. 相当于对随机时间序列以最佳线性预测外推后估计信号的功率谱密度；4. 估计的功率谱密度是有限数量自相关函数值条件下最平坦的功率谱密度；5. 功率谱估计等效于对输入时间序列进行最佳白化处理。35第五章第五章 随机信号的功率谱估计随机信号的功率谱估计u 功率谱估计的经典和现代方法u AR模型法的功率谱估

28、计u AR模型法的主要性质u Yule-Walker方程的方程的Levinson Durbin求解算法求解算法u 格型滤波器u AR模型参数提取算法u 噪声对AR谱估计的影响u ARMA和MA模型法简介u 白噪声中正弦波频率的估计36关于Yule Walker方程组的求解： 求解AR(p)模型的Yule - Walker方程组的常规方法运算量大（p3级）； 计算结果仅为p阶模型的参数； Yule - Walker的系数矩阵很有规律性，可以构造迭代算法进行递推求解，以减少计算量。设已经迭代求得了AP(k)模型的参数：2,1,2,kkk kkaaa需要用AP(k)模型的参数求出AP(k1)模型的参

29、数：21,11,21,11kkkkkaaa372,1,1(0)(1)( )(1)(0)( -1)0( )( -1)(0)0kkk kRRR kaRRR kaR kR kR也就是说，已知k步的参数满足（条件方程组条件方程组）：需要求解（k+1）步目标方程组目标方程组的参数：211,11,1,11(0)(1)( )(1)(1)(0)( -1)( )0( )( -1)(0)(1)0(1)( )(1)(0)0kkkkkkRRR kR kaRRR kR kaR kR kRRaR kR kRR 38,00(1- ) , 1kkk ikiDa R kia式中，为此，由条件方程组构建扩大方程组扩大方程组：2,

30、1,(0)(1)( )(1)1(1)(0)( -1)( )0( )( -1)(0)(1)0(1)( )(1)(0)0kkk kkRRR kR kRRR kR kaR kR kRRaR kR kRRD 将扩大方程组与目标方程组比较，发现模型参数缺少一个自由度，故利用Toeplitz性质再由扩大方程组构造一个预备方程组，为增加自由度创造一个条件：39,12(0)(1)( )(1)00(1)(0)( -1)( )0( )( -1)(0)(1)(1)( )(1)(0)1kk kkkDRRR kR kRRR kR kaR kR kRRaR kR kRR 预备方程组：预备方程组：目标解由扩大方程组的解和预

31、备方程组的解线性表示，即构成模型系数的迭代关系式：1,1,1,111,1,1,11111001kkkik ikk kkkkk kkkk kikaaaaaaaaaa (1,2, )ik反射系数40111222221221120000(1)kkkkkkkkkkkkkkkDDDD显然，方程组间的右边列矢量也满足对应的线性关系，由此得到模型参数 的迭代关系式，即：2k2,1(0)( )kkk iiRa R i式中，小结：(1)Y-W方程组求解的L-D算法过程是一个从AR(1)到AR(p)的迭代过程；(2)迭代计算需要 或 可获得p组模型参数(包括系数和噪声方差)；(3)利用模型参数可估计对应阶AR模型

32、的时间序列功率谱。( )R m( )R m41设已经得到AR(k)模型的参数，迭代求AR(k+1)模型参数的步骤共6步：102,10,000120222101,1,101,111 .(0)( );2 .(1),1;3 .;4 .(1);5 .(1,2, );6 .()kkkk iikkk ikikkkkkkik ikk kikkkRa R iDa R kiaDaaaika 由系数迭代关系的最后一个方程42L-D迭代算法的实现流程图：(1)R k ( )R k(1)R k (2)R(1)R( )R k(1)R k (2)R(1)R11,1ka,1ka,2ka,2ka,1k ka,1k ka,k

33、ka,k ka1k1k1k1k21k2k12(1)k1,kka1,1kka1,2ka1,1ka21k11,1kkka (0)RkD43222212,022,12,220,01,02,031( )(0)1(1)(2)42(2)1( )( )( )-11xkkx nRRRARA zSzA zL Dazaaaaa例：已知平稳随机序列的自相关值：、，先采用模型的参数估计方法对其功率谱进行估计：，试用算法递推出模型参数：、。（）解：先给出L-D算法的6个递推公式：1020,100022211201,1,101,111 .(0)( ); 2 .(1);3 .; 4 .(1);5 .(1,2, );6 .k

34、kkkk ikk iiikkkkkkik ikk kikkkRa R iDa R kiDaaaika 4412000,02220110201,112111,011,12201103(0)1; (1)437; (1)4163 ; ()4171; (2)(6611)11kRDaRDaakDa Ra RD 当时，有：当时，有：2222212,222,11,122,01,13(1)71671; 7; ;(7) 1aaaaa 45第五章第五章 随机信号的功率谱估计随机信号的功率谱估计u 功率谱估计的经典和现代方法u AR模型法的功率谱估计u AR模型法的主要性质u Yule-Walker方程的Levin

35、son Durbin求解算法u 格型滤波器格型滤波器u AR模型参数提取算法u 噪声对AR谱估计的影响u ARMA和MA模型法简介u 白噪声中正弦波频率的估计46(1)( )(2) (1)(3)kAR pkkAR模型随机序列的功率谱估计过程与最佳线性预测等效；阶预测误差与格型滤阶预测误差之间存在递推关系；各阶反射系数与波器设计的出发各阶模型点：一一对应。,0101( )( )( ) ( )()( )()1 ( )( ) ( ) kkkk ikk ikiikkEzA z X zx nx nenx na x niena x nia 定义 ：时间序列前向预测值和前向预测误差：,1,0002( )()

36、( ) ()() ( )()( ()1kkkk iikkik ikk ik kRikiikx nx nkenx nka x nkiena x nkiEzAz Xax nia 定义 ：时间序列后向预测值和后向预测误差：( ) z47,00,0( )(0,1,)( ) ( )(1)(0,1,)( )( ) ( )kk ikkikk ikiRk ik k ikkkk imRik k iiAR kaikzA zA za zaaaikzA zAzAzaza 模型的系数序列为，其 变换为多项式：的倒序列为，其 变换为的倒序多项式：0,11,0 ()()k mk mm kkkmkk mkmzzazzA z

37、1,1,11,1,1111(1,2, ) kik ikk kikkkiiikik ikk kiiiiaaaikazazaz 4811,11,1,11111(1)1,1( ) 1 ()kkkkkkiiikik ikk kiiiikikikiAzkakazazazazza 11,1,111( )( )(1)11Rkkkkkiiikk kiiiAzzkkAzzazz 1(1)11111 ( )( )() ( )( )(1) kRkkkRkkkAzzAzAzA zz Az(1)11111( )111) ()() ( ()kkRkkRkkkkz AzAzkkkkkRkkkzz AA zz AzzzA z

38、zzAA zz ( )(2)z49111111111111(1)(2)( )( )1 ( )( )( )( )( )( )1 ( )(kkkkkRRkkkkkkkAzA zAzAzx nzX zEzEzEzzzzEz由和有：等式两边同时乘以的 变换，有1111)( )( )1 1( )(1)kkkkkkzzenenenen等式两边进行 反变换，得前向、后向预测误差的迭代公式：5000000000 0( )( )1 ( )( )( )( )( )( ) ( )( )( )( )( )( ) RRkA zAzEzA z X zX zenx nEzAz X zX zenx n迭代递推初值：设，则或

39、( ) ( )ppenen递推求解格型滤波器的递推流程如下图：( )pen( )pen12p1z1z1z( )x n5111111111( )( )11( )(1) ( )=( )+(1)( ( )(1)( )()kkkkkkkkkkkkkkpenenenenenenenenenenenx n由，有于是，若需由前向预测误差倒递推得到时间序列反向递推( )AR p，相当于建模，则仅需将上图格型滤波器结构上半部分反向，并将减法运算变成加法运算，如下图：( )pen( )pen12p1z1z1z( )x n52,1,(1)(2)1(3)(4)( )( ), (1), ()1,( )( )( )kkk

40、kkkkkkenx n x nx nkaaFIRA zennnex格型滤波器由各级的反射系数唯一描述；格型滤波器工作稳定的充分必要条件是；格型滤波器各级间互不耦合，均为最佳滤波，级联后为全程最佳；前向预测误差是序列激励单位脉冲响应为的预测误差滤波器产生的响应；或是 时刻几个重要事项以前的数据预测产明生的误差：说：,1,1(5)( )( ), (1), (),1( )( )()( )()()(6)( )()()kkRk kkk kkkenx n x nx nkaaaFIRAzennx nx nnkx nkenx nkx nk；后向预测误差是序列激励单位脉冲响应为的预测误差滤波器产生的响应；或是时

41、刻以后的数据预测产生的误差：。前向和后向预测处理的示意图如下：,1,2,2,1 k kk kk kkkaaaaa,1,2,2,1, kkk kk kk kaaaaa 1 1 -2 2 nknnknknn)( )(x nx n预测)()(x nkx nk预测53AR模型法的稳定性和功率谱估计界：稳定条件：u 一般自相关函数没有误差时能自动满足u H(z)的极点都在单位圆内u 12 22 p2 0u |k| 1，k=1, 2 p谱估计的界： Burg证明了AR谱的动态范围满足：111 |1 |(0)( )(0)1 |1 |ppiiiiiiRSR下界上界u 任何一个反射系数i接近于1时，上界变大，下

42、界变小。u 具有大反射系数模值的AR过程，其谱一定具有尖锐的峰。54第五章第五章 随机信号的功率谱估计随机信号的功率谱估计u 功率谱估计的经典和现代方法u AR模型法的功率谱估计u AR模型法的主要性质u Yule-Walker方程的Levinson Durbin求解算法u 格型滤波器u AR模型参数提取算法模型参数提取算法u 噪声对AR谱估计的影响u ARMA和MA模型法简介u 白噪声中正弦波频率的估计55问题的提出：问题的提出： Levinson - Durbin算法要先估计出几个自相关函数值，而自相关函数的估计是有偏估计，这将导致AR模型的估计精度降低，是否有其他办法提高AR模型的估计精

43、度？ Levinson - Durbin算法要先算自相关函数值，可否直接利用随机采样样本提取AR模型参数的方法？问题的定义：怎样根据x0, x1 xN1 这一随机采样样本估计随机时间序列的AR(p)模型参数？,1,1,1,Tpp pp paaaa,1( )( - )( )pp kkx nax n kn 56主要思想：主要思想：u AR模型法谱估计与最佳线性预测滤波等效。u AR模型的参数与线性误差预测滤波器的单位脉冲响应相同：E(z) = A(z) X(z)。u 对于平稳个态历经随机时间序列，可用时间平均代替集平均。u AR模型参数的提取可化为以下优化问题：2, 1, 2 ( ) kpakpM

44、inEen或2,1, 2 ( ) kpakpMinEen,1,( )(1) ( )(1)()( ) 1TTpp px nnpxx nx nx npAR paaa设为无限长时间序列，且 时刻的个样值列矢量为：模型的系数矢量为：57实际上，由于输入时间序列的样值有限，其自相关函数矩阵只能通过取样自相关（时间平均）进行估计，并按误差平方时间平均最小准则进行模型参数估计（提取）。构造的方式决定不同的参数提取算法： Y-W法（相关法）； 协方差法； Burg法（前两种方法的结合）。,022,1( )()( ) ( ) =1pTTpp iiTTTpTTTpTp pp penax nia xx aE ena

45、 E xxaaaE ena E xxaaaaaaRR 最佳线性前向预测误差：前向预测误差均方值：同理，对后向预测误差，有：这里，。58Yule - Walker法（自相关法）：法（自相关法）：误差平方时间平均最小准则：-1-1022, 1, 2 , 1, 2 01 ( ) ( )kkNpNpnnppakpakpMinenMinenN或Yule - Walker法（自相关法）计算 的示意图( )pen00(0)x0(1)x np ( )x n(1)x N 000(1)x n()x np11pappa,1p pa11pappa,1p pa11pappa,1p pa前补0后补0第一个数据x(0)进入

46、滤波器产生(0)pe最后一个数据x(N-1)留在滤波器产生(-1)peNp,1,2,( )(1)(1,)( (0), (1), (1)()pppp penpaaaNxxx NNp是长度为的预测误差滤波器单位脉冲响应与长度为 的数据序列的卷积结果，长度为。卷积过程中需要对有限长序列进行前、后补零，也等效为无限长序列加窗。如下图：59估计结果：-1-12,0000-1,000-10( )()() () ()1( - ) () () NpNppppp ip jnnijppNpp ip jijnNpnenax niax njax ni x njaR ijx ni x njN 取样自相关值：,00 (

47、- )ppTp ip jijNaR ij aNaaR取样自相关矩阵,00( )( - ) , 1ppp ipienax n ia602(1),1,2,min,1,2,2() ( , )(1,2,;0,1,2,) 1( (0)(1)(2)( )=ppTppp pppp pNaYWaR i jip jpaaaaN RaRRaR paa RR00R令取样自相关矩阵表示的方程式中，另外，因22min) 1NN，于是-1-10-1-()()0( )(1) (1)11() () ()( ) () 1 ()( )( ) () 1 ( ) (NpNp in i mnmiNpkjij kmkjkjkj kpR

48、kppR ijx ni x njx m x mijNNR ijR kx m x mkNx m x mN R是由取样自相关值组成的阶方阵。因为：-10), 0-1NkmkkN取样自相关序列61相关法的性能：相关法的性能：u 取样自相关矩阵是正定的，能保证系统的稳定性；u 求时相当于对时间序列进行了加窗处理（前后均补0），因此，估计精度不高。2, 1, 2 - 1 ( )kpakpNn pMinen误差平方时间平均最小准则：未对数据两端加未对数据两端加0(未加窗未加窗)，对观测时间以外数据不做任何假设。，对观测时间以外数据不做任何假设。(0)x(1)x p( )x n(1)x Np(2)x N (

49、1)x N ()x np(1)x np11pappa,1p pa11pappa,1p pa11pappa,1p pa(1)x( )x p(1)x n()x Np用协方差法计算 的原理图pe (n)协方差法：协方差法：6212,00-1 ( )()()( - )(1) (1)1 ( - )() () ppNpp ip jn pijNn penaR ij aNp aaR ijppR ijx ni x njNpRR误差平方的时间平均：为协方差估计值构成的矩阵：2min1-2NaNpaaR0R0性能特点：性能特点： 求 时没有对随机信号进行加窗处理，因此估计精度较高； 自相关矩阵不一定正定的，不能够保

50、证系统的稳定性。63例：试根据信号的4个取样值x(n)=2,4,1,3，分别用自相关法和协方差法估计AR(1)模型参数。例4.1（教材P143）解：(1)(1)自相关法自相关法412111,00( ) , ( )()inienena x ni( ):x n31421( )en11,1a1(0)2e11,1a11,1(1)42ea11,1a11,1(2)14ea 11,1a11,1(3)3ea11,1a11,1(4)3ea（右端添（右端添0）（左端添（左端添0）41101,11,1( )2( )0nenenaa令，有4111111111101111( )( )2 42)4(14)(3)900.5

51、nenenaaaaaa ，（64(2)协方差法协方差法( ):x n31421( )en111a，1(0)2e111a，111(1)42ea，111a，111(2)14ea ，111a，111(3)3ea，111a，111(4)3ea，31111111( )2( )0nenenaa，令111111112 42)4(14)(3)00.714aaaa ，（3211( )nen(不外推添不外推添0)(不外推添不外推添0)65101 ( )( ) ()(01) ( )2,4,1 ,3(0 )2; (1)4; (2)1; (3) 3Nknx nR kx n x nkkNNxxxx 另外一种方协方差法法，

52、直接通过和求解：以相关法为例请进行练习自相关序。列，协方差序列2011,1,11 ( )( ) ()(03)41(0)(2 24 4 +1 1 3 3)7.54 1(1)(2 44 1+1 3)3.754 (1)(0)0. 0 5 nR kx n x nkkRRRRaa 2min1,1 (0)(1)7.53.75 0.55. 6254 5.62522 5 .RRa66Burg法：法：概述：概述：u自相关法：自相关法：计算效率高，能保证预测滤波器是最小相位，但对数据两端添加了0（加窗处理），估计精度下降，短数据时性能下降。u 协方差法：协方差法：计算效率高，未加窗，估计精度高，但估计自相关矩阵不

53、一定正定，导致预测滤波器不是最小相位系统，可能不稳定。u Burg法法的思想：一方面希望利用已知数据段以外的未知数据（但不做主观臆测）；另一方面使预测误差滤波器是最小相位的。不直接估计AR模型参数，先估计反射系数，再利用Levinson递推算法由反射系数求得AR模型参数。67误差平方时间平均最小准则：22, -11, 2 ( )( ) kppaNkpn pMinenenBurg法前后向预测误差产生原理图(0)x(1)x p(1)x Np(2)x N (1)x N 1,1pa,p pa,1p pa(1)x( )x p()x Np1,1pa,p pa,1p pa,p pa,1p pa1,1pa,p

54、 pa,1p pa1,1pa(-1)peN(-1)peN( )pep( )pep参数计算方法参数计算方法：利用Levinson - Durbin迭代算法以及格形滤波器预测误差的迭代关系式，在已知p-1阶模型参数的情况下，只要求出p即可求出p阶模型参数。68-1-11( )( )1( )( -1)ppppppenenenen1-1-12( )(1)( )( )0Nppppn ppen enen en11-11221-12( )(1)( )(1) Nppn ppNppn pen enenen 1p显然，k反射系数 的迭代递推关系。69,-1,-1, -1,0,-1-100220 , 1, 2 1,

55、(1,)1( )( ) 1( )(1)( )( )( ), 011p ipippp ipp ppppppppaaaipaaenenenenenenx nnNxN 由反射系数的递推关系得模系数预测误型参数的递推关差初始条：件系：100,0( )1Nnna70 设已知有限数据序列x(n),n=0,1,N-1，可按下步骤计算模型参数系数，并在此基础上计算功率谱。1. 置置k=0, 计算初值：计算初值：2. k=k+1, 计算反射系数：计算反射系数：00122000,0( )( )( ) (01); 1 ( ) 1Nnenenx nnNxnNa1-11-122112( )(1)( )(1) Nkkn

56、kkNkkn ken enenen71 5. 计算计算k阶预测误差功率：阶预测误差功率：2221(1)kkk6. 回到步骤回到步骤(2)-(5), 进行下一次迭代。进行下一次迭代。4. 计算前、后向预测误差：计算前、后向预测误差：1-11( )( )1( )(1)kkkkkkenenenen3.计算滤波器系数：计算滤波器系数：,1,1,1,0,(1,.,1),1,k ikikkk ikk kkaaaikaa 721-11122112( )(1)( )(1)()Nppn ppNpppn ppen enenenwnwn算法修正：算法修正：u如果处理数据来自AR过程，则可获得精确的结果，同时系统的稳

57、定性也有保证。u如果处理的是正弦信号会遇到一些困难，例如：谱线分裂、谱峰位置受相位影响大等。u为了减小相位的影响，可对反射系数估计公式进行如下修正：其中：wp(n)为某一非负的窗函数。73第五章第五章 随机信号的功率谱估计随机信号的功率谱估计u 功率谱估计的经典和现代方法u AR模型法的功率谱估计u AR模型法的主要性质u Yule-Walker方程的Levinson Durbin求解算法u 格型滤波器u AR模型参数提取算法u 噪声对噪声对AR谱估计的影响谱估计的影响u ARMA和MA模型法简介u 白噪声中正弦波频率的估计74问题的提出问题的提出：uAR谱估计对观测噪声比较敏感：n噪声会使谱

58、峰展宽，分辨率下降。n噪声会使谱峰偏离正确位置。u在信噪比低的情况下，AR谱估计不会优于周期图方法。u研究怎样减小噪声对AR谱估计的恶化影响？假设x(n)是一个AR(p)过程， 为与x(n)不相关且方差为 的白噪声， x(n)与 混合相加后得到y(n)，即：( )n2( )n ( )( )( )y nx nn于是，2221211( ) ()( )( ) ()( ) ()xxyA z A zSzA z A zA z A z噪声使噪声使AR(p)过程变为过程变为ARMA(p, p)过程！过程！75减小噪声对减小噪声对AR谱估计影响的方法：谱估计影响的方法： 使用与实际序列相符的ARMA(p, p)

59、模型进行谱估计方法； 对数据进行滤波，使用维纳滤波器进行波形估计以减小噪声； 补偿自相关函数或反射系数估计中噪声的影响：2( )( )( )yxR kR kk 采用髙阶AR模型：ARMA(p,p)模型可用AR()模型描述；-1( )( )1( )kkkA zC zc zB z 2-1-1( )/( ) ()/()( )yA zB zSzA zB z22-122-12-1-1( ) ()( ) ( ) ()( ) ()( ) ()xyxA z A zSzA z A zB z B zA z A z令：76AR(p)模型谱估计总结：模型谱估计总结：AR谱估计由于对随机采样数据进行了外推，因此得到了比

60、经典谱估计方法高的分辨率，同时能够较好地调和偏差和方差之间的矛盾。几种常用的现代谱估计方法，例如：最大熵谱估计、线性预测谱估计和AR谱估计等，本质上都是线性外推的思路，因此是等效的。AR模型可以用三种方法等效地表达：自相关函数值、AR模型参数以及反射系数。在自相关法、协方差法和Burg法三种AR模型参数提取方法中，自相关法的性能稍差些，其他两种方法的性能相近。AR谱估计方法对噪声比较敏感，可以从四个方面减小噪声对AR谱估计结果的影响。77N个样值个样值x(0), x(1)x(N-1)功率谱密度 自相关函数 相关图周期图AR模型p+1个参量自相关函数的值Y-W法协方差法修正协方差法反射系数Bur