医学信号处理现代谱估计

《医学信号处理现代谱估计》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

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