功率谱估计浅谈讲解

功率谱估计浅谈
摘要：介绍了几种常用的经典功率谱估计与现代功率谱估计的方法原理，并利用Matlab对随机信号进行功率谱估计，对两种方法做出比较，分别给出其优缺点。

关键词：功率谱；功率谱估计；经典功率谱估计；现代功率谱估计
前言
功率谱估计是从频率分析随机信号的一种方法,一般分成两大类:一类是经典谱估计;另一类是现代谱估计。

由于经典谱估计中将数据工作区以外的未知数据假设为零，这相当于数据加窗，导致分辨率降低和谱估计不稳定。

现代谱估计则不再简单地将观察区外的未知数据假设为零，而是先将信号的观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。

周期图、自相关法及其改进方法(Welch)为经典(非参数)谱估计方法, 其以相关和傅里叶变换为基础,对于长数据记录较适用,但无法根本解决频率分辨率低和谱估计稳定性的问题,特别是在数据记录很短的情况下,这一问题尤其突出。

以随机过程的参数模型为基础的现代参数法功率谱估计具有更高的频率分辨率和更好的适应性,可实现信号检测或信噪分离,对语音、声纳雷达、电磁波及地震波等信号处理具有重要意义,并广泛应用于通信、自动控制、地球物理等领域。

在现代参数法功率谱估计方法中,比较有效且实用的是AR模型法，Burg谱估计法,现代谱估计避免了计算相关,对短数据具有更强的适应性,从而弥补了经典谱估计法的不足,但其也有一些自身的缺陷。

下面就给出这两类谱估计的简单原理介绍与方法实现。

经典谱估计法
经典法是基于传统的傅里叶变换。

本文主要介绍一种方法：周期图法。

周期图法
由于对信号做功率谱估计，需要用计算机实现，如果是连续信号，则需要变换为离散信号。

下面讨论离散随机信号序列的功率谱问题。

连续时间随机信号的功率谱密度与自相关函数是一对傅里叶变换对，即：
()()j x x S R e d +∞
-Ω-∞Ω=⎰τττ
若()x R m 是()x R Ω的抽样序列，由序列的傅里叶变化的关系，可得
()()j j n x x m S e R m e ω
ω∞-=-∞=∑
即()j x S e ω与()x R m 也是一对傅里叶变换对。

显然，由序列傅里叶的频谱特性可知()j x S e ω是以2π为周期的。

而实际计算只能从离散随机信号序列x(n)的有限长(长度为N)的数据来对()j x S e ω与()x R m 进行估计。

设有限长离散序列为x(n),则：
1()[()*()]()[()]N N N x N N j x x R m x m x m N S e DFT R m =
-=ω 由DFT 的下列卷积特性：
若()[()]j N X e DFT X m ω=，则：
*()[()]j N X e DFT X m ω=-
从而：
1[()][()][()]N x N N DTFT R m DFT x m DFT x m N
=
- 即 211()()*()()N j j j j x S e X e X e X e N N
ωωωω== 综上所述，先用FFT 求出随机离散信号N 点的DFT ，再计算幅频特性的平方，然后除以N ，即得出该随机信号的功率谱估计。

由于这种估计方法在把()x R τ离散化的同时，使其功率谱周期化，故称之为“周期图法”，也称为经典谱估计法。

周期图法进行谱估计，是有偏估计，由于卷积的计算过程会导致功率谱真实值的尖峰附近产生泄漏，相对地平滑了尖峰值，因此造成谱估计的失真。

另外，当N →∞时，功率谱估计的方差不为零，所以不是一致性估计。

并且功率谱估计在ω等于2/N π整数倍的各数字频率点互不相关。

其谱估计的波动比较显著，特别是当N 越大、2/N π越小时，波动越明显。

但如果N 取得太小，又会造成分辨率的下降。

图1. 原始信号1
图2 原始信号1的功率谱估计图3. 原始信号2
图4. 原始信号2的功率谱估计
图5. 平均周期图法(4*256)
图6. 平均周期图法(重叠一半)
ππ，其中，图1所示的信号为sin(2**50*)2sin(2**120*)(1,)
xn t t randn N
=++
randn是正态分布随机数组，N为256，t是从0到1，dt为1/256。

图2为该信号的功率谱估计。

图2所示的信号为
ππ，其中，randn是正态分布随=++
xn t t randn N
sin(2**50*)2sin(2**120*)(1,)
机数组，N为1024，t是从0到1，dt为1/1024。

图4为该信号的功率谱估计。

图5是将图2所示的信号分为四段，每段的范围分别为
(1,256),(257,512),(513,768),(768,1024).每一道都没有重叠。

然后对分段分别作傅里叶变换，再把功率谱加起来做平均，得到图5。

图6是将图2所示的信号分为六段，分别为(1:256)，(129:384)，(257:512)，(385:640)，(513:768)，(641:896)，(769:1024)。

每两段之间都重合一半。

图1和图3相比，图1较为平滑，相应的，图1的功率谱也比较平滑。

图5和图6比，图6较为平滑，这是因为图6的谱是六段的平均。

对信号加入窗函数的话，功率谱的变化也是很明显的。

图7. 加入矩形窗原始信号和512点、1024点功率谱
图8.Bartlett平均周期图法
现代谱估计法
现代参数法功率谱估计方法中,比较有效且实用的是AR 模型法，Burg 谱估计法，在本文中介绍的是AR 模型法。

AR 模型法
经典谱的主要缺点是频率分辨率低。

这是由于周期图法在计算中把观测到的有限长的N 个数据以外的数据认为是零，这显然与事实不符。

如果把已观测到的数据估计出一白噪声激励，就不必认为N 个以外的数据全为零，就有可能克服经典谱估计的缺点。

一个实际中的随机过程总是可以用以下模型很好的表示： 11()1p
i i i p k
n k b z H z a z -=-==
+∑∑
当除0b 外的所有i b 均为零时的形式称为p 阶自回归模型即AR 模型，又称为全极点
模型。

当方差为2σ的白噪声通过AR 模型时，输出的功率谱密度为：
221()1xx p j k k P a e ω
ω
σω-==+∑
若已知参数12,,......p a a a 及2ωσ,就可以得到信号的功率谱估计。

它们之间是
Yule-Walker 方程。

解这个方程是一个复杂的数学问题，这里不做讨论。

图9. 原始信号3
图10. 自相关函数的无偏估计
图11. AR 模型求出的功率谱
图7所示的信号长度为155个点，长度较周期图法里的信号短，信号为sin(2**20*)sin(2**30*)(1/20)*(())X n n randn size n ππ=++，其中，n 为
0~155/100，间隔为1/100。

图8为自相关函数的无偏估计。

图9为功率谱。

从图9可以看出，这种方法具有极高的分辨能力。

只是在20和30处有两个峰值，在其它地方的值为零。

将信号改变成以下信号：
sin(2**1*)sin(2**2*)(1/20)*(())X F n F n randn size n =++ππ
则图11变成下图：
输出误差功率 0.0356 输出阶数P 2 （F1=20，F2=22）
输出误差功率 0.0485 输出阶数P 11 （F1=20，F2=25）
输出误差功率 0.0403 输出阶数P 9 （F1=20，F2=23）试验一下发现两个峰值频率为20和23的时候开始可以分辨出来。

与经典谱估计方法进行一个对比（基于第一个信号）：
相比与前面几种方法得到的结果来看，相差非常大，不仅仅是分辨率的提高，其余无效噪音或者说扰动都是非常小的。

结论
通过实验可以直接看出以下特性：
1)经典功率谱估计的方差大，谱分辨率差。

采用经典的傅里叶变换及窗口截断，对长序列有良好估计。

2)现代谱的分辨率较高。

这是由于在时域的开窗，使得在频域发生“频谱泄漏”，即功率谱的主瓣能量泄漏到旁瓣中，导致弱信号的主瓣被强信号的旁瓣所湮没，造成谱的模糊。

3)平均周期图法的收敛性较好，曲线平滑，但是功率谱主瓣较宽，分辨率低。

这是由于对随机序列的分段处理引起了长度有限所带来的 Gibbs 现象而造成的。

参考文献
1.刘志刚,李录明,赵冬梅. 现代谱估计法及应用效果[J]. 石油地球物理勘
探,2009,S1:5-9+167+9.
2.樊剑,刘铁,胡亮. 基于现代时频分析技术的地震波时变谱估计[J]. 振动与
冲击,2007,11:79-82+98+184.
3.蔡希玲,刘学伟,吕英梅,曹锡娜. 统计F-X谱估计方法及应用[J]. 勘探地球
物理进展,2008,03:181-186+163.
4.李刚. 宽带信号空间谱估计算法研究[D].哈尔滨工程大学,2007.
5.姚武川,姚天任. 经典谱估计方法的MATLAB分析[J]. 华中理工大学学
报,2000,04:45-47.
6.王凤瑛,张丽丽. 功率谱估计及其MATLAB仿真[J]. 微计算机信
息,2006,31:287-289.
7.王福杰,潘宏侠. MATLAB中几种功率谱估计函数的比较分析与选择[J]. 电子
产品可靠性与环境试验,2009,06:28-31.
附录：谱估计的Matlab实现
程序1：经典法谱估计
clf;
Fs=1000;
N=256;Nfft=256;%数据的长度和FFT所用的数据长度
n=0:N-1;t=n/Fs;%采用的时间序列
xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);
Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列
subplot(2,1,1),plot(f,Pxx);%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('周期图 N=256');grid on;
程序2：经典谱加窗分析
Fs=1000;%采样频率
n=0:1/Fs:1;%产生含有噪音的信号
xn=sin(2*pi*50*n)+2*sin(2*pi*120*n)+randn(size(n));
plot(n,xn);
window=boxcar(length(xn));%矩形窗
nfft=512;
[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法
figure(1)
plot(f,10*log10(Pxx));
window=boxcar(length(xn));%矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法
figure(2)
plot(f,10*log10(Pxx));
程序3：Bartlett平均周期图法
fs=1000;
n=0:1/fs:1;
xn=sin(2*pi*50*n)+2*sin(2*pi*120*n)+randn(size(n));
nfft=1024;
window=hamming(nfft);
noverlap=0;
p=0.9;
[Pxx,Pxxc]=psd(xn,nfft,fs,window,noverlap,p);
index=0:round(nfft/2-1);
k=index*fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxx(index+1));
figure(1)
plot(k,plot_Pxx)
figure(2)
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
程序4：现代谱估计法
clear all
%%%%%%%%%%%%%%%%%%%%%%%% 产生信号
FS=100; %设采样频率为100，则根据F1/FS=0.2，F2/FS=0.3 OR 0.25 ，可以得到F1，F2便于计算；
N=155; %数据长度改变数据长度会导致分辨率的变化；
n=0:1/FS:N/FS;
F1=20; %第一个SIN信号的频率；
F2=30; %第二个SIN信号的频率，取25或者30；X=sin(2*pi*F1*n)+sin(2*pi*F2*n)+(1/20)*randn(size(n)) ; %产生信号，由信噪比为10dB推出噪声功率；信号长度从X（1）到X（N+1）；
%%%%%%%%%%%%%%%%%%%%%%%%% 产生自相关函数数组
for m=1:N+1 %初始化R（m），R（m）用来存放自相关函数；由于R（0）在MATLAB里无效，所以从1开始到N+1；
R(m)=0;
end
for m=1:N+1 %做自相关函数的无偏估计；
S=0;
for n=1:N+2-m
H=X(n)*X(n+m-1);
S=S+H;
end
R(m)=S/N;
end %估计完毕
subplot(3,1,1);plot(X);subplot(3,1,2);plot(R); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Levinson算法主体部分
for k=1:N+1 %初始化后面算法中要用到的数组，其中a(k,k)用来存放AR模型参数，FPE(k)是最终预测误差准则函数，
a(k,k)=0; % U2(k)是AR模型中的另一参数即误差功率；
FPE(k)=0;
U2(k+1)=0;
end
U2(1)=R(1);
a(1,1)=-R(2)/R(1); %由Levinson算法，计算一阶模型参数
a11；
U2(2)=(1-(abs(a(1,1)))^2)*U2(1); %由Levinson算法，计算一阶模型参数误差功率；
FPE(1)=U2(2)*(N+2)/N; %计算一阶模型的最终预测误差准则函数；
S=0;
for k=2:N %由Levinson梯归公式计算K阶模型参数和FPE函数；
for l=1:k-1
M=a(k-1,l)*R(k-l+1);
S=S+M;
end
a(k,k)=-(R(k+1)+S)/U2(k);
for i=1:k-1
a(k,i)=a(k-1,i)+a(k,k)*a(k-1,k-i);
end
U2(k+1)=(1-(abs(a(k,k)))^2)*U2(k);
FPE(k)=U2(k+1)*(N+k+1)/(N-k+1);
S=0;
end %梯归计算完毕；
%% 确定阶数P
min=FPE(1); %求出使得FPE函数取最小值的阶数P；
for k=2:N
if FPE(k)<min
min=FPE(k);
p=k;
end
end
%p=2 %为了调整效果可以在这里自行指定阶数；disp('输出模型参数a');
for k=1:p
disp(a(p,k));
end
disp('输出误差功率');
disp(U2(p+1));
%%%%%%%%%%%%%%%%% AR模型参数确定后计算出功率谱
Z=0;
W=0:0.01:pi; %功率谱以2*pi为周期，又信号为实信号，只需输出0到PI即可；
for k=1:p
Z=Z+(a(p,k).*exp(-j*k*W));
end
PXX=U2(p+1)./((abs(1+Z)).^2); %得到功率谱函数；
F=W*FS/(pi*2); %将角频率坐标换算成HZ坐标，便于观察；subplot(3,1,3);plot(F,abs(PXX));
disp('输出阶数P')
disp(p);
grid on;。