功率谱估计：提取被淹没在噪声中的有用信号

weixin

　　功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。下面对谱估计的发展过程做简要回顾：

　　英国科学家牛顿最早给出了“谱”的概念。后来，1822年，法国工程师傅立叶提出了著名的傅立叶谐波分析理论。该理论至今依然是进行信号分析和信号处理的理论基础。

　　傅立叶级数提出后，首先在人们观测自然界中的周期现象时得到应用。19世纪末，Schuster提出用傅立叶级数的幅度平方作为函数中功率的度量，并将其命名为“周期图”(periodogram)。这是经典谱估计的最早提法，这种提法至今仍然被沿用，只不过现在是用快速傅立叶变换(FFT)来计算离散傅立叶变换(DFT)，用DFT的幅度平方作为信号中功率的度量。


　　周期图较差的方差性能促使人们研究另外的分析方法。1927年，Yule提出用线性回归方程来模拟一个时间序列。Yule的工作实际上成了现代谱估计中最重要的方法——参数模型法谱估计的基础。

　　Walker利用Yule的分析方法研究了衰减正弦时间序列，得出Yule-Walker方程，可以说，Yule和Walker都是开拓自回归模型的先锋。

　　1930年，著名控制理论专家Wiener在他的著作中首次精确定义了一个随机过程的自相关函数及功率谱密度，并把谱分析建立在随机过程统计特征的基础上，即，“功率谱密度是随机过程二阶统计量自相关函数的傅立叶变换”，这就是Wiener-Khintchine定理。该定理把功率谱密度定义为频率的连续函数，而不再像以前定义为离散的谐波频率的函数。

　　1949年，Tukey根据Wiener-Khintchine定理提出了对有限长数据进行谱估计的自相关法，即利用有限长数据估计自相关函数，再对该自相关函数球傅立叶变换，从而得到谱的估计。1958年， Blackman和Tukey在出版的有关经典谱估计的专著中讨论了自相关谱估计法，所以自相关法又叫BT法。 周期图法和自相关法都可用快速傅立叶变换算法来实现，且物理概念明确，因而仍是目前较常用的谱估计方法。

　　1948年，Bartlett首次提出了用自回归模型系数计算功率谱。自回归模型和线性预测都用到了1911年提出的Toeplitz矩阵结构，Levinson曾根据该矩阵的特点于1947年提出了解Yule-Walker的快速计算方法。这些工作为现代谱估计的发展打下了良好的理论基础。

　　1965年，Cooley和Tukey提出的FFT算法，也促进了谱估计的迅速发展。

　　现代谱估计主要是针对经典谱估计的分辨率差和方差性能不好的问题而提出的。现代谱估计从方法上大致可分为参数模型谱估计和非参数模型谱估计两种，前者有AR模型、MA模型、ARMA模型、PRONY指数模型等;后者有最小方差方法、多分量的MUSIC方法等。

　　周期运动在功率谱中对应尖锋，混沌的特征是谱中出现"噪声背景"和宽锋。它是研究系统从分岔走向混沌的重要方法。在很多实际问题中(尤其是对非线性电路的研究)常常只给出观测到的离散的时间序列X1, X2, X3,...Xn，那么如何从这些时间序列中提取前述的四种吸引子(零维不动点、一维极限环、二维环面、奇怪吸引子)的不同状态的信息呢?

　　我们可以运用数学上已经严格证明的结论，即拟合。我们将N个采样值加上周期条件Xn+i=Xi，则自关联函数(即离散卷积)为然后对Cj完成离散傅氏变换，计算傅氏系数。 Pk说明第k个频率分量对Xi的贡献，这就是功率谱的定义。当采用快速傅氏变换算法后，可直接由Xi作快速傅氏变换，得到系数 然后计算，由许多组{Xi}得一批{Pk'}，求平均后即趋近前面定义的功率谱Pk。

　　从功率谱上，四种吸引子是容易区分的，混沌的功率谱表现为"噪声背景"及宽锋。考虑到实际计算中，数据只能取有限个，谱也总以有限分辨度表示出来，从物理实验和数值计算的角度看，一个周期十分长的解和一个混沌解是难于区分的，这也正是功率谱研究的主要弊端。

　　信号的频谱分析是研究信号特性的重要手段之一，通常是求其功率谱来进行频谱分析。功率谱反映了随机信号各频率成份功率能量的分布情况，可以揭示信号中隐含的周期性及靠得很近的谱峰等有用信息，在许多领域都发挥了重要作用。然而，实际应用中的平稳随机信号通常是有限长的，只能根据有限长信号估计原信号的真实功率谱，这就是功率谱估计。

　　功率谱估计分为经典谱估计和现代谱估计。经典谱估计是将数据工作区外的未知数据假设为零，相当于数据加窗，主要方法有相关法和周期图法;现代谱估计是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱，主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的，应用最广的是AR参数模型。本文就应用上述两种方法，在MATLAB中对一个语音信号进行功率谱估计，然后进行比较。

　　经典功率谱估计

　　1、相关法
       是利用维纳-辛钦定理该方法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。

　　MATLAB程序如下：
　　[xn,Fs,bits]=wavread('1.wav');
　　n=0:1/Fs:1;
　　nfft=512;
　　cxn=xcorr(xn,'unbiased');
　　CXk=fft(cxn,nfft);
　　Pxx=abs(CXk);

　　index=0:round(nfft/2-1);
　　k=index*Fs/nfft;
　　a=log(10);
　　b=log(Pxx(index+1));
　　c=b/a;
　　plot_Pxx=10*c;

　　figure(1);
　　plot(k,plot_Pxx);
　　xlabel('frequeney(hz)');
　　ylabel('power spectraldensity(Db-Hz)');
　　title('recorrelation psd estimate');

　　2. 周期图法
　　周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得x(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。

　　MATLAB程序如下：
　　[y,Fs,bits]=wavread('1.wav');
　　N=512;
　　Nfft=1024;
　　n=0:N-1;
　　XF=fft(y,Nfft);
　　Pxx=abs(XF).^2/length(n);
　　index=0:round(Nfft/2-1);
　　f=index*Fs/Nfft;
　　plot(f,10*log(Pxx(index+1))),
　　grid

　　现代功率谱估计

　　在经典谱估计进行功率谱估计时，由于将所有在窗口外的数据都视为0，这便使得谱估计的质量下降。现代谱估计是利用待研究信号的先验知识，对信号在窗口外的数据作出某种较合理的假设，以达到提高谱估计质量的目的。现代功率谱估计主要是利用白噪声输入参数模型之后得到输出序列，当改变系统参数时得到的序列也不同，这样当改变参数使得输出与已知有限序列相同或者近似时就可以利用下面的公式求得其功率谱。

　　可以说现代谱估计实际上是对模型参数的估计。现代谱估计的参数模型有自回归滑动平均(ARMA)模型、自回归(AR)模型、滑动平均(MA)模型，Wold分解定理阐明了三者之间的关系：任何有限方差的ARMA或MA模型的平稳随机过程可以用无限阶的AR模型表示，任何有限方差的ARMA或MA模型的平稳随机过程可以用无限阶的MA模型表示。但是由于只有AR模型参数估计是一组线性方程，而实际的物理系统往往是全极点系统，因而AR应用最广。

　　AR模型：利用AR模型进行功率谱估计须通过L-D递推算法由Yule-Walker方程求得AR的参数：a1，a2，…ap，σp2，或者利用线性预测滤波器系统函数与AR模型系统函数的关系，通过Burg递推算法求得AR模型参数。具体推导就不再给出。在Matlab仿真中可调用pburg函数直接画出基于burg算法的功率谱估计的曲线图。

　　MATLAB程序如下：
　　N=512;
　　Nfft=1024;
　　[xn,fs,bits]=wavread('1.wav');
　　%xn=<13921*2> fs=22050;bits=16;
　　xn=xn(:,1);
　　xn1=xn(1:N);
　　order1=50;
　　figure
　　pburg(xn1,order1,Nfft,fs);

　　结语

　　(1) 功率谱估计中的相关函数法和周期图法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。

　　(2) 现代谱估计方法曲线明显比经典谱估计方法光滑，说明其处理结果的方差比经典谱估计方法处理的结果小，证明了现代谱估计方法利用参数模型对窗口外的数据的假设是比较合理有效的，能达到提高谱估计质量的目的，这也是现代谱估计优于经典谱估计的主要原因。

　　本文转载自新浪阿彬的博客，原文为Geoinformatics的博文