原文出自:http://www.cnblogs.com/jacklu/p/5140913.html

功率谱估计在分析平稳各态遍历随机信号频率成分领域被广泛使用,并且已被成功应用到雷达信号处理、故障诊断等实际工程中。本文给出了经典功率谱估计的几类方法,并通过Matlab的实验仿真对经典功率谱估计方法性能进行了分析,最后说明了经典功率谱估计法的局限性和造成这种局限性的原因。

1.引言

给定一个标准的正弦信号,我们可以通过傅里叶变换来分析它的频率成分。然而,实际工程应用中,由于存在着各种干扰、噪声,我们得到的信号往往不是理想的,如图1-1这种信号,具有不确定性,幅度不能预知,非周期,但往往服从一定的统计特性,这种信号叫作随机信号。需要注意的是,本文所说的随机信号是指平稳各态遍历的随机信号,关于非平稳随机信号的分析方法本文不予讨论。

经典功率谱估计及Matlab仿真-LMLPHP

图1-1 一种随机信号时域形式

对于图1-1的随机信号,我们可以通过功率谱来分析它的频率成分,如图1-2所示为图1-1随机信号的功率谱。实际过程中,我们只能获得随机信号的一些离散数据点(假设为N个),本文将讨论如何利用这N个数据点,来得到一个"非精确"的功率谱来对真实随机信号的功率谱进行估计,并讨论如何更好的估计,即在下一章要讲述的几个经典的功率谱估计法。

经典功率谱估计及Matlab仿真-LMLPHP经典功率谱估计及Matlab仿真-LMLPHP

图1-2 上图所示的随机信号功率谱

2.经典功率谱估计法

上一章我们已经知道功率谱估计法是通过利用已经获得的N个数据点,来得到一个"非精确"的功率谱对真实随机信号的功率谱进行估计,所以在给出具体的方法之前,如何来评价我们得出的这个"非精确"的功率谱的好坏呢?

评价功率谱性能好坏的标准有很多,本文只给出两个影响最大的标准:分辨率和方差。分辨率即功率谱上能够区分的最小相邻频率成分,分辨率越高,我们观察信号的频率成分越清晰;方差大小则反映到功率谱波动性的大小,如果方差太大,功率谱波动性大,则很容易造成有用的频率成分被噪声淹没。所以,我们希望得到的这个"非精确"的功率谱,分辨率越高越好,方差越小越好。

同时,我们给出概率论与数理统计中所学的一致估计和非一致估计的概念,假定真实信号的功率谱为经典功率谱估计及Matlab仿真-LMLPHP,估计得到的"非精确"功率谱为经典功率谱估计及Matlab仿真-LMLPHP,如果满足公式(2-1),则称为一致估计。

经典功率谱估计及Matlab仿真-LMLPHP                                         (2-1)

2.1周期图法

2.1.1周期图法原理

我们已知N个离散的数据点经典功率谱估计及Matlab仿真-LMLPHP,对这些数据点进行傅里叶变换,得到式(2-2):经典功率谱估计及Matlab仿真-LMLPHP (2-2)

再对(2-2)式取模的平方,除以N,即得到一个"非精确"的谱,如式(2-3)所示,这就是周期图法的原理。

经典功率谱估计及Matlab仿真-LMLPHP (2-3)

2.1.2周期图法性能(Matlab仿真)

上一小节我们已经给出了周期图法的原理。本节将通过Matlab仿真给出数据点数N对功率谱性能好坏的影响,正如上文所述,将通过对所得功率谱的分辨率和方差两方面进行分析。

我们在Matlab中通过三个正弦函数和白噪声叠加,构造了一个随机信号。其数学形式如式(2-4),其中频率经典功率谱估计及Matlab仿真-LMLPHP经典功率谱估计及Matlab仿真-LMLPHP经典功率谱估计及Matlab仿真-LMLPHP,幅值经典功率谱估计及Matlab仿真-LMLPHP经典功率谱估计及Matlab仿真-LMLPHP经典功率谱估计及Matlab仿真-LMLPHP,相位经典功率谱估计及Matlab仿真-LMLPHP为相互独立在经典功率谱估计及Matlab仿真-LMLPHP上服从均匀分布随机相位,经典功率谱估计及Matlab仿真-LMLPHP为均值为0,方差为1的实值高斯白噪声,采样频率为1000。信号的时域形式如图2-1所示。

经典功率谱估计及Matlab仿真-LMLPHP (2-4)

经典功率谱估计及Matlab仿真-LMLPHP

图2-1 实验所用的随机信号

当数据点数N分别为128、256、512和1024时,得到的功率谱分别如图2-2、图2-3、图2-4和图2-5所示。分辨率能够直观的通过功率谱图形看出,方差的数值由表2-1给出。

经典功率谱估计及Matlab仿真-LMLPHP

图2-2 N=128时周期图法得到的功率谱

经典功率谱估计及Matlab仿真-LMLPHP

图2-3 N=256时周期图法得到的功率谱

经典功率谱估计及Matlab仿真-LMLPHP

图2-4 N=512时周期图法得到的功率谱

经典功率谱估计及Matlab仿真-LMLPHP

图2-5 N=1024时周期图法得到的功率谱

表2-1 不同N值得到功率谱的方差值

N

128

256

512

1024

方差

92.7108

130.9109

160.9187

483.5894

通过上面实验结果的比较,我们很容易发现,周期图法得到的功率谱随着数据点数N的增大,分辨率变大、方差变也大。

2.1.3平均周期图法

周期图法得到的功率谱与我们所期望的"分辨率大、方差小"是矛盾的。为了进一步降低方差,将N个观测样本数据点经典功率谱估计及Matlab仿真-LMLPHP分为L段,每段数据长度为M, 分别对每段数据求周期图功率谱估计,然后求平均值,这种方法称平均周期图法。

那么这种方法会如何改善方差呢?我们给出证明:经典功率谱估计及Matlab仿真-LMLPHP (2-5)

其中:

经典功率谱估计及Matlab仿真-LMLPHP

经典功率谱估计及Matlab仿真-LMLPHP

(2-6)

由式(2-5)我们可以看出,平均周期图法将原来的方差变为原来的经典功率谱估计及Matlab仿真-LMLPHP,L为分段数。

2.1.4平均周期图法性能(Matlab仿真)

当数据点数N为1024,分段数分别为8、4、2时,平均周期图法得到的功率谱分别如图2-6、图2-7、图2-8所示。分辨率能够直观的通过功率谱图形看出,方差的数值由表2-2给出。

经典功率谱估计及Matlab仿真-LMLPHP

图2-6 L=8时平均周期图法得到的功率谱

经典功率谱估计及Matlab仿真-LMLPHP

图2-7 L=4时平均周期图法得到的功率谱

经典功率谱估计及Matlab仿真-LMLPHP

图2-8 L=2时平均周期图法得到的功率谱

表2-2 不同L值得到功率谱的方差值

L

8

4

2

1

方差

96.3756

190.9647

400.6464

483.5894

L=1时,平均周期图法退化为周期图法。通过上面实验结果的比较,我们很容易发现,平均周期图法得到的功率谱随着分段数L变大,方差变小,但分辨率变小。

当观测样本序列数据个数N固定时,要降低方差需要增加分段数L。当N不大时分段长度M取值较小,则功率谱分辨率降低到较低的水平。若分段数L固定时,增加分辨率需要分段长度M,则需要采集到更长的检测数据序列。实际中恰恰是检测样本序列长度不足。

2.1.5修正的平均周期图法

上一节已经提到实际中检测样本序列长度是有限的。对现有数据长度N,如果能获得更多的段数分割,将会得到更小的方差。允许数据段间有重叠部分,来得到更多的段数。对段间重叠长度的选取,最简单是取为段长度M的一半。由式(2-5)可知更多的段数可以进一步降低方差。

数据截断的过程中相当于数据加矩形窗,矩形窗幅度较大的旁瓣会造成"频谱泄漏"。我们分段时采取的窗函数更为多样(三角窗,海明窗等), 以减小截断数据(加矩形窗)窗函数带来的影响

2.1.6修正的平均周期图法性能(Matlab仿真)

利用修正平均周期图法,分别使用矩形窗、Blackman窗和Hamming窗得到的功率谱如图2-9所示。

经典功率谱估计及Matlab仿真-LMLPHP经典功率谱估计及Matlab仿真-LMLPHP经典功率谱估计及Matlab仿真-LMLPHP

图2-9 不同窗函数的修正平均周期图法得到的功率谱

可以发现,矩形窗的分辨率最高,但是方差也最大,这是由于矩形窗频谱主瓣最窄,分辨率因此最高,旁瓣也高,导致频谱泄漏最严重,方差最大。

2.1.7总结

周期图法获得的功率谱随着样本点数越多,分辨率越大、方差越大;平均周期图法以牺牲分辨率来进一步改善方差;修正的平均周期图法允许段的重叠来进一步增大分段数、或者分段数相同,每段样本点数变多。无论是哪种方法都没有彻底结局方差与分辨率之间的矛盾。

2.2相关功率谱估计法-BT法

正如我们之前介绍的,要提高功率谱估计的分辨率,必须增加数据序列的长度N,但是较长的数据序列,由噪声引起的随机性得到更为充分的体现-较大的方差。事实上,当N无穷大时,方差为一非零常数。即周期图法无法实现功率谱的一致估计。而这节讲述的相关功率谱估计法(下文称作BT法),是一致估计。

2.2.1 BT法的原理

维纳辛钦定理指出,随机信号的相关函数与它的功率谱是一对傅里叶变换对。BT法就是基于这个原理。先由观测数据经典功率谱估计及Matlab仿真-LMLPHP估计出自相关函数,然后求自相关函数的傅立叶变换,以此变换作为对功率谱的估计,也称为间接法。BT法要求信号长度N以外的信号为零,这也造成BT法的局限性。

经典功率谱估计及Matlab仿真-LMLPHP的自相关函数经典功率谱估计及Matlab仿真-LMLPHP定义如式(2-7)所示,得到的功率谱记为经典功率谱估计及Matlab仿真-LMLPHP,则BT法可以表述为式(2-8)。

经典功率谱估计及Matlab仿真-LMLPHP 经典功率谱估计及Matlab仿真-LMLPHP (2-7)

经典功率谱估计及Matlab仿真-LMLPHP (2-8)

2.2.2 BT法的性能(Matlab仿真)

数据点数N分别为128、256、512和1024的BT法,得到的功率谱如图2-10、图2-11、图2-12和图2-13所示。

经典功率谱估计及Matlab仿真-LMLPHP

图2-10 N=128时,BT法得到的功率谱

经典功率谱估计及Matlab仿真-LMLPHP

图2-11 N=256时,BT法得到的功率谱

经典功率谱估计及Matlab仿真-LMLPHP

图2-12 N=512时,BT法得到的功率谱

经典功率谱估计及Matlab仿真-LMLPHP

图2-13 N=1024时,BT法得到的功率谱

由上面实验可以发现,M随着N的增大而增大时,分辨率提高,方差变大。BT法仍然没有解决分辨率与方差之间的矛盾,但是BT法得到的功率谱当N为无穷大时,方差会趋向于零,即为一致估计。

2.2.3 周期图法与BT法的关系

相关函数经典功率谱估计及Matlab仿真-LMLPHP可以写成如式(2-9)的卷积形式:经典功率谱估计及Matlab仿真-LMLPHP (2-9)

设序列经典功率谱估计及Matlab仿真-LMLPHP的傅立叶变换为经典功率谱估计及Matlab仿真-LMLPHP,则当M=N-1时,功率谱的估计可表示为式(2-10)的形式。可以看出周期图法可以看作BT法在取M=N-1时的特例。

经典功率谱估计及Matlab仿真-LMLPHP(2-10)

结 论

本文通过Matlab仿真,以一个具体的随机信号为例,简单介绍了周期图法、平均周期图法、修正的平均周期图法以及BT法的基本原理,并对这些方法的性能进行分析。可以看出,无论是周期图法及其改进算法还是BT法都没有从根本上解决分辨率与方差的矛盾。经典功率谱估计是利用傅里叶变换估计功率谱,而我们之前分析随机信号不满足傅里叶变换的条件,所以经典功率谱估计方法不得不从无限长数据点截取有限长数据点,加入限制条件(周期图法实际上假定N点外数据周期重复、BT法假定N点外数据为零)来"强制"作傅里叶变换,这也是造成它局限性的原因。

参考资料

[1] 朱哲,钟宏伟. 非平稳随机信号分析处理方法研究[J] 安徽电子信息技术学院学报 2008.6:28-28

[2] 皇甫堪. 现代数字信号处理[M]. 电子工业出版社

部分matlab程序代码:

周期图法:(by宋同学)

 Fs=;
f1=;
f2=;
f3=;
N=;
Nfft=N;
n=:N-;
t=n/Fs;%采用的时间序列
xn=cos(*pi*f1*t)+1.5*cos(*pi*f2*t)+cos(*pi*f3*t)+1.5*randn(size(n));
figure;
plot(n,xn);grid on;title('时域信号');
P=*log10(abs(fft(xn,Nfft).^)/N);%Fourier振幅谱平方的平均值,并转化为dB
f=(:length(P)-)*Fs/length(P);%给出频率序列
figure;
plot(f(:N/),P(:N/));grid on;title('功率谱(dB图)');ylabel('功率谱/dB');
xlabel('频率/Hz');
Pxx=abs(fft(xn,Nfft).^)/N;%Fourier振幅谱平方的平均值,并转化为dB
f=(:length(Pxx)-)*Fs/length(Pxx);%给出频率序列
figure;
plot(f(:N/),Pxx(:N/));%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱');
title('周期图 N=128');grid on;
std(Pxx)^ N=;
Nfft=N;
n=:N-;
t=n/Fs;%采用的时间序列
xn=cos(*pi*f1*t)+1.5*cos(*pi*f2*t)+cos(*pi*f3*t)+1.5*randn(size(n));
figure;
plot(n,xn);grid on;title('时域信号');
P=*log10(abs(fft(xn,Nfft).^)/N);%Fourier振幅谱平方的平均值,并转化为dB
f=(:length(P)-)*Fs/length(P);%给出频率序列
figure;
plot(f(:N/),P(:N/));grid on;title('功率谱(dB图)');ylabel('功率谱/dB');
xlabel('频率/Hz');
Pxx=abs(fft(xn,Nfft).^)/N;%Fourier振幅谱平方的平均值,并转化为dB
f=(:length(Pxx)-)*Fs/length(Pxx);%给出频率序列
figure;
plot(f(:N/),Pxx(:N/));%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱');
title('周期图 N=256');grid on;
std(Pxx)^ N=;
Nfft=N;
n=:N-;
t=n/Fs;%采用的时间序列
xn=cos(*pi*f1*t)+1.5*cos(*pi*f2*t)+cos(*pi*f3*t)+1.5*randn(size(n));
figure;
plot(n,xn);grid on;title('时域信号');
P=*log10(abs(fft(xn,Nfft).^)/N);%Fourier振幅谱平方的平均值,并转化为dB
f=(:length(P)-)*Fs/length(P);%给出频率序列
figure;
plot(f(:N/),P(:N/));grid on;title('功率谱(dB图)');ylabel('功率谱/dB');
xlabel('频率/Hz');
Pxx=abs(fft(xn,Nfft).^)/N;%Fourier振幅谱平方的平均值,并转化为dB
f=(:length(Pxx)-)*Fs/length(Pxx);%给出频率序列
figure;
plot(f(:N/),Pxx(:N/));%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱');
title('周期图 N=512');grid on;
std(Pxx)^ N=;
Nfft=N;
n=:N-;
t=n/Fs;%采用的时间序列
xn=cos(*pi*f1*t)+1.5*cos(*pi*f2*t)+cos(*pi*f3*t)+1.5*randn(size(n));
figure;
plot(n(:),xn(:));grid on;title('时域信号');
P=*log10(abs(fft(xn,Nfft).^)/N);%Fourier振幅谱平方的平均值,并转化为dB
f=(:length(P)-)*Fs/length(P);%给出频率序列
figure;
plot(f(:N/),P(:N/));grid on;title('功率谱(dB图)');ylabel('功率谱/dB');
xlabel('频率/Hz');
Pxx=abs(fft(xn,Nfft).^)/N;%Fourier振幅谱平方的平均值,并转化为dB
f=(:length(Pxx)-)*Fs/length(Pxx);%给出频率序列
figure;
plot(f(:N/),Pxx(:N/));%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱');
title('周期图 N=1024');grid on;
std(Pxx)^

平均周期图法

 clear;
Fs=;
f1=;
f2=;
f3=;
n=:/Fs:;
xn=cos(*pi*f1*n)+1.5*cos(*pi*f2*n)+cos(*pi*f3*n)+1.5*randn(size(n)); N=;Nsec=;%数据的长度和FFT所用的数据长度
Pxx1=abs(fft(xn(:),Nsec).^)/Nsec; %第一段功率谱
Pxx2=abs(fft(xn(:),Nsec).^)/Nsec;%第二段功率谱
Pxx=*log10((Pxx1+Pxx2)/);%Fourier振幅谱平方的平均值,并转化为dB
(std((Pxx1+Pxx2)/))^
f=(:length(Pxx)-)*Fs/length(Pxx);%给出频率序列
figure;
plot(f(:Nsec/),Pxx(:Nsec/));%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('N=2*512');grid on; N=;Nsec=;%数据的长度和FFT所用的数据长度
Pxx1=abs(fft(xn(:),Nsec).^)/Nsec; %第一段功率谱
Pxx2=abs(fft(xn(:),Nsec).^)/Nsec;%第二段功率谱
Pxx3=abs(fft(xn(:),Nsec).^)/Nsec;%第三段功率谱
Pxx4=abs(fft(xn(:),Nsec).^)/Nsec;%第四段功率谱
Pxx=*log10((Pxx1+Pxx2+Pxx3+Pxx4)/);%Fourier振幅谱平方的平均值,并转化为dB
std((Pxx1+Pxx2+Pxx3+Pxx4)/)^
f=(:length(Pxx)-)*Fs/length(Pxx);%给出频率序列
figure;
plot(f(:Nsec/),Pxx(:Nsec/));%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('N=4*256');grid on; N=;Nsec=;%数据的长度和FFT所用的数据长度
Pxx1=abs(fft(xn(:),Nsec).^)/Nsec; %第一段功率谱
Pxx2=abs(fft(xn(:),Nsec).^)/Nsec;%第二段功率谱
Pxx3=abs(fft(xn(:),Nsec).^)/Nsec;%第三段功率谱
Pxx4=abs(fft(xn(:),Nsec).^)/Nsec;%第四段功率谱
Pxx5=abs(fft(xn(:),Nsec).^)/Nsec; %第五段功率谱
Pxx6=abs(fft(xn(:),Nsec).^)/Nsec;%第六段功率谱
Pxx7=abs(fft(xn(:),Nsec).^)/Nsec;%第七段功率谱
Pxx8=abs(fft(xn(:),Nsec).^)/Nsec;%第八段功率谱
Pxx=*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8)/);%Fourier振幅谱平方的平均值,并转化为dB
std((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8)/)^ f=(:length(Pxx)-)*Fs/length(Pxx);%给出频率序列
figure;
plot(f(:Nsec/),Pxx(:Nsec/));%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('N=8*128');grid on;

修正的平均周期图法(by周同学)

 clear;%求1024点功率谱以及方差
Fs=;
f1=;
f2=;
f3=;
n=:/Fs:;
xn=cos(*pi*f1*n)+1.5*cos(*pi*f2*n)+cos(*pi*f3*n)+1.5*randn(size(n));
M=;N=;%数据总点数1024,每段长度M
window=hamming(M);
Pxxtotal=;
L=N/(M/)-;
for i=::(L-)
m=abs(fft(window'.*(xn((M/2+M/2*i-M+1):(M/2+M/2*i))),M).^2)/M;
Pxxtotal=Pxxtotal+m;
end
window=hamming(Fs-(N-M+1)+1);
mend=abs(fft(window'.*(xn((N-M+):Fs)),M).^)/M;
Pxxtotal=(Pxxtotal+mend)/L;
Pxx=*log10((Pxxtotal));%Fourier振幅谱平方的平均值,并转化为dB
f=(:length(Pxx)-)*Fs/length(Pxx);%给出频率序列
w1=fft(window,N);
w10=abs(fftshift(w1)); plot(f(:M/),Pxx(:M/));%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
grid on;
B=var(Pxxtotal)

BT法(by阳同学)

 Fs=;
n=:/Fs:;
x=cos(*pi**n)+1.5*cos(*pi**n)+cos(*pi**n)+1.5*randn(size(n));
nfft=;
ncxk=*nfft/;
cxn=xcorr(x,'unbiased');
CXk=fft(cxn,ncxk);
Pxx=abs(CXk);
index=:round(ncxk/-);
k=index*Fs/ncxk;
C=Pxx(index+);
P=(Pxx(index+));
plot(k,P);grid on
var(C)
title('BT法功率谱估计 N=1024');
xlabel('频率 Hz');
ylabel('功率');
04-16 11:13