古典法功率谱估计.docx

1、古典法功率谱估计古典法功率谱估计一、 信号的产生（一）信号组成在本实验中，需要事先产生待估计的信号，为了使实验结果较为明显，我产生了由两个不同频率的正弦信号（频率差相对较大）和加性高斯白噪声组成的信号。（二）程序xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN信号figure(1)plot(n,xn);title(a) 两个正弦信号与白噪声叠加的时域波形)（三） 信号波形二、相关法功率谱估计（一）算法原理简介此方法以相关函数为媒介来计算功率谱，所以又叫间接法。它是1958年由Blackman 和T

2、ukey提出。这种方法的具体步骤是： 第一步：从无限长随机序列x(n)中截取长度N的有限长序列 第二步：由N长序列求（2M-1）点的自相关函数序列。 第三步：由相关函数的傅式变换求功率谱。 以上过程中经历了两次截断，一次是将x(n)截成N长，称为加数据窗，一次是将想x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计。一般取MN，因为只有当M较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。因此，在FFT问世之前，相关法是最常用的谱估计方法。当FFT问世后，情况有所变化。因为截断后的)(nxN可视作能量信号，由相关卷积定理可

3、得这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。我们可对上式两边取（2N-1）点DFT，则有于是将时域卷积变为频域乘积，用快速相关求自相关函数估值的完整方案如下： 1. 对N长)(nxN的充（N-1）个零，成为（2N-1）长的。 2. 求（2N-1）点的FFT，得3、求 4、 求（2N-1）点的IFFT（二）运算简要框图 输出矩形窗截断相关法谱估计运算简要框图图中快速相关的输出时从-（N-1）到（N-1）的2N-1点，加窗后截取的是-（M-1）到（M-1）的 ，最后做（2M-1）点FFT，即可得到结果。（三）程序示例程序的主要思路就是按照运算框图一步一步进行计算，下面附程序并进行简

4、要解释：N=512,n=0:N-1; %N是FFT的变换区间xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号Xk=fft(xn,1024); %进行2N-1点FFT,系统会自动补0Sk=abs(Xk).*(abs(Xk)./N; %取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计Rn=ifft(Sk);Sk1=fft(Rn,512);figure(2)subplot(2,1,1);plot(n/N,Sk1);ylabel(Sk)title(b) 相关法估计功率谱密度)Sk2=10*

5、log(Sk1); %对估计出的Sk取对数，使画出的图更加突出特点subplot(2,1,2);plot(n/N,Sk2);ylabel(10log(PSD)（四）结果分析 下面是程序运行后的结果从上图中我们可以较为明显的看到信号中有两个频率分量，一个在0.2处，一个在0.4处，与产生的信号相一致。但是我们不难看出，估计出的功率谱谱线非常不平坦，有很多起伏。三、 周期图法谱估计（一）算法原理简介周期图法又称直接法。它是从随机信号x(n)中截取N长的一段，把它视为能量有限x(n)真实功率谱的估计的抽样。其具体步骤如下：第一步：由获得的N点数据构成有限长序列直接求傅里叶变换，得频谱。第二步：取频谱

6、幅度的平方，并除以N，以此作为对x(n)真实功率谱的估计。事实上，周期图法谱估计与自相关法谱估计的差异只是估计自相关函数的方法不同。（二）运算简要框图 矩形窗（长度N）截断X(n) 图中用FFT来代替傅里叶变换（三）程序示例N=512,n=0:N-1;xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号Xk1=fft(xn,512); %进行N点FFTSk3=abs(Xk1).*(abs(Xk1)./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计Sk4=10*log(Sk3);

7、 %对估计出的Sk取对数，使画出的图更加突出特点figure(3)subplot(2,1,1);plot(n/N,Sk3);title(c) 周期图法估计功率谱密度)ylabel(Sk)subplot(2,1,2);plot(n/N,Sk4);ylabel(10log(PSD)（四）结果分析下面是程序运行后的结果从上图中我们同样可以看到信号中有两个频率分量，一个在0.2处，一个在0.4处，与产生的信号相一致。但是我们不难看出，估计出的功率谱谱线与相关法功率谱估计一样非常不平坦，有很多起伏。四、 Bartlett法功率谱估计（一）算法原理简介当我们用相关法或者周期图法对信号的功率谱进行估计时，都

8、不是对的一致估计，主要原因是方差大。于是就产生了周期图法的改进。改进的主要途径是平滑和平均。平滑是用一个适当的窗函数与计算的功率谱进行卷积，是谱线平滑。这种方法的出的谱估计是无偏的，方差也小，但分辨率下降。平均就是将截取的数据段再分成L个小段，分别计算功率谱后取功率谱的平均。因为L个平均的方差比随机变量的单独方差小L倍，所以当L趋于无穷时，L个平均的方差趋于零，可以达到一致估计的目的。（二）运算简要框图 矩形窗截断X(n) 分成L小段 输出（三）程序示例%L=2时bartlett法N=256,n=0:255;x1n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*ra

9、ndn(size(n);%产生加有均值为0，方差为1的AWGN的信号的前半段Xk1=fft(x1n,N); %进行N点FFTSk5=abs(Xk1).2./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计n=256:511;x2n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号后半段Xk2=fft(x2n,N); %进行N点FFTSk6=abs(Xk2).2./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计Sk7=(Sk5+Sk6)/2; %相加求平均Sk8=1

10、0*log(Sk7);n=0:255;figure(4)subplot(2,1,1);plot(n/N,Sk7);title(d) Bartlett法估计功率谱密度 L=2)ylabel(Sk)subplot(2,1,2);plot(n/N,Sk8);ylabel(10log(PSD)%L=4时bartlett法N=128,n=0:127;x1n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号的1/4段Xk1=fft(x1n,N);%进行N点FFTSk1=abs(Xk1).2./N;%取频谱幅度

11、的平方，并除以N，以此作为对xn真实功率谱的估计n=128:255;x2n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号的1/4段Xk2=fft(x2n,N);%进行N点FFTSk2=abs(Xk2).*(abs(Xk2)./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计N=128,n=256:383;x3n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号的1/4段Xk3=ff

12、t(x3n,N);%进行N点FFTSk3=abs(Xk3).2./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计n=384:511;x4n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号的1/4段Xk4=fft(x4n,N);%进行N点FFTSk4=abs(Xk4).*(abs(Xk4)./N;%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计Sk5=(Sk1+Sk2+Sk3+Sk4)/4; %相加求平均Sk6=10*log(Sk5);n=0:127;figure(5

13、)subplot(2,1,1);plot(n/N,Sk5);title(e) Bartlett法估计功率谱密度 L=4)ylabel(Sk)subplot(2,1,2);plot(n/N,Sk6);ylabel(10log(PSD)（四）结果分析下面是程序运行后的结果上图分别是L=2和L=4时用bartlett法进行信号功率谱估计的波形。从上图中我们同样可以看到信号中有两个频率分量，一个在0.2处，一个在0.4处，与产生的信号相一致。但是我们不难看出，估计出的功率谱谱线与之前相关法功率谱估计和周期图法功率谱估计相比，波形相对平坦了一些。L=2和L=4时用bartlett法进行信号功率谱估计的波

14、形相比我们可以很明显的看出L大的那个波形更加平坦，这与之前在算法原理中介绍的一样，L越大，平均后的方差就越小，越能达到一致估计的目的。五、 Welch法功率谱估计（一）算法原理简介现在比较常用的功率谱估计改进方法是Welch法，又叫加权交叠平均法。这种方法以加窗（加权）求取平滑，以分段重叠求得平均，因此集平均与平滑的优点于一体，同时也不可避免的带有两者的缺点。其主要步骤如下：第一步：将N长的数据段分成L个小段，每小段M点，相邻小段见交叠M/2点。第二步：对个小段加同样的品挂窗后求傅里叶变换第三步：求个小段功率谱的平均，得这里 （二）运算简要框图 矩形窗截断 X（n） 分成L小段 输出（三）程序

15、示例N=128;n=0:127;x1n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号一部分wn=hanning(128);L=7;U=61.7942;%U是128长窗函数的能量,那个函数不会写，这个数是自己加出来的。x11n=x1n.*wn;Xk1=fft(x11n,128);Sk1=abs(Xk1).2.*1/U;N=128;n=64:191;x2n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);x22n=x2n.*wn;Xk2=fft

16、(x22n,128);Sk2=abs(Xk2).2.*(1/U);N=128,n=128:255;x3n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);x33n=x3n.*wn;Xk3=fft(x33n,128);Sk3=abs(Xk3).2.*(1/U);N=128,n=192:319;x4n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);x44n=x4n.*wn;Xk4=fft(x44n,128);Sk4=abs(Xk4).2.*(1/U);N=128,n=256:383;x5n

17、=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);x55n=x5n.*wn;Xk5=fft(x55n,128);Sk5=abs(Xk5).2.*(1/U);N=128,n=320:447;x6n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);x66n=x6n.*wn;Xk6=fft(x66n,128);Sk6=abs(Xk6).2.*(1/U);N=128,n=384:511;x7n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);

18、x77n=x7n.*wn;Xk7=fft(x77n,128);Sk7=abs(Xk7).2.*(1/U);n=0:127;Sk=(Sk1+Sk2+Sk3+Sk4+Sk5+Sk6+Sk7)./L;figure(6)title(f) Welch法（汉宁窗，L=7,64点交叠）)subplot(2,1,1);plot(n/N,Sk);ylabel(Sk)Skk=10*log(Sk);subplot(2,1,2);plot(n/N,Skk);ylabel(10log(PSD)（四）结果分析下面是程序运行后的结果上图分别是L=7，64点交叠用welch法进行信号功率谱估计的波形。从上图中我们同样可以看到信号中有两个频率分量，一个在0.2处，一个在0.4处，与产生的信号相一致。并且我们不难看出，估计出的功率谱谱线与之前上午功率谱估计相比，波形是最平坦的。可见这种改进方法是比较准确的。六、各波形直观比较 由以上六个图可以很清晰的看出不同功率谱估计的方法的优劣。