《7-医学信号处理现代谱估计.ppt》由会员分享，可在线阅读，更多相关《7-医学信号处理现代谱估计.ppt（86页珍藏版）》请在三一办公上搜索。
1、1,第七章功率谱估计的现代方法 现 代 谱 估 计,2,经典谱估计以傅立叶变换为基础，具有计算效率高的优点，但是由于将未观测数据认为0和数据加窗，而具有频率分辨率低、旁瓣泄漏等严重的缺陷。现代谱估计与经典谱估计不同，它以参数模型为基础，能够得到小方差和高分辨率，特别是数据长度很短的情况，更具优势。,7.1 概述,3,现代谱估计法的基本思想：,处理步骤：1 确定或选择一个合适的模型依赖于对所研究随机过程进行理论分析和实验研究；2 根据观测数据估计模型参数涉及各种算法的研究；3 由模型参数计算功率谱。关键 1、模型选择问题（AR，MA，ARMA）2、参数确定方法（导致产生了各种算法）,4,7.2
2、自回归模型(AR)谱估计,数字系统的数学模型:有理分式传递函数的模型如下图：,式中ak为自回归系数，称为AR系数；bk为滑动平均系数，称为MA系数。,模型传递函数为：,5,有理分式传递函数的模型的差分方程为：,令a0=1有：,w(n)为高斯白噪声，,6,求功率谱的实质变为确定系统参数的问题,模型的功率谱密度：即系统输出功率谱和输入功率谱之间关系为（假定h(n)为实序列）：,7,如果除b0外其它的MA系数都等于0，即,AR（p）模型,全极点模型,8,自回归模型,9,如果除a0外其它的AR系数都等于0，即,全零点模型,MA（q）模型,10,如果除a0=1和b0=1外其它的AR系数和MA系数都不全等
3、于0，即,称为ARMA（p,q）模型，即极点-零点模型。,11,到底选择什么模型？三种模型之间关系如何？Word分解定理,Wold分解定理：任何一个有限方差的平稳ARMA过程可以分为完全随机的部分和确定的部分。推论：任何有限方差的ARMA或MA平稳过程可以用一个无限阶的AR模型表示；同样，任何ARMA或AR模型可以用一个无限阶的MA模型表示。因 此，如果在这三个模型中选了一个与信号不匹配的模型，利用高的阶数仍然可以得到对信号的好的逼近。,12,结论：由于对AR模型参数的估计，得到的是线性方程。故AR模型比ARMA以及MA模型有计算上的优点，即只需解一组线性方程，而ARMA或MA模型一般需要解一
4、组非线性方程。同时，实际的物理系统往往是全极点系统。AR模型得到了深入的研究和广泛的应用。,13,已知：自相关函数,要求：AR模型的阶数p，以及p个AR 参数a(i),激励源方差,Yule-Walker方程,7.3 AR模型的Yule-Walker方程,14,7.3.1 Yule-Walker方程的推导,1.对 进行求逆z变换 2.直接由模型差分方程推导，把模型的差分方程代入x(n)的自相关函数,15,如何根据自相关函数确定系统参数,16,可见，AR模型输出信号的自相关函数具有递推性质，即：,Yule-Walker方程（Y-W方程）,17,选择m0的前P个方程并写成单一正规矩阵的形式为：,以上
5、利用了自相关函数的偶对称性。Y-W方程表明：只要已知输出平稳随机信号的自相关函数，就能求出AR模型中的参数ak，并且需要的观测数据较少。,18,AR模型谱估计,N个样值x(0),x(1)x(N),自相关函数R(0),R(1).R(N),AR模型参数和a1,a2，ap激励源方差,功率谱密度,Y-W方程,19,Yule-Walker方程的求解,1、采用高斯消元法，解线性方程组常用方法，运算量数量级为p的三次方。2、用Levinson-Durbin算法，Y-W方程的高效解法，即按阶次进行递推运算量数量级为p的二次方。,7.3.2 Levinson-Durbin算法,20,Levinson-Durbi
6、n递推算法：算法的关键就是要推导出由第K阶AR模型的参数计算第k+1阶AR模型AR(k+1)参数的迭代计算公式。,首先以AR(0)和AR(1)模型参数作为初始条件，计算AR(2)模型参数，然后根据这些参数计算AR(3)模型参数，等等，一直到计算出AR(p)模型参数为止。,21,22,23,递推公式为：,其中akk称为反射系数,将所估计的模型参数代入即可计算功率谱估计值：,24,AR模型参数和a1,a2，ap激励源方差,功率谱密度,AR模型谱估计,25,给定初始值和AR模型的阶数p，可按照L-D算法流程进行估计，流程终止规则为 或,MATLAB里有专门实现L-D算法的函数可估计AR模型参数：a
7、E=aryule(x,p),a为模型参数，E为噪声方差。分析：AR模型的稳定性；L-D算法的收敛性。,26,AR模型谱估计的L-D算法流程,、给定N个观察数据xN(n),n=0,1,N-1;、由xN(n)估计自相关函数值，m=0,1,p;、利用L-D递推算法，根据 计算AR模型参数的估计值。首先令p=1，按下式计算a11和,然后，使p=p+1,按下式计算app,api,、重复以上递推过程，直到满足p=m或者。、代入 计算公式估计功率谱。,27,例7-1、已知实数据序列,的自相关为：,用Levinson-Durbin递推算法求AR模型的参量：,解：,28,29,一、AR模型的稳定性具有下面性质：
8、,H(z)全部的极点在单位圆内自相关矩阵正定激励信号方差随阶次增加而递减,7.3.3 确定AR模型的阶,30,阶太低，功率谱平滑的太厉害，平滑后的谱分辨不出真实谱中的两个峰；阶太高，可以提高谱估计的分辨率，但会出现许多虚假谱峰。,真实谱,虚假谱峰,二、有关AR模型的阶的问题：,31,所以，估计一个AR(p）过程，选取AR(k）阶数要求 kp,但k不能太大。如果估计精确的话，kp时，AR(p)模型参数估计为：,AR模型谱估计方法，既要估计模型参数，又要估计模型的阶，在这样复杂的情况下，如何评价各种谱估计的性能，目前尚无定论。,32,三、确定AR模型的阶的方法 一般的观察方法，简单而直观,不断增加
9、阶数，观察预测误差功率，下降到最小时，对应的阶选为模型的阶；不断增加阶数，观察各阶模型预测误差序列的周期图，最接近平坦时对应于最佳的阶；,33,1、FPE(最终预测误差）,N为观测数据长度，为拟合残差方差，随阶增加而减小，FPE的最小值对应的阶数为最后确定的阶。,四、确定AR模型的阶的方法根据误差准则确定,34,2、Akaike（AIC)信息准则,AIC(i)=最小所对应的阶。,i为模型的阶，为模型误差，随着阶的增加而减小，而式中第二项随阶次增加而增加。AIC定义式有一个最小值。适用于AR模型。,35,此外，还有CAT等准则。,通过实验发现：在将这些准则用于估计AR模型的阶，对于实际数据，所得
10、到的谱估计结果常常无太大区别。对于短数据，以上准则都不理想。在实际应用中，应该参照实验结果对模型的阶加以适当调整。,36,7.4 线性预测谱估计,假设x(n)是一个N阶AR过程，现在时刻x(n)的值 可以由过去N个时刻的取样值的加权来预测，加权系数为-ak，那么,N阶线性预测器：可看作用序列x(n-N),x(n-N-1),x(n-1)激励一个冲击响应为-ak的线性时不变系统的输出值。,x(n-N),x(n-N-1),x(n-1),-ak,预测误差为：,预测误差功率为：,38,确定系数ak的一个原则是使预测误差功率最小。根据这一原则推导出的预测器系数-ak与x(n)的自相关序列Rxx(m)之间的
11、关系为：,将两个关系式写成矩阵展开式分别为：,40,将(1)和(2)两个关系式合并为一个式子：,41,将(3)写成矩阵展开形式为：,可以看出：N阶线性预测器的系数ak与AR模型中的AR系数相等，预测误差概率最小值Pmin与AR模型中的输入噪声方差 相等。所以，线性预测谱估计与AR谱估计是等效的。,42,熵是信息量的一种量度，也是不确定性的一种量度。信息量与事件发生概率之间有类似于反比例的关系，信息量与概率之间存在对数关系。复合事件的信息量等于各独立事件信息量之和。对于事件A有：,7.5 最大熵谱估计（MESE）Maximum Entropy Spectral Estimation）,7.5.1
12、 按最大熵谱外推自相关函数,43,N个符号组成信号系统传递消息，每个符号出现的概率为pi，接收到第i个符号的信息量为I(i)，消息中总的平均信息量为：,这个平均信息量称为具有符号i和概率pi的信号系统的熵。对于随机过程，应该用联合概率密度函数来定义熵：x0,x1,xN为随机过程的N+1个取样值。,44,对于零均值的高斯平稳随机过程则有：,其中X=x0,x1,xN为由N+1个取样值构成行矩阵，XH为X的共轭转置矩阵。,45,detR(N)是行列式的值。于是有下列式子：,46,均值为0的高斯平稳随机过程的熵的表达式，它是R(N)的函数。,47,最大熵谱估计：为了使得H取得最大值，应当使detR(N
13、)取最大值。根据外推或预测方法，求出使detR(N+1)取最大值的Rxx(N+1)：,由,得到：,48,结论：,上述方程为Rxx(N+1)的一元一次方程，可由已知或估计的N+1个自相关值Rxx(0)、Rxx(1)、Rxx(N)求出Rxx(N+1)。以此类推。可以证明这种按最大熵外推自相关函数的结果与AR模型是等价的。所以，上式实质为Yulerwalker方程。,49,实际上，假定在线性预测谱估计中，用外推的方法得到了Rxx(N+1)，即：,7.5.2 MESE与AR谱估计等效,50,使联立方程有非零解的充分必要条件是系数行列式等于0：,结论 因此：最大熵谱估计与AR模型谱估计和线性预测谱估计是
14、等效的。此外，还可以证明：AR谱估计等效于最佳白化处理。,51,7.6.1 预测误差格型滤波器,7.6 预测误差格型滤波器及伯格(Burg)递推算法,已知n个观测数据x(1),x(2),x(n-1)，利用p阶线性预测滤波器估计x(n)为估计误差为：,52,代入有：,53,上式中,因此可得：,54,格型前向预测误差滤波器,55,格型预测误差滤波器传递函数为,同时有：,56,也就是相当于：,57,后向预测误差,58,格型后向预测误差滤波器传递函数为,又因：,因此：,59,7.6.2 Burg递推算法Kp的确定,根据信号的有限个取样值估计AR模型参数的方法自相关法、协方差法和Burg递推法。自相关法
15、和协方差法都是直接估计AR参数，而Burg法是先估计反射系数，然后利用Levinson-Durbin递推算法由反射系数求得AR参数。Burg法可以保证：,60,算法准则是前向均方误差和后向均方误差之和最小。,如果用前向预测方法以均方误差最小为准则确定Kp，用 表示为：,令,61,如果用后向预测方法以均方误差最小为准则确定Kp，用 表示为：,令,62,Burg算法准则是前向均方误差和后向均方误差之和最小，令,63,于是得到：,对于平稳随机序列，集合平均可用时间均值来代替，故上式为：,通过观察发现，下式成立：,64,对于长度为N的有限长序列，有下面一系列递推关系式。当p=1时：,65,于是有：,6
16、6,于是,其中,67,继续代入求得K3，b3(n)，e3(n)，以此类推。便可以由x(n)求得各阶Kp以及前向与后向误差及其各个akk。,68,Burg法估计AR(p)模型参数的具体步骤为：1、确定初始条件：2、按照公式计算Kp。3、按照公式计算ep(n)和bp(n)。4、计算均方误差：5、p=p+1。6、重复第25步，直至满足条件为止。,69,例5-2、设N=5的数据记录为x(0)=1,x(1)=2,x(2)=3,x(3)=4,x(4)=5,AR模型的阶次p=3,试用相关函数法确定AR参量及预测值.,解：先由数据求自相关函数式：,70,用Levinson-Durbin递推算法求AR模型的参量
17、分别是：,71,根据所得的AR(3)参量，预测值：,若使用的是二阶线性预测器，有例51所得的结果，则,可分别由前向与后向预测得到如下：,72,例5-3、设仍利用例52中的记录数据，试用伯格法求AR(2)的参量。解：用上述递推公式，i=1时：,e1(n)和b1(n),73,p=2时：,若使用此二阶线性预测，可得：,74,算法比较 Levinson-Durbin Burg算法 真正样值,25.9090 3.0650,25.5403 0.17415,1.2700 0.8549 1.0000,2.8983 4.5825 5.0000,显然，伯格算法要比莱文森德宾算法优越得多。短数据！,75,比较Wel
18、ch方法和Burg方法在噪声信号的功率谱估计中的效果。为高斯型白噪声。,例子：现代谱估计和经典谱估计方法的比较。,利用MATLAB中的pburg和pwelch函数分别用Burg方法和Welch方法对上述噪声信号进行功率谱估计并比较结果。,76,function=burgwelchpsd()fs=1000;t=0:1/fs:1;xn=sin(2*pi*80*t)+2*sin(2*pi*140*t)+randn(size(t);plot(xn);pw,f=pwelch(xn,fs,twosided);pb1,f=pburg(xn,17,fs,twosided);pb2,f=pburg(xn,13,
19、fs,twosided);figurew=10*log10(pw)10*log10(pb1)10*log10(pb2);plot(f,w);gridxlabel(frequency(Hz);ylabel(amplitude(dB);axis(0 200-50 0);legend(welch方法,Burg方法高阶,Burg方法低阶);,77,信号：,78,很明显，Burg方法比Welch方法更光滑。但是，当AR模型阶数降低时，谱峰的频移越来越明显，频率分辨率降低。,79,7.7 AR模型谱估计存在的问题,7.7.1 谱线分裂,由正弦信号叠加噪声构成的随机信号，在下列四种情况下容易出现谱线分裂的现
20、象，即谱线频率偏移或出现两个靠得很近的谱峰。高信噪比；正弦信号分量的初始相位是/4的奇数倍；数据长度为正弦分量的1/4周期的奇数倍；AR模型参数的数目与数据的个数相比的百分比较大，即二者大小可比拟。,80,对于Burg算法，谱线分裂是由于第一个反射系数K的计算误差引起的，K1的估计并没有使预测误差功率最小。改善措施：1、用解析信号代替实值信号，克服信号相位的影响；2、调整修正反射系数，以使预测误差功率达到最小。,81,7.7.2 附加噪声使分辨率下降,AR谱估计方法对观测噪声比较敏感，从而限制了其应用范围。噪声使谱峰展宽，导致分辨率下降，使谱峰偏离正确位置。原因是：AR谱估计假设的全极点模型在
21、有观测噪声时，不再成立。,82,设x(n)是p阶AR过程，有观测噪声v(n)存在时，成为y(n)，y(n)=x(n)+v(n)。若v(n)为与x(n)不相关的、方差为 的白噪声，则：其功率谱为 分别为w(n)、v(n)的方差。,83,令且则有于是看到，由于噪声的存在，使得AR模型变成一个ARMA过程。,84,减小噪声对AR谱估计影响的措施：1、补偿自相关函数或反射系数估计中噪声的影响。对于Burg算法，检查反射系数是否小于1。2、对数据进行滤波减小噪声。3、采用ARMA谱估计方法。4、采用高阶AR模型。AR谱估计分辨率为：分辨率随阶的增加而增加，但是考虑到虚假谱峰的问题，模型阶数最高不应该超过数据点数的一半。,85,作业1、已知序列x(n)=4.684,7.247,8.423,8.650,8.640,8.392是由模型 x(n)=1.70 x(n-1)-0.72x(n-2)+u(n)产生的，这里u(n)是均值为0、方差为1的白噪声。试用Burg算法求AR(2)的模型参数,并画出二阶格型预测误差滤波器结构。2、已知某自回归过程的五个观测值为1，1，1，1，1利用L-D算法求一阶和二阶反射系数；求该自回归过程的功率谱估计。,86,实验2,例题中增加周期图法和Bartlett法谱估计，然后将5种方法的功率谱估计画在图上，并分析比较。,
