基于Matlab功率谱密度估计方法

基于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

由此可知，滤波器输入白噪声序列的输出信号的功率谱或自相关可以确定滤波器的频率特性。

（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

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 ，

式中, P () 为信号功率谱密度函数, r () 为信号的自相关函数。根据式 ( 1) 可知, 只要想办法得到信号的自相关函数, 经过傅里叶变换后即可得到信号的功率谱密度函数。

1.3 平稳信号实现流程图 在Matlab 仿真软件的环境中, 信号可以通过接口输入仿真界面或建模实现信号。产生信号后, 对信号进行抽样, 然后进行自相关运算[ 6], 再将自相关运算结果进行傅里叶变换即可得到信号的功率谱密度。流程如图1所示。

图1 实现功率谱密度估计得流程图

2 建模、编程实例及结果分析

2 1 建模、编程实例

根据功率谱估计原理及实现流程图, 可以通过一些技巧性的处理在Matlab 中实现对无表式平稳信号的功率谱的估计。以下通过实例展示如何通过建模及编程实现之, 然后通过仿真结果与理论计算结果的对比, 验证这种方法的正确性。图 2是产生 密勒码的模型图(建模时，命名为Milermdl) 。根据奈 奎斯特( Nyquist) 抽样定理, 抽样频率必须是信号频的2倍以上才可以在抽样后的信号中恢复原信号。在本例的模型中, 为了进一步减小仿真的误差, 采样频率设置为信号频率的5倍。另外, 为了使输出的样值在进行自相关运算时的结果更准确, 抽样样值的点数应尽量多取一些。本例中输出至工作空间的抽样样值, 即进行自相关运算的样

值设置为3000点。

图2 产生密勒码及输出抽样样值的模型图

图2中的模型实现了图1中的信号产生及抽样, 图1中信号的自相关运算、傅里叶变换及画图通过编程实现。程序见下面的 Miler.m 。

Miler .m

% 自相关运算

clear

M lag= 50;

sim( Miler. m dl)’; % 执行模型的仿真

x corrmile r= xcorr( Miler, Mlag, biased )';

Ts= 1/5, % 设置抽样时间间隔

N =505, % 设置 FFT 的点数;

% 对密勒码的自相关输出进行 FFT 变换

fftm iler= fft ( xcorr miler, N )

fT = 1 / ( N* T s) .* [ 0;N-1] ; 设置归一化频率样点

M fftmiler= 2* T s. * fftm iler, 将 FFT 值换算为 FT 值

plot ( fT , abs ( M fftm ile r) , .' ,' MarkerSize, '14) ;

hold on;

% 计算并打印理论计算结果

Tf = 0. 001 :0 .01: 2 .0;

w = pi* fT;

pf= 2* ( 23- 2* cos( w ) - 22* cos( 2* w ) - 12* cos( 3* w ) + 5* cos( 4* w ) + 12*cos( 5* w ) + 2*cos( 6* w ) - 8* cos( 7* w ) + 2 *cos( 8 *w ) ). / ( 2 * ( w* w ) * ( 17+ 8 *cos( 8* w ) ) );

plot( fT, p f) ;

axis( [ 0 2 – 0. 1 3] ); x label( f T)'; y label( Spectral density )'; legend( 密勒码仿真结果 , '密勒码理论结果 )';

运行程序后得到密勒码的功率谱密度如图3所示。另根据 Hecht& Guida所导出的密勒码的功率谱密度的表达式(3)，编程画出理论的密勒码的功率谱密度如图4所示，仿真结果与理论结果对比如图5。图3~ 5都是单边谱, 因此在程序中也进行了相应的换算处理。图 3~ 5的纵坐标都是功率谱密值，为了方便作图对比，这3个图中的横坐标都归一化处理为fT (其中f 是频率，T 是密勒码的一个码元的持续时间) 。

(f ) =1/2( 17+ 8cos8 )*( 23- 2cos - 22cos 2- 12cos 3+5cos 4 + 12cos 5+2cos 6+ - 8cos 7+ 2cos 8) (= fT ) (3)

2 2 结果分析

从图5中可以看出, 由建模及编程实现的密勒码的功率谱密度的仿真结果在主峰及副峰值处均能非常好地与理论值拟合，虽然在功率谱变化平坦的地方，由于信号的数字化及FFT 变换的误差有一些波动，但仿真结果仍然具有相当高的准确性，符合实际应用的要求。

图3 仿真的密勒码功率谱

图4 理论的密勒码功率谱

图5 仿真结果与理论结果对比

3结语

经典功率谱估计方法实现简单, 不需要任何先验的模型假设, 但存在估计的可靠性和分辨率之间的矛盾。本文在维纳 - 辛钦 (W iener Khintch in e) 定理的基础上, 实现了一种基于Matlab 的无表达式信号的功率谱密度估计方法, 通过建模及编程对比密勒码的功率谱的仿真值及理论值, 显示出这种方法能相当准确地获得 信号的功率谱密 度。由于只需 要输入信号到Matlab 的接口, 或者建模产生相关的信号, 即可方便地得到信号的功率谱密度样值, 因此该方法具有一定的实用价值。

参考文献:

[1] PROAK IS J G. Digital Communications[M ]. 4 thed Beijing Publishing House of Electronics Industry 2008.

[2] 邵玉斌. Matlab /Simulink通信系统建模与 仿真实例分析 [M ]. 北京: 清华大学出 版社, 2008m.

[3] 赵树杰, 赵建勋. 信号检测与估计理论 [ M ]. 北京: 清华大学出版社, 2005.

[4] 胡广书. 数字信号处理理论、算法与实现 [ M ]. 2版. 北京: 清华清华大学出版社, 2003.

[5] HECHT M, GUIDA A. Delay Modulation[ C ]. P roc IEEE, 1969( 57): 1314- 1316..

[6] TITSW ORTH i R C, WELCH L R. Power Spectra l of Signals Modulated by Rando and Pseudo rando Sequences[ C ]. JPLmm Tech. Rep. , 1961( 32) : 140.

[7] 王晓燕, 朱志峰, 方世良. 非合作水 声低载噪比直接序列扩频信号的载频估计方法 [ J]. 电子与信息学报, 2009, 31( 10): 2438- 2442.

[8]张朋友, 吕明；近代信号处理理论与犯法. 北京：国防工业出版社, 2005

[9]张旭东, 陆明泉. 离散信号处理；北京：清华大学出版社, 2005

[10]张贤；达现代信号处理(第二版) ； 北京：清华大学出版社,2002；

[11]汪源源；现代信号处理理论和方法；上海：复旦大学出版社, 2002；