作业：谱分析.doc

随机信号的两种功率谱估计对比 姓 名：董兴蒙 学 号：S11010190 班 级：地研11-62011年12月1日随机信号的两种功率谱估计对比摘要：功率谱估计是信号处理领域的重要问题之一。功率谱估计利用以观测到的一定数量样本数据估计一个平稳随机信号的功率谱密度，因其能够分析信号的功率谱密度，因其能够分析信号的能量随频率变化的分布特性，在许多实际应用中功率谱的分析与估计已变得越来越重要。 功率谱估计可以分为经典谱估计方法与现代谱估计方法。本文首先阐述了两种谱估计方法的原理，然后用一随机信号采用matlab进行了编程实现，通过各种谱估计方法的结果，进一步阐述了不同谱估计方法的优点和实用性。关键词：经典谱估计、现代谱估计、周期图法、Welch法、最大熵法、AR模型法第一章 经典谱估计与现代谱估计的原理概述1、1 引言功率谱密度研究二阶平稳随机过程特征，揭示随机过程中与隐含的周期及相邻的谱峰等有用信息；用有限长的N个样本数据点估计该平稳随机过程的功率谱密度称为功率谱估计。此种估计是建立在时间平均的方法之上，并假定具有遍历性。功率谱估计，简称谱估计，涉及信号与系统、随机信号分析、概率统计、随机过程、线性代数等基础学科，广泛应用于雷达、声纳、通信、地质勘探、天文、生物医学工程等众多领域。其内容和方法不断更新，是非常活跃的研究领域。谱估计的方法大体分为参数法和非参数法，重点介绍常用的几种。经典谱估计线性非参数化方法：周期图法，相关图法等，采用经典的傅氏变换及窗口截断，对长序列有良好的估计；现代谱估计非线性参数化方法：最大似然估计等，对短序列有良好估计；1、2 经典谱估计经典谱估计是采用经典的傅里叶变换及窗口截断对信号功率谱进行估计。其中最简单的就是周期图法，又分为直接法与间接法。直接法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换（即频谱），然后取频谱与其共轭的乘积得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计；间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。经典功率谱估计对长序列有良好估计。对于周期图法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此我们需要对方法进行改进。我们采用了改进的周期图法bartlett法和Welch法。bartlett法将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并在周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。welch法改善了估计谱曲线的光滑性，大大提高谱估计的分辨率。1、3 现代谱估计现代谱估计主要是针对经典谱估计的分辨率差和方差性能不好的问题而提出的。现代谱估计从方法上大致可分为参数模型谱估计和非参数模型谱估计两种，前者有AR模型、MA模型、ARMA模型、PRONY指数模型等；后者有最小方差方法、多分量的MUSIC方法等。常用的一种现代谱估计方法是最大熵谱估计。它是把信息熵的概念收入信号处理中，有时又称为时序谱分析方法。这是一种自相关函数外推的方法，在分析过程中，没有固定的窗函数。在每一步外推过程中，使估计的相关函数包含过程的信息最多，即要求在过程的熵达到最大的条件，确定未知的自相关函数值，借以达到谱估计的逼真和稳定程度最好的目的。也就是采用谱熵为最大的准则来估计功率谱。同传统的谱估计方法相比，最大熵谱（MEM谱）不仅没有传统谱受到数据加窗这一致命弱点带来的一系列缺陷，而且由于它是连续谱，从理论上讲，谱光滑，谱峰陡峭，频率分辨力无限高，不存在传统谱是离散谱带来频率分辨误差这一问题，对传统谱而言，样本数越短，一方面是窗影响越加严重，另一方面是频率分辨误差成比例增大，以致在短样本时传统谱无法应用，传统谱对许多瞬变过程的研究无能为力，而MEM谱不存在这些问题，所以MEM谱适用于短数据缓慢变化过程的谱估计，MEM谱的峰值可靠性差。第二章 两种谱估计的编程举例对比2、1 原始随机信号 采用含有噪声序列的随机信号,其中信号的采样频率为Fs=1000Hz, xn=cos(2*pi*50*n)+3*cos(2*pi*125*n)+randn(size(n)，对两种谱估计方法进行了对比。 2-1 原始随机信号 2-2 原始随机信号（连续函数形式）Matlab代码示例：clear;Fs=1000; %采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*50*n)+3*cos(2*pi*125*n)+randn(size(n);%plot(n,xn);stem(n,xn);title(原始信号);ylabel(xn);xlabel(n)grid;2.2经典谱估计方法2.2.1直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。定义：长度为N的实平稳随机信号序列的周期图为： 其中，由于的DFT有周期性，所以也有周期性，是有偏估计。（1）先计算N个数据的Fourier变换（即频谱）：（2）然后取频谱和其共轭的乘积，得到功率谱； 2-3 周期图法谱估计Matlab代码示例：clear;Fs=1000; %采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*50*n)+3*cos(2*pi*125*n)+randn(size(n);window=boxcar(length(xn); %矩形窗nfft=1024;Pxx,f=periodogram(xn,window,nfft,Fs); %直接法plot(f,10*log10(Pxx);title(周期图法谱估计);xlabel(frequency);ylabel(Power Spectrum (dB)grid;2.2.2间接法（BT法） 间接法又称相关图法。由维纳定理：先估计出有限长信号x(n)的自相关函数，即： ，它可以看成是与的褶积运算结果除以N，即：，两边取Z变换，求的DFT，得到的功率谱估计： ;由于功率谱是自相关函数间接求出的，称间接法。为直接法，令，适用于FFT来计算。先根据N个样本数据估计x(n)的样本自相关函数： 得到功率谱： （2）由于在计算（1）（2）两式的Fourier变换时，分别将x(n)和视作周期函数。由直接法和间接法估计的功率谱常称为周期图。周期图方法估计的功率谱为有偏估计，为减少其偏差，通常需要加窗函数对周期图进行平滑。加窗函数有两种不同的方法：一种是直接将窗函数C(n)直接加在样本函数，得到功率谱称为修正周期图，定义为：是窗函数的Fourier变换。另一种方法是将窗函数加在样本自相关函数，得到的功率谱称为周期图平滑。BT法其功率谱定义为： 直接给数据加窗函数称为数据窗，而加给自相关函数的窗函数称为滞后窗，其Fourier变换则称作谱窗。加窗函数虽然能够减少周期图的偏差，改善功率谱曲线的光滑性，但作为参数化谱估计，周期图具有分辨率低的因有缺陷不能适应高分辨率功率谱估计的需要。2-4相关图法谱估计Matlab代码示例：%间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n) %的功率谱估计。clear;Fs=1000; %采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*50*n)+3*cos(2*pi*125*n)+randn(size(n); nfft=1024;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);plot(k,plot_Pxx);title(相关图法谱估计);xlabel(frequency);ylabel(Power Spectrum (dB)grid;2.2.3直接法和间接法的比较直接法进行谱估计，是有偏估计，由于卷积的运算过程会导致功率谱真实值的尖峰附近产生泄漏，相对地平滑了尖峰值，因此造成谱估计的失真。另外，当时，功率谱估计的方差不为零，所以不是一致性估计。并且功率谱估计在等于整数倍的各数字频率点互不相关。其谱估计的波动比较显著，特别是当N 越大、 越小时，波动越明显。但如果N 取得太小，又会造成分辨率的下降。2.2.4直接法的改进对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。1. Bartlett法Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。 2-5 改进的直接法Bartlett法谱估计Matlab代码示例：clear；Fs=1000;n=0:1/Fs:1;xn=cos(2*pi*50*n)+3*cos(2*pi*125*n)+randn(size(n);nfft=1024;window=boxcar(length(n); %矩形窗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);pause;title(改进的周期图法-bartlett法);xlabel(frequency);ylabel(Power Spectrum (dB)grid;figure(2)plot(k,plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc);2. Welch法Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。2-6改进的直接法Welch法谱估计Matlab代码示例：clear;Fs=1000; %采样频率.nfft=1024;n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*50*n)+3*cos(2*pi*125*n)+randn(size(n);window=rectwin(100); %矩形窗window1=hamming(100); %海明窗window2=blackman(100); %blackman窗noverlap=20; %数据无重叠range=half; %频率间隔为0 Fs/2，只计算一半的频率Pxx3,f3=pwelch(xn,window,noverlap,nfft,Fs,range);Pxx4,f4=pwelch(xn,window1,noverlap,nfft,Fs,range);Pxx5,f5=pwelch(xn,window2,noverlap,nfft,Fs,range);plot_Pxx3=10*log10(Pxx3);plot_Pxx4=10*log10(Pxx4);plot_Pxx5=10*log10(Pxx5);subplot(3,1,1);plot(f3,plot_Pxx3);title(Welch法-矩形窗);xlabel(normalized frequencies);ylabel(Power Spectrum (dB)grid;subplot(3,1,2);plot(f4,plot_Pxx4);title(Welch法-海明窗);xlabel(normalized frequencies);ylabel(Power Spectrum (dB)grid;subplot(3,1,3);plot(f5,plot_Pxx5);title(Welch法-blackman窗);xlabel(normalized frequencies);ylabel(Power Spectrum (dB);grid;2.3常用的现代功率谱估计方法2.3.1 AR模型AR模型又称为自回归模型，它是一个全相点的模型，可用差分方程如下表示： （1）式中是均值为零，方差为的白噪声序列，p为AR模型的阶数，为p阶模型的参数。就以序列可以看作是白噪声通过AR模型系统的输出。可得AR模型系统的转移函数数学表达式： （2）而得到AR模型功率谱估计的计算方式：其中为白噪声序列的方差。从该式可以看出，要利用AR模型进行功率谱估计，必须先求得AR模型的参数及白噪声的方差，对（1）两边同乘以，并求其均值，可得到（2）所表示的平稳信号序列的自相关函数：写成矩阵形式：实时傅氏变换（STFT）用specgram函数tfrstf函数。 2 -7 AR模型算法Matlab代码示例clear;Fs=1000; %采样频率% nfft=1024;n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*50*n)+3*cos(2*pi*125*n)+randn(size(n);%最大熵谱估计Maxmum Entropy Methodorder=20;Nfft=256; Pxx7,f7=pyulear(xn,order,Nfft,Fs,onesided); plot(f7,10*log10(Pxx7); xlabel(Frequency (Hz); ylabel(Power Spectrum (dB); title(AR模型法); Grid;2.3.1 最大熵谱估计 就是根据已知数据信息，在不进行任何新的假设(不增加任何虚假信息)的情况下，合理地预测未知延迟离散时间上的相关函数。即在根据已知信息外推相关函数时，每一步都保持未知事件的不确定性或熵为最大。 可见熵是消息源发出每个消息的平均信息量。 对于高斯分布的随机变量，布卡乔夫证明了其熵和自协方差矩阵间存在关系： 当时间序列为零均值时，熵和自相关函数之间存在关系 ： 当过程为无限长时，用熵率作为信息的度量： 。时间序列功率谱密度和熵率的关系： 。 时间序列的频率范围是-fc ，fc 。 若已知自相关函数Rx(m)的前2M+1个序列值，则选择未知自相关函数要使：， 。 从而可以外推出Rx(M+1)。并依此类推得到其它自相关函数值。于是功率谱 。若选择，可以得到，。令，整理后得到，。令，。最大熵谱估计：。一、最大熵谱估计 2-8 最大熵谱估计Matlab代码示例clear;Fs=1000; %采样频率% nfft=1024;n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*50*n)+3*cos(2*pi*125*n)+randn(size(n);%最大熵谱估计Maxmum Entropy Methodorder=20;Nfft=256;Pxx6,f6=pmem(xn,order,Nfft,Fs);plot(f6,10*log10(Pxx6);xlabel(Frequency (Hz);ylabel(Power Spectrum (dB);title(最大熵谱估计);grid;二、最大熵谱估计的Burg法 2-9 最大熵谱估计的Burg法Matlab代码示例 clear; Fs=1000; %采样频率% nfft=1024; n=0:1/Fs:1; %产生含有噪声的序列 xn=cos(2*pi*50*n)+3*cos(2*pi*125*n)+randn(size(n); %最大熵谱估计Maxmum Entropy Method order=20; Nfft=256;Pxx8,f8=pburg(xn,order,Nfft,onesided);plot(f8*5