基于Matlab功率谱密度估计方法
摘要
在实际情况下, 许多平稳信号无法导出数学表达式, 要准确获取这些信号的功率谱密度存在一定的困难。根据维纳-辛钦 (Wiener Khintchine)定理,提出一种基于Matlab编程实现这类信号的功率谱密度的估计方法。通过仿真实验表明, 该方法简单易行,准确性较高。
关键词：平稳信号；功率谱密度；估计方法；
Estimation Method for Power Spectral Density of Stationary
Signals
Abstract
In practice m any mathematical expressions can no t be derived from stationary signals as there are: some difficulties in getting the power spectral density of the signals According to Wiener Khintchine theorem, the paper suggested an estimation method for power spectral density of the signals based on Matlab programming Simulation results show that the method is simple and comparatively accurate.iim. Keywords：Stationary signal; Power spectral density; estuation method;
前言
信号的功率谱密度在通信系统的设计、信号传输、信号分析及信号处理等方面有很重要的作用。频谱分析是信号在频率域上的重要手段，他反映了信号的频率成分以及分布情况。信号频谱估计是信号分析的重要手段，也是信号综合的前提
在实际情况下, 许多信号很难导出闭合的数学表达式, 例如密勒码 ( Miler code)、双极性不归零码 ( NRZI) 等信号至今在时域中还没有数学表达式, 虽然Hecht& Guida在1969年导出了密勒码的功率谱密度的表达式并做出了图形，但推导过程及得出的表达式非常复杂。因此在进行理论计算或推导时, 常常需要将实际使用的信号理想化地近似为标准的正弦信号或有数学表达式的信号, 才能推导出相应的功率谱密度表达式。
而这种近似由于与实际信号有差别, 因此不能代表实际信号的情况, 且有时根据这种近似推导得出的结果可能存在错误。本文首先介绍多种功率谱密度相关方法的MATLAB实现，并根据维纳 - 辛钦定理, 提出一种基于Matlab编程实现信号的功率谱密度的估计方法：在不知道信号的数学表达式的情况下, 通过建模及编程运算就可以实现该信号的功率谱密度的样值输出及画出图形。
1 功率谱密度估计方法的实现原理及流程
1.1功率谱密度相关方法的MATLAB实现
1.1 1.Welch法估计及其MATLAB函数
Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计(PSD如上例中的下半部分的求法)和两个信号序列的互功率谱估计(CSD)。
MATLAB信号处理工具箱函数提供了专门的函数PSD和CSD自动实现Welch法估计，而不需要自己编程。
(1) 函数psd利用Welch法估计一个信号自功率谱密度，函数调用格式为：
[Pxx[,f]]=psd(x[,Nfft,Fs,window,Noverlap,’dflag’])
式中，x为信号序列；Nfft为采用的FFT长度。这一值决定了功率谱估计速度，当Nfft采用2的幂时，程序采用快速算法；Fs为采样频率；Window定义窗函数和x分段序列的长度。窗函数长度必须小于或等于Nfft，否则会给出错误信息；Noverlap为分段序列重叠的采样点数(长度)，它应小于Nfft;dflag为去除信号趋势分量的选择项：’linear’,去除线性趋势分量，’mean’去除均值分量，’none’不做去除趋势处理。Pxx为信号x的自功率谱密度估计。f为返回的频率向量，它和Pxx对应，并且有相同长度。
在psd函数调用格式中，缺省值为：Nfft=min(256,length(x)),Fs=2Hz, window=hanning(Nfft),noverlap=0. 若x是实序列，函数psd仅计算频率为正的功率
注意程序前半部分中频率向量f的创建方法。它与函数psd的输出Pxx长度的关系如下：若x为实序列，当Nfft为奇数时，f=(0:(Nfft+1)/2-1)/Nfft;当Nfft为偶数时，f=(0:Nfft/2)/Nfft。
函数还有一种缺省返回值的调用格式，用于直接绘制信号序列x的功率谱估计曲线。
函数还可以计算带有置信区间的功率谱估计，调用格式为：
[Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,Noverlap,p)
式中，p为置信区间，0<=p<=1。
由此可知，滤波器输入白噪声序列的输出信号的功率谱或自相关可以确定滤波器的频率特性。
(2)函数csd利用welch法估计两个信号的互功率谱密度，函数调用格式为：
[Pxy[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,’dflag’)
[Pxy,Pxyc[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,p)
这里，x,y为两个信号序列；Pxy为x,y的互功率谱估计；其他参数的意义同自功率谱函数psd。
可以看到，两个白噪声信号的互功率谱(上图)杂乱无章，看不出周期成分，大部分功率谱在-5dB以下。然而白噪声与带有噪声的周期信号的功率谱在其周期(频率为1000Hz)处有一峰值，清楚地表明了周期信号的周期或频率。因此，利用未知信号与白噪声信号的互功率谱也可以检测未知信号中所含有的频率成分。
1.1.2 多窗口法
多窗口法(Multitaper method,简称MTM法)利用多个正交窗口(Tapers)获得各自独立的近似功率谱估计，然后综合这些估计得到一个序列的功率谱估计。相对于普通的周期图法，这种功率谱估计具有更大的自由度，并在估计精度和估计波动方面均有较好的效果。普通的功率谱估计只利用单一窗口，因此在序列始端和末端均会丢失相关信息，而且无法找回。而MTM法估计增加窗口使得丢失的信息尽量减少。
MTM法简单地采用一个参数：时间带宽积(Time-bandwidth product)NW，这个参数用以定义计算功率谱所用窗的数目，为2*NW-1。NW越大，功率谱计算次数越多，时间域分辨率越高，而频率域分辨率降低，使得功率谱估计的波动减小。随着NW增大，每次估计中谱泄漏增多，总功率谱估计的偏差增大。对于每一个数据组，通常有一个最优的NW使得在估计偏差和估计波动两方面求得折中，这需要在程序中反复调试来获得。
MATLAB信号处理工具箱中函数PMTM就是采用MTM法估计功率谱密度。函数调用格式为：
[Pxx[,f]]=pmtm(x[,nw,Nfft,Fs])
式中，x为信号序列；nw为时间带宽积，缺省值为4。通常可取2，5/2，3，7/2；Nfft为FFT长度；Fs为采样频率。
上面的函数还可以通过无返回值而绘出置信区间，如pmtm(x, nw, Nfft, Fs,’option’,p)绘制带置信区间的功率谱密度估计曲线，0<=p<=1。 1.1.3最大熵法(Maxmum entropy method，MEM法)
如上所述，周期图法功率谱估计需要对信号序列“截断”或加窗处理，其结果是使估计的功率谱密度为信号序列真实谱和窗谱的卷积，导致误差的产生。 最大熵功率谱估计的目的是最大限度地保留截断后丢失的“窗口”以外信号的信息，使估计谱的熵最大。主要方法是以已知的自相关序列rxx(0),rxx(1),?,rxx(p)为基础，外推自相关序列rxx(p+1), rxx(p+2),?，保证信息熵最大。
最大熵功率谱估计法假定随机过程是平稳高斯过程，可以证明，随机信号的最大熵谱与AR自回归(全极点滤波器)模型谱是等价的。
MATLAB信号处理工具箱提供最大熵功率谱估计函数pmem，其调用格式为：
[Pxx,f,a]=pmem(x,p,Nfft,Fs,’xcorr’)
式中，x为输入信号序列或输入相关矩阵；p为全极点滤波器阶次；a为全极点滤波器模型系数向量；’xcorr’是把x认为是相关矩阵。
比较最大熵功率谱估计(MEM)和改进的平均周期图功率谱估计，可见，MEM法估计的功率谱曲线较光滑。在这一方法中，MEM法选定全极点滤波器的阶数取得越大，能够获得的窗口外的信息越多，但计算量也越大，需要根据情况折中考虑。
1.1.4 多信号分类法
MATLAB信号处理工具箱还提供另一种功率谱估计函数pmusic。该函数执行多信号分类法(multiple signal classification，Music法)。将数据自相关矩阵看成由信号自相关矩阵和噪声自相关矩阵两部分组成，即数据自相关矩阵R包含有两个子空间信息：信号子空间和噪声子空间。这样，矩阵特征值向量(Eigen vector)也可分为两个子空间：信号子空间和噪声子空间。为了求得功率谱估计，函数pmusic计算信号子空间和噪声子空间的特征值向量函数，使得在周期信号频率处函数值最大，功率谱估计出现峰值，而在其他频率处函数值最小。其调用格式为：
[Pxx,f,a]=pmusic(x,p[[,thresh],Nfft,Fs,window,Noverlap])
式中,x为输入信号的向量或矩阵；p为信号子空间维数；thresh为阈值，其他参数的意义与函数psd相同。
1.2平稳信号实现原理
根据维纳 - 辛钦定理, 平稳信号的自相函数和功率谱密度互为傅里叶变换对, 表达式如下：
P()=r()=
d ,d，
<<< ( 1) <<< ( 2)
式中, P () 为信号功率谱密度函数, r () 为信号的自相关函数。根据式 ( 1) 可知, 只要想办法得到信号的自相关函数, 经过傅里叶变换后即可得到信号的功率谱密度函数。
