MATLAB实验

周期图法估计平稳随机信号功率谱的MATLAB 实现

分析：使用周期图法估计功率谱的MATLAB 函数为psd()（已被pwelch()代替了，实验中使用后者）：

welch 方法为改进的功率谱法，是用相干平均对周期图法功率谱估计进行二次处理，改善功率谱的估计方差，welch 方法对取样序列采取重叠分段方法，MATLAB 实现函数为pwelch()：

[pxx,w]=pwelch(x,window,noverlap,nfft)

[pxx,f]=pwelch(x,window,noverlap,nfft,fs)

各参数的含义和其余谱估计方法相似，window 默认为hamming 窗，使用前者为产生在w 下的谱密度，w 的单位为rad/sample ，实信号w 范围为[0,π]复信号为[0,2π]；后者fs 为采样频率，获得单位Hz 下的谱密度，f 的范围为[0,f/2]或

[0,fs]，判断与w 相同。无输出参数时，MATLAB 画出相应的频谱图。

假设()()()()n v n n n x ++=ππ4.0cos 35.0cos

其中，v(n)为均值0方差1的随机噪声信号，用welch 法估计其功率谱密度。在使用pwelch 时一定要注意w 或f 与pxx 的对应关系，当采样频率fs 改变时，w 及f 要做出相

应的调整。结果如下图六.1（对应程序为shiyan6.m）。

实验中可以看出，随分段长度的增加，即时域窗函数长度的增加，功率谱估计的方差特性变差，为改变周期图法功率谱估计方差特性差的缺点，Thomson提出了功率谱的多窗估计法（Multitaper Method），图六.2为用多窗估计法对上例中的x(n)做功率谱估计（对应程序仍为shiyan6.m）。

MATLAB中使用多窗估计法估计功率谱的函数为pmtm()：

[pxx,pxxc,w]=pmtm(x,nw,nfft,fs)

pxx为功率谱估计，nw为离散扁长球体序列（Discrete Prolate Spheroidal Sequences, DPSS, or Slepian Sequences），默认值为4，其余参数的含义与以上相同。实验中可看出，随窗函数的增加，方差特性得到进一步改善，但是频率分辨率下降。

00.10.20.30.40.50.60.7

0.80.91-20-10010

20Normalized Frequency(rad/sample)

P o w e r /F r e q u e n c y (d B /r a d /s a m p l e )00.10.20.30.40.50.60.70.80.91-20-10010

20Normalized Frequency(rad/sample)P o w e r /F r e q u e n c y (d B /r a d /s a m p l e )Power Spectral Density Estimate via Welch(modified periodogram method)

图.1

00.10.20.30.40.50.60.70.80.91-20-10010

20Normalized Frequency(rad/sample)P o w e r /F r e q u e n c y (d B /r a d /s a m p l e )Power Spectral Density Estimate via Thomson Multitaper 00.10.20.30.40.50.60.7

0.80.91-20-10010

20Normalized Frequency(rad/sample)

P o w e r /F r e q u e n c y (d B /r a d /s a m p l e )

图.2

附，源程序:

%%------------------------------------------------------------------------

%%功能：使用MA TLAB 对平稳随机信号以改进的功率谱法（welch 法）以及Thomson

%% 多窗估计法（Thomson Multitaper Method ）做功率谱估计 %%------------------------------------------------------------------------

fs=1;

n1=01fs1000;

n=n1(,1end-1);

ln=length(n);

vn=randn(1,ln);

xn=cos(0.35pin)+cos(0.4pin)+vn;

%%

%%以下用修正的周期图法welch法对平稳随机信号xn作功率谱分析%%Figuer 1 示意随分段长度（窗函数时域长度）增加，方差特性变差

figure(1);

[pxx,w]=pwelch(xn,128,64);

subplot(2,1,1);

plot(wpifs,10log10(abs(pxx)));

%横轴需根据采样频率sample作相应调整

%无输出参数时绘图的横轴为radsample，对应采样率为1时绘出正确图形，否则需调整横轴

axis([0 1 -20 20])

xlabel('Normalized Frequency(radsample)');

ylabel('PowerFrequency(dBradsample)');

title('Power Spectral Density Estimate via Welch(modified periodogram