功率谱估计仿真实验

功率谱估计仿真实验

选题条件：对于给定的一个信号()()()t t f t f t x ?ππ++=212sin 2)2sin(，其中1f =50Hz ，

2f =100Hz ，()t ?为白噪声，采样频率Fs 为1000Hz ，对其进行功率谱估 计。

仿真目标：采用多种方法对该指定信号进行功率谱估计，计算其功率谱密度，比较

各种估计方法的优劣。

设计思路：本仿真实验采用经典谱估计中的周期图法对给定信号进行谱估计。但是

由于其自身的缺陷，使得频率分辨率较低。为了不断满足需要，找到恰 当的估计法，实验使依次使用了周期图法的改进型方法如分段周期图法、 窗函数法以及修正的周期图法进行功率谱估计，对四种方法得出的谱估 计波形进行比较分析，得出估计效果最好的基于周期图法的谱估计方法。

仿真指标：频率分辨率、估计量的方差、频谱光滑度

平台说明：本实验采用MATLAB7.0仿真软件，基于WINDOWS-XP 系统。Matlab 是

一个集数值分析、矩阵运算、信号处理和图形显示于一体的工程分析处理软件。它提供的部分算法函数为功率谱估计提供了一条可行的方便途径，如PSD 和CSD 可以自动实现Welch 法估计，而不需要自己编程。但是较为有限，大部分需要自己编写相应的M 文件来实现。

实现方法： 一、周期图法

周期图法是直接将信号的采样数据()n x 进行傅立叶变换求功率谱密度估计。假设有限长随机信号序列()n x ，将它的功率谱按定义写出如下：

()()???

?????+=∑-=-∞→2121lim N N n n

j N j xx e n x N E e P ωω 如果忽略上式中求统计平均的运算，观测数据为：()n x 10-≤≤N n ，便得到了周期图法的定义：

()()2

10

^

1n

j N n j xx

e n x N e P ωω--=∑=， 式中的绝对值符号内的部分可以用FFT 计算，这样就可得到周期图法的计算框图如下所示：

()

ω

j xx e ^

图1 周期图法计算功率谱框图

采用周期图法时，可以分取不同的信号长度256、512和1024，分别进行功率谱

估计，并进行观察分析。仿真程序如下：

clf

Fs=1000;

N=256;

Nfft=256;

n=0:N-1;

t=n/Fs;

xn=sin(2*pi*50*t)+2*sin(2*pi*100*t)+randn(1,N); Pxx=10*log10(abs(fft(xn,Nfft).^2)/(N+1));

f=(0:length(Pxx)-1)*Fs/length(Pxx);

subplot(211)

plot(f,Pxx)

xlabel('Frequency (Hz)');

ylabel('Power spectrum (dB)');

title('Periodogram N=256')

grid

程序运行结果如下图所示：

a.N=256

b.N=512

c.N=1024

图2 周期图法功率谱估计N 分别为256、512、1024

从图2可以看出，在频率50Hz 和100Hz 处，功率谱有两个峰值，说明信号中含有50Hz 和100Hz 的周期成分，这点与实际信号相吻合。功率谱密度在很大范围波动，随着信号取样点数由256增加为1024，摆动的幅度并未减小，只是摆动的频率加快，功率谱估计效果并没有什么改进。

用有限长样本序列的周期图法来表示随机序列的功率谱虽然只是一种估计或近似，不可避免地存在误差，为了减小误差，使功率谱估计更加平滑，可以采用以下方法进行改进。 二、平均周期图法

将信号序列()n x ，10-≤≤N n ，分成互不重叠的L 个小段，每个小段有m 个采样值，则Lm=N 。对每小段信号序列进行功率谱估计，第i 组的周期图用下式表示：

()()210

1

∑-=-=

M n n

j i

i e

n x M

I ωω。然后求他们的平均值作为整个序列()n x 的功率谱估计，公

式如下：()

()ωω

∑==L

i i j xx I L e P 1

^

1

算法框图如下：

图3 分段周期图法框图

本仿真实验中可以自行设计分段数分别为2、4、8段，只需将仿真代码中的分段数进行调整即可实现。仿真程序设计如下（分四段）： clf

Fs=1000; N=1024; Nsec=256; n=0:N-1;

t=n/Fs;

xn=sin(2*pi*50*t)+2*sin(2*pi*100*t)+randn(1,N); pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;

pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;

pxx3=abs(fft(xn(513:768),Nsec).^2)/Nsec;

pxx4=abs(fft(xn(769:1024),Nsec).^2)/Nsec;

Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);

f=(0:length(Pxx)-1)*Fs/length(Pxx);

subplot(211)

plot(f,Pxx)

xlabel('Frequency(Hz)');

ylabel('Power Spectrum(dB)');

title('Averaged Periodogram(no overlap)N=2*512') grid

程序运行结果如图4所示：

a．分段数L=2

b.分段数L=4

c.分段数L=8

图4 分段平均周期图法功率谱估计

图4中，分别采用了不同的分段数2、4、8，从图中可以清楚地看到，随着分段数的增加功率谱曲线越来越平滑，功率谱估计值在0dB 附近摆动的幅度越来越小。

但是由于数据量N=1024是个定值，段数加大，每一段的数据量必然减少，因此估计量方差减小了，使偏移加大，分辨率降低。因此，估计量的方差和分辨率是一对矛盾，它们的效果可以互换，可以根据实际情况适当地选择L 和M 。如果对分辨率要求不高，可以取L 大些；反之，只好将M 的值取得大些。

图4与图2相比，谱估计效果有了明显改善。 三、窗函数法

窗函数法是使用一种适当的功率谱窗函数()ωj e W 与周期图进行卷积，来达到使周期图平滑的目的，如下式所示：

()

()()

()

θθπθωππωd e

W I e j N j xx P --?=21^

式中，()()n

j N N m N e

m xx I r ωω----=∑=

1

)

1(^

，

()m xx

r ^

是有偏自相关函数。周期图和频谱函数

卷积得到功率谱，等效于在频域对周期图进行修正，使周期图通过一个线性非频变

系统，滤除掉周期图中的快变成分。计算框图如下：

图5 窗函数法框图

仿真程序设计如下： clf

Fs=1000;

N=1024; Nsec=256; n=0:N-1; t=n/Fs;

w=hanning(256)';

xn=sin(2*pi*50*t)+2*sin(2*pi*100*t)+randn(1,N); pxx1=abs(fft(w.*xn(1:1024),Nsec).^2)/norm(w)^2; subplot(211) plot(f,Pxx)

Pxx=10*log10(conv(pxx1，w)); xlabel('Frequency(Hz)');

ylabel('Power Spectrum (dB)');

title('Averaged Modified Periodogram (no overlop)N=4*256') grid

程序运行结果如图6所示：

图6 窗函数法功率谱估计

分析：与图2、图4相比，图6的功率谱曲线更加光滑，主瓣宽度比较宽，估计误差变小了，但是偏移加大了，使分辨率降低。这点可以从窗函数基本知识可以得到，采用合适窗函数对信号进行处理可以减少频谱泄漏，同时可增加频峰宽度。分辨率和估计方差两者之间的矛盾还是比较明显。为了折中两者之间的矛盾，可以采用修正的周期图求平均法。

四、Welch 修正的周期图求平均法

Welch 算法是由Welch 提出的修正周期图法，是经典谱估计中获得的一项有效的算法。Welch 算法谱估计采取数据分段加窗处理再求平均的办法，把窗函数加到每一个数据段上，求出每一段的周期图，形成修正的周期图，再对每一个修正的周期

图进行平均。第i 段的修正周期图为()()()210

1

∑-=-=

M n j i

i e

n n x U

I ω

ωωi=1，2，3…M-1

式中()n M

U M n ∑-==

10

21

ω，将每一段的修正的周期图之间看成互不相关，最后的功率谱

估计为()

()ωω

∑==L

i i j xx I L e P 1

^

1

Welch 法谱估计流程图如下图所示：

图7 Welch 修正的周期图法框图

仿真程序设计如下：

clf

N=1024;

Nfft=256; Fs=1000; n=0:N-1; t=n/Fs;

window=hanning(256); noverlap=128; dflag='none';

xn=sin(2*pi*50*t)+2*sin(2*pi*100*t)+randn(1,N); Pxx=psd(xn,Nfft,Fs,window,noverlap,dflag); f=(0:Nfft/2)*Fs/Nfft; subplot(211)

plot(f,10*log(Pxx))

xlabel('Frequency (Hz)');

ylabel('Power Spectrum(dB)'); title('PSD----Welch Wethod') grid

程序运行结果如图8所示：

图8 修正的周期图求平均法功率谱估计

从图7可以看出，谱波形更加光滑，摆动幅度较小。由于加了hanning窗，大大压低了旁瓣宽度，使得低电平信号清晰可见，但由于主瓣宽度加宽，功率谱波峰变宽了，从而降低了信号的分辨率。与前几种估计的波形相比，Welch修正的周期图法所得到的标准方差比其他几种周期图法要小，这说明经过分段、加窗后方差也会减小，从而说明经过加窗平滑方法后的周期图估计也越来越正确。

五．结论：

通过仿真实验的波形可以直观地看出以下特性：

（1）平均周期图法、窗函数法以及修正的周期图法的收敛性较好，曲线较周期图法更为光滑，估计的结果方差较小。但是功率谱主瓣较宽，分

辨率较低。这是由于对随机序列的分段处理引起了长度有限所带来的

问题，由于只有N个观测数据，观测不到的信号被认为是0。对于N以外

的数据，信号仍有较大的相关性，这样估计出的功率谱就会有很大的

偏差。

（2）窗函数法和修正的周期图法与周期图法和平均周期图法相比，谱估值比较平滑，但是分辨率较差。其原因是给每一段序列用适当的窗函数

加权后，在得到平滑的估计结果的同时，使功率谱的主瓣变宽，因此

分辨率有所下降。

由于本人对谱估计的基础知识掌握的不是很扎实，对各种估计方法的原理掌握的不够透彻，在进行仿真实验时，对可能发生的各种情况估计不周，对一些关键参数的作用还不是很了解，种种原因导致实验进行的不顺利，对实验结果的分析也没有做到完全正确，可能与实际存在一些出入，希望在以后的学习生活中能够不断进步，深层次地把握各种功率谱估计方法之精髓。

总之，本实验中所涉及的这几种传统的功率谱估计方法无论采用哪一种改进方法，总是以减少分辨率为代价，换取估计方差的减少，提高分辨率的方法无法从根本上解决。