第5章 功率谱估计.ppt

第 章功率谱估计 5 1经典谱估计5 2自回归模型法5 3最大熵谱估计5 4AR模型参数的求解 频域分析又称谱分析 对于确定性信号 可直接对信号进行傅立叶变换求得其幅度频谱 数字处理中可用快速傅立叶变换 FFT 求得 但是对于一个无始无终的平稳随机信号 它的能量是无限的 其傅立叶变换是不存在的 因而不能求得这种信号的频谱 所以 一个随机信号的频谱在数学上是不存在的 但它的功率谱是存在的 因此我们可以用功率谱来表征一个随机过程的谱特性 根据维纳 辛钦定理 广义平稳随机过程的功率谱是自相关函数的傅立叶变换 因此对于一个观察到的随机信号 重要的是确定它的功率谱密度函数 PSD 和自相关函数 在实际应用中 可以利用的观察数据往往是有限的 所以要准确计算功率谱是不可能的 我们只能通过一个好的估计来得到 至于怎样得到好的估计 这就是这一章里我们要研究的内容 5 1经典谱估计经典谱估计方法实质上就是传统的傅立叶分析法 包含有BTPSD估计法和周期图法 5 1 1BTPSD估计法对于均值为零的平稳随机信号 其功率谱密度函数与自相关函数是一对傅立叶变换对 即 5 1 BTPSD估计法是1958年由Blackman与Tukey提出的 它首先是通过 4 46 对自相关函数进行估计 然后对进行傅立叶变换得到功率谱估计值 即 5 2 上式中为功率谱密度函数 简写成PSD 由于这种方法是将功率谱用有限个自相关函数值的傅立叶变换代替无限个自相关函数值的傅立叶变换求得的 这相当于将无限序列乘上了一个矩形窗函数 这样必将使谱分辨率大大降低 5 1 2周期图法对于平稳随机信号 根据各态历经假设 集合的平均可以用时间的平均代替 于是有 5 3 代入式 5 1 得 令上式可写成 5 4 实际上 式 5 4 在时是不可能收敛的 这是因为对于无限时域的随机信号 它的傅里叶变换是不存在的 对于随机信号的有限个样本序列 由式 5 4 可得到功率谱密度的一个估计为 5 5 这里 是的离散傅立叶变换 即 5 6 显然 是周期的 由式 5 5 所得到的功率谱估计称为周期图 并用表示 即 5 7 周期图法都是用获得的N个数据对随机信号进行谱估计可利用FFT进行计算 因而有计算效率高的优点 在谱分辨率要求不高的地方可用这种方法进行谱估计 但它又一个突出的缺点就是谱分辨率低 因为周期图法隐含着对无限长数据在时域加了一个长度为N的矩形窗 时域中与矩形窗函数的相乘对应于频域中与矩形窗频谱相卷积 所以估计谱就相当于真实谱与矩形窗频谱相卷积的结果 矩形窗频谱为 5 8 图5 1矩形窗频谱的幅度函数它的频谱图如图5 1所示 得到的功率谱估计是它与真实功率谱的卷积 由于它与函数比较有二方面的差别 一是主瓣不是无限窄 二是有旁瓣 因此卷积的结果必然造成失真 由于主瓣不是无限窄的 与主瓣卷积后使功率向附近频域扩散 使信号模糊 降低了分辨率 主瓣愈宽分辨率愈差 如图5 2所示 图5 2 a 是真实谱的两个峰 b 矩形窗谱与真实谱得卷积结果 c 是两个峰离得比较远的情况 此时原峰还能分辨 d 是两个峰离得比较近的情况 此时无法分辨出两个峰的位置 由于矩形窗谱存在旁瓣 也将产生两个结果 其一是PSD主瓣内的能量 一是功率谱主瓣内的能量 泄漏 到旁瓣使谱估计的方差增大 二是与旁瓣卷积后得到的功率谱完全属于干扰 严重情况下 强信号与旁瓣的卷积可能大于弱信号与主瓣的卷积 使弱信号淹没在强信号的干扰中而无法检测出来 图5 2谱分辨率示意图 a 真实谱的两个峰 b 矩形窗谱与真实谱的卷及结果 c 原峰离得较远的情况 d 原峰离得较近的情况 对于BTPSD谱估计法 由于它是按式 5 2 将功率谱用有限个自相关函数的傅立叶变换得到的 这相当于将序列乘上了一个矩形窗函数 因此也存在周期图法同样的缺点 为了减少泄漏和提高谱估计的分辨率 改善窗函数的形状是必要的 但是发现 所有能降低旁瓣的窗口函数都是以主瓣的增宽为代价的 反之亦然 这两个缺点只能互换而不能同时改善 因此用经典法无法克服分辨率低的缺点 可以证明 周期图法是不满足一直估计的条件 即当时 的方差趋于 所以必须对周期图法做进一步的改进 使其满足一致估计的条件 改进周期图的主要方法是平滑或平均 5 1 3巴特利特 Bartlett 平均周期图的方法这种方法是将截取的数据段再分成L个小段 分别计算周期图后取周期图的平均 因为L个平均的方差比随机变量单独的方差小L倍 将序列分成L段 每段有M个点 因而 第段可写成 5 9 第段的周期图为谱估计可定义为L段周期图的平均 即 5 10 可以证明 式 5 10 进行功率谱估计满足一致估计的条件 实际上 对于一个固定的纪录长度N 周期图分段的数目L愈大 则M愈小 估计方差减小 但分辨率降低 估计的方差和分辨率是一对矛盾 它们的效果可以互换 可以根据实际情况适当的选择L和M 如果对分辨率要求不高 可以取L大些 最好将N取大些 分辨率和估计误差都能适当满足要求 5 1 4窗口处理法平滑周期图这种方法是用一适当的功率谱窗函数与周期图进行卷积 来达到使周期图平滑的目的 将式 5 2 重写 5 11 为了减少谱估计的方差 可用窗函数对自相关函数进行加权 此时谱估计公式为 5 12 所以 有 5 13 式中 5 14 也就是说 与分别是和的傅氏反变换 并设序列长 我们知道 是的实 偶 非负函数 为了使是一个实 偶 非负函数 应是一个偶序列 并且满足条件三角窗函数是满足这个条件的 但哈明窗 汉宁窗并不满足这个条件 虽然这两个窗函数能够提供较好的频率分辨率以及较低的旁瓣 但会产生负的功率谱估计 利用窗函数法可平滑周期图 减少估计误差 但是估计偏差加大了 使分辨率降低 现在比较常用的改进方法是Welch法 又叫加权交叠平均法 简记为WOSA法 这种方法先将N长的数据段分成L个小段 每小段M点 相邻小段间交叠点 于是段数为 5 15 然后对各小段加同样的平滑窗后求傅立叶变换 得各小段周期图 5 16 求各小段周期图的平均 得功率谱为 5 17 这里 代表窗函数的平均功率 所以是M长窗函数的能量 这种的优点是无论什么样的窗函数均可使谱估计非负 由于在分段时 可使各段之间有重叠 结果会使估计方差减小 总之 经典谱估计方法总是以减少分辨率为代价 换取估计方差的减少 提高分辨率的问题无法从根本上解决 经典谱估计隐含着数据窗以外的序列值为零的假设 显然这是不合理的 如何利用有限的数据记录 尽可能得到PSD的良好估计是现代谱估计重点研究的内容 5 2自回归模型法任何具有有理功率谱密度的随机信号都可以看成是一白噪声激励一因果稳定的可逆线性网络所形成 自回归模型法谱估计就是由观察获得的随机序列估计的参数 然后根据4 3节中式 4 26 来估计功率谱 即实平稳随机信号的功率谱可表示为 5 18 其中 是白噪声的功率谱 为常数 式 4 2 说明平稳随机信号的功率谱可以用系统的参数来表示 因此 谱估计的问题就转化为模型参数的估计问题 只要估计出模型的参数 就可以得到随机信号的功率谱 由白噪声产生随机序列的模型可以用一个线性差分方程来表示 即 5 19 将上式进行Z变换 得 5 20 其中 5 21 白噪声的功率谱为 输出的功率谱为 5 22 如果能确定与各及就可求得 设且 则式 5 19 称为 5 23 式 5 23 称为p阶自回归模型 简称AR Autoregressive 模型 AR模型的传递函数为 5 24 AR模型的只有极点 没有除原点以外的零点 所以又称全极点模型 式 5 22 成为 5 25 此时 只要我们能求得和所有 就可求得 下面我们就来推导这些AR模型参数与的自相关函数之间的关系 按定义将式 5 23 代入上式 得 5 26 由于只与相关而与无关 故有代入式 5 26 得 5 27 即 5 28 将分别代入上式并写成矩阵形式 得 5 29 上式就是Yule Walker方程 令 5 30 称为自相关矩阵 它有三个性质 它是一个Hermitian矩阵 即 为的转置矩阵 它是一个Toeplitz矩阵 即沿任一对角线上的元素相等 它是一个正定矩阵 即它是特征值权大于零的实对称阵 由Yule Walker方程式 5 28 解出AR模型参数以及 就可得到AR模型的传递函数为 5 31 其功率谱密度为 5 32 5 3最大熵谱估计BTPSD经典谱估计方法用已知的有限个 N个 自相关函数序列的估值求得功率谱密度 它将此有限个估值以外的自相关序列的数据认为是零 因而不能得到好的分辨率 Burg提出的最大熵谱估计 MaximunEntropySpectralEstimation 方法是将有限长度的自相关序列以外的数据用外推法求得 而不是将它们认为是零 假定已知自相关序列 通过合理的外推求得 并要保证外推后的自相关矩阵是非负定的 一般有无限多种可能的外推方法 都能得到比较合适的自相关序列 Berg证明了如果外推后的自相关序列所对应的时间序列具有最大熵 那么这种外推方法才是最合理的 为了帮助理解这种方法的合理性 先简单介绍一下熵的概念 在香农 Shannon 信息论中 离散型随机变量X的熵定义为 5 33 式中表示这一事件发生的概率 是信息量的定义 信息量是解除事件不确定性所需信息的度量 上式表明 熵是平均不确定性的度量 它在数值上等于平均信息量 对于必然事件 由于其发生的概率为1 故其信息量等于零 对应的熵等于零 越是小概率事件 其信息量越大 对应的熵值越大 根据上面的概念不难想象 如果外推后的自相关序列所对应的时间序列具有最大熵 则意味着在具有已知的个自相关值的所有时间序列中 该时间序列的不确定性最大 将是最随机或最不可预测的 统计学认为 最大商是最合理 最自然 最无主观性的假定 这正是选择最大熵准则外推自相关序列的合理性所在 可以求得N维高斯分布信号的熵为 5 34 式中代表矩阵R的行列式 要使H最大 就要求最大 如果已知现欲求得 由于自相关矩阵必是正定的 故矩阵的行列式大于零 即 5 35 为了求得最大熵 要求最大 为此用对上式微分 使微分结果为零 即求得使最大的 满足下列方程 5 36 由此式可解出 于是 用类似的方法可求得 以此类推 每步都按最大熵的原则外推一个自相关序列值 可以外推任意多个而不必认为它们是零 这就是最大熵谱估计的基本思想 可以证明这种最大熵外推自相关函数的结果与AR模型是等价的 使熵最大的功率谱为 5 37 其中与由下列矩阵方程确定 5 38 上式与式 5 28 表示的Yule Walker方程是一致的 因此这就直接说明了最大熵谱估计与AR模型是等价的 5 4AR模型参数的求解由以上的讨论可以看出 最大熵谱估计与AR模型谱估计是等价的 他们都可归结为求解Yule Walker方程中的各AR系数的问题 但是直接从Yule Walker方程求解参数需要作求逆矩阵的计算 当阶P数较大时 运算量很大 并且每当模型阶数增加一阶 矩阵增大一维 需要全部重新计算 所以必须要有对Yule Walker方程提供一个高效率的算法 5 4 1Levinson Durbin递推算法Levinson Durbin递推算法又称自相关法 它是解Yule Walker方程的一种快速高效的算法 此算法是按下列递推法进行的 依次求得附加的的第一个下标是指AR模型的阶数 最后P阶的解即是所要求的解 递推算法以一阶AR模型 求一阶参数及 开始 按式 5 38 一阶AR模型的Yule Walker矩阵方程应为从这个矩阵方程解出及 分别为 5 39 5 40 再从二阶AR模型的矩阵方程 解得分别为 5 41 5 42 5 43 以此类推的递推公式 5 44 5 45 5 46 当我们从式 5 39 和 5 40 得到初始的及的数据以及各 就可按式 5 44 5 45 5 46 以此递推出各阶的 从式 5 46 可知 一般来讲阶数预先是不知道的 当我们递推到第k阶 满足所允许的值 就可选阶数 5 4 2Burg递推算法用Levinson Durbin递推算法求解Yule Walker矩阵方程中的AR模型系数虽然可以简化计算 但需要知道自相关序列 实际上自相关序列只能从时间序列的有限个数据得到它的估计值 当时间序列短时 估计误差很大 这将对AR模型系数的计算引入很大的误差 导致功率谱估计出现谱线分裂与谱峰频率偏移等现象 Burg于1967年提出最大熵谱估计 后来又在另一篇文章里提出直接由时间序列估计AR模型系数的方法 被人们称为Burg递推算法 这种方法与预测误差格型滤波器有密切关系 它是在Levinson关系式 5 45 的约束下 用使前向与后向预测误差能量之和为最小的方法来求得各的值 按线性预测理论 的估计值可用的各过去值的加权之和表示 即 5 47 用表示估计误差 5 48 称为线性预测器的前向误差 因为是由以前的各数据加权之和得到的 再由式 5 45 计算 即 5 49 又令 这里的称为部分相关系数 PARCOR 或反射系数 将这些关系式代入式 5 48 得 5 50 这里所以 5 51 这里 5 52 是由以后的各数据加权之和得到的 故为后向预测误差 将式 5 49 代入式 5 50 可得 5 53 5 54 当时有 5 55 式 5 53 与式 5 54 为前向及后向预测误差的递推公式 Burg递推算法的优点是不需要估计自相关函数 可以从已知的序列求得参数 如果按前向均方误差最小的准则确定并用表示 则按式 5 54 可得令即所以 5 56 如果按后向均方误差最小的准则确定并用表示 同理可得 5 57 Burg算法是以前向均方误差与后向均方误差之和最小的准则求得 令可得 5 58 对于平稳随机过程 集合平均可用时间平均代替 因此上式可写成 5 59 不难看出 5 60 如果已知为有限长序列 当时 则按式 5 59 可得 5 61 而按式 5 50 及 5 53 可得 5 62 5 63 将及代入式 5 58 可得再代入式 5 50 及 5 53 又可得及 再将及代入式 5 58 又可得 将与及代入式 5 50 及 5 53 又可求得及 以此类推 可直接从时间序列求得各阶的以及前向与后向误差及 将各代入式 5 49 又可求得各 于是将这些求得的AR系数代入式 5 52 其中 也可用式 5 46 求出 即可求得功率谱密度的估计值 Burg递推算法的优点是求得的AR模型保证稳定 不需要估算自相关函数 尤其实在短时间序列时能得到较好的估计 具有较高的谱分辨率 图5 3所示的是AR模型法谱估计的Levinson Durbin递推算法和Burg递推算法对白噪声中的两个正弦信号进行谱估计的结果 两个正弦信号对采样频率的归一化频率分别为和 数据窗采用矩形窗 长度取 图5 3参数模型法谱估计 N 512 a 两个正弦信号与白噪声叠加的时域波形 b Levinson Durbin算法 c Burg算法 由于Burg递推算法仍然受到Levinson Durbin递推算法的约束 也就是说AR模型系数除外 其它系数均由Levinson关系式递推求得 而Levinson关系式在方程组得系数矩阵 自相关矩阵 为Toeplize阵的条件下推出的 实际上只有无始无终得平稳随机序列所对应的自相关矩阵才具有这种性质 因此 Burg递推算法不能完全克服Levin