1、1,第7讲 功率谱估计 (教材 第14章）,2,7.1 基于有限长信号观察的功率谱估计,数据序列的有限记录长度是功率谱估计质量的主要限制:,当处理统计平稳信号时，数据记录越长，可从数据提取的信号估计就越好。,当信号统计是非平稳的，则要选择仅可能短且仍能解析不同信号分量谱特征的数据记录。,3,7.1 基于有限长信号观察的功率谱估计,能量密度谱计算,有限长数据序列确定性信号的谱计算,（1）直接方法,（2）间接方法,4,7.1 基于有限长信号观察的功率谱估计,将序列x(n)的长度限制到N点等价于对x(n)乘上一个矩形窗,旁瓣泄露问题：一般要求提供的窗谱W(f)比X(f)相对要窄。然而，即使W(f)比
2、X(f) 窄，X(f)与W(f)旁瓣的卷积产生了 旁瓣能量，在此频带处的真实信号谱X(f)=0，这种旁瓣能量称为泄露。,5,7.1 基于有限长信号观察的功率谱估计,X(f)与长度为N=61的矩形窗卷积,例：,6,7.1 基于有限长信号观察的功率谱估计,为了减少旁瓣泄露，一般选择具有低旁瓣的窗。如选择平滑的时域截止特性代替矩形窗陡峭的截止特性，但却使频谱特性X(f)的平滑特性或展宽特性得到增加。,7,随机信号的自相关和功率谱估计：周期图,对于该随机过程的单个实现，从中估计该过程的功率密度谱,广义平稳随机过程不具有有限能量，因此不能进行傅立叶变换。这类信号具有有限平均功率，因此表征为功率密度谱。,
3、时间平均自相关函数：,如果平均随机过程的均值和自相关函数是各态遍历的：,8,随机过程的单个实现样本估计功率密度谱：周期图,7.1 基于有限长信号观察的功率谱估计,rxx(m)是统计自相关函数xx(m)的一致估计,当N时， varrxx(m)收敛于零,9,7.1 基于有限长信号观察的功率谱估计,7.1 基于有限长信号观察的功率谱估计,估计谱的均值是真实谱的平滑版，受损于由有限数据点导致的 相同的谱泄露。,估计的自相关rxx(m)是真实自相关函数xx(m)的一致估计。 但其傅立叶变换Pxx(f)，即周期图不是真实功率谱密度 的一致估计，而是渐近无偏估计。但对于一个有限长序 列， Pxx(f)均值包
4、含了偏移，因此真实功率谱密度产生 了失真。,7.2 功率谱估计的非参数化方法,对于一个有限长序列， Pxx(f)均值包含了偏移，因此真实功率谱密度产生了失真。估计谱受损于巴特利窗的平滑效应和具体的泄露。,1. Bartlett方法：平均周期图,首先，N点序列被划分为K个不重叠段，每段的长度为M,然后，对每一段计算周期图,最后，每K段的周期图进行平均得到Bartlett功率谱估计,7.2 功率谱估计的非参数化方法,2. Welch方法：平均修正的周期图,首先，数据段重叠分段,然后，先对数据段进行开窗，再计算周期图,最后，每K段的周期图进行平均得到Bartlett功率谱估计,不同的窗函数将生成不同
5、的方差。若为三角形窗,13,7.2 功率谱估计的非参数化方法,3. Blackman和Tukey方法：平滑周期图,首先对采样自相关序列开窗，然后进行傅立叶变换生成功率谱估计,14,7.3 功率谱估计的参数化方法,非参数功率谱估计方法相对简单，并且使用FFT算法容易计算。然而这些方法需要可用的长数据记录，以得到很多应用所必需的频率分辨率。而且，这些方法受损于有限长数据记录所固有的由开窗带来的谱泄露效应。,非参数功率谱估计方法存在两种假设：假设自相关估计rxx(m)在mN时是零值；假设数据是周期的，周期为N。,信号生成模型可用观察数据估计到的一些参数来构造。从该模型和估计参数，可计算模型所包含的功
6、率谱密度。,7.3 功率谱估计的参数化方法,以参数模型为基础的谱估计方法一般按照下列三个步骤：,（1） 为被估计的随机过程确定或选择一个合理的模型 一般将数据序列x(n)模型化为一个线性系统的输出，该线性系统由有理系统函数形式所表征,16,（2） 根据已知观测数据估计模型的参数。通常模型参数比观测数据的数据量少很多。,（3） 用估计得到的模型参数计算功率谱。,7.3 功率谱估计的参数化方法,给定数据序列，0nN-1，则可估计ak和bk,假定输入序列w(n)是一个零均值白噪声序列,ww (m)=w2(m),7.3 功率谱估计的参数化方法,h(n) H(z)=B(z)/A(z),w(n),x(n)
7、 X(z),w(n)数学期望为0,方差为2的白噪声序列,x(n)的功率谱为xx(ej)：,18,自相关和模型参数之间的关系,19,Yule-Walker方程,对于AR过程：,自相关和模型参数之间的关系,20,自相关和模型参数之间的关系,对于MA过程：,21,AR模型参数的Yule-Walker求解方法,自相关估计的偏移形式：,22,AR模型阶数的选择,使用AR模型最重要的方面是阶数p 的选择。如果选择一个太地阶的模型，将会得到一个高度平滑谱。如果p 选择得太高，则会冒险在谱中引入假的低峰。,（1）最后预测误差(FPE)准则，选择阶数以最小化性能指标,模型阶数选择的准则：,（2）Akaike信息
8、准则（AIC) ，通过选择阶数以最小化,Sper(e j)称为周期图，它是Sxx (ej)的渐近无偏估计,随机过程的功率谱估计,周期图及其估计质量,假设有一个随机过程：x(n)=Asin(0n+)+v(n)（ v(n)白噪声）,x(n)的真实功率谱:,Sxx(e j)=v2+(A2)/2 (-0)+(+0),x(n)的周期图的期望值,ESper(ej)=v2+A2/4 WB(e j(-0)+WB(e j(+0),由于WB(e j)的主瓣不是无限窄的，导致正弦信号中的功率扩展到带宽约为4/N的整个主瓣范围内，使本来是一根谱线的正弦信号的功率谱变成了与滞后窗傅立叶变换主瓣形状相同的功率谱。,由于W
9、B(e j)有许多旁瓣，使与正弦信号功率谱相卷积的结果，在k0+-2k/N等频率点上形成其他的谱峰。,周期图及其估计质量,周期图作为功率谱的估计，不仅会产生偏差，而且由于滞后窗频率特性主瓣的平滑作用，限制了周期图分辨率中任何两个频率相近的窄带成分的能力或频率分辨力。,由于WB(e j)的主瓣宽度随数据记录的减少而增加，因此，对于一定的数据长度N， WB(e j)的主瓣宽度是一定的。这样，周期图能够分辨两个频率相近的正弦信号的能力就是一定的。通常把这种频率分辨率用WB(e j)的主瓣宽度来度量。 = 0.89(2/N)，因此，周期图的频率分辨率为：,例：为使周期图的频率分辨率不大于0.05，数据
10、记录长度应为多少？,周期图及其估计质量,周期图的计算是基于傅立叶变换，在利用FFT做谱分析应考虑的参数选择,(1)采样频率应满足采样定理,(2) 数据长度（N）,(3)采样信号的持续时间 tp=NT,采样间隔,频域分辨率和谱图表示,x(t)=cos(2*1.25*103t)+0.5cos(2*1*103t)+0.25cos(2*0.8*103t),已知一连续信号为,试用DFT计算其幅度谱，并且与原信号进行比较,(2)取采样频率和采样间隔分别为：fs=5kHz, T=0.2ms,(3)N=tp/T=20ms/0.2ms=100,(1)确定该信号的周期：20ms,周期图及其估计质量,高密度频谱和高
11、分辨率频谱的比较。设信号为x(n)=cos(0.48n)+cos(0.52n)利用有限长序列的FFT来分析下列情况的幅度谱：,(2)采集数据长度N=10，但补90个零，做100点的DFT，画出幅度谱。,(1)采集数据长度N=10，即0n9，做10点的DFT，画出幅度谱。,(3)采集数据长度100点，同样画出幅度谱。,t=1:0.1:10; signalD=2*t+5*sin(2*pi*0.5*t); noiseD=randn(size(signalD); measureD=signalD+noiseD; testsig=measureD; s=spectrum.periodogram; Hps
12、d=psd(s,testsig); figure; subplot(1,2,1),plot(t, testsig); title(time series); subplot(1,2,2); plot(Hpsd);,平稳随机过程的功率谱估计,signalD的周期图,32,白噪声序列及其功率谱,功率谱求解程序,34,非平稳随机过程分析,单道脑电信号时序图,35,非平稳随机过程分析,脑电信号直方图分析,36,非平稳随机过程分析,脑电信号的自相关函数,37,非平稳随机过程分析,对脑电信号进行AR参数建模，利用FPE准则，对AR建模阶次p取520之间进行判断。,选择p=11,38,非平稳随机过程分析,脑
13、电实测序列与AR模型的预测输出,39,脑电序列及AR模型输出信号的功率谱,40,非平稳随机过程脑电信号分析,clear; close all; clc; load eeg8s; testsig = eeg8s; Figure,subplot(1,2,1),plot(testsig); xlim(1,1100) subplot(1,2,2), histfit(testsig); Figure, parcorr(testsig,50); %parcorr( )求解自相关函数 %确定AR模型阶数 k=1; for p=5:20m=armax(testsig, p 0); %armax函数确定模型参数
14、并进行预测FPEvalue(k) = m. EstimationInfo.FPE; pvalue(k)=p; k=k+1; end,41,非平稳随机过程脑电信号分析,j= find(FPEvalue=min(FPEvalue); Porder=pvalue(j); figure, plot(pvalue, FPEvalue, -.o, pvalue(j), FPEvalue(j), *); xlabel(p-value); ylabel(FPE value); title (find the minimum FPE); m = armax(testsig, Porder, 0); d1=m.a
15、; p1= m.NoiseVariance; H1,w1= freqz(sqrt(p1),d1); s= spectrum.periodogram; Hpsd = psd (s, testsig); %用周期图估计信号的功率谱 Figure , plot(Hpsd); hold on,42,非平稳随机过程脑电信号分析,hp = plot(w1/pi, 20*log10(2*abs(H1)/(2*pi),r); Set(hp, LineWidth,2); xlabel(Normalized frequency (timespi rad/sample) ylabel(One-sided PSD (
16、dB/rad/sample) legend(PSD estimate of x, PSD of model output); m=armax(testsig, Porder,0); ye= predict( m, testsig,1 ); %模型预测输出 figure, stem (testsig ,ye); xlabel(Sample time); ylabel(Signal time); legend(Original autoregressive signal, signal estimate from linear predictor); P2 = norm (testsig-ye),
17、2)2/(length(testsig)-1);,应用实例,心电信号的计算机分析,心电信号的计算机分析,分析步骤： 心电数据采集：500Hz采样频率 心电信号预处理：滤除干扰（基线漂移、50Hz、肌电） 特征点检测：P、QRS、ST、T波 自动诊断：心律失常分析与波形分类,QRS波形检测算法：,经典的QRS波检测算法包括三部分 ： 线性滤波； 非线性变换； 决策规则。,线性滤波一般采用中心频率在1025Hz之间带宽为510Hz的带通滤波器，用于减除ECG信号中的非QRS波频率成分，提高信噪比。,非线性变换的目的是将每个QRS波信号变换为单向正波峰。,决策规则一般用峰值检测器或自适应阈值检测器来
18、检测QRS波。,由于心电信号是非平稳随机信号，利用经典的QRS波检测算法往往受到以下两个因素的制约：, QRS复合波的信号带宽对于不同的病人乃至同一病人的不同心搏均有所不同。,因此，人们一直致力于采用新的信号处理方法来分析QRS波，常用的有人工神经网络算法和小波变换算法。, 各种噪声与QRS波的通带相互交叠。,心电信号的计算机分析,实时采集的心电信号x(n) 如下图所示：,实时采集的心电信号x(n)通过上限截止频率15Hz的三阶巴特沃兹低通滤波器，滤除高频干扰，得y1(n)；,采用高斯函数一阶导数导出的小波，对y1(n)进行尺度S=22的小波变换，突出信号特征点，消除基线漂移，得y2(n)；,对y2(n)计算差分，取绝对值，并进行三点平滑，得y3(n)；,对y3(n)再计算差分，取绝对值、平滑，得y4(n)；,对y3(n)与y4(n)逐点求和，再平滑，得y5(n)；,进行阈值检测：取连续前3秒内的待检测信号y5(n)的振幅P（可自适应调整），设检测阈值d1=0.25P，d2=0.15P，对超过d1的信号再降低阈值以d2作双向检测，大于d2则赋值为1，得QRS模板信号y6(n)，并记录每个模板区内y5(n)的峰值时刻和峰值；,
