|||
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-5 图6-4最低频处。横轴放大16倍,纵轴放大1.2倍。
%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法估计功率谱。分段无重叠
由于数据太长,无法在一张图片中显示整个结果,只能分别取最低频最高频段各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))
运行,得:
什么也没有。所以在个人电脑上是无法使用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。
其中绿线nw=2,红线nw=3,蓝线nw=4。
图中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))
运行,得:
觉得还是很难确定nw到底取多大为好。改程序中权重算法因子为unity(一致权重的线性组合算法)、eigen(特征值权重的线性组合算法),运行,得:
发现上面三张图中,曲面形状非常一致。注:上面三张图及图6-7~图6-10的纵坐标应乘10才是分贝值,因为程序中“surf(log10(M))”应为“surf(10*log10(M))”。懒得再截图了,说明一下。
回到Spectrum Viewer窗口,取定nw=3,Nfft也不变,分别取Weights='adapt'、'unity'、'eigen',得到3条曲线。令其在同一窗口显示,如下图:
可以看出,最低频段处,3条曲线基本上重合了,后面的绿线、蓝线基本上都被前面的红线挡住。
最高频段处,绿线与红线曲线基本上重合了,蓝线在某些部位不与红线、绿线重合。
所以,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))
运行,得:
可以看出,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
(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))
运行,得:
上图数据太密集,图形太暗。
M(1:2:240,:)=[];%删去一部分数据
surf(10*log10(M))
图6-25(1) 用pburg法估计TW461512_0功率谱,搜索阶数参数Order
可以看出,低频段,阶数达到8、9000以上,高频段,阶数达到5、6000以上后,谱峰形状与位置都很稳定。这也符合我在用welch方法估计功率谱时所说的“我认为一种好的估计方法,当需要选择的参数落在一个大范围后,其估计结果应该对参数的变化是“不敏感”的”的特点。
贴几张Spectrum Viewer窗口谱图:
从波形上来看,总的感觉,参数估计法比非参数估计法效果就是要好得多。
(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-25也极其相似。再贴几张Spectrum Viewer窗口谱图:
图6-32 图6-30最高频段处,横轴放大128倍,纵轴放大为2倍
感觉上面三图与图6-26、图6-27及图6-28也极其相似。将图6-27与图6-31、图6-28与图6-32中的曲线分别放在同一窗口中进行比较。见下二图:
绿线Welch法谱估计参数是hanning窗,窗宽Nwind=25600,窗函数重叠点数Overlap=19200。
可以看到,蓝线与红线重叠程度非常高。绿线作为一种非参数估计方法,其曲线形状变化步调与蓝线、红线能够保持这么高,已经很不简单了。
上面就目前最常用的8种谱估计方法,对体温信号序列TW461512_0作了初步的功率谱估计。对这些谱估计的结果作进一步的分析处理,就放到以后的篇目中再做。
(本文首发于:http://blog.sina.com.cn/s/blog_6ad0d3de0100otj3.html
首发时间:2011-01-21 11:11:50)
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 07:17
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社