基于Welch算法的经典功率谱估计的Matlab分析

摘 要:从经典功率谱估计周期图法原理入手,从理论上分析了其存在的局限性,借助Welch算法对其进行修正。依靠Matlab强大的数值分析和信号处理能力,进行实验仿真,比较不同的窗函数,不同的数据长度对Welch法谱估计质量的影响,并分析了造成这些影响的原因。

关键词:功率谱估计;周期图法;Welch算法;Matlab

中图分类号:TP911 文献标识码:A

文章编号:1004-373X(2010)03-007-03

Matlab Simulation Analysis of Power Spectrum Estimation Based on Welch Method

YI Xin,QU Aihua

(Naval Command College,Nanjing,211800,China)

Abstract:The paper mainly introduces the principles of Periodogram method of classical PSD estimation,analyzes the deficiency of Periodogram method in theory,and makes use of Welch to amend Perodogram method.By the use of simulation in Matlab,the impacts of different window function and different lenghth of data on estimation quality of Welch are discussed and the reasons of the impacts are analyzed.

Keywords:power spectrum estimation;periodogram method;Welch method;Matlab

0 引 言

随机信号在时间上是无限的,在样本上是无穷多,因此随机信号的能量是无限的,它应是功率信号。功率信号不满足傅里叶变换的绝对可积条件,因此严格意义上其傅里叶变换是不存在的。因此,对随机信号的频域分析,不再是简单的频谱,而是功率谱。

功率谱估计(PSD)是利用有限长的数据估计信号的功率谱,它涉及信号系统、随机信号分析、概率统计、随机过程、矩阵代数等一系列的基础学科,广泛应用于雷达、声纳、通信、地质勘探、天文、生物医学工程等众多领域[1]。功率谱估计可以分为经典谱估计和现代谱估计[2],本文主要研究经典谱估计中的Welch算法。目前经典谱估计中主要采用两种方法进行估计,分别是周期图法和自相关法[3],而Welch算法就是在周期图法的基础上进行改进得到的。

1 周期图法

Schuster于1899年首先提出了周期图这一概念,因为它是直接由傅里叶变换得到的,又称之为直接法。该方法是把随机信号x(n)的N点观察数据xN(n)视为一能量有限信号,直接取xN(n)的傅里叶变换,得xN(ω),然后再取其幅值的平方,并除以N,作为对x(n)真实的功率谱P(ejω)У墓兰,即:

(ω)=1NXN(ω)2=1N∑N-1n=0x(n)e-jωn2

(1)

由此可见,式(1)在n∞时,既不存在均值,也不存在极限,它只能看作是对真实谱做均值运算时的一个样本。缺少了统计平均,在记录的信号序列长度一定的条件下,要保证足够高的谱分辨率,谱估计的方差会很大,谱的正确性会很差。因此周期图的方差性能不好,而且,当数据长度N太大时,谱曲线呈现较大的起伏;当数据长度N太小时,谱的分辨率又不好。据此,直接法进行谱估计不满足一致性估计条件。因此,必须对周期图进行改进。

2 直接法谱估计的改进

为了改进周期图法的估计性能,常用的方法有两种:一是平均,就是对同一信号做多次周期图估计后再平均,在一定程度上弥补上述所缺的求均值运算。该方法事实上就是经典谱估计的间接法,本文不做讨论。┒是平滑,就是用适当的窗函数对谱进行平滑。其思想是把一长度为N的数据xN(n)分成L段,每段的长度为M,分别求每一段的功率谱,然后加以平均,以达到所希望的目的。

若对分段的数据加矩形窗[1,4],则第i段的数据变为:

xiN(n)=xN[n+(i-1)M]d1[n+(i-1)M],

0≤n≤M-1,1≤i≤L

(2)

由此得到修正后的周期图,即平均周期图:

P(ω)=1ML∑Li=1∑M-1n=0xiNe-jωn2

(3)

此方法就是Bartlett法,它很好地改善了直接法的方差特性,但是它是以牺牲偏差和分辨率为代价的。

Welch法是对Bartlett法的改进。主要改进在两个方面:一是在对xN(n)Х侄问,允许每段数据存在部分的交叠;二是每一段的数据窗口可以不是矩形窗口。这样可以改善由于矩形窗所造成的分辨率较差的影响。然后按照Bartlett法求每一段的功率谱,并对结果进行归一化,从而得到进一步修正的周期图,即:

(ω)=1MU∑Li=1∑M-1n=0xiNd2(n)e-jωn2

(4)

式中:U为归一化因子;d2(n)是数据窗口。

因为Welch算法各段允许交叠,从而增大了段数L,这样可以更好地改善方差特性。但是,数据的交叠又减小了每一段的不相关性,使方差的减小不会达到理论计算的程度。另外,选择合适的窗函数可以减小频谱的泄漏,改善分辨率。

3 Matlab仿真

为了分析Welch算法的性能,利用Matlab中提供的pwelch[5]函数实现Welch算法的功率谱估计。采样频率为1 000,采样点数为1 000,FFT点数为256,交叠数为20,窗函数采用矩形窗、海明窗和blackman窗[6-8],仿真程序如下:

Fs=1000;

n=0:1/Fs:1;

xn=2sin(2*pi*100*n)+4*cos(2*pi*150*n)+randn(size(n));

nfft=256;

[Pxx,f]=pwelch(xn,rectwin(100),20,nfft,Fs,′half′);%矩形窗

[Pxx1,f]=pwelch(xn,hamming(100),20,nfft,Fs,′half′);%海明窗

[Pxx2,f]=pwelch(xn,blackman(100),20,nfft,Fs,′half′);%blackman窗

xpsd=10*log10(Pxx);

xpsd1=10*log10(Pxx1);

xpsd2=10*log10(Pxx2);

figure(1)

plot(f,xpsd);

figure(2)

plot(f,xpsd1);

figure(3)

plot(f,xpsd2);

如图1所示,由矩形窗处理的谱估计的主瓣宽度最窄,分辨率最好,但是其旁瓣比其他窗函数的旁瓣要高,因此其正弦谱线附近的旁瓣泄漏比较严重,而且其起伏性较大,所以其方差特性最差[9]。由blackman窗和海明窗处理谱估计的主瓣宽度最宽,因此其分辨率相对较差,但其旁瓣较小,大大改善了由矩形窗处理的谱估计旁瓣较大所产生的谱失真。究其原因,选择不同的窗函数其主瓣宽度不一样,造成谱估计的分辨率也不相同;另外,选择不同的窗函数旁瓣的衰减速度也不相同,因而谱估计旁瓣的泄漏程度也不一样。

图1 长数据条件下不同窗函数的Welch算法谱估计

改变采样点数为100,交叠数为10,窗长度为50,其他参数不变。由于保持采样频率,降低采样点数,所以要对信号采取抽取,这里使用的是FIR滤波器,采用的是凯撒窗,其中对于矩形窗使用的是24阶,海明窗和blackman窗采用的是50阶。采样点数的大大降低,造成采样后的信号已经和原信号存在较大的误差,这时谱变得较平滑,但是各个窗的主瓣宽度变得更宽,造成分辨率降低。换句话说,减小数据长度,造成窗函数主瓣的宽度将变宽,从而减小谱的起伏。同时,由于信号不够尖锐,在功率较小的情况下就会淹没在噪声中。较大的边瓣掩盖了功率谱中的较弱的成分,或者产生假的峰值,这些情况在图2中已经出现。

图2 短数据长度条件下不同窗函数的Welch算法估计

4 结 语

综上所述,尽管Welch算法在短数据情况下存在局限性,但是由于其原理简单、实现容易,目前得到了广

泛的应用[10]。在实际的应用中,Welch算法对周期图

的平滑和平均是和窗函数的使用紧密相关联的。平滑和平均主要是用来改善周期图的方差性能,但往往又减小了分辨率和增大了偏差。不存在任何窗函数能够综合地改善估计谱的方差、偏差和分辨率等性能。因此,应该根据不同的信号、不同的处理目的合理地选择窗函数。

参考文献

[1]宋宁,关华.经典功率谱估计及其仿真[J].现代电子技术,2008,31(11):159-161.

[2]冯磊.经典功率谱估计与现代功率谱估计的对比[J].商业文化,2009(5):239-243.

[3]胡广书.数字信号处理[M].北京:清华大学出版社,2007.

[4]姚武川,姚天任.经典谱估计方法的Matlab分析[J].华中理工大学学报,2000,28(4):45-47.

[5]宁长春,陈天禄,索郎桑姆,等.数字信号处理中常用的Matlab工具箱函数简介[J].科技,2007(12):75-78.

[6]魏鑫,张平.周期图法功率谱估计中的窗函数分析[J].现代电子技术,2005,28(3):14-15.

[7]邵玉斌.Matlab/Simulink通信系统建模与仿真实例分析[M].北京:清华大学出版社,2008.

[8]范瑜,邬正义.功率谱估计的Welch方法中的窗函数研究[J].常熟高专学报,2000,14(7):36-39.

[9]瞿海雁,李鹂,钱小凌.如何在Matlab中优化基本周期图法对随机信号进行的功率谱估计[J].首都师范大学学报:自然科学版,2006,27(5):33-36.

[10]罗敏,刘嵩.基于Welch算法的功率谱估计的实现[J].北京工商大学学报:自然科学版,2007,25(3):58-59.