c语言编译功率谱密度函数,科学网—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-7 TW461512_0的Welch法估计功率谱。分段无重叠