中医现代化研究分享 http://blog.sciencenet.cn/u/baishp 当你们都还在想象时,我已经在路上了。这注定是一场一个人的战争吗?

博文

6、功率谱密度函数估计

已有 9816 次阅读 2012-11-27 22:05 |个人分类:斤斤计较|系统分类:论文交流| 生物钟, 中医现代化, 生物医学信号, 功率谱估计, 时间医学

 
   功率谱密度函数估计,在随机信号处理中具有极其重要的意义。不管是为了目前增加对信号基本属性的了解,还是为了以后对信号作进一步分析处理,现在对敝人各生理信号作一下功率谱估计,都是很有必要的。

    MATLAB信号处理工具sptool中,有8种已经很成熟的功率谱估计方法(FFT法,Welch法,Covariance法,Mod.Covariance法,MTM法,MUSIC法,Burg法,Yule AR法)可供调用,非常方便。

    将体温信号序列TW461512_0导入sptool,并Create(加载)于spectra栏,即打开功率谱处理、观察窗口Spectrum Viewer。

    (1)、先看看快速傅立叶变换FFT方法,也就是所谓的经典周期图法吧。取FFT执行点数Nfft=2^16=65536。在2的n次幂中,它比TW461512_0的长度55380长,且最接近55380。采样频率Fs取Fs=Nfft=65536(次/单位采样时间,我亦称之为“Hz”,将“1次/单位采样时间”的采样频率提高Nfft倍,是为了使频率轴都用正整数表示。

        图6-1 TW461512_0快速傅立叶FFT法估计的功率谱图
 

    此图中左标尺所标示的谱线位置,为频率f=5461Hz处。其波形周期T=Fs/f=65536Hz/5461Hz=12.0007单位采样时间,即1天。其小数点后面的误差为数字化误差。此谱线右边的谱线为其倍频谱线。

 

        图6-2 FFT法功率谱最低频处。频率轴放大128倍,幅值轴放大4倍。
 

    波峰狭窄尖锐,幅值相近,噪声谱与信号谱混杂在一起,很难分辨样本信号谱峰。

         图6-3 FFT法功率谱最高频处。放大倍数同前。

    谱峰形状与最低频处一样。

 

    (2)、来看看使用Welch(经典的韦尔奇)方法估计功率谱。此方法需要选择的参数比较多。除了FFT执行点数Nfft外,还需要确定窗函数Window及其长度Nwind、分段重叠点数Overlap。对于没有足够的先验知识的人来说,要一次选对各个参数是很难的。先比较随意地选定一组参数:Window=hanning,Nwind=4096,Overlap=1024,作一功率谱图如下:


        图6-4 TW461512_0经典welch法估计的功率谱图

 

        图6-5 图6-4最低频处。横轴放大16倍,纵轴放大1.2倍。

 


        图6-6 图6-4最高频处。横轴放大16倍,纵轴无放大。

    直觉上感到此谱图各波形太过平缓,波峰幅值相近,难以确定信号波峰及其频率点。现在以窗函数宽度及其重叠的比例为变量,作搜索式计算。编程如下:

 

%welch法搜索参数
L=length(TW461512_0);%=55380;
m0=512;%窗函数最小宽度
c=floor(L/(2*m0));%信号序列最大分段数
for i=1:4
M=zeros(300,c);
for j=1:c
m=j*m0;%窗宽按每次增加一个最小窗宽数m0搜索
n=m*(i-1)/4;%数据分段重叠比例
Pxx=PWELCH(TW461512_0,hanning(m),n,65536,65536);
M(1:149,j)=Pxx(1:149);%数据太长,无法全部清晰作图显示,只能

%取最低频与最高频段各149点
M(150:151,j)=ones(2,1);%人工隔板
M(152:300,j)=Pxx(length(Pxx)-148:length(Pxx));
end
M=M';%横坐标表示频率,纵坐标表示窗宽
figure(i)
surf(log10(M))
end

 

    运行,得: 

        图6-7 TW461512_0的Welch法估计功率谱。分段无重叠

 


        图6-8 TW461512_0的Welch法估计功率谱。分段重叠比例为1/4

        图6-9 TW461512_0的Welch法估计功率谱。分段重叠比例为2/4

        图6-10 TW461512_0的Welch法估计功率谱。分段重叠比例为3/4
 

    由于数据太长,无法在一张图片中显示整个结果,只能分别取最低频最高频段各149Hz合在一起显示。中间的“隔板”是另加的。

    上面四张图片,分别是信号序列分段无重叠、重叠1/4、2/4、3/4时的谱图。其纵坐标表示窗函数的宽度(以512点为一个单位)。我认为谱估计效果最好的应该是最后一张图,即信号序列分段重叠3/4的时候。因为它在整个纵轴方向上,随着窗宽增加,波峰略为变高变窄,而位置(频率点)都很少改变,而不像前面三张图,各波形波峰高度及其频率点都随着窗函数宽度的改变而有所改变,尤其是窗函数宽度比较小的时候。

    我认为一种好的估计方法,当需要选择的参数落在一个大范围后,其估计结果应该对参数的变化是“不敏感”的。原因很简单:真理只有一个。如果反之,估计结果对参数的变化是“敏感”的,参数稍微一变化,估计结果就跟着很快变化,那么我们最终到底要选哪一个参数?如果随便选一个,那岂不就象买彩票一样?科学要追求确定性的结果,而不是中彩。而这里的参数必定是由人来选定的。如果参数可以由一套算法来计算、确定,那就不用设参数,而设置一个函数好了。

    根据上面的分析,本信号序列采用Welch方法估计功率谱时,窗宽重叠率须取1/2以上。至于窗宽,取一个中间偏大的值如512×(35~54)都可(第一行的窗宽为512点,最后一行窗宽为512×54=27648点)。谱图图片暂就不贴了,等到在同一窗口比较各方法谱估计曲线时一并贴出。

 

    (3)Covariance(协方差)法。此法需选择参数,除了快速傅立叶变换执行点数Nfft,就是阶数Order了。取Nfft=65536,任取几个阶数Order运行,发现此法非常耗用计算机资源。当Order>500时即Out of memory(超出内存)了。

    我的内存为2×1G。要增加内存,牵涉到整机升级的问题,不是换两根内存条就可以了的。据说可以设置虚拟内存,但我设置之后,也没见运算能力的明显提高。不知是何缘故。

    自己编程,看看Order<=500时的谱估计情况:

 

%Covariance法搜索阶数的计算

M=zeros(300,10);
for i=1:10
order=50*i;
Pxx=pcov(TW461512_0,order,65536,65536);
M(:,i)=Pxx(1:300);
end
M=M';
surf(log10(M))

 

    运行,得:

       图6-11 Covariance法搜索阶数的计算 

   什么也没有。所以在个人电脑上是无法使用Covariance(协方差)法的。

 

    (4)、Mod.Covariance(改进的协方差)法。在Spectrum Viewer窗口中

略一运行,便知此法比Covariance法更耗电脑资源。pass。

 

    (5)、MTM(Multitaper多窗口)法。此方法需要选择的参数有时间-带宽乘积nw、fft执行点数Nfft、采样频率Fs、权重算法Weights。

    取Nfft=Fs=65536,Weights='adapt'(Thomson的非线性自适应组合算法),各取nw=2、3、4,得到3条谱估计曲线。使3条曲线在同一窗口显示,见图6-12。


    图6-12 时间-带宽乘积nw不同,其余参数相同时,TW461512_0的功率谱估计曲线。

其中绿线nw=2,红线nw=3,蓝线nw=4。

 


        图6-13 图6-12最低频段处。横轴放大128倍,纵轴放大(为原来)2倍。

 


        图6-14 图6-12最高频段处。横轴放大128倍,纵轴放大(为原来)2倍。
 

    图中3条曲线在大的走势上虽然保持了一致,但在更细节的形状上差异也不小,而且很多波峰形状也不明显,幅值接近,使人很难判断真实的谱线。

    此方法中,2*nw-1表示采用的窗口数。nw的典型取值有2,5/2,3,7/2,4。但我发现在上述取值范围内,谱估计函数Pxx=pmtm(x,nw,Nfft,Fs,Weights)对所有实数nw都是有意义的。仿照我前面搜索welch法参数的方法,作搜索nw的计算。程序如下:

 

%搜索MTM谱估计法中时间-带宽乘积参数nw

Nfft=65536;fft执行点数
Fs=65536;%采样频率
M=zeros(300,30);
for i=16:45
nw=0.1*i;%时间-带宽乘积
Pxx=pmtm(TW461512_0,nw,Nfft,Fs,'adapt');
M(1:149,i)=Pxx(1:149);%最低频段149HZ
M(150:151,i)=ones(2,1);%人造隔板
M(152:300,i)=Pxx(length(Pxx)-148:length(Pxx));%最高频段149HZ
end
M=M';
surf(log10(M))

    运行,得:


        图6-15 TW461512_0的pmtm法估计功率谱。权重算法为adapt

 

    觉得还是很难确定nw到底取多大为好。改程序中权重算法因子为unity(一致权重的线性组合算法)、eigen(特征值权重的线性组合算法),运行,得:

        图6-16 TW461512_0的pmtm法估计功率谱。权重算法为unity。


        图6-17 TW461512_0的pmtm法估计功率谱。权重算法为eigen。

 

    发现上面三张图中,曲面形状非常一致。注:上面三张图及图6-7~图6-10的纵坐标应乘10才是分贝值,因为程序中“surf(log10(M))”应为“surf(10*log10(M))”。懒得再截图了,说明一下。

    回到Spectrum Viewer窗口,取定nw=3,Nfft也不变,分别取Weights='adapt'、'unity'、'eigen',得到3条曲线。令其在同一窗口显示,如下图:


    图6-18 Weights不同,其余参数相同时,TW461512_0的功率谱估计曲线。最低频段处,横轴放大128倍,纵轴放大为2倍。其中绿线Weights='eigen',红线Weights='unity',蓝线Weights='adapt'

 

    可以看出,最低频段处,3条曲线基本上重合了,后面的绿线、蓝线基本上都被前面的红线挡住。


        图6-19 TW461512_0的功率谱估计曲线。最高频段处。其余情况同图6-18


    最高频段处,绿线与红线曲线基本上重合了,蓝线在某些部位不与红线、绿线重合。

    所以,MTM法中参数Weights的选择是不敏感的,选哪一个都差不多。

  

 (6)、MUSIC(多信号分类)法。该法需要选择的参数比较多,有信号子空间数Signal Dim.、阈值、fft长度Nfft、窗函数、窗函数宽度Nwind、窗函数重叠点数Overlap。随选几组参数,在Spectrum Viewer窗口试运行,即知此方法非常占用电脑资源。信号子空间Signal Dim.参数达到4000、加窗宽度Nwind达到7680(=512*15)时便“Out of memory”了。

    现在我觉得,Matlab中Sptool的Spectrum Viewer窗口功率谱估计,最适合的是那些做固定项目工程应用的人。他们对谱估计相关的参数选择都已经很熟悉,只是由于样本数据的变化,而在此窗口查看一下谱估计的结果。对于做科研工作的人来说,大多数是第一次接触到新项目,谱估计参数选择不熟悉。此窗口只能作辅助性的研究工具。

    准备看看在Signal Dim.<=4000、Nwind<=7680时的功率谱估计情况。先固定窗函数及其宽度、窗函数重叠点数,对Signal Dim.作搜索计算。程序在傍晚时开始运行,结果运行到第二天临晨还没有结束,只好强行中断。虽然感觉再要不了多久就会结束,但也不是十分有把握。如果还需要运行几十、几百……个小时,那不麻烦了。第二天先试循环少量几次,对整个程序的运行时间有一个大概的评估后才开始正式的运行。所以对有循环语句的程序,这个工作是少不了的,除非非常有经验。

    将原程序略作改动,分段运行:

 

pack%整理工作空间
Nfft=65536;fft执行点数
Fs=65536;%采样频率
d=100;%每段程序循环搜索次数

k=20;%信号子空间数每次增加数
Dim0=2000;

M_2000_20_4000=zeros(300,d);
m=5120;%窗函数宽度
n=m/2;%数据分段重叠比例
for i=1:d

Dim=Dim0+i*k;
Pxx=pmusic(TW461512_0,Dim,Nfft,Fs,hanning(m),n);
M_2000_20_4000(1:150,i)=Pxx(1:150);
M_2000_20_4000(151:300,i)=Pxx(length(Pxx)-149:length(Pxx));
end
M_2000_20_4000=M_2000_20_4000';

 

    运行完毕后,将其中d=100; k=20;Dim0=2000;分别改成d=100、100、75;

k=10、7、4;Dim0=1000、300、0;M_2000_20_4000改成M_1000_10_2000、M_300_7_1000、M_0_4_300,分4段运行。4段程序大约运行了15、6个小时。

    最后:

M=[M_0_4_300;M_300_7_1000;M_1000_10_2000;M_2000_20_4000];

surf(10*log10(M))

    运行,得:


    图6-20 TW461512 _ 0的pmusic法估计功率谱,对信号子空间Signal Dim.的搜索。hanning窗宽5120点,重叠2560点


    可以看出,Signal Dim.太小的时候,基本上没有什么谱峰。即便在Signal Dim.最大值4000处,谱峰也很平缓。这是低频段的情况。在高频段则始终没有出现谱峰。

    总而言之,本人现在的电脑配置是没法使用MUSIC方法估计本样本数据的功率谱的。最后再贴几张不同参数下的功率谱估计图吧:

        图6-21 Dim=4000,Nwind=15*512=7680

 

        图6-22 Dim=4000,Nwind=7*512=3584

 

        图6-23 Dim=3200,Nwind=2*512=1024

 


        图6-24 Dim=1000,Nwind=2*512=1024 

 

    (7)、Burg法。前面6种估计方法都是属于非参数估计方法,Burg法及最后的Yule AR法属于参数估计方法。在Spectrum Viewer窗口中,Burg法需要选择的参数,除了fft长度Nfft,就只有阶数Order。试选择几个Order数运行,知道其可运行范围很大(Order=20000都可以)。反正我对于Order的具体选择也没有经验,那就编程作搜索式运算吧。

 

%用pburg法估计TW461512_0功率谱,搜索阶数参数Order
Nfft=65536;
M=zeros(Nfft/2+1,240);
for i=1:240
Order=50*i;
Pxx=pburg(TW461512_0,Order,Nfft);
M(:,i)=Pxx;
end
P_TW4615_b50_50_12000=M;%保存
M(151:(Nfft/2+1)-150,:)=[];
M=M';
surf(10*log10(M))

 

    运行,得:


        图6-25 用pburg法估计TW461512_0功率谱,搜索阶数参数Order
 

    上图数据太密集,图形太暗。

 

M(1:2:240,:)=[];%删去一部分数据
surf(10*log10(M))

 

        图6-25(1) 用pburg法估计TW461512_0功率谱,搜索阶数参数Order

 

    可以看出,低频段,阶数达到8、9000以上,高频段,阶数达到5、6000以上后,谱峰形状与位置都很稳定。这也符合我在用welch方法估计功率谱时所说的“我认为一种好的估计方法,当需要选择的参数落在一个大范围后,其估计结果应该对参数的变化是“不敏感”的”的特点。

    贴几张Spectrum Viewer窗口谱图:


        图6-26 用pburg法估计TW461512_0功率谱,阶数Order=10000
 

        图6-27 图6-26最低频段处,横轴放大128倍,纵轴放大为2倍

 


        图6-28 图6-26最高频段处,横轴放大128倍,纵轴放大为2倍

 

    从波形上来看,总的感觉,参数估计法比非参数估计法效果就是要好得多。 

 

    (8)、Yule AR法。此法与Burg法极其相似,在Spectrum Viewer窗口中需选择的参数完全一样。搜索阶数Order的程序只需将其谱估计函数“pburg”换成“pyulear”就可以了。

 

%用pyulear法估计TW461512 _ 0功率谱,搜索阶数参数Order。
Nfft=65536;
M=zeros(Nfft/2+1,120);
for i=1:120
Order=100*i;
Pxx=pyulear(TW461512_0,Order,Nfft);
M(:,i)=Pxx;
end
P_TW4615_y100_100_12000=M;%保存
M(151:(Nfft/2+1)-150,:)=[];
M=M';
surf(10*log10(M))

 

    运行,得:


        图6-29 用pyulear法估计TW461512_0功率谱,搜索阶数参数Order

 

    第一感觉,就是此图与图6-25也极其相似。再贴几张Spectrum Viewer窗口谱图:


        图6-30 用pyulear法估计TW461512_0功率谱,阶数Order=10000
 

        图6-31 图6-30最低频段处,横轴放大128倍,纵轴放大为2倍


        图6-32 图6-30最高频段处,横轴放大128倍,纵轴放大为2倍

    感觉上面三图与图6-26、图6-27及图6-28也极其相似。将图6-27与图6-31、图6-28与图6-32中的曲线分别放在同一窗口中进行比较。见下二图:


        图6-33 最低频段处。蓝线是Burg法、红线是Yule AR法、绿线是Welch法

 


        图6-34 最高频段处。其余情况同上图

 

    绿线Welch法谱估计参数是hanning窗,窗宽Nwind=25600,窗函数重叠点数Overlap=19200。

    可以看到,蓝线与红线重叠程度非常高。绿线作为一种非参数估计方法,其曲线形状变化步调与蓝线、红线能够保持这么高,已经很不简单了。

    上面就目前最常用的8种谱估计方法,对体温信号序列TW461512_0作了初步的功率谱估计。对这些谱估计的结果作进一步的分析处理,就放到以后的篇目中再做。

 

(本文首发于:http://blog.sina.com.cn/s/blog_6ad0d3de0100otj3.html

首发时间:2011-01-21 11:11:50)

 

 


https://blog.sciencenet.cn/blog-825323-636888.html

上一篇:5、互相关函数估计
下一篇:7、主要参考书目
收藏 IP: 14.153.189.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 07:36

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部