1、第6章.功率谱的估计6.1 经典法 101NXiimxN21201NXiXixmN 101,0,1,10,NmiimXix xmNRmNmN两种经典谱估计方法 1.直接法（周期图法）2.间接法（BT法）21XNGXN Nj kXXkNGRk e 101,0,1,1NkXiikiRkx xkNN21()lim(,)2XTTGXT 频率分辨率 10,1,10nNd nothers *NXXD假设两个正弦信号之和 频谱泄漏 11,Mj kXXkMGRkk eMN Nj kXXkNGRk e 101,0,1,1NkXiikiRkx xkNN 1,0,mmM MNkMothers经典谱估计的改进窗函数法
2、 平均法 由概率论可知，对L个具有相同的均值和方差的独立随机变量，新随机变量的均值不变，方差减小了L倍。解决矛盾的方法：数据交叠平均法的矛盾 要减小方差，需要增加段数L。每一段的数据M不能太少，否则谱峰将展宽，偏倚变大，从而分辨率会变差。所以段数L不能太大，即方差减少不多。Welch法01002003004005006007008009001000-30-20-100102030Frequency(Hz)Power Spectrum(dB)Periodogram N=25601002003004005006007008009001000-30-20-10010203040Frequency(H
3、z)Power Spectrum(dB)Periodogram N=102401002003004005006007008009001000-10-50510152025Frequency(Hz)Power Spectrum(dB)Averaged Periodogram(no overlap)N=4*25601002003004005006007008009001000-50510152025Frequency(Hz)Power Spectrum(dB)Averaged Periodogram(half overlap)N=1024经典谱估计的说明 经典谱估计，都用FFT快速计算 谱的分辨率
4、较低 由于不可避免有窗函数的影响，使得谱在窗口主瓣内的功率向边瓣部分“泄漏”方差性能不好，不是一致估计，且N增大时谱曲线起伏加剧 周期图的平均和窗函数的使用紧密相关。平均的目的是改善方差性能，但往往会减小分辨率。谱估计的实际问题 数据采样率 每段数据的长度 数据总长度 数据预处理：滤除直流分量和周期分量(市电干扰)2scffsfMf自相关和谱估计的应用 检测混有周期性确定信号的随机信号(48g)相关测速(fla)故障诊断(48h1-1)各阶固有频率的识别（48a)振型分析（48b）机械系统和基础振动传递特性的分析(48c)结构与设备的振动监测与故障诊断 查找电机噪声源 监视机器的工作状态或作故
5、障诊断 查找各种振动源和噪声源6.2 6.2 参数模型功率谱估计参数模型功率谱估计三种参数模型三种参数模型ARAR模型参数的估计模型参数的估计参数模型功率谱估计参数模型功率谱估计 假定所研究的过程x(n)是由一个输入序列w(n)激励一个线性系统H(z)的输出 由已知的x(n)，或其自相关函数Rx(m)估计H(z)的参数 由H(z)的参数来估计x(n)的功率谱 图1随机信号的参数模型)(nw)(nx)(zH)(zW)(zX)(nh 10pqkkkkx na x nkb w nk 0kx nh k w nk B zH zA z 1001pkkkqkkkkkA za zB zb zH zh k z
6、22222jjjxjB eP eH eA e两边取z变换w(n)为白噪声方差为21 1三种参数模型三种参数模型1.1 MAMA（滑动平均）模型（滑动平均）模型 随机信号 由当前的激励 和若干次过去的激励 线性组合产生：(7-1)该模型的系统函数是:(7-2)表示系统阶数，系统函数只有零点，没有极点，所以该系统一定是稳定的系统，也称为全零点模型，用MA（）来表示。)(nx)(nw)(knwqkkknwbnx0)()(0()qkkkH zb zqq1.2 AR（自回归模型（自回归模型)随机信号 由本身的若干次过去值 和当前的激励值 线性组合产生：(7-3)该模型的系统函数是：(7-4)是系统阶数，
7、系统函数中只有极点，无零点，也称为全极点模型，系统由于极点的原因，要考虑到系统的稳定性，因而要注意极点的分布位置，用AR（）来表示。)(nx)(knx)(nwpkkknxanwnx1)()()(pkkkzazH111)(pp1.3 ARMA(自回归滑动平均）模型自回归滑动平均）模型 ARMA是AR与MA模型的结合：(7-5)该模型的系统函数是：(7-6)它既有零点又有极点，所以也称极零点模型，要考虑极零点的分布位置，保证系统的稳定，用ARMR（，）表示。pkkqkkknxaknwbnx10)()()(pkkkqkkkzazbzH101)(pq2.1 AR模型参数和自相关函数的关系模型参数和自相
8、关函数的关系根据式(7-3)：对该式两边同时乘以 ，然后求均值：(7-7)pkkknxanwnx1)()()()(mnx)()()()()()(1pkkmnxknxamnxnwEmnxnxE2 AR模型参数的估计模型参数的估计 因为自相关函数：所以自相关函数呈现偶对称，(7-7)式化为：(7-8)系统的单位脉冲响应 是因果的，所以输出的平稳随机信号和输入的白噪声之间的互相关函数有下列推导：(7-9)(7-10)()()xxR mRm1()()()pxxwkxkR mRma R mk)(nh)()()(mnwnxEmRxw0)()()()()(kknwkhnhnwnx (7-11)所以 (7-1
9、2)带入式(7-8)得到：(7-13)()()()()()()(00mnwknwEkhmnwknwkhEmRkkxw2200()()()()()wwwkkh k Rmkh kmkhm 00)(0)(2mmmhmRwxw121()0()()()0pkxkxpkxwka RmkmRma Rmkhmm由于 由z变换的定义：因而显然，AR模型输出信号的自相关函数具有递推的性质，即：(7-14)上式就是著名的YuleWalker（Y-W）方程，将上式变换：(7-15)pkkkzazH111)(1)0(h1()()0pxkxkR ma R mkm 10()()0pxkxkR ma R mkm lim0zH
10、 zh 从（7-13）求得输入的白噪声方差为：(7-16)将(7-15)和(7-16)结合，把该式的下标简化并写成矩阵的形式，可以写成单一的正规矩阵方程：(7-17)21(0)()0pwxkxkRa Rkm001)0()1()()1()0()1()()1()0(21wpaaRpRpRpRRRpRRR 【例7-1】已知自回归信号模型AR(3)为：式中 是具有方差 =1的平稳白噪声，求 a.自相关序列 ，m0，1，2，3，4，5。b.用a求出的自相关序列来估计AR(3)的参数 ，以及输入白噪声的方差 大小。c.利用给出的AR模型，用计算机仿真给出32点观测值 ,用观测值的自相关序列直接来估计AR(
11、3)的参数 以及输入白噪声的 。)()3(241)2(249)1(2414)(nwnxnxnxnx)(nw2w()xR m ka 2w)(nx ka 2w 32点观测值 0.4282 1.1454 1.5597 1.8994 1.6854 2.3075 2.4679 1.9790 1.6063 1.2804 -0.2083 0.0577 0.0206 0.3572 1.6572 0.7488 1.6666 1.9830 2.6914 1.2521 1.8691 1.6855 0.6242 0.1763 1.3490 0.6955 1.2941 1.0475 0.4319 0.0312 0.58
12、02 -0.6177 解：a.已知的是模型参数 ，14/24 9/24，1/24，来求自相关序列 。)(nx ka1a2a3a()xR m 利用式（7-17），把 代入，利用自相关函数的偶对称，得到一个44的的矩阵：001)0()1()()1()0()1()()1()0(21wpaaRpRpRpRRRpRRR ka000124/124/924/141)0()1()2()3()1()0()1()2()2()1()0()1()3()2()1()0(RRRRRRRRRRRRRRRR解线性方程组得：R(0)=4.9377 R(1)4.3287 R(2)4.1964 R(3)3.8654利用式（7-14
13、）可以求出R(4)，R(5)3.6481，3.4027当然还可以求出无穷多的自相关序列值。1()()0pxkxkR ma R mkm 31(4)(4)xkxkRa Rk 31(5)(5)xkxkRa Rk b.已知自相关序列值，来估计3阶AR模型的参数 以及 利用式（7-17）得到矩阵：(7-18)解线性方程组得到：14/24，9/24，1/24，1 ka 2w0001)0()1()2()3()1()0()1()2()2()1()0()1()3()2()1()0(2321waaaRRRRRRRRRRRRRRRR1 a2 a3 a2wc.利用给出的32点观测值，先求自相关序列（按照上节的样本自相
14、关定义 计算）由于偶对称只给出m0，1，231的 1.9271 1.6618 1.5381 1.3545 1.1349 0.9060 0.8673 0.7520 0.7637 0.8058 0.8497 0.8761 0.9608 0.8859 0.7868 0.7445 0.6830 0.5808 0.5622 0.5134 0.4301 0.3998 0.3050 0.2550 0.1997 0.1282 0.0637 0.0329 -0.0015 -0.0089 -0.0143 -0.008311()nxii miR mx xn()xR m 把头4个相关序列值代入矩阵（7-18）求得估计
15、值：0.6984 ,0.2748 0.0915，0.4678 与真实AR模型参数误差为：0.1151，0.1002，0.0498，1 a2 a3 a2w1e2e3e)(nw)(nx)(zH)(zW)(zX)(nh 10pqkkkkx na x nkb w nk 1001pkkkqkkkkkA za zB zb zH zh k z 22222jjjxjB eP eH eA e差分方程功率谱pkkknxanwnx1)()()(pkkkzazH111)(211(0)(1)()(1)(0)(1)0()(1)(0)0wpRRR paRRR paR pR pRAR模型（全极点模型）转移函数AR模型的正则方
16、程2.2 YW方程的解法方程的解法L-D算法算法mkmknxkanx1)()()()(ka式中 ，k1，2，m，代表m阶预测器的预测系数线性预测(7-19)若序列的模型已知而用过去观测的数据来推求现在和将来的数据称为前向预测器已知x(n)在n时刻之前的m个数据x(n-m),x(n-m+1),x(n-1),预测n时刻的值x(n)把 看成是系统的输出，看成是系统的输入，得到系统函数：)(nemkmknxkanxnxnxne1)()()()()()()(ne)(nxmkkmzkazXzE1)(1)()(预测出来的结果与真实的结果存在预测误差(7-20)假如mp，且预测系数和AR模型参数相同(7-21
17、)即有 ，即前向预测误差系统中的输入为 ，输出为预测误差 等于白噪声。前向预测误差系统对观测信号起了白化的作用。AR模型和前向预测误差系统的系统函数互为倒数，所以求AR模型参数就可以通过求预测误差系统的预测系数来实现。)(nw)(nx)(zH)(zW)(zX)(nh)(ne)(nx)(1zH)(zE)(zX)()(nenw)(nx)(ne图2预测误差系统和AR模型)()()()(212mkmknxkanxEneE111(0)2()()()()()mmmxmxmmxkklRak R kal ak R lk1()()()1,2,mxmxkR la k R l klm预测误差均方值最小均方估计的正交
18、原理：欲选取最佳的系数向量，应是数据向量和误差正交欲选取最佳的系数向量，应是数据向量和误差正交(7-22)(7-23)0,1,2,E x nl e nlm21()(0)()()mmxmxkEenRak R k21()(0)()ppxkxkE e nRa R kpmkaamk)()(ne)(nw22)(wpneE最小均方误差或(7-24)(7-25)(7-26)也就是p阶预测器的预测系数等于p阶AR模型的参数，由于 ，所以最小均方预测误差等于白噪声方差，即 E x n e n 特点：每一阶次参数的计算是从低一阶次的模型参数推算出来的，既可减少工作量又便于寻找最佳的阶数值，满足精度时就停止递推。2
19、1()(0)()()mmxmxkEe nRak R k2,wkpa 2111,wa 22221,2,waa 2(1),ppwpaap估计AR模型参数，也即要估计参数L-D算法的基本思想根据L-D递推算法是模型阶数逐渐加大的一种算法1()()0pxkxkR ma R mkm 先计算阶次m1时的预测系数然后计算m2时的预测系数一直计算到mp阶时的)2()1()1()1()0()0()1()1(111RaRERaR)1(1)0(21121aREw)0()2()1()1()2()1()2()0()1()1(2222RaRaRRaRaR2212121212(1)(2)(1)(1)(1)(2)(1)(0)
20、(2)(1)(1)(2)(1)(2)(1)(1)(2)(0)(0)RaRaaaaRRaRRaaaRaRRm1m=2：把（7-27）的（1）代入上式得到：(7-27)(7-28)1()()()1,2,mxmxkR lak R lklm 21()(0)()ppxkxkE e nRa R k)2()1()1()2()1()1()0()1()1()2()2(11112ERaRRaRRaRa)2()2()1()1()0(222RaRaRE)0()2()1()1()2()1()1()0(2222RaRaaRaR)0()2()1()2(1)1()0(2222RaRaaR)0()2()1()2(1)2(1)1
21、()0(22221RaRaaaR122122)2(1)1()1()0()2(1 EaRaRa根据(7-25)，估计的方差为：把(7-28)的（1）代人：(7-29)其中 称为反射系数，从上式知道整个迭代过程需要已知自相关函数，给定初始值 ，以及AR模型的阶数 ，就可以按照下图所示流程图进行估计。)0(0RE 1)0(0ap()mam这样递推下去可得到预测系数和均方误差估计的通式：(7-30)11111122211()()()()(1)()()()()(2)1()(0)1()(3)mmmmmmkmmmmwmmmkkakakam amkR mak R mkamEEamERak 开始10),(Nnn
22、x输入()xRm计 算000)01xmERa初 始 化(（）,1mmmamE计 算，1.1)(mkkampm(),1.pakkp结束是否1 mm【例7-3】已知自回归信号模型AR（3）为：式中 是具有方差 =1的平稳白噪声，利用给出的AR模型，用计算机仿真给出32点观测值 0.4282 1.1454 1.5597 1.8994 1.6854 2.3075 2.4679 1.9790 1.6063 1.2804 -0.2083 0.0577 0.0206 0.3572 1.6572 0.7488 1.6666 1.9830 2.6914 1.2521 1.8691 1.6855 0.6242 0
23、.1763 1.3490 0.6955 1.2941 1.0475 0.4319 0.0312 0.5802 -0.6177，用L-D算法来估计AR(3)的参数 以及输入白噪声的方差 )()3(241)2(249)1(2414)(nwnxnxnxnx)(nw2w)(nx ka 2w 解：步骤1 利用给出的32点观测值，先求自相关序列，由于偶对称，只给出m0，1，231的 1.9271 1.6618 1.5381 1.3545 1.1349 0.9060 0.8673 0.7520 0.7637 0.8058 0.8497 0.8761 0.9608 0.8859 0.7868 0.7445 0
24、.6830 0.5808 0.5622 0.5134 0.4301 0.3998 0.3050 0.2550 0.1997 0.1282 0.0637 0.0329 -0.0015 -0.0089 -0.0143 -0.0083)(mRxx0(0)1.9271xxER10a10211(1)1.6618(1)0.86231.9271(0)1(1)0.4942RaEERa 步骤2，初始化：步骤3，根据式（7-30）计算，m1：1212122212(2)(1)(1)(2)0.2127(1)(1)1(2)0.67891(2)0.4718RaRaEaaaEEa2232323232322323(3)(1)
25、(2)(2)(1)(3)0.0914(1)(1)(3)(2)0.6983(2)(2)(3)(1)0.27481(3)0.4679RaRaRaEaaaaaaaaEEa m2：m3：3p6983.01a2748.02a0914.03a4679.032 Ew因而当 时，估计到的AR模型参数为：估计的输入信号的方差为3 以参数模型做谱估计下图中随机信号x的谱估计为：2j2wjx)e(H)e(P)(nw)(nx)(zH)(zW)(zX)(nhp1k2kj)p(k2wjxp1kk)p(k|ea1|)e(Pza11)z(HAR代入建模为例以以参数模型做谱估计的特点 AR谱的平滑特性00.10.20.30.4
26、0.50.60.70.80.9-34-32-30-28-26-24-22-20-18-16-14Frequency(kHz)Power/frequency(dB/Hz)Power Spectral Density Estimate via Modified Covariance050100150200250300350400450500-70-60-50-40-30-20-100Frequency(Hz)Power Spectral Density(dB/Hz)Periodogram PSD Estimate AR谱的分辨率0102030405060708090100-70-60-50-40-
27、30-20-100Frequency(Hz)Power Spectral Density(dB/Hz)Periodogram PSD Estimate0102030405060708090100-50-40-30-20-1001020Frequency(Hz)Power/frequency(dB/Hz)Power Spectral Density Estimate via Modified Covariance4 信号建模的本质信号建模的思路：假定所研究的过程x(n)是由一个输入序列w(n)激励一个线性系统H(z)的输出问题？任意地给定一个平稳过程x(n)，均可由一个白噪序列w(n)激励一个线
28、性系统H(z)来精确地产生001)0()1()()1()0()1()()1()0(21wpaaRpRpRpRRRpRRR对信号x(n)建立参数模型，并不是要求模型的输出x(n)在时域等于x(n),而是要求它们在某一阶次上的统计特性相同“准确建模”的定义：设平稳随机过程x(n)存在 r 阶模型，使得模型的输出x(n)在 r 阶统计特性上和x(n)的同价统计特性相一致，则把x(n)称为在 r 阶统计意义上可准确建模的随机过程5 关于线性预测的进一步讨论mkmknxkanx1)()()(12mffmkffffxnak x nkenx nxnEE en 前向预测x(n-m),x(n-m+1),x(m-
29、1)min1(0)()()mffxmxkERak R k1()()()1,2,mfxmxkR lak R lklm min1(0)()()mbbxmxkERak R k 12mbbmkbbbffxnmak x nmkenmenx nmxnmEE en 1()()()1,2,mbxmxkR lak R lklm 1mbbmkxnak x nk 后向预测x(n+1),x(n+2),x(n+m)2fbminminEE bfmmakak由此可得 前、后向预测误差和反射系数在不同阶次的递推关系：111111ffbmmm mbb*fmmm menenk enenenk en mmkam反射系数 11111
30、1bfmmmfbmmcov en,enkvar envar en 前后预测误差之间的相关系数6 AR模型系数的求解 n的取值范围 在令预测误差功率最小时，是单独使前向预测误差功率最小，还是使前后向预测误差功率之和最小 是先估计自相关函数再求解还是直接由数据递推求解。6.1自相关法 由自相关函数，利用L-D递推算法所得到的AR模型的参数等效于前向预测器的系数，AR模型白噪声的方差等效于前向预测的最小预测误差功率。6.2 Burg算法 令前后向预测误差功率之和最小12fbfbEEE 预测误差功率的求和范围 121211Nffppn pNbbppn pEenNpEenNp 111111ffbmmm
31、mbb*fmmm menenk enenenk en0fbmEk 111112211211Nfb*mmn mmNNfbmmn mn men enkenen 1121()()()()1mmmmmmmmmakakk amkamkEkE令 ，可得 Burg算法的步骤 初始条件 计算 m=1时参数 计算 m=2时参数 重复上述过程，直到m=p 002201fbNneex nx nN1k11a,E211fbe,ek和 22212a,a,E 习习 题题1.对于一个随机信号，可以对它进行频谱分析，叙述AR谱法和周期图法相比的优点。2.设已知 用L-D算法为此信号估计p1，2，3阶AR模型的系数和激励白噪的功率。2,1,0),4.0(113)8.0(1114)(mmRmms
