古典法功率谱估计

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

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

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

4、n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号*k=fft(*n,1024); %进展2N-1点FFT,系统会自动补0Sk=abs(*k).*(abs(*k)./N; %取频谱幅度的平方，并除以N，以此作为对*n真实功率谱的估计Rn=ifft(Sk);Sk1=fft(Rn,512);figure(2)subplot(2,1,1);plot(n/N,Sk1);ylabel(Sk)title(b) 相关法估计功率谱密度)Sk2=10*log(Sk1); %对估计出的Sk取对数，使画出的图更加突出

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

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

7、ot(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个平均的方差趋于零，可以到达一致估计的目的。二运算简要框图矩形窗截断对求得结果进展平均周期图法对每段进展计算*(n)分成L小段输出三程序例如%L=2时bartlett法N=256,n=0:255;*1n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AW

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

10、ubplot(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;*1n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);%产生加有均值为0，方差为1的AWGN的信号的1/4段*k1=fft(*1n,N);%进展N点FFTSk1=abs(*k1).2./N;%取频谱幅度的平方，并除以N，以此作为对*n真实功率谱的估计n=128

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

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

13、;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法进展信号功率谱估计的波形相比我们可以很明显的看出L大的那个波形更加平坦，这与之前在算

14、法原理中介绍的一样，L越大，平均后的方差就越小，越能到达一致估计的目的。五、 Welch法功率谱估计一算法原理简介现在比较常用的功率谱估计改进方法是Welch法，又叫加权交叠平均法。这种方法以加窗加权求取平滑，以分段重叠求得平均，因此集平均与平滑的优点于一体，同时也不可防止的带有两者的缺点。其主要步骤如下：第一步：将N长的数据段分成L个小段，每小段M点，相邻小段见交叠M/2点。第二步：对个小段加同样的品挂窗后求傅里叶变换第三步：求个小段功率谱的平均，得这里二运算简要框图矩形窗截断相加求平均交叠进展功率谱估计 *n分成L小段输出三程序例如N=128;n=0:127;*1n=2*cos(2*pi*

15、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长窗函数的能量,那个函数不会写，这个数是自己加出来的。*11n=*1n.*wn;*k1=fft(*11n,128);Sk1=abs(*k1).2.*1/U;N=128;n=64:191;*2n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);*22n=*2n.*wn;*k2=fft(*22n,128);Sk2=abs(*k2).2.*(1/U

16、);N=128,n=128:255;*3n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);*33n=*3n.*wn;*k3=fft(*33n,128);Sk3=abs(*k3).2.*(1/U);N=128,n=192:319;*4n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);*44n=*4n.*wn;*k4=fft(*44n,128);Sk4=abs(*k4).2.*(1/U);N=128,n=256:383;*5n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.

17、4*n)+2*randn(size(n);*55n=*5n.*wn;*k5=fft(*55n,128);Sk5=abs(*k5).2.*(1/U);N=128,n=320:447;*6n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);*66n=*6n.*wn;*k6=fft(*66n,128);Sk6=abs(*k6).2.*(1/U);N=128,n=384:511;*7n=2*cos(2*pi*0.2*n)+ cos(2*pi*0.4*n)+2*randn(size(n);*77n=*7n.*wn;*k7=fft(*77n,128);Sk7=abs(*k7).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);plo