数字信号处理-时域离散随机信号处理(丁玉美)第4章.ppt

第四章功率谱估计 4 1引言4 2经典谱估计4 3现代谱估计中的参数建模4 4AR模型谱估计的性质4 5AR谱估计的方法4 6最大熵谱估计与最大似然谱估计4 7特征分解法谱估计 4 1引言 我们知道 对信号和系统进行分析研究 处理有两类方法 一类是在时域进行 前面我们学习的维纳 卡尔曼滤波和自适应滤波都属于这种方法 本章则是在频率域进行研究的另一类方法 这两类方法都是信号处理的重要方法 对确定性信号傅里叶变换是在频率域分析研究的理论基础 但对于随机信号 其傅里叶变换并不存在 因此转向研究它的功率谱 按照Weiner Khintchine定理 信号的功率谱和其自相关函数服从一对傅里叶变换关系 公式如下 4 1 1 4 1 2 4 1 3 4 1 1 式被称做功率谱的定义 对于平稳随机信号 服从各态历经定理 集合平均可以用时间平均代替 由 4 1 1 式还可以推出功率谱的另一个定义 推导如下 将 4 1 3 式中的集合平均用时间平均代替 得到 4 1 4 将 4 1 4 式代入 4 1 1 式 得到 令l n m 则 4 1 5 上式中x n 是观测数据 Pxx ej 是随机变量 必须对Pxx ej 取统计平均值 得到 4 1 6 上式被认为是功率谱的另一定义 4 1 1 式表明功率谱是无限多个自相关函数的函数 但观测数据只有有限个 只能得到有限个自相关函数 按照 4 1 6 式求功率谱 也需要无限个观测数据 因此根据有限个样本数据 分析计算随机序列的真正功率谱 是求功率谱的中心问题 毫无疑问 这是一个功率谱的估计问题 在第一章已介绍了统计估计的一般估计准则 主要有偏移 估计量方差和估计量的均方误差 有效性 这里不再重复 下面直接用它们分析估计质量 现代谱估计以信号模型为基础 图4 1 1表示的是x n 的信号模型 输入白噪声w n 均值为0 方差为 2w x n 的功率谱由下式计算 4 1 7 如果由观测数据能够估计出信号模型的参数 信号的功率谱可以按照 4 1 7 式计算出来 这样 估计功率谱的问题变成了由观测数据估计信号模型参数的问题 模型有很多种类 例如AR模型 MA 模型等等 针对不同的情况 合适地选择模型 功率谱估计质量比较经典谱估计的估计质量有很大的提高 遗憾的是 尚无任何理论能指导我们选择一个合适的模型 我们只能根据功率谱的一些先验知识 或者说一些重要的谱特性 选择模型 图4 1 1平稳随机序列的信号模型 4 2经典谱估计 4 2 1BT法 BT法是先估计自相关函数 然后按照 4 1 1 式进行傅里叶变换得到功率谱 设对随机信号x n 只观测到一段样本数据 n 0 1 2 N 1 关于如何根据这一段样本数据估计自相关函数 第一章已经作了详细介绍 结果是共有两种估计方法 即有偏自相关函数估计和无偏自相关函数估计 有偏自相关函数估计的误差相对较小 这种估计是一种渐近一致估计 将该估计公式重写如下 4 2 1 对上式进行傅里叶变换 得到BT法的功率估计值为 4 2 2 为了减少谱估计的方差 经常用窗函数w m 对自相关函数进行加权 此时谱估计公式为 4 2 3 式中 M 1 m M 1 其它 M N 4 2 4 有时称 4 2 3 式为加权协方差谱估计 它要求加窗后的功率谱仍是非负的 这样窗函数w m 的选择必须满足一个原则 即它的傅里叶变换必须是非负的 例如巴特利特窗就满足这一条件 为了采用FFT计算 4 2 3 式 设FFT的变换域为 0 L 1 必须将求和域 M 1 M 1 移到 0 L 1 功率谱的计算公式如下 k 0 1 2 L 1 4 2 7 按照 4 2 1 式估计自相关函数 我们已经证明这是渐近一致估计 但经过傅里叶变换得到功率谱的估计 功率谱估计却不一定仍是渐近一致估计 可以证明它是非一致估计 是一种不好的估计方法 下面我们将证明 BT法中用有偏自相关函数进行估计时 它和用周期图法估计功率谱是等价的 因此BT 法估计质量和周期图法的估计质量是一样的 4 2 2周期图法 将功率谱的另一定义 4 1 6 式重写如下 如果忽略上式中求统计平均的运算 观测数据为 x n 0 n N 1 便得到周期图法的定义 4 2 8 图4 2 1用周期图法计算功率谱框图 1 周期图与BT法的等价关系 周期图法的功率谱估计公式用 4 2 8 式表示 下面由该公式出发推导它们的等价关系 令m k n 即k m n 则 上式中的方括号部分正是有偏自相关函数的计算公式 因此得到 因此证明了利用有偏自相关函数的BT法和周期图法的等价关系 2 周期图法谱估计质量分析1 周期图的偏移 已知自相关函数的估计值 m N 1 N N 1 0 1 2 N 1 按照 4 2 2 式求功率谱的统计平均值 得到 有偏自相关函数统计平均值已由第一章 1 3 30 式确定 将该式代入上式 得到 4 2 9 式中 4 2 10 4 2 9 式中 两序列乘积的傅里叶变换 在频域服从卷积关系 得到 4 2 11 式中 4 2 12 WB ej 称为三角谱窗函数 4 2 11 式表明 周期图的统计平均值等于它的真值卷积三角谱窗函数 因此周期图是有偏估计 但当N 时 wB m 1 三角谱窗函数趋近于 函数 周期图的统计平均值趋于它的真值 因此周期图属于渐近无偏估计 2 周期图的方差 由于周期图的方差的精确表示式很繁冗 为分析简单起见 通常假设x n 是实的零均值的正态白噪声信号 方差是 x2 即功率谱是常数 x2 其周期图用IN 表示 N表示观测数据的长度 按照周期图的定义 周期图表示为 下面先求周期图的均值 再求其均方值 式中 4 2 13 上式说明周期图是无偏估计 但前面已推导出周期图是有偏估计 一般情况 这里由于对信号作了实白噪声的假设 才有无偏估计的结果 在求均方值时 先求两个频率 1和 2处的均方值 最后令 1 2 利用正态白噪声 多元正态随机变量的多阶矩公式 有 将上式代入周期图的均方值公式中 得到 4 2 14 将 1 2代入上式 得到 4 2 15 显然 当N趋于无限大时 周期图的方差并不趋于0 而趋于功率谱真值的平方 即 4 2 16 这里无论怎样选择N 周期图的方差总是和 4x同一个数量级 我们知道 信号的功率谱真值是 2x 说明周期图的方差很大 周期图的均方误差也是非常大 用这种方法估计的功率谱在 2x附近起伏很大 故周期图是非一致估计 是一种很差的功率谱估计方法 为了进一步说明数据长度N对功率谱估计的影响 下面求两个频率处的协方差函数 将 4 2 13 式和 4 2 14 式代入上式 得到 令 1 2 k N 2 2 l N 式中k l均是整数 得到 4 2 17 图4 2 2白噪声的周期图 4 2 3经典谱估计方法改进 1 平均周期图法平均周期图法是基于这样的思想 对一个随机变量进行观测 得到L组独立记录数据 用每一组数据求其均值 然后将L个均值加起来求平均 这样得到的均值 其方差将是用一组数据得到的均值的方差的1 L 假设随机信号x n 的观测数据区间为 0 n M 1 共进行了L次独立观测 得到L组记录数据 每一组记录数据用xi n i 1 2 3 L表示 第i组的周期图用下式表示 4 2 18 将得到的L个周期图进行平均 作为信号x n 的功率谱估计 公式如下 4 2 19 为了分析偏移 对上式求统计平均 得到 4 2 20 4 2 21 周期图的统计平均值已经求出 如 4 2 11 4 2 12 式所示 重写如下 上式表明 平均周期图仍然是有偏估计 偏移和每一段的数据个数M有关 由于M N 平均周期图的偏移比周期图的偏移大 表现在三角谱窗主瓣的宽度比周期图主瓣的宽度宽 由于三角谱窗主瓣的宽度变宽 分辨率更加降低 因此也可以说 偏移的大小反映分辨率的低与高 按照 4 2 19 式求方差 由于是L次独立观测 L个周期图相互独立 因此平均周期图的方差为 4 2 22 即平均周期图的估计方差是周期图的方差的1 L 显然 是以分辨率的降低换取了估计方差的减少 当然 估计的均方误差也减少 图4 2 3平均周期图法 2 窗函数法 这种方法是用一适当的功率谱窗函数W ej 与周期图进行卷积 来达到使周期图平滑的目的的 4 2 23 式中 是有偏自相关函数 M 1 n M 1 4 2 24 那么 4 2 25 将 4 2 25 式和 4 2 3 式进行对比 它们是一样的 说明周期图的窗函数法就是前面提到的BT法的加权协方差谱估计 在 4 2 23 式中 周期图和谱窗函数卷积得到功率谱 等效于在频域对周期图进行修正 使周期图通过一个线性非频变系统 滤除掉周期图中的快变成分 谱窗函数需具有低通特性 对 4 2 25 式求统计平均 得到 将 1 3 30 式和 1 3 32 式代入上式 得到 4 2 26 式中 上式表明 周期图的窗函数法仍然是有偏估计 其偏移和wB m w m 两个窗函数有关 如果w m 窗的宽度比较窄 M比N小得多 这样 m N 则wB m 1 4 2 26 式可以近似写成下式 4 2 27 3 修正的周期图求平均法这种方法和平均周期图法一样 首先把数据长度为N的信号x n 分成L段 每一段数据长度为M N LM 然后把窗函数w n 加到每一个数据段上 求出每一段的周期图 形成修正的周期图 再对每一个修正的周期图进行平均 第i段的修正周期图为 i 1 2 3 M 1 4 2 28 式中 4 2 29 同样 将每一段的修正的周期图之间近似看成互不相关 最后功率谱估计为 4 2 30 对上式求统计平均 得到 4 2 31 式中 4 2 32 U称为归一化因子 4 2 31 式的证明可参考文献 2 为使功率谱估计渐近无偏 4 2 28 式和 4 2 32 式中的归一化因子是不可缺少的 这种在计算周期图之前 先对各数据段加窗函数的方法 使平均周期图方法的估计方差减少 当然分辨率同样减少 但这种方法对窗函数没有限制 此外 分段时 相邻的两段可以有重叠 进一步使方差减少 可以重叠50 6 总之 传统的功率谱估计方法无论采取哪一种改进方法 总是以减少分辨率为代价 换取估计方差的减少 提高分辨率的问题无法根本解决 对于由白噪声和正弦信号或者窄带信号组成的信号 在计算周期图之前 给数据加窗是有利的 如果不加数据窗 相近的低电平信号可能被高电平信号的旁瓣淹没掉 在图4 2 4 a 中 0 12处的正弦信号的旁瓣几乎掩盖了在 0的信号 经过给数据加哈明窗处理后 功率谱如图4 2 4 b 所示 由于加了哈明窗 大大压低了旁瓣 使低电平信号清晰可见 但由于主瓣加宽 功率谱波峰变宽了 却降低了信号的分辨率 一般来说 两个等幅的正弦信号的频率相隔很近 可以不加数据窗 频率间隔应该大于2 N才能分辨 图4 2 4利用数据窗减少窄带过程周期图的旁瓣 a 没有数据窗 b 加哈明数据窗 4 3现代谱估计中的参数建模 在第一章中我们介绍了用参数模型来描述随机信号的方法 如果能确定信号x n 的信号模型 根据信号观测数据求出模型参数 系统函数用H z 表示 模型输入白噪声方差为 w2 信号的功率谱用下式求出 按照这种求功率谱的思路 功率谱估计可分成三个步骤 1 选择合适的信号模型 2 根据x n 有限的观测数据 或者它的有限个自相关函数估计值 估计模型的参数 3 计算模型输出功率谱 4 3 1模型选择选择模型的主要考虑是模型能够表示谱峰 谱谷和滚降的能力 对于具有尖峰的谱 应该选用具有极点的模型 如AR和ARMA模型 对于具有平坦的谱峰和深谷的信号 可以选用MA模型 既有极点也有零点的谱应选用ARMA模型 相对地说 ARMA模型适用范围较宽 对于滚降太快的谱 没有一种模型可以准确地表示功率谱 可以选用高阶的AR模型近似表示 如果选择不合适 例如 选用MA模型去估计具有尖峰的功率谱 估计效果会很差 图4 3 1表示的是用MA模型估计二阶AR信号功率谱的例子 图中 a b c 中的AR信号谱峰较平坦 用二阶MA模型功率谱拟和真实谱时 差别较大 随着阶数的提高 估计的谱愈来愈近似于真实的谱 但是对于 d e f 中AR信号谱峰很窄 再用MA信号模型拟和时 直到MA模型阶数提高到10阶 其效果仍很差 实际信号中一般都含有和信号不相关的噪声 对这一类带有噪声的信号 如果信号是AR模型 由于有噪声需要用ARMA模型 如果用AR模型 则需要阶数更高 5 选择模型的另一个考虑是尽量减少模型参数 当然这和选择模型是否合适有关系 前面介绍过三种信号模型均有普遍应用价值 但模型选择合适 功率谱估计和实际的谱拟合得好 如果不合适 只有提高阶数才能得到较近似的谱 这样需要估计的参数增多 同样也会降低谱估计的质量 因此应该在选择模型合适的基础上 尽量减少模型的参数 4 3 2模型参数和自相关函数之间的关系根据观测数据如何来确定模型的参数对各种功率谱估计方法不尽相同 这里先介绍模型参数和信号自相关函数之间的关系 这些关系在功率谱估计中起着很重要的作用 假设模型的差分方程和系统函数分别用下式表示 4 3 1 4 3 2 图4 3 1对AR 2 信号的模型选择 图4 3 1对AR 2 信号的模型选择 1 ARMA模型的系数和信号自相关函数之间的关系 在第一章我们已推导出系统输出功率谱与输入功率谱之间的关系 用 1 4 7 式表示 重写如下 式中 Pxx z 和Pww z 分别表示模型输出信号x n 和输入信号w n 功率谱的Z变换形式 H z 表示模型的系统函数 将上式两边同乘以A z 得到 4 3 3 将上式进行Z反变换 公式左边的反变换为 公式右边的反变换为 式中 代入上式 得到 这样得到 由于模型H z 是因果的 h n 0 n 0 可以得到 m 0 1 2 3 q m q 1 此式即是ARMA模型输出自相关函数与模型参数之间的关系式 如果能由信号的观测数据估计出信号的自相关函数 可以按照 4 3 4 式求出ARMA模型的参数 再由 1 4 8 式求出信号的功率谱 但我们看到 在公式中由于项的存在 模型输出自相关函数和模型参数之间的关系是非线性的 从而增加了估计功率谱的困难 但是当m q时 4 3 4 式却是一个线性方程 用矩阵方程表示如下 上式共有p个方程 可以用该方程首先计算出AR部分的p个系数hA i i 1 2 3 p 然后代入 4 3 4 式 设法求出MA部分的系数 2 AR模型的系数和信号自相关函数之间的关系 AR模型的系统函数为 H z 1 A z 相当于ARMA模型中B z 1的情况 这样在公式中 因为h n 是因果性的 因此m 0时h m 0 将上面公式代入到 4 3 4 式中 得到 m 1 m 0 4 3 6 也可以将上式中m 1的情况写成矩阵形式 4 3 7 4 3 6 式的矩阵形式如下式 4 3 8 或者用模型参数表示 4 3 9 令 4 3 10 称为自相关矩阵 它满足 H表示共轭转置 是一个埃尔米特 Hermitian 矩阵 且沿任一对角线的元素相等 是一个托布列斯 Toeplitz 矩阵 也是正定矩阵 上面推导出的 4 3 6 式或 4 3 9 式确定了AR模型参数 包括模型输入噪声方差 和信号自相关函数之间的关系 我们注意到这是一个线性方程 如果能够由信号的观测数据求出其自相关函数 可以按照 4 3 7 式 通过解一组线性方程得到模型参数 相对ARMA模型 这是AR模型的优点 3 MA模型的系数和信号自相关函数之间的关系 MA模型的系统函数H z B z 相当于ARMA模型中A z 1 hA n n 的情况 此时h n hB n 由 4 3 4 式得到 MA模型系数和信号自相关函数的关系为 4 3 11 上式表明 MA模型的参数和信号自相关函数之间也是非线性关系 上面我们分别推导了三种信号模型的参数和信号自相关函数之间的关系 这些关系式为我们提供了一种估计功率谱的方法 即首先根据信号观测数据估计信号自相关函数 再按照所选择信号模型 解上面相应的方程 求出模型参数 最后按照下式求出信号的功率谱 4 3 12 4 4AR模型谱估计的性质 4 4 1AR模型的线性预测在第二章中 我们已推导出维纳线性一步预测器系数和信号自相关函数之间的关系式 称为Yule Walker方程 重写如下 4 4 1 式中 e n 表示线性一步预测误差 其公式为 4 4 2 E e2 n min表示e n 的均方差最小值 api i 1 2 3 p 表示预测器的系数 它和线性预测器单位脉冲响应h n 之间差一符号 即 对 4 4 2 式进行Z变换 得到 4 4 3 令He z E z X z 由上式 得到 4 4 4 称He z 为线性一步预测误差滤波器 其作用是将信号x n 转换成预测误差e n 如图4 4 1所示 一般认为e n 具有白噪声的性质 因此He z 也称为白化滤波器 图4 4 1预测误差滤波器 我们知道AR模型的系统函数为 4 4 5 对比 4 4 4 4 4 5 两式 当api ai i 1 2 3 p 时 He z 和H z 互为逆滤波器 He z 1 H z 基于以上分析 也可以将AR模型定义为 4 4 6 式中 是基于信号前p个样本的最佳一步线性预测 公式为 w n 是模型输入白噪声 将AR模型参数和信号自相关函数之间的关系式 4 3 9 和 4 4 1 式进行对比 得到 api ai w n e n E e2 n min 2w 信号自相关函数和它的AR模型参数之间的关系服从Yule Walker方程 注意 只有当AR模型的阶数与线性预测器的阶数相同时 以上结论才是正确的 由于AR模型具有这种特性 因而AR模型法也称为线性预测AR模型法 AR模型与线性预测之间的关系 可以被用来解卷积 假设信号s n 通过一个AR系统 系统单位脉冲响应为h n 响应是x n 即 x n s n h n 如果s n 具有白噪声性质 可以利用AR模型与预测滤波器之间的关系对上式进行解卷积 得到s n 信号 方法是 先由x n 的观测数据估计AR模型的参数 得到AR模型的系统函数H z 再让x n 通过H z 的逆滤波器H 1 z 便得到信号s n 4 4 2预测误差滤波器的最小相位特性我们知道 AR模型必须因果稳定 即极点均在单位圆内 才能保证信号x n 是平稳随机信号 于是AR模型H z 和预测误差滤波器He z 互为逆滤波 那么He z 应为最小相位系统 但是由解Yule Walker方程得到AR模型的参数 其极点不一定在单位圆内 下面将证明当最佳P阶线性预测系数与AR模型参数相同时 由此得到的极点保证在单位圆内 AR滤波器稳定 预测误差滤波器He z 或者A z 是最小相位系统 这里自相关矩阵是正定的 解 Yule Walker方程得到的是最佳线性预测滤波器的系数 此时预测误差滤波器输出功率Pe达到最小 用Pemin表示 即 4 4 7 式中 ai是最佳预测系数 下面先用反证法证明A z 的全部零点不在单位圆外部 即全部零点在单位圆上或者单位圆内部 设A z 的第i个零点zi在单位圆外部 即 zi 1 用1 z i代替zi 这时A z 的幅度函数不变 按照 4 4 7 式计算出的预测误差滤波器输出功率Pe应不变 仍是最小的 但是下面将用公式证明预测误差滤波器输出功率Pe不是最小的 这一矛盾的结论只能说明A z 的零点不可能在单位圆外部 4 4 8 式中 将 4 4 8 式代入 4 4 7 式 得到 4 4 9 式中 因为 因此 4 4 10 4 4 3AR模型隐含自相关函数延拓特性AR模型的自相关函数和模型系数之间的关系服从Yule Walker方程 重写如下 上式中 对于m 1的情况 公式本身就是一个递推方程 如果已由观测数据计算出p 1个自相关函数 用 m 0 1 2 p表示 对于m p的情况 可以用该公式外推得到 公式如下 4 4 11 上式中 系数hA l 需用 4 3 6 式求出 因此AR模型隐含着自相关函数外推的特性 我们知道 经典谱估计BT法中 自相关函数只能限于由观测数据计算出的有限个自相关函数 其它的认为是0 造成了谱估计分辨率低 模糊 也正是AR模型具有自相关函数外推特性 使它具有高分辨率的优点 4 5AR谱估计的方法 4 5 1自相关法 列文森 Levenson 递推法 自相关法的出发点是选择AR模型的参数使预测误差功率最小 预测误差功率为 4 5 1a 假设信号x n 的数据区在0 n N 1范围 有p个预测系数 N个数据经过冲激响应为api i 0 1 2 p 的滤波器 输出预测误差e n 的长度为N P 因此应用下式计算 4 5 1b 显然 e n 的长度长于数据的长度 上式中数据x n 的两端需补充零点 这相当于无穷长的信号经过加窗处理 得到长度为N的数据 用 4 5 1b 式对系数api的实部和虚部求微分的方法使预测误差功率最小 得到 4 5 2 式中自相关函数采用有偏自相关估计 即 m 0 1 2 p m p 1 p 2 1 4 5 3 对比 4 3 7 式 4 5 2 式就是已推导出的Yule Walker方程 因此自相关法也是基于解Yule Walker方程的一种方法 首先由信号的观测数据估计出其自相关函数 再解该方程 得到模型参数 便可求出信号的功率谱 因此该方法也称为 Yule Walker法 但是直接解该方程 需要计算逆在矩阵 不方便 在第三章自适应滤波器中 曾介绍了基于Yule Walker方程中自相关矩阵的性质 导出Levenson Durbin递推法 这是一种高效的解方程方法 下面把已推出的 Levenson CD 2 Durbin递推法简化重写如下 4 5 4 i 1 2 3 k 1 4 5 5 4 5 6 由k 1开始递推 递推到k p 依次得到 a11 21 a21 a22 22 ap1 ap2 app 2p AR模型的各个系数以及模型输入白噪声方差求出后 信号功率谱用下式计算 4 5 7 4 5 6 式表明 说明随着阶数增加 预测误差功率将减少或者不变 为此要求 akk 1 akk称为反射系数 另外 递推公式提供了一种确定模型阶数的实验方法 如模型的阶数不知道 由低阶开始递推 当递推到M阶时 预测误差满足允许的值 停止递推 选AR模型的阶数为M 这种递推法效率高 且当阶数变化时 无需从头计算 如果知道信号的N个观测数据 x n 0 n N 1 利用列文森递推法计算功率谱的计算流程图如图4 5 1所示 图中采用有偏自相关估计 4 5 3 式计算 2是要求的方差 图4 5 1利用列文森递推法计算功率谱的流程图 4 5 2协方差法与修正协方差法1 协方差法这种方法和自相关法一样 仍然利用使预测误差功率最小的方法求模型参数 但由观测数据求预测误差功率的公式如下式 4 5 8 将该式对比自相关法中求预测误差功率的公式 4 5 1b 不同的是求和限不同 该公式中使用的观测数据均已得到 不需要在数据两端补充零点 因此比较自相关法去掉了加窗处理的不合理假设 为求得模型参数仍然应用复梯度法使 4 5 8 式达到最小 公式如下 4 5 9 式中 4 5 10 白噪声的方差为 4 5 11 由观测数据x n n 0 1 2 N 1 利用上面三个公式可以求出模型的参数 api i 1 2 3 p 2w 按照定义 4 5 9 式中的cxx j k 可以称作协方差函数 它有两个变量 因此也适合于非平稳随机信号 4 5 9 式中的协方差矩阵是埃尔米特 Hermitian 矩阵 cxx k j c xx j k 是半正定的 这种方法近似于自相关法 一些实验结果说明它的分辨率优于自相关法 5 另外对于纯正弦信号数据 可以有效地估计正弦信号的频率 进一步的内容可参考文献 14 15 2 修正协方差法修正协方差法使用前向和后向预测误差平均值最小的方法 估计AR模型的参数 进而估计信号的功率谱 该方法最初由 Nuttall 在1976年提出 称为前向 后向法 同年 Ulrych 和Ctayton也独立提出 称之为最小二乘法 因此该法也称为前后向线性预测误差功率最小的最小二乘法 前面我们已推导出信号的前向和后向预测 分别重写如下 4 5 12 4 5 13 式中apk 是AR模型的参数 最小预测误差平均功率是模型输入白噪声的方差 即 p 2w 前后向预测误差平均功率为 4 5 14 式中前向和后向预测误差功率 pe pb分别用下式表示 4 5 15 4 5 16 和协方差法一样 上式中仅对那些用到的观测数据预测误差求和 为了使预测误差平均功率最小 求 p对apk k 1 2 3 p 的实部和虚部的微分 或者用复梯度法求 得到 4 5 17 经过简化 得到 令 4 5 18 将上式写成矩阵形式 4 5 19 白噪声的方差估计值为 4 5 20 例4 5 1已知信号的四个观察数据为x n x 0 x 1 x 2 x 3 2 4 1 3 分别用自相关法和协方差法估计 AR 1 模型参数 解 1 自相关法 2 协方差法 4 5 3伯格 Burg 递推法 设信号x n 观测数据区间为 0 n N 1 前向 后向预测误差功率分别用 p e和 p b表示 预测误差平均功率用 p表示 公式分别为 4 5 21 4 5 22 4 5 23 下面解释 4 5 21 4 5 22 式中的求和限 我们已经知道 前向 后向预测误差公式分别为 4 5 24 4 5 25 在上面两个公式中 信号项的自变量最大的是n 最小的是n p 为了保证计算范围不超出给定的数据范围 在 4 5 21 和 4 5 22 式中 选择求和范围为 p n N 1 为求预测误差平均功率 p最小时的反射系数kp 令 4 5 26 将第三章中已推导出的前向 后向预测误差递推公式重写如下 4 5 27 4 5 28 将上面两个公式带入 4 5 21 4 5 22 4 5 23 式中 得到 利用复梯度法 求 p对kp的实部和虚部的微分 并令其结果等于0 得到 解出kp 得到 4 5 29 上式就是利用伯格递推法求第p个反射系数的公式 将伯格递推法求AR模型参数的递推公式总结如下 4 5 30 4 5 31 4 5 32 4 5 33 4 5 34 对于p阶AR模型的输入白噪声方差 2w p 利用伯格递推法求AR模型参数的流程图如图4 5 2所示 图4 5 2伯格递推法流程图 4 5 4关于AR模型阶次的选择在AR模型谱估计中 模型阶次的选择也是一个关键问题 一般模型的最好选择是先验未知的 实际中需预先选定模型阶次 如果是纯P阶AR信号 选择模型阶次k P时 将产生对谱的平滑作用 降低谱的分辨率 如图4 5 3所示 图中 AR信号P 4 选择模型的阶次k 2 产生的平滑作用使两个峰变成一个峰 分辨率明显降低 如果选择k P 且假定观测的数据没有误差 没有干扰 估计的参数应是 图4 5 3AR模型阶次太小时的平滑作用 因此 对于纯P阶AR信号 应选择阶次k P 一般对于纯AR信号 只要记录数据的长度不十分短 模型阶次的估计结果比较满意 21 但是如果是白噪声中的AR信号 观测数据有误差或者信号中含有白噪声 此时选择ARMA模型合适 如选择了AR模型 其阶次应加大 较低的阶次会使谱估计产生偏移 降低分辨率 当然 这也和信噪比有关 信噪比愈低 平滑作用愈严重 愈需要高的阶次 因此信噪比低应选高的阶次 一般来说 阶次愈高 分辨率愈高 但阶次太高 会使估计误差加大 谱峰分裂 因此 对于白噪中的AR信号 其阶次的选择应折衷考虑 5 下面介绍三种模型阶次的估计方法 为了确定模型的阶次 所提出的各种方法都是基于估计出的预测误差功率 即选择模型的阶次 使预测误差功率极小化 这是因为所讨论的AR模型参数估计方法 预测误差功率都是随模型的阶次的增加而减少 或者不变的 但是不能仅仅考虑预测误差功率的减少 也要考虑因模型阶次增高引起模型参数加多 引起估计误差的加大 根据这一原理 阿凯克提出两种估计方法 第一种方法称为最终预测误差 FPE 准则 用下式计算后再取最小 估计模型阶次 4 5 40 式中 表示k阶AR模型的白噪声方差 预测误差功率 的估计值 N是观测数据的长度 k是模型阶次 上式中 随k增加而减少 但 N k N k 却随k增加而加大 第二种方法应用更广泛 称为阿凯克信息论准则 4 5 41 选择使上式最小的k值作为模型的阶次 但这种方法有过高估计模型阶次的趋势 例如真实的阶次是2 而AIC k 最小值对应的是4 为改进 提出用klnN代替公式中的2k的方法 22 23 AIC和FPE估计阶次的性能是相近的 但对于短记录数据 建议使用AIC准则 24 对于长记录 N 两种方法得到相同的模型阶次 第三种方法由Parzen在1976年提出 称为自回归传递函数准则 CAT 25 定义是 4 5 42 式中 同样 模型的阶次取使上式最小的k值 以上三种估计模型阶次的方法性能类似 26 还有其它的模型阶次估计方法 可参考文献 27 29 4 6最大熵谱估计与最大似然谱估计 4 6 1最大熵谱估计 1 利用最大熵的原则外推自相关函数按照Shannon对熵的定义 当随机变量X取离散值时 熵的定义为 4 6 1 式中pi是出现状态i的概率 当X取连续值时 熵的定义为 4 6 2 式中 p x 是X的概率密度函数 对于离散随机序列 概率密度函数用联合概率密度函数代替 显然 熵代表一种不确定性 最大熵代表最大的不确定性 或者说最大的随机性 下面我们研究对于有限的自相关函数值不作任何改变 对于未知自相关函数用最大熵原则外推 即不作任何附加条件的外推方法 假设x n 是零均值正态分布的平稳随机序列 它的N维高斯概率密度函数为 式中 按照 4 6 2 式 x n 信号的熵为 4 6 3 式中det Rxx N 表示矩阵Rxx N 的行列式 由上式表明为使熵最大 要求det Rxx N 最大 若已知N 1个自相关函数值rxx 0 rxx 1 rxx N 下面用最大熵方法外推rxx N 1 设rxx N 1 确实是信号自相关函数的第N 2个值 根据自相关函数的性质 由N 2个自相关函数组成的矩阵为 4 6 4 它必须是非负定的矩阵 即 4 6 5 将行列式展开 det Rxx N 1 是rxx N 1 的二次函数 该二次函数系数的符号是 1 1 N 2 1 1 N 1 1 且det Rxx N 1 对rxx N 1 的二次导数是 2det Rxx N 1 它是负值 负值表示det Rxx N 1 对rxx N 1 的一次导数是减函数 det Rxx N 1 作为rxx N 1 的函数 凹口向下 那么只有一个最大值 为选择rxx N 1 使det Rxx N 1 最大 解下列方程 4 6 6 用数学归纳法 得到 4 6 7 上式是rxx N 1 的一次函数 可以解出rxx N 1 继续再将rxx N 1 代入Rxx N 2 和det Rxx N 2 中 求det Rxx N 2 对rxx N 2 的最大值 得到rxx N 2 以此类推 可推出任意多个其它自相关函数值 而不必假设它们为零 这就是最大熵谱估计的基本思想 2 最大熵谱估计与AR模型谱估计的等价性 我们已经知道AR模型信号自相关函数与模型参数服从 Yule Walker方程 即 将m 1的情况写成矩阵形式 式中 ai hA i i 1 2 3 N ai是AR模型系数 解该方程 可以得到模型系数ai 即 4 6 8 4 6 9 4 6 10 4 6 11 在 4 3 6 式中 令m N 1 得到 4 6 12 将以上求出的系数a1 a2 aN代入上式 求出rxx N 1 而最大熵外推自相关函数的公式是 4 6 7 式 按照该公式的最后一行展开 得到 4 6 13 上式即是最大熵外推自相关函数的公式 对比 4 6 12 式 两公式完全一样 证明了AR模型功率谱估计和最大熵谱估计的等价性 这里最大熵外推自相关函数等价于已知N 1个自相关函数 匹配一个N阶AR信号模型的系数 一旦通过解Yule Walker方程 解出模型参数 最大熵谱估计用下式计算信号功率谱 4 6 14 4 6 2最大似然谱估计 最小方差谱估计最大似然谱估计是用一个FIR滤波器实现 该滤波器对所关心频率的正弦信号 可以无失真地通过 而对于其它频率的信号 让其频响尽可能地小 亦即将它们尽可能地滤除 此时 滤波器输出的均方值 就作为信号的功率谱估计 设实信号用x n 表示 FIR滤波器系统函数用A z 表示 输出y n 为 4 6 15 式中 输出信号的均方值为 4 6 16 上式中T表示转置 H表示共轭转置 Rp E XXT 是Toeplith自相关矩阵 为求 必须先求FIR滤波器的系数 求这些系数的原则是 在所关心频率 i处 信号x n 无失真地通过 即在 i处的传输函数为1 式中 4 6 17 另外一个原则是在 i附近的频率分量尽量衰减掉 即 i处 滤波器输出y n 的均方差最小 即 4 6 16 式最小 此时作为信号x n 的功率谱估计 因此 最大似然谱估计称为最小方差谱估计更为合适 但由于习惯也可以仍称为最大似然谱估计 在以上原则下 使方差最小的滤波器系数和分别为 30 31 应该指出 此时并不是真正意义上的信号功率谱 只是描述了信号功率谱的相对强度 下面分析最小方差谱估计与AR模型谱估计之间的关系 4 6 20 4 6 21 4 7特征分解法谱估计 4 7 1正弦波用退化AR模型表示无论是实正弦波还是复正弦波 都可以用一个退化AR模型表示 设P个实正弦波组成的信号用下式表示 4 7 1 式中 初相位 i是在区间 均匀分布的随机变量 首先分析下面的三角恒等式 令x n sin n 则上式变为 4 7 2 将上式进行Z变换 得到 4 7 3 这样 4 7 2 式的特征多项式为 4 7 4 上式的两个根分别是 z1 ej z2 e j 它们共轭成对 且模为1 由这两个根可以确定正弦波的频率 对比AR模型的系统函数 可以把正弦波信号用一个特殊的AR 2 模型表示 括弧中的2表示模型是二阶的 该AR模型的激励白噪声方差趋于0 极点趋于单位圆 通常称为退化的AR模型 这一模型系数有两个 即2cos 和1 4 7 2 式是模型的差分方程 对于P个实正弦波 特征多项式是 4 7 5 上式是z 1的2P阶多项式 可以表示为 4 7 6 注意上式中的系数ak k 1 2 3 2P 必须保证它的根共轭成对 考虑到根共轭成对 也可表示为 4 7 7 这样由 4 7 6 式 P个正弦波组合的模型用下面2P阶差分方程描述 4 7 8 对于复正弦波情况 P个复正弦波组成的信号是 4 7 9 用一个退化的AR p 模型表示的差分方程为 4 7 10 其特征多项式为 4 7 11 其根为 1 i P 注意这里的根不是共轭成对出现的 总结以上P个正弦波组合是一个退化的AR 2P 过程 独立参量个数为P个 P个复正弦波的组合是退化的AR P 过程 独立参量个数仍为P个 实正弦过程相应的退化AR过程的阶数比复正弦情况的阶数高1倍 4 7 2白噪声中正弦波组合用一特殊的ARMA模型表示白噪声中正弦波组合的信号为 4 7 12 式中 w n 为白噪声 且 将 4 7 12 式中x n 的用AR 2P 表示 即将 4 7 8 式带入 4 7 12 式中 得到 4 7 13 将 4 7 12 式中的n用n i代替 x n i y n i w n i 再将上式带入 4 7 13 式 得到 4 7 14 上式可以看成一个特殊的ARM