功率谱估计及其MATLAB仿真詹红艳（201121070630控制理论与控制工程）摘要:从介绍功率谱的估计原理入手分析了经典谱估计和现代谱估计两类估计方法的原理、各自特点及在Matlab中的实现方法。
关键词:功率谱估计;周期图法;AR参数法;MatlabPower Spectrum Density Estimation and the simulation inMatlabZhan Hongyan（201121070630Control theory and control engineering）Abstract:Mainly introduces the principles of classical PSD estimation and modern PSD estimation,discusses the characteristics of the methods of realization in Matlab.Moreover,It gives an example of each part in realization using Matlab functions.Keywords:PSDPstimation,Periodogram method,AR Parameter method,Matlab1引言现代信号分析中,对于常见的具有各态历经的平稳随机信号,不可能用清楚的数学关系式来描述,但可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD)。
它是数字信号处理的重要研究内容之一。
功率谱估计可以分为经典功率谱估计(非参数估计)和现代功率谱估计(参数估计)。
功率谱估计在实际工程中有重要应用价值,如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。
Matlab是MathWorks公司于1982年推出的一套高性能的数值计算和可视化软件,人称矩阵实验室,它集数值分析、矩阵运算、信号处理和图形显示于一体,构成了一个方便的、界面友好的用户环境,成为目前极为流行的工程数学分析软件。
也为数字信号处理进行理论学习、工程设计分析提供了相当便捷的途径。
本文的仿真实验中,全部在Matlab6.5环境下调试通过;随机序列由频率不同的正弦信号加高斯白噪声组成。
2经典功率谱估计经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。
经典功率谱估计方法分为:相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法,其中周期图法应用较多,具有代表性。
1.1相关函数法(BT法)该方法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。
当延迟与数据长度相比很小时,可以有良好的估计精度。
Matlab代码示例1:Fs=500;%采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));nfft=512;cxn=xcorr(xn,'unbiased');%计CXk=fft(cxn,nfft);Pxx=abs(CXk);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));figure(1)%plot(k,plot_Pxx);1.2周期图法(periodogram)周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列,直接计算x(n)的离散傅立叶变换,得X(k),然后再取其幅值的平方,并除以N,作为序列(n)真实功率谱的估计。
Matlab代码示例2:Fs=600;%采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+0.1*randn(size(n));window=boxcar(length(xn));%矩形窗nfft=512;[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法plot(f,10*log10(Pxx));window=boxcar(length(xn));%矩形窗nfft=1024;[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法 figure(1)plot(f,10*log10(Pxx));1.3平均周期图法和平滑平均周期图法对于周期图的功率谱估计,当数据长度N 太大时,谱曲线起伏加剧,若N 太小,谱的分辨率又不好,因此需要改进。
两种改进的估计法是平均周期图法和平滑平均周期图法。
Bartlett 法:Bartlett 平均周期图的方法是将N 点的有限长序列x(n)分段求周期图再平均。
Matlab 代码示例3: fs=600; n=0:1/fs:1;xn=cos(2*pi*20*n)+3*cos(2*pi*90*n)+randn(size(n)); nfft=512;window=hamming(nfft);%矩形窗 noverlap=0;%数据无重叠 p=0.9;%置信概率[Pxx,Pxxc]=psd(xn,nfft,fs,window,noverlap,p); index=0:round(nfft/2-1); k=index*fs/nfft;plot_Pxx=10*log10(Pxx(index+1)); plot_Pxxc=10*log10(Pxxc(index+1)); figure(1)plot(k,plot_Pxx); figure(2)plot(k,[plot_Pxxplot_Pxx-plot_Pxxcplot_Pxx+plot_Pxxc]);Welch 法:Welch 法对Bartlett 法进行了两方面的修正,一是选择适当的窗函数w(n),并在周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。
二是在分段时,可使各段之间有重叠,这样会使方差减小。
Matlab 代码示例4: Fs=600; n=0:1/Fs:1;xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n)); nfft=512;window=boxcar(100);%矩形窗window1=hamming(100);%海明窗window2=blackman(100);%blackman 窗 noverlap=20;%数据无重叠range='half';%频率间隔为[0 Fs/2],计算一半的频率 [Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range); [Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range); [Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range); plot_Pxx=10*log10(Pxx); plot_Pxx1=10*log10(Pxx1); plot_Pxx2=10*log10(Pxx2); figure(1)plot(f,plot_Pxx); figure(2)plot(f,plot_Pxx1); figure(3);plot(f,plot_Pxx2);3 现代功率谱估计现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。
主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。
主要方法有最大嫡谱分析法(AR 模型法)、Pisarenko 谐波分解法、Prony 提取极点法、Prony 谱线分解法以及Capon 最大似然法等。
其中AR 模型应用较多,具有代表性。
常用的模型 有ARMA 模型、AR 模型、MA 模型。
2.1模型的公式表达ARMA 模型功率谱数学式子表达为:222()|1|/|1|ppjwk jwk k jwk x k lk lP e b ea e σ--===++∑∑其中σ2是激励白噪声的方差,Px(ej ω)为功率谱密度,ak 和bk 为模型参数。
如果ARMA 模型参数b1,b2......bq 全为0,就演化为AR 模型:22()/|1|pjwk jwk x k lP e a e σ-==+∑如果ARMA 模型参数a1,a2......aq 全为0,就演化为MA 模型:22()/|1|pjwk jwk x k lP e b e σ-==+∑在实际中,AR 模型的参数估计比较简单,对其有充分的研究,而对于ARMA 模型,其参数比较复杂,对其算法的研究和改进还在完善中。
2.2 AR 模型功率谱估计AR 模型功率谱估计又称为自回归模型,它是一个全极点的模型,要利用AR 模型进行功率谱估计须通过levinson_dubin 递推算法由Yule-Walker 方程求得AR 的参数:σ2α1α2...αp 。
在Matlab 仿真中可调用Pburg 函数直接画出基于burg 算法的功率谱估计的曲线图。
Matlab 代码示例5:用周期图法求出的功率谱曲线和burg 算法求出的AR 功率谱曲线(p=50) fs=200; n=0:1/fs:1;xn=cos(2*pi*40*n)+cos(2*pi*41*n)+3*cos(2*pi*90*n)+ 0.1*randn(size(n));window=boxcar(length(xn)); nfft=512;[pxx,f]=periodogram(xn,window,nfft,fs); subplot(121)plot(f,10*log10(pxx)) xlabel('frequency(hz)');ylabel('power spectral density(Db/Hz)'); title('periodogrampsd estimate'); order1=50; range='half'; magunits='db'; subplot(122)pburg(xn,order1,nfft,fs,range))周期图法求出的功率谱曲线和burg 算法求出的 AR 功率谱曲线(p=50)如图所示:4 常见谱估计法的比较通过实验仿真可以直观地看出以下特性:(1)功率谱估计中的相关函数法和周期图法所得到的结果是一致的,其特点是离散性大,曲线粗糙,方差较大,但是分辨率较高。