《【精选】谱估计基础及仿真分析》由会员分享，可在线阅读，更多相关《【精选】谱估计基础及仿真分析（9页珍藏版）》请在金锄头文库上搜索。
1、谱估计基础及仿真分析0 引言 现代信号分析中 ,对于常见的具有各态历经的平稳随机信号,不可能用清楚的数学关系式来描述,但可以利用给定的 N 个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD) 。它是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计(非参数估计)和现代功率谱估计(参数估计) 。功率谱估计在实际工程中有重要应用价值,如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。一 相关的数学基础1.1 概率论：1.1.1
2、多维高斯分布高斯分布公式： (1)2()1()xpxe准则为： 内 68.26%； ：95.44% ； ：99.74%。3,3向量的形式的公式为： (2)11/21/2()exp()()| TxxNpx CC其中 (TxxxCE； (3)xijiijjEE1.1.2 mont-carlo 仿真为指数分布， 瑞利分布，还有正交的高斯分布等。1ln,(0)yu2lnxyu1.2 随即过程平稳随即过程：多维联合概率密度和时间起点无关，狭义的平稳。数字特征： ()tExt相关函数： (4)1212(,)()xrxt协方差： (5)122()xxCttt广义的平稳： ， (6)1221(,)()(xxr
3、rtr,()C各态历经性：时间平均代替集平均。平稳的功率谱：维纳辛钦定理(7)()()jxxPred(8)12jr能量有限；有限不连续点等条件才能傅立叶变换成立。下面给出几种随即过程：白噪声 ；只跟自己相关。2(0)()xxrP随相正弦波； (9)020()cos(),(,2)()()xxtAftrP高斯随机过程完备表达式。101/21/2(.)exp()()(| TN xNp CC平稳白高斯的 (10)2210x I1.3 数理统计1.3.1 参数估计无偏性： E有效性：CRB ： (11)21var()ln(;)px一致性： lim(|)0Np1.3.2 方法1 数字特征法(12)1Nkx
4、(13)221()k2 最大似然估计使似然函数最大的估计方法：(14)(;)0px(15)ln从理论出发，渐近无偏；渐近有效；渐近高斯；不变性。最大似然估计高斯分布的均值和方差：(16)1Nkx，有偏的，渐近无偏的。221()k3 最小二乘估计已知样本模型和样本解方程： xH准则： ，可以看作广义逆。21|()Hx4 线性最小均方误差估计线性样本组合求10*Nix(17)2|minE1xCr二：经典谱估计傅立叶变换 ()()jtXxed(18)12jtt时间离散化：DTFT(19)2()()1jnnjXxexd离散 DFT(20)210()()NjknNnjkkXxex功率谱估计经典的功率谱估
5、计是将数据工作区外的未知数据假设为零，相当于数据加窗。经典功率谱估计方法分为：相关函数法、周期图法以及改进的周期图估计法即平均周期法和平滑周期图法；2.1 相关函数法该方法先由序列估计出自相关函数 R，然后对其进行傅立叶变换，得到功率谱估计。当延时和数据长度相比很小的时候可以有良好的估计精度。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(CX
6、k);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1);figure(1)plot(k,plot_Pxx);0 50 100 150 200 250-10-505101520253035 估估估估估估估估估估图 1 相关函数功率谱估计2.2 周期图法把随机序列 x(n)的 N 个观测数据视为一能量有限的序列，直接计算 x(n)的离散傅立叶变换，得 X(k)，然后再取其幅值的平方，并除以 N，作为序列 x(n)真实功率谱的估计。Fs=600;%采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos
7、(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);0 50 100 150 200 250 300-70-60-50-40-30-20-
8、10010 估估估估估估估估估图 2 周期图法功率谱估计2.3 平均周期图法和平滑平均周期图法对于周期图的功率谱估计，当数据长度 N 太大时，谱曲线起伏加剧，若 N 太小，谱的分辨率又不好，因此需要改进。两种改进的估计法是平均周期图法和平滑平均周期图法。Bartlett 法：Bartlett 平均周期图的方法是将 N 点的有限长序列 x(n)分段求周期图再平均。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;%数据无重叠
9、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_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc);0 50 100 150 200 250 300-30-20-100102030 估估估估估估_Bart
10、lett估_估估估估估0 50 100 150 200 250 300-50-40-30-20-1001020304050 估估估估估估_Bartlett估_估估估估估图 3 平均周期图法功率谱估计Welch 法： Welch 法对 Bartlett 法进行了两方面的修正：一是选择适当的窗函数 w(n)，并在周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时,可使各段之间有重叠，这样会使方差减小。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=
11、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=
12、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);0 50 100 150 200 250 300-35-30-25-20-15-10-50 估估估估估估_Welch估_估估估估估0 50 100 150 200 250 300-35-30-25-20-15-10-50 估估估估估估_Welch估_估估估估估0 50 100 150 200 250 300-40-35-30-25-20-15-10-50 估
13、估估估估估_Welch估_估估估估估图 3 Welch 法功率谱估计三：参数模型谱估计现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。主要方法有最大嫡谱分析法(AR 模型法) 、Pisarenko 谐波分解法、Prony 提取极点法、Prony 谱线分解法以及 Capon 最大似然法等。其中 AR 模型应用较多,具有代表性。常用的模型有 ARMA 模型、AR 模型、MA 模型。用周期图法求出的功率谱曲线和 burg 算法求出的 AR 功率谱曲线(p=50)fs=200;n=0:1/f
14、s: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(periodogram psd estimate);order1=50;range=half;magunits=db;su
15、bplot(122)pburg(xn,order1,nfft,fs,range);0 50 100-60-50-40-30-20-10010frequency(hz)power spectral density(Db/Hz)periodogram psd estimate0 20 40 60 80 100-60-50-40-30-20-10010Frequency (Hz)Power/frequency (dB/Hz)Power Spectral Density Estimate via Burg图4 Burg法估计通过实验仿真可以直观地看出以下特性：(1)功率谱估计中的相关函数法和周期图法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。(2)平均周期图法和平滑平均周期图法的收敛性较好，曲线平滑，估计的结果方差较小，但是功率谱主瓣较宽，分辨率低。这是由于对随机序列的分段处理引起了长度有限所带来的Gibbs现象而造成的。(3)平滑平均周期图法与平均周期图法相比，谱估值比较平滑，但是分辨率较差。其原因是给每一段序列用适当的窗口函数加权后，在得到平滑的估计结果的同时，使功率谱的主瓣变宽，因此分辨率有所下降
