生物医学信号处理5ppt课件

1 第六章功率谱估计的经典方法 功率谱就是指功率谱密度 powerspectraldensity PSD 简称谱 2 研究二阶平稳随机过程特征 功率谱密度 揭示随机过程中所隐含的周期及相邻的谱峰等有用信息 要用有限长的N个样本数据去估计该平稳随机过程的功率谱密度 谱估计的方法 此种估计是建立在时间平均的方法之上 并假定具有遍历性 3 经典谱估计 线性 非参数化方法 周期图法 相关图法等 采用经典的傅里叶变换及窗口截断 对长序列有良好估计 现代谱估计 非线性 参数化方法 最大似然估计 最大熵法 AR模型法 预测滤波器法 ARMA模型等 对短序列的估计精度高 与经典法相互补充 是融合经典变换理论 统计估计理论 系统辨识 信息论 时间序列分析及计算方法等理论与技术 新学科 应用广泛 发展迅速 4 经典谱估计 现代谱估计 完整的谱估计理论 长数据序列 必需应用经典谱估计 性能良好短数据序列 现代谱估计 估计精度分辨度大大高于经典法 5 6 1概述 6 功率谱密度当x t 为广义平稳过程时 其能量通常是无限的 则需研究其功率在频域上的分布 即功率谱密度 对于平稳随机过程 谱分析是采用自相关函数 自相关函数是滞后积的数学期望 Wiener Kinchine定理将自相关函数与功率谱密度联系起来 7 功率谱的物理意义及定义 Wiener Kinchine维纳 钦辛定理 是以0对称 周期为2 反变换为 自相关函数是下列随机过程的数学期望 8 6 2自相关函数和谱密度的估计 已知随机序列 的一段样本数据 利用这段数据估计自相关函数的方法有 为遍历性零均值广义平稳随机过程 无偏估计 对实数序列 该估计的特点 尽管是无偏估计 但不是一致估计 方差很大 不是一种好的估计方法 9 有偏估计 该估计的特点 尽管求平均时 只用N去除不合理 但其服从渐进一致估计的原则 比无偏估计法误差小 因此工程上常用该方法计算 对于实数序列 10 利用自相关函数计算功率谱的实质 周期图 间接法 相关图法 11 即加窗截断为有限长序列 则有 12 功率谱的估计可写成 为数据序列 的离散傅氏变换DFT 直接法 周期图法 13 估计方差为1 长度N 32的白噪声序列的自相关函数和周期图 实现程序为 editransig m然后编辑ransig m文件为运行结果 14 function m ms sq cc cs ransig N 32 x randn 1 N disp 均值 m mean x disp 均方值 ms x x N disp 方差 xstd std x 1 sq xstd xstd cc xcorr x unbiased figure 1 plot x xlabel n ylabel x n title origalsignal figure 2 plot cc xlabel n ylabel cor title auto correlation cs 10 log10 abs fft x 2 N 1 figure 3 f 0 length cs 1 length cs plot f cs xlabel Frequency ylabel PowerSpectrum dB title Periodogram grid 15 均值m 0 2326均方值ms 0 8507方差sq 0 7966 16 17 6 3相关图法 CorrelogramMethod 根据Wiener Kinchine定理 先估计出有限长信号x n 的自相关函数 即 易见自相关函数是偶函数 其长度为2N 1 实际计算时 由于x n 只有N个观测值 则对于每个延迟 m 可用的数据只有N 1 m 个 故可将上式改写成 18 第二步 求 的DFT 得到x n 的功率谱估计 由于功率谱是经自相关函数间接求出的 故称间接法 平稳随机过程的自相关函数应为 19 现在 来讨论相关图法的性能 即 估计偏差 分析相关图法谱估计的性能 周期图是随机变量集合的函数 须从统计观点讨论周期图的估计性能 即上述收敛问题 这是由于 20 所以估计的偏差为 故可得出结论 m 0 为有偏估计 当N m 有限 偏差 0 故为渐进无偏估计 注意到 当 m 接近于N时 滞后积序列ym n 非常短 用来计算平均值的样值数目非常少 得到的将远离 21 周期图不是R m 的傅立叶变换 为有偏估计 22 估计方差 推导过程略 l n k 显然 当N N m l 注意到 当 m 接近于N时 的方差将变得很大 式不再成立 23 m越大 分辨率越高 但自相关的偏差及方差也相应增大 通常取m N 10 N 5 较好 自相关函数的无偏一致估计为 但是这个估计的方差永远大于前面的有偏一致估计 而且用它计算出的功率谱出现负值 与功率谱非负的性质不符 所以工程上常用前者估计 24 6 4周期图法 PeriodogramMethod 定义 长度为N的实平稳随机信号序列 的周期图为 式中 25 分析周期图法谱估计的性能 下面分析偏差和方差 26 不难看出 是个三角函数 两个矩形函数的卷积 被称为Bartlett窗函数 用表示 27 自相关函数的真值 28 Bartlett窗函数及其傅立叶变换示意图 1 29 30 当时 为渐进无偏估计 31 方差的讨论 假定x是方差为的白噪声随机过程 由于周期图的方差与随机过程的4阶矩有关 计算一般随机过程的方差非常困难 但可以计算高斯白噪声的周期图的方差 其结果对一般随机信号有重要参考价值 32 33 34 矩分解定理 35 36 37 38 39 40 41 说明 方差不等于零 周期图的方差不趋近于零 不满足一致估计的条件 周期图不是功率谱的一致估计 方差较大 不是功率谱的好估计 使得取一次估计得到的周期图的一致性和稳定性差 说明此时周期图的方差为常数 与记录数据长度无关 即周期图的方差不会随记录数据的长度的增加而下降 42 高斯白噪声序列的周期图 3个周期图的平均值都近似等于1 但是其方差并不随数据记录长度增加而下降 43 结论 周期图法不是功率谱的一致估计 而且周期图都在真实功率谱附近随机起伏 包括主瓣平滑作用和旁瓣泄漏作用 但是 增加信号的数据点数 数据量增大 相当于滞后窗加宽 相应的傅立叶变换的主瓣变窄 即功率扩散的频率范围变窄 同时 方差不会降低 随机起伏现象仍然存在 且数据点越多 随机起伏越密集 44 随机信号 取N 64和256 观测周期图以及周期图的平均 讨论 点数增大 滞后窗加宽 其傅立叶变换的主瓣变窄 功率扩散的频率范围变窄 45 function sinvnperiod clf case1 N 64x n sin 2 pi f n p v n N 64or256 n 0 N 1 f 0 2 wn1 randn 1 N wn2 randn 1 N wn3 randn 1 N wn4 randn 1 N wn5 randn 1 N xn1 sin 2 pi f n wn1 xn2 sin 2 pi f n wn2 xn3 sin 2 pi f n wn3 xn4 sin 2 pi f n wn4 xn5 sin 2 pi f n wn5 px1 10 log10 abs fft xn1 2 N 1 px2 10 log10 abs fft xn2 2 N 1 px3 10 log10 abs fft xn3 2 N 1 px4 10 log10 abs fft xn4 2 N 1 px5 10 log10 abs fft xn5 2 N 1 figure 1 x 0 length px1 1 length px1 plot x px1 x px2 x px3 x px4 x px5 pxmean 10 log10 px1 px2 px3 px4 px5 5 figure 2 plot x pxmean editsinvnperiod m 46 47 48 49 50 经典谱估计的缺点 加窗造成主瓣平滑作用 使得频率分辨率低 加窗效应 相关图法主观认为未观测数据都等于0 造成频谱能量的泄漏 周期图法不是功率谱的一致估计 而且周期图都在真实功率谱附近随机起伏 周期图法假设数据是以N为周期的周期性延拓 把不真实的信息加于随机过程之上 限制了频率分辨率和谱估计的质量指标 对短时间序列误差太大 51 一些改进方法的思路 将长度为N的序列分k段进行谱估计 分别求出每一段的功率谱 然后加以平均 甚至可以每段数据有部分重叠 得平均周期图 这时 如果各段数据相互独立 则所得估计的方差为原来不分段时得1 k 达到一致估计 缺点是点数减少 分辨率下降 这样会造成有偏的 不用矩形窗 采用合适的窗函数来消除矩形窗旁瓣带来的谱失真 但同时频谱宽度增大 即改进的周期图法 以上两者结合为平均改进周期图法 即韦尔奇法 经典谱估计法 1 2 3 52 例如 为了说明经平均后的周期图作为功率谱估计的实际效果 设有一零均值高斯分布的随机过程 其功率谱密度为 53 与的特性 54 平均后的周期图 每段取16个数据 55 经典谱估计法总结 实质 根据测得的数据 通过计算FFT 傅立叶变换 获得相应的谱密度信息 特点 计算简单 效率高 缺点 由于将测得数据以外的数据均看作零 实质相当于乘上了一个矩形窗口 因此频率分辨率低 频谱能量向旁瓣泄漏 且不是功率谱的一致估计 即使进行改进 也不能从根本上改善周期图的性能 只适用于数据较长和频率分辨率要求不高的场合 56 6 5改进的周期图法 分析 滞后窗的加窗效应源自无限长序列取样时加的数据窗 周期图的期望值被主瓣平滑的程度和功率向旁瓣泄露的多少 取决于加在信号x上的数据窗w n 的类型 当w n 为矩形窗 主瓣窄 但是旁瓣泄漏也最严重 所以 可以选择合适的窗 比如海明窗 以减弱旁瓣泄漏 使弱窄带信号也显现出来 从而提高频率分辨率 但是旁瓣的下降是以主瓣的变宽为代价的 57 修正的周期图 加窗周期图法 定义 N为数据w n 的长度 U为数据窗的平均能量 修正周期图是功率谱的渐进无偏估计 但不是一致估计 58 汉宁窗函数 海明窗函数 59 Bartlett法 周期图的平均 由上式分析知 如果得到周期图的期望的一致估计 即为功率谱密度的一致估计 而一个随机变量的一组互不相关的观测数据的算术平均是该随机变量的均值的一致估计 那么 如果对一个随机信号的若干互不相关的取样序列的周期图进行平均 得到的周期图的平均将是该随机信号的功率谱的一致估计 60 定义 xi n 是x n 的K个互不相关的取样序列 每个序列长度为L 则xi n 的周期图为则这些周期图的平均为 性能分析略 将近无偏估计 且在K和L都趋于无穷大时为一致估计 61 性能分析略 将近无偏估计 且在K和L都趋于无穷大时为一致估计 Bartlett法 对于长度为N的一次观测x n 将它分成长度为L 互不重叠的K段序列 N KL 每个子序列为 则Bartlett周期图为 62 例子 为高斯白噪声 用平均周期图法和Bartlett法求信号的功率谱密度估计 63 N 1024 Ns 256 pxxo 10 log10 abs fft xn 2 N 1 周期图无重叠分段平均法pxx1 abs fft xn 1 256 Ns 2 Ns pxx2 abs fft xn 257 512 Ns 2 Ns pxx3 abs fft xn 513 768 Ns 2 Ns pxx4 abs fft xn 769 1024 Ns 2 Ns pxxb 10 log10 pxx1 pxx2 pxx3 pxx4 4 64 65 66 分析结果 与周期图谱估计比较 无重叠分段平均周期图法 Bartlett法估计曲线均较为平坦 但是波峰加宽 频率分辨率下降 67 Welch法 对Bartlett法进行两点修正 让子序列xi n 有部分重叠 对每个子序列都加数据窗w n 设子序列为xi n 长度为L 相邻两个子序列有L D点重叠 于是第i个子序列为 K个子序列总长度为N L D K 1 若D L 则N KL 即为Bartlett法 性能分析略 将近无偏估计 估计方差比Bartlett法大 68 pxxo 10 log10 abs fft xn 2 N 1 周期图无重叠分段加窗平均周期图w hanning 256 pxx1 abs fft w xn 1 256 Ns 2 norm w 2 pxx2 abs fft w xn 257 512 Ns 2 norm w 2 pxx3 abs fft w xn 513 768 Ns 2 norm w 2 pxx4 abs fft w xn 769 1024 Ns 2 norm w 2 pxxp 10 log10 pxx1 pxx2 pxx3 pxx4 4 重叠分段加窗平均周期图法 即Welch法pxx1 abs fft w xn 1 256 Ns 2 Ns pxx2 abs fft w xn 129 384 Ns 2 norm w 2 pxx3 abs fft w xn 257 512 Ns 2 norm w 2 pxx4 abs fft w xn 385 640 Ns 2 norm w 2 pxx5 abs fft w xn 513 768 Ns 2 norm w 2 pxx6 abs fft w xn 641 896 Ns 2 norm w 2 pxx7 abs fft w xn 769 1024 Ns 2 norm w 2 pxxw 10 log10 pxx1 pxx2 pxx3 pxx4 pxx5 pxx6 pxx7