2011年 8月 15日第 34卷第 16期现代电子技术M odern Electro nics T echniqueA ug. 2011V ol. 34N o. 16基于 Matlab 实现现代功率谱估计王春兴(山东师范大学物理与电子科学学院 , 山东济南 250014摘要 :功率谱估计可以分为经典谱估计和现代谱估计。
现代谱的估计可建立 A R 模型对离散信号进行谱估计、建立 M A 模型和 A RM A 模型进行谱估计。
基于 M atlab 对三种模型进行仿真 , 并对结果进行了分析。
结果显示 , 三种模型对现代谱的获得是有效的 , 并得到较好的谱估计。
关键词 :P SE; 现代功率谱估计 ; AR 模型法 ; A RM A中图分类号 :T N911-34; G202 文献标识码 :A 文章编号 :1004-373X (2011 16-0065-03Modern Power Spectrum Estimation Based on MatlabW AN G Chun -x ing(Colleg e o f Physics and Elect ro nics, Shando ng No rm al U niversity , Jinan 250014, Chi naAbstract :Po wer spectr um estimation can be divided into classical spectr al estimat ion and modern spectr al estimation. M odern spectr al estimation model can establish AR mo del, M A mo del and ARM A model fo r discr ete sig nals to per for m spec -t ralestimatio n. T hese t hr ee models can be simulated based o n M atlab, and the r esults ar e analy zed. T he r esult s sho w that the three models of mo der n spect rum are valid, and can get better spectrum estimatio n.Keywords :PSE; mo der n pow er spect rum est imatio n; A R model method; A RM A收稿日期 :2011-03-26基金项目 :国家自然科学基金项目资助 (10874103随机信号在时域上是无限长的 , 在测量样本上也是无穷多的 , 因此随机信号的能量是无限的 , 应用功率信号来描述。
然而 , 功率信号不满足傅里叶变换的狄里赫莉绝对可积的条件 , 因此严格意义上随机信号的傅里叶变换是不存在的。
因此 , 要实现随机信号的频域分析 , 不能简单从频谱的概念出发进行研究 , 而是功率谱 [1-2]。
信号的功率谱密度描述随机信号的功率在频域随频率的分布。
利用给定的 N 个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。
谱估计方法分为两大类 :经典谱估计和现代谱估计。
作为经典谱估计 , 其主要缺陷 [3]是描述功率谱波动的数字特征方差性能较差 , 频率分辨率低 ; 方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算。
分辨率低的原因是在周期图法中 , 假定工作区域以外的数据全为零 , 而对相关函数法中 , 假定延迟窗以外的自相关函数全为零。
这是不符合实际情况的 , 因而产生了较差的频率分辨率。
而现代谱估计的目标都是旨在改善谱估计的分辨率。
1 现代功率谱估计的参数模型如果将平稳的随机过程 x (n 看成是一个输入序列u(n (白噪声过程 , 其方差为 2激励线性系统 H (z 的输出 , 那么由已知的 x (n 或其自相关函数 r(m 就可估计 H (z 的参数 , 再由 H (z 的参数估计 x (n 的功率谱。
即当输入一线性系统 H (z 时 , 其输出 x (n 的功率谱可用系统 H (z 的参数来估计。
按照该思路 , x (n 的功率谱可表示为 :P X (e j= 21+qk=1b k e-j k2/1+pk=1a k e-j k2(1式中 : 2是激励白噪声的方差 ; P X (e j为功率谱密度 , a k 和 b k 为模型参数。
式 (1 也称为 ARM A 模型。
如果式 (1 参数 b 1, b 2, , b q 全为 0, 就变为 AR 模型 , 它是一个全极点模型 [4], 其含义是该模型现在的输出是现在的输入和过去 p 个输出的加权和。
该情况下的功率谱密度为 :P X (e j=21+k=1ake-j k2(2如果 ARMA 模型参数 a 1, a 2, a p 全为 0, 则变为 M A 模型 , 它是一个全零点模型。
P X (e j= 21+qk=1b k e -j k2(3AR 模型的正则方程是一组线性方程 , 而 M A 和ARM A 模型是非线性方程。
2 基于 AR 模型的功率谱AR 模型又称为自回归模型 , 它是一个全极点模型 , 它当前输出是现在输入和过去输入的加权和 :x (n =-pk=1a k x (n -k +u(n(4式中 :u(n 为白噪声信号 ; p 为 AR 模型的阶数。
对式 (4 求 z 变换 , 便得到系统函数 :H (z =A (z =1/(1+ pk=1a k z -k(5 由随机信号通过线性系统理论知输出序列的功率谱 :P X (e j= 2/1+pk=1a k e-j k2(6式中 : 2为白噪声序列的方差 , 因此进行功率谱估计 , 需要知道 AR 模型的参数 a k (k =1, 2, , p 及 2, 这些由 AR 模型的正则方程 [5]决定。
AR 模型功率谱估计是一个全极点的模型 , 要利用 AR 模型进行功率谱估计须通过 levinson _dubin递推算法由正则方程求得 AR 的参数 : 2, a 1, a 2, , a p 。
在Matlab 仿真中可调用 Pburg 函数直接画出基于 bur g 算法的功率谱估计的曲线图如图 1所示。
用周期图法求出的功率谱曲线和 burg 算法求出的 AR 功率谱曲线(p =50 , 其程序如下 :fs=200;n=0:1/fs:1;xn=cos(2*pi*40*n +cos(2*pi*41*n +34*co s(2*pi+90*n +0. 1*randn(size(n ; window =bo xcar (1eng th(x n ; nfft=512;[pxx , f]=perio do gr am(x n, w indow , nfft, fs ; subplot(121 ;plot(f, 10+log 10(pxx ; xlabel( ffequency(hz ;ylabel( pow er spectr al density(Db/H z ; title( perio do gr am PSE estimate ; or der l=50; range= half ; magunits= db ; subplot(122 ;pburg (x n, or der l, nfft, fs, r ang e;图 1 基于 AR 模型的功率谱3 基于 MA 模型的功率谱基于 M A 模型算法原理是由 N 点数据 x (n 建立一个 p 阶的 AR 模型 , 从而求出 p 阶 AR 系数 a, 再利用AR 系数建立线性预测 , 等效 q 阶 AR 模型 (p q , 再利用求解 AR 模型的出 b, 从而实现 M A 功率谱估计 [6], 如图 2所示。
根据上述原理建立函数 function[b , epb]=my ma(x , ip, iq , 其中 ip 和 iq 分别是 2次 AR 模型的阶次 , b 是得到最终 MA 模型的参数 , 程序如下 :A 1A 14=[1, a4 ]B14=fliplr (conv (fliplr (B1 , fliplr (A 14 ; y 24=filter (B14, A 1, randn(1,N ; %.*[zer os(1, 200 , ones(1, 256 ];[A ma4, Ema4]=arburg (y24, 32 , B1b4=arburg (A ma4, 4%---求功率谱 ---%w=linspace(0, pi, 512 ; %H1=fr eqz(B1, A 1, w H14=f reqz(b4, A 14, w ; %Ps1=abs(H 1. ^2;Py 14=abs(H14 . ^2;%if P y14>200%PP y14=200; %elseif P y14<200%PP y14=P y14; %endSPy 14=SPy14+P y14;V SPy14=V SP y14+abs(Py 14 . ^2;fig ure(4plot(w. /(2*pi , P s1, w. /(2*pi , Py 14 ;leg end( 真实功率谱 , 20次 ARM A(4, 4 的估计图 ; hold on;R8=zer os(16, 8 ; %A RM A (8, 8 的 R r8=zer os(16, 1 ; %A RM A (8, 8 的 r for i=1:16r8(i, 1 =-Ry (264+i ; for j=1:8R8(i, j =Ry(264+i-j ; End End R8r8a8=inv (R8 *R8 *R8 *r8%利用最小二乘法得到的 a的估计参数图 2 基于 M A 模型的功率谱及均值4 基于 ARMA 的功率谱该算法原理是事先估计 AR 参数 , 同时对已知数据 x (n 用 FIR 滤波器 , 将其输出近似 MA (q 过程 , 并求66现代电子技术2011年第 34卷解 MA(q 参数的方法求出 b , 从而实现 ARM A 模型的参数估计。
根据上述原理建立函数 function[b , epb]=myarm a(x , ip, iq, M 。
其中 , M 是延迟 , 程序中首先求出 x 的延时 M 的自卷积 aa, 通过赋值给 r , b , 然后通过最小平方估计 a =-r /b 求出AR, 进一步通过 MA 模型来估计参数 , 程序通过求出 AR 参数来计算最终ARM A 模型参数 [7-9], 如图 3所示。