医学信号处理现代谱估计应用

第七章功率谱估计的现代方法 现代谱估计 经典谱估计以傅立叶变换为基础 具有计算效率高的优点 但是由于将未观测数据认为0和数据加窗 而具有频率分辨率低 旁瓣泄漏等严重的缺陷 现代谱估计与经典谱估计不同 它以参数模型为基础 能够得到小方差和高分辨率 特别是数据长度很短的情况 更具优势 7 1概述 现代谱估计法的基本思想 处理步骤 1确定或选择一个合适的模型 依赖于对所研究随机过程进行理论分析和实验研究 2根据观测数据估计模型参数 涉及各种算法的研究 3由模型参数计算功率谱 关键1 模型选择问题 AR MA ARMA 2 参数确定方法 导致产生了各种算法 7 2自回归模型 AR 谱估计 数字系统的数学模型 有理分式传递函数的模型如下图 式中ak为自回归系数 称为AR系数 bk为滑动平均系数 称为MA系数 模型传递函数为 有理分式传递函数的模型的差分方程为 令a0 1有 w n 为高斯白噪声 求功率谱的实质变为确定系统参数的问题 模型的功率谱密度 即系统输出功率谱和输入功率谱之间关系为 假定h n 为实序列 如果除b0外其它的MA系数都等于0 即 AR p 模型 全极点模型 自回归模型 如果除a0外其它的AR系数都等于0 即 全零点模型 MA q 模型 如果除a0 1和b0 1外其它的AR系数和MA系数都不全等于0 即 称为ARMA p q 模型 即极点 零点模型 到底选择什么模型 三种模型之间关系如何 Word分解定理 Wold分解定理 任何一个有限方差的平稳ARMA过程可以分为完全随机的部分和确定的部分 推论 任何有限方差的ARMA或MA平稳过程可以用一个无限阶的AR模型表示 同样 任何ARMA或AR模型可以用一个无限阶的MA模型表示 因此 如果在这三个模型中选了一个与信号不匹配的模型 利用高的阶数仍然可以得到对信号的好的逼近 结论 由于对AR模型参数的估计 得到的是线性方程 故AR模型比ARMA以及MA模型有计算上的优点 即只需解一组线性方程 而ARMA或MA模型一般需要解一组非线性方程 同时 实际的物理系统往往是全极点系统 AR模型得到了深入的研究和广泛的应用 已知 自相关函数 要求 AR模型的阶数p 以及p个AR参数a i 激励源方差 Yule Walker方程 7 3AR模型的Yule Walker方程 7 3 1Yule Walker方程的推导 1 对进行求逆z变换2 直接由模型差分方程推导 把模型的差分方程代入x n 的自相关函数 如何根据自相关函数确定系统参数 可见 AR模型输出信号的自相关函数具有递推性质 即 Yule Walker方程 Y W方程 选择m 0的前P个方程并写成单一正规矩阵的形式为 以上利用了自相关函数的偶对称性 Y W方程表明 只要已知输出平稳随机信号的自相关函数 就能求出AR模型中的参数 ak 并且需要的观测数据较少 AR模型谱估计 N个样值x 0 x 1 x N 自相关函数R 0 R 1 R N AR模型参数和a1 a2 ap激励源方差 功率谱密度 Y W方程 Yule Walker方程的求解 1 采用高斯消元法 解线性方程组常用方法 运算量数量级为p的三次方 2 用Levinson Durbin算法 Y W方程的高效解法 即按阶次进行递推运算量数量级为p的二次方 7 3 2Levinson Durbin算法 Levinson Durbin递推算法 算法的关键就是要推导出由第K阶AR模型的参数计算第k 1阶AR模型AR k 1 参数的迭代计算公式 首先以AR 0 和AR 1 模型参数作为初始条件 计算AR 2 模型参数 然后根据这些参数计算AR 3 模型参数 等等 一直到计算出AR p 模型参数为止 递推公式为 其中akk称为反射系数 将所估计的模型参数代入即可计算功率谱估计值 AR模型参数和a1 a2 ap激励源方差 功率谱密度 AR模型谱估计 给定初始值和AR模型的阶数p 可按照L D算法流程进行估计 流程终止规则为或 MATLAB里有专门实现L D算法的函数可估计AR模型参数 aE aryule x p a为模型参数 E为噪声方差 分析 AR模型的稳定性 L D算法的收敛性 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或者 代入计算公式估计功率谱 例7 1 已知实数据序列 的自相关为 用Levinson Durbin递推算法求AR模型的参量 解 一 AR模型的稳定性具有下面性质 H z 全部的极点在单位圆内自相关矩阵正定激励信号方差随阶次增加而递减 7 3 3确定AR模型的阶 阶太低 功率谱平滑的太厉害 平滑后的谱分辨不出真实谱中的两个峰 阶太高 可以提高谱估计的分辨率 但会出现许多虚假谱峰 真实谱 虚假谱峰 二 有关AR模型的阶的问题 所以 估计一个AR p 过程 选取AR k 阶数要求k p 但k不能太大 如果估计精确的话 k p时 AR p 模型参数估计为 AR模型谱估计方法 既要估计模型参数 又要估计模型的阶 在这样复杂的情况下 如何评价各种谱估计的性能 目前尚无定论 三 确定AR模型的阶的方法 一般的观察方法 简单而直观 不断增加阶数 观察预测误差功率 下降到最小时 对应的阶选为模型的阶 不断增加阶数 观察各阶模型预测误差序列的周期图 最接近平坦时对应于最佳的阶 1 FPE 最终预测误差 N为观测数据长度 为拟合残差方差 随阶增加而减小 FPE的最小值对应的阶数为最后确定的阶 四 确定AR模型的阶的方法 根据误差准则确定 2 Akaike AIC 信息准则 AIC i 最小所对应的阶 i为模型的阶 为模型误差 随着阶的增加而减小 而式中第二项随阶次增加而增加 AIC定义式有一个最小值 适用于AR模型 此外 还有CAT等准则 通过实验发现 在将这些准则用于估计AR模型的阶 对于实际数据 所得到的谱估计结果常常无太大区别 对于短数据 以上准则都不理想 在实际应用中 应该参照实验结果对模型的阶加以适当调整 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 预测误差为 预测误差功率为 确定系数ak的一个原则是使预测误差功率最小 根据这一原则推导出的预测器系数 ak与x n 的自相关序列Rxx m 之间的关系为 将两个关系式写成矩阵展开式分别为 将 1 和 2 两个关系式合并为一个式子 将 3 写成矩阵展开形式为 可以看出 N阶线性预测器的系数ak与AR模型中的AR系数相等 预测误差概率最小值Pmin与AR模型中的输入噪声方差相等 所以 线性预测谱估计与AR谱估计是等效的 熵 是信息量的一种量度 也是不确定性的一种量度 信息量与事件发生概率之间有类似于反比例的关系 信息量与概率之间存在对数关系 复合事件的信息量等于各独立事件信息量之和 对于事件A有 7 5最大熵谱估计 MESE MaximumEntropySpectralEstimation 7 5 1按最大熵谱外推自相关函数 N个符号组成信号系统传递消息 每个符号出现的概率为pi 接收到第i个符号的信息量为I i 消息中总的平均信息量为 这个平均信息量称为具有符号i和概率pi的信号系统的熵 对于随机过程 应该用联合概率密度函数来定义熵 x0 x1 xN 为随机过程的N 1个取样值 对于零均值的高斯平稳随机过程则有 其中X x0 x1 xN 为由N 1个取样值构成行矩阵 XH为X的共轭转置矩阵 det R N 是行列式的值 于是有下列式子 均值为0的高斯平稳随机过程的熵的表达式 它是R N 的函数 最大熵谱估计 为了使得H取得最大值 应当使det R N 取最大值 根据外推或预测方法 求出使det R N 1 取最大值的Rxx N 1 由 得到 结论 上述方程为Rxx N 1 的一元一次方程 可由已知或估计的N 1个自相关值Rxx 0 Rxx 1 Rxx N 求出Rxx N 1 以此类推 可以证明这种按最大熵外推自相关函数的结果与AR模型是等价的 所以 上式实质为Yuler walker方程 实际上 假定在线性预测谱估计中 用外推的方法得到了Rxx N 1 即 7 5 2MESE与AR谱估计等效 使联立方程有非零解的充分必要条件是系数行列式等于0 结论因此 最大熵谱估计与AR模型谱估计和线性预测谱估计是等效的 此外 还可以证明 AR谱估计等效于最佳白化处理 7 6 1预测误差格型滤波器 7 6预测误差格型滤波器及伯格 Burg 递推算法 已知n个观测数据x 1 x 2 x n 1 利用p阶线性预测滤波器估计x n 为估计误差为 代入有 上式中 因此可得 格型前向预测误差滤波器 格型预测误差滤波器传递函数为 同时有 也就是相当于 后向预测误差 格型后向预测误差滤波器传递函数为 又因 因此 7 6 2Burg递推算法 Kp的确定 根据信号的有限个取样值估计AR模型参数的方法 自相关法 协方差法和Burg递推法 自相关法和协方差法都是直接估计AR参数 而Burg法是先估计反射系数 然后利用Levinson Durbin递推算法由反射系数求得AR参数 Burg法可以保证 算法准则是前向均方误差和后向均方误差之和最小 如果用前向预测方法以均方误差最小为准则确定Kp 用表示为 令 如果用后向预测方法以均方误差最小为准则确定Kp 用表示为 令 Burg算法准则是前向均方误差和后向均方误差之和最小 令 于是得到 对于平稳随机序列 集合平均可用时间均值来代替 故上式为 通过观察发现 下式成立 对于长度为N的有限长序列 有下面一系列递推关系式 当p 1时 于是有 于是 其中 继续代入求得K3 b3 n e3 n 以此类推 便可以由x n 求得各阶Kp以及前向与后向误差及其各个akk Burg法估计AR p 模型参数的具体步骤为 1 确定初始条件 2 按照公式计算Kp 3 按照公式计算ep n 和bp n 4 计算均方误差 5 p p 1 6 重复第2 5步 直至满足条件为止 例5 2 设N 5的数据记录为x 0 1 x 1 2 x 2 3 x 3 4 x 4 5 AR模型的阶次p 3 试用相关函数法确定AR参量及预测值 解 先由数据求自相关函数式 用Levinson Durbin递推算法求AR模型的参量分别是 根据所得的AR 3 参量 预测值 若使用的是二阶线性预测器 有例5 1所得的结果 则 可分别由前向与后向预测得到如下 例5 3 设仍利用例5 2中的记录数据 试用伯格法求AR 2 的参量 解 用上述递推公式 i 1时 e1 n 和b1 n p 2时 若使用此二阶线性预测 可得 算法比较Levinson DurbinBurg算法真正样值 25 90903 0650 25 54030 17415 1 27000 85491 0000 2 89834 58255 0000 显然 伯格算法要比莱文森 德宾算法优越得多 短数据 比较Welch方法和Burg方法在噪声信号的功率谱估计中的效果 为高斯型白噪声 例子 现代谱估计和经典谱估计方法的比较 利用MATLAB中的pburg和pwelch函数分别用Burg方法和Welch方法对上述噪声信号进行功率谱估计并比较结果 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方法低阶 信号 很明显 Burg方法比Welch方法更光滑 但是 当AR模型阶数降低时 谱峰的频移越来越明显 频率分辨率降低 7 7AR模型谱估计存在的问题 7 7 1谱线分裂 由正弦信号叠加噪声构成的随机信号 在下列四种情况下容易出现谱线分裂的现象 即谱线频率偏移或出现两个靠得很近的谱峰 高信噪比 正弦信号分量的初始相位是 4的奇数倍 数据长度为正弦分量的1 4周期的奇数倍 AR模型参数的数目与数据的个数相比的百分比较大 即二者大小可比拟 对于Burg算法 谱线分裂是由于第一个反射系数K的计算误差引起的 K1的估计并没有使预测误差功率最小 改善措施 1 用解析信号代替实值信号 克服信号相位的影响 2 调整修正反射系数 以使预测误差功率达到最小 7 7 2附加噪声使分辨率下降 A