原文出自：http://www.cnblogs.com/jacklu/p/5140913.html
功率谱估计在分析平稳各态遍历随机信号频率成分领域被广泛使用，并且已被成功应用到雷达信号处理、故障诊断等实际工程中。本文给出了经典功率谱估计的几类方法，并通过Matlab的实验仿真对经典功率谱估计方法性能进行了分析，最后说明了经典功率谱估计法的局限性和造成这种局限性的原因。
1.引言
给定一个标准的正弦信号，我们可以通过傅里叶变换来分析它的频率成分。然而，实际工程应用中，由于存在着各种干扰、噪声，我们得到的信号往往不是理想的，如图1-1这种信号，具有不确定性，幅度不能预知，非周期，但往往服从一定的统计特性，这种信号叫作随机信号。需要注意的是，本文所说的随机信号是指平稳各态遍历的随机信号，关于非平稳随机信号的分析方法[1]本文不予讨论。
图1-1 一种随机信号时域形式
对于图1-1的随机信号，我们可以通过功率谱来分析它的频率成分，如图1-2所示为图1-1随机信号的功率谱。实际过程中，我们只能获得随机信号的一些离散数据点（假设为N个），本文将讨论如何利用这N个数据点，来得到一个"非精确"的功率谱来对真实随机信号的功率谱进行估计，并讨论如何更好的估计，即在下一章要讲述的几个经典的功率谱估计法。
图1-2 上图所示的随机信号功率谱
2.经典功率谱估计法
上一章我们已经知道功率谱估计法是通过利用已经获得的N个数据点，来得到一个"非精确"的功率谱对真实随机信号的功率谱进行估计，所以在给出具体的方法之前，如何来评价我们得出的这个"非精确"的功率谱的好坏呢？
评价功率谱性能好坏的标准有很多，本文只给出两个影响最大的标准：分辨率和方差。分辨率即功率谱上能够区分的最小相邻频率成分，分辨率越高，我们观察信号的频率成分越清晰；方差大小则反映到功率谱波动性的大小，如果方差太大，功率谱波动性大，则很容易造成有用的频率成分被噪声淹没。所以，我们希望得到的这个"非精确"的功率谱，分辨率越高越好，方差越小越好。
同时，我们给出概率论与数理统计中所学的一致估计和非一致估计的概念，假定真实信号的功率谱为
，估计得到的"非精确"功率谱为
，如果满足公式（2-1），则称为一致估计。
（2-1）
2.1周期图法
2.1.1周期图法原理
我们已知N个离散的数据点
，对这些数据点进行傅里叶变换，得到式（2-2）：
(2-2)
再对（2-2）式取模的平方，除以N，即得到一个"非精确"的谱，如式（2-3）所示，这就是周期图法的原理。
（2-3）
2.1.2周期图法性能（Matlab仿真）
上一小节我们已经给出了周期图法的原理。本节将通过Matlab仿真给出数据点数N对功率谱性能好坏的影响，正如上文所述，将通过对所得功率谱的分辨率和方差两方面进行分析。
我们在Matlab中通过三个正弦函数和白噪声叠加，构造了一个随机信号。其数学形式如式（2-4），其中频率
、
、
，幅值
、
、
，相位
为相互独立在
上服从均匀分布随机相位，
为均值为0，方差为1的实值高斯白噪声，采样频率为1000。信号的时域形式如图2-1所示。
（2-4）
图2-1 实验所用的随机信号
当数据点数N分别为128、256、512和1024时，得到的功率谱分别如图2-2、图2-3、图2-4和图2-5所示。分辨率能够直观的通过功率谱图形看出，方差的数值由表2-1给出。
图2-2 N=128时周期图法得到的功率谱
图2-3 N=256时周期图法得到的功率谱
图2-4 N=512时周期图法得到的功率谱
图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个观测样本数据点
分为L段，每段数据长度为M, 分别对每段数据求周期图功率谱估计，然后求平均值，这种方法称平均周期图法。
那么这种方法会如何改善方差呢？我们给出证明：
（2-5）
其中：
（2-6）
由式（2-5）我们可以看出，平均周期图法将原来的方差变为原来的
，L为分段数。
2.1.4平均周期图法性能（Matlab仿真）
当数据点数N为1024，分段数分别为8、4、2时，平均周期图法得到的功率谱分别如图2-6、图2-7、图2-8所示。分辨率能够直观的通过功率谱图形看出，方差的数值由表2-2给出。
图2-6 L=8时平均周期图法得到的功率谱
图2-7 L=4时平均周期图法得到的功率谱
图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]
2.1.6修正的平均周期图法性能（Matlab仿真）
利用修正平均周期图法，分别使用矩形窗、Blackman窗和Hamming窗得到的功率谱如图2-9所示。
图2-9 不同窗函数的修正平均周期图法得到的功率谱
可以发现，矩形窗的分辨率最高，但是方差也最大，这是由于矩形窗频谱主瓣最窄，分辨率因此最高，旁瓣也高，导致频谱泄漏最严重，方差最大。
2.1.7总结
周期图法获得的功率谱随着样本点数越多，分辨率越大、方差越大；平均周期图法以牺牲分辨率来进一步改善方差；修正的平均周期图法允许段的重叠来进一步增大分段数、或者分段数相同，每段样本点数变多。无论是哪种方法都没有彻底结局方差与分辨率之间的矛盾。
2.2相关功率谱估计法-BT法
正如我们之前介绍的，要提高功率谱估计的分辨率，必须增加数据序列的长度N，但是较长的数据序列，由噪声引起的随机性得到更为充分的体现－较大的方差。事实上，当N无穷大时，方差为一非零常数。即周期图法无法实现功率谱的一致估计。而这节讲述的相关功率谱估计法（下文称作BT法），是一致估计。
2.2.1 BT法的原理
维纳辛钦定理指出，随机信号的相关函数与它的功率谱是一对傅里叶变换对。BT法就是基于这个原理。先由观测数据
估计出自相关函数，然后求自相关函数的傅立叶变换，以此变换作为对功率谱的估计，也称为间接法。BT法要求信号长度N以外的信号为零，这也造成BT法的局限性。
的自相关函数
定义如式（2-7）所示，得到的功率谱记为
，则BT法可以表述为式（2-8）。
（2-7）
（2-8）
2.2.2 BT法的性能（Matlab仿真）
数据点数N分别为128、256、512和1024的BT法，得到的功率谱如图2-10、图2-11、图2-12和图2-13所示。
图2-10 N=128时，BT法得到的功率谱
图2-11 N=256时，BT法得到的功率谱
图2-12 N=512时，BT法得到的功率谱
图2-13 N=1024时，BT法得到的功率谱
由上面实验可以发现，M随着N的增大而增大时，分辨率提高，方差变大。BT法仍然没有解决分辨率与方差之间的矛盾，但是BT法得到的功率谱当N为无穷大时，方差会趋向于零，即为一致估计[2]。
2.2.3 周期图法与BT法的关系
相关函数
可以写成如式（2-9）的卷积形式：
(2-9）
设序列
的傅立叶变换为
，则当M=N-1时，功率谱的估计可表示为式（2-10）的形式。可以看出周期图法可以看作BT法在取M=N-1时的特例。
（2-10）
结 论
本文通过Matlab仿真，以一个具体的随机信号为例，简单介绍了周期图法、平均周期图法、修正的平均周期图法以及BT法的基本原理，并对这些方法的性能进行分析。可以看出，无论是周期图法及其改进算法还是BT法都没有从根本上解决分辨率与方差的矛盾。经典功率谱估计是利用傅里叶变换估计功率谱，而我们之前分析随机信号不满足傅里叶变换的条件，所以经典功率谱估计方法不得不从无限长数据点截取有限长数据点，加入限制条件（周期图法实际上假定N点外数据周期重复、BT法假定N点外数据为零）来"强制"作傅里叶变换，这也是造成它局限性的原因。
参考资料
[1] 朱哲,钟宏伟. 非平稳随机信号分析处理方法研究[J] 安徽电子信息技术学院学报 2008.6:28-28
[2] 皇甫堪. 现代数字信号处理[M]. 电子工业出版社
部分matlab程序代码：
周期图法：（by宋同学）
1 Fs=1000; 2 f1=50; 3 f2=125; 4 f3=135; 5 N=128; 6 Nfft=N; 7 n=0:N-1; 8 t=n/Fs;%采用的时间序列 9 xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n)); 10 figure; 11 plot(n,xn);grid on;title('时域信号'); 12 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB 13 f=(0:length(P)-1)*Fs/length(P);%给出频率序列 14 figure; 15 plot(f(1:N/2),P(1:N/2));grid on;title('功率谱（dB图）');ylabel('功率谱/dB'); 16 xlabel('频率/Hz'); 17 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier振幅谱平方的平均值，并转化为dB 18 f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 19 figure; 20 plot(f(1:N/2),Pxx(1:N/2));%绘制功率谱曲线 21 xlabel('频率/Hz');ylabel('功率谱'); 22 title('周期图 N=128');grid on; 23 std(Pxx)^2 24 25 N=256; 26 Nfft=N; 27 n=0:N-1; 28 t=n/Fs;%采用的时间序列 29 xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n)); 30 figure; 31 plot(n,xn);grid on;title('时域信号'); 32 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB 33 f=(0:length(P)-1)*Fs/length(P);%给出频率序列 34 figure; 35 plot(f(1:N/2),P(1:N/2));grid on;title('功率谱（dB图）');ylabel('功率谱/dB'); 36 xlabel('频率/Hz'); 37 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier振幅谱平方的平均值，并转化为dB 38 f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 39 figure; 40 plot(f(1:N/2),Pxx(1:N/2));%绘制功率谱曲线 41 xlabel('频率/Hz');ylabel('功率谱'); 42 title('周期图 N=256');grid on; 43 std(Pxx)^2 44 45 N=512; 46 Nfft=N; 47 n=0:N-1; 48 t=n/Fs;%采用的时间序列 49 xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n)); 50 figure; 51 plot(n,xn);grid on;title('时域信号'); 52 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB 53 f=(0:length(P)-1)*Fs/length(P);%给出频率序列 54 figure; 55 plot(f(1:N/2),P(1:N/2));grid on;title('功率谱（dB图）');ylabel('功率谱/dB'); 56 xlabel('频率/Hz'); 57 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier振幅谱平方的平均值，并转化为dB 58 f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 59 figure; 60 plot(f(1:N/2),Pxx(1:N/2));%绘制功率谱曲线 61 xlabel('频率/Hz');ylabel('功率谱'); 62 title('周期图 N=512');grid on; 63 std(Pxx)^2 64 65 N=1024; 66 Nfft=N; 67 n=0:N-1; 68 t=n/Fs;%采用的时间序列 69 xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+1.5*randn(size(n)); 70 figure; 71 plot(n(1:1000),xn(1:1000));grid on;title('时域信号'); 72 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB 73 f=(0:length(P)-1)*Fs/length(P);%给出频率序列 74 figure; 75 plot(f(1:N/2),P(1:N/2));grid on;title('功率谱（dB图）');ylabel('功率谱/dB'); 76 xlabel('频率/Hz'); 77 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier振幅谱平方的平均值，并转化为dB 78 f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 79 figure; 80 plot(f(1:N/2),Pxx(1:N/2));%绘制功率谱曲线 81 xlabel('频率/Hz');ylabel('功率谱'); 82 title('周期图 N=1024');grid on; 83 std(Pxx)^2
平均周期图法
1 clear; 2 Fs=1000; 3 f1=50; 4 f2=125; 5 f3=135; 6 n=0:1/Fs:1; 7 xn=cos(2*pi*f1*n)+1.5*cos(2*pi*f2*n)+cos(2*pi*f3*n)+1.5*randn(size(n)); 8 9 N=1024;Nsec=512;%数据的长度和FFT所用的数据长度 10 Pxx1=abs(fft(xn(1:512),Nsec).^2)/Nsec; %第一段功率谱 11 Pxx2=abs(fft(xn(513:1000),Nsec).^2)/Nsec;%第二段功率谱 12 Pxx=10*log10((Pxx1+Pxx2)/2);%Fourier振幅谱平方的平均值，并转化为dB 13 (std((Pxx1+Pxx2)/2))^2 14 f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 15 figure; 16 plot(f(1:Nsec/2),Pxx(1:Nsec/2));%绘制功率谱曲线 17 xlabel('频率/Hz');ylabel('功率谱/dB'); 18 title('N=2*512');grid on; 19 20 N=1024;Nsec=256;%数据的长度和FFT所用的数据长度 21 Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec; %第一段功率谱 22 Pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;%第二段功率谱 23 Pxx3=abs(fft(xn(513:768),Nsec).^2)/Nsec;%第三段功率谱 24 Pxx4=abs(fft(xn(769:1000),Nsec).^2)/Nsec;%第四段功率谱 25 Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);%Fourier振幅谱平方的平均值，并转化为dB 26 std((Pxx1+Pxx2+Pxx3+Pxx4)/4)^2 27 f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 28 figure; 29 plot(f(1:Nsec/2),Pxx(1:Nsec/2));%绘制功率谱曲线 30 xlabel('频率/Hz');ylabel('功率谱/dB'); 31 title('N=4*256');grid on; 32 33 N=1024;Nsec=128;%数据的长度和FFT所用的数据长度 34 Pxx1=abs(fft(xn(1:128),Nsec).^2)/Nsec; %第一段功率谱 35 Pxx2=abs(fft(xn(129:256),Nsec).^2)/Nsec;%第二段功率谱 36 Pxx3=abs(fft(xn(257:384),Nsec).^2)/Nsec;%第三段功率谱 37 Pxx4=abs(fft(xn(385:512),Nsec).^2)/Nsec;%第四段功率谱 38 Pxx5=abs(fft(xn(513:640),Nsec).^2)/Nsec; %第五段功率谱 39 Pxx6=abs(fft(xn(641:768),Nsec).^2)/Nsec;%第六段功率谱 40 Pxx7=abs(fft(xn(769:896),Nsec).^2)/Nsec;%第七段功率谱 41 Pxx8=abs(fft(xn(897:1000),Nsec).^2)/Nsec;%第八段功率谱 42 Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8)/8);%Fourier振幅谱平方的平均值，并转化为dB 43 std((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8)/8)^2 44 45 f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 46 figure; 47 plot(f(1:Nsec/2),Pxx(1:Nsec/2));%绘制功率谱曲线 48 xlabel('频率/Hz');ylabel('功率谱/dB'); 49 title('N=8*128');grid on;
修正的平均周期图法（by周同学）
1 clear;%求1024点功率谱以及方差 2 Fs=1000; 3 f1=50; 4 f2=125; 5 f3=135; 6 n=0:1/Fs:1; 7 xn=cos(2*pi*f1*n)+1.5*cos(2*pi*f2*n)+cos(2*pi*f3*n)+1.5*randn(size(n)); 8 M=512;N=1024;%数据总点数1024，每段长度M 9 window=hamming(M); 10 Pxxtotal=0; 11 L=N/(M/2)-1; 12 for i=1:1:(L-1) 13 m=abs(fft(window'.*(xn((M/2+M/2*i-M+1):(M/2+M/2*i))),M).^2)/M; 14 Pxxtotal=Pxxtotal+m; 15 end 16 window=hamming(Fs-(N-M+1)+1); 17 mend=abs(fft(window'.*(xn((N-M+1):Fs)),M).^2)/M; 18 Pxxtotal=(Pxxtotal+mend)/L; 19 Pxx=10*log10((Pxxtotal));%Fourier振幅谱平方的平均值，并转化为dB 20 f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 21 w1=fft(window,N); 22 w10=abs(fftshift(w1)); 23 24 plot(f(1:M/2),Pxx(1:M/2));%绘制功率谱曲线 25 xlabel('频率/Hz');ylabel('功率谱/dB'); 26 grid on; 27 B=var(Pxxtotal)
BT法（by阳同学）
1 Fs=1000; 2 n=0:1/Fs:1; 3 x=cos(2*pi*50*n)+1.5*cos(2*pi*125*n)+cos(2*pi*135*n)+1.5*randn(size(n)); 4 nfft=1024; 5 ncxk=3*nfft/4; 6 cxn=xcorr(x,'unbiased'); 7 CXk=fft(cxn,ncxk); 8 Pxx=abs(CXk); 9 index=0:round(ncxk/2-1); 10 k=index*Fs/ncxk; 11 C=Pxx(index+1); 12 P=(Pxx(index+1)); 13 plot(k,P);grid on 14 var(C) 15 title('BT法功率谱估计 N=1024'); 16 xlabel('频率 Hz'); 17 ylabel('功率');
