谱估计基础及仿真分析.doc

谱估计基础及仿真分析0引言 现代信号分析中,对于常见的具有各态历经的平稳随机信号,不可能用清楚的数学关系式来描述,但可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD)。它是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计(非参数估计)和现代功率谱估计(参数估计)。功率谱估计在实际工程中有重要应用价值,如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。一 相关的数学基础1.1概率论：1.1.1多维高斯分布高斯分布公式：(1)准则为：内68.26%；：95.44%；：99.74%。向量的形式的公式为：(2)其中；(3)1.1.2 mont-carlo仿真为指数分布，瑞利分布，还有正交的高斯分布等。1.2 随即过程平稳随即过程：多维联合概率密度和时间起点无关，狭义的平稳。数字特征：相关函数：(4)协方差：(5)广义的平稳：，(6)各态历经性：时间平均代替集平均。平稳的功率谱：维纳辛钦定理(7)(8)能量有限；有限不连续点等条件才能傅立叶变换成立。下面给出几种随即过程：白噪声；只跟自己相关。随相正弦波；(9)高斯随机过程完备表达式。平稳白高斯的(10)1.3 数理统计1.3.1 参数估计无偏性：有效性：CRB：(11)一致性：1.3.2 方法1 数字特征法(12)(13)2 最大似然估计使似然函数最大的估计方法：(14)(15)从理论出发，渐近无偏；渐近有效；渐近高斯；不变性。最大似然估计高斯分布的均值和方差：(16)，有偏的，渐近无偏的。3 最小二乘估计已知样本模型和样本解方程：准则：，可以看作广义逆。4 线性最小均方误差估计 线性样本组合求(17)二：经典谱估计傅立叶变换(18)时间离散化：DTFT(19)离散DFT(20)功率谱估计经典的功率谱估计是将数据工作区外的未知数据假设为零，相当于数据加窗。经典功率谱估计方法分为：相关函数法、周期图法以及改进的周期图估计法即平均周期法和平滑周期图法；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(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.2周期图法把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。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);图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;%数据无重叠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);图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=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 Welch法功率谱估计三：参数模型谱估计现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。主要方法有最大嫡谱分析法(AR模型法)、Pisarenko谐波分解法、Prony提取极点法、Prony谱线分解法以及Capon最大似然法等。其中AR模型应用较多,具有代表性。常用的模型有ARMA模型、AR模型、MA模型。用周期图法求出的功率谱曲线和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(periodogram psd estimate);order1=50;range=half;magunits=db;subplot(122)pburg(xn,order1,nfft,fs,range);图4 Burg法估计通过实验仿真可以直观地看出以下特性：(1)功率谱估计中的相关函数法和周期图法所得到的结果是一致的，其特点是离散性大，曲线粗糙，方差较大，但是分辨率较高。(2)平均周期图法和平滑平均周期图法的收敛性较好，曲线平滑，估计的结果方差较小，但是功率谱主瓣较宽，分辨率低。这是由于对随机序列的分段处理引起了长度有限所带来的Gibbs现象而造成的。(3)平滑平均周期图法与平均周期图法相比，谱估值比较平滑，但是分辨率较差。其原因是给每一段序列用适当的窗口函数加权后，在得到平滑的估计结果的同时，使功率谱的主瓣变宽，因此分辨率有所下降。经典功率谱估计的分辨率反比于有效信号的长度，但现代谱估计的分辨率可以不受此限制。这是因为对于给定的N点有限长序列x(n)，虽然其估计出的自相关函数也是有限长的，但是现代谱估计的一些隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度，不像经典谱估计那样受窗函数的影响。因而现代谱的分别率比较高,而且现代谱线要平滑得多，从上图可以清楚看出。四 结论近十几年来，随着计算机的广泛应用，各种新算法不断提出，谱估计得到了迅速的发展,出现了许多功率谱估计的实现方法，也有很