日志

6、功率谱密度函数估计（一）

已有 1774 次阅读2011-4-14 16:19 |个人分类:研究|

6、功率谱密度函数估计（一）

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

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、功率谱密度函数估计（一）

图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、功率谱密度函数估计（一）

图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