古典法功率谱估计.docx
- 文档编号:2279279
- 上传时间:2022-10-28
- 格式:DOCX
- 页数:16
- 大小:190.53KB
古典法功率谱估计.docx
《古典法功率谱估计.docx》由会员分享,可在线阅读,更多相关《古典法功率谱估计.docx(16页珍藏版)》请在冰豆网上搜索。
古典法功率谱估计
古典法功率谱估计
古典法功率谱估计
一、信号的产生
(一)信号组成
在本实验中,需要事先产生待估计的信号,为了使实验结果较为明显,我产生了由两个不同频率的正弦信号(频率差相对较大)和加性高斯白噪声组成的信号。
(二)程序
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)两个正弦信号与白噪声叠加的时域波形')
(三)信号波形
(二)运算简要框图
输出
矩形窗截断
相关法谱估计运算简要框图
图中快速相关的输出时从-(N-1)到(N-1)的2N-1点,加窗后截取的是-(M-1)到(M-1)的,最后做(2M-1)点FFT,即可得到结果。
(三)程序示例
程序的主要思路就是按照运算框图一步一步进行计算,下面附程序并进行简要解释:
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,系统会自动补0
Sk=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*log(Sk1);%对估计出的Sk取对数,使画出的图更加突出特点
subplot(2,1,2);
plot(n/N,Sk2);
ylabel('10log(PSD)')
(四)结果分析
下面是程序运行后的结果
从上图中我们可以较为明显的看到信号中有两个频率分量,一个在0.2处,一个在0.4处,与产生的信号相一致。
但是我们不难看出,估计出的功率谱谱线非常不平坦,有很多起伏。
三、周期图法谱估计
(一)算法原理简介
周期图法又称直接法。
它是从随机信号x(n)中截取N长的一段,把它视为能量有限x(n)真实功率谱的估计的抽样。
其具体步骤如下:
第一步:
由获得的N点数据构成有限长序列直接求傅里叶变换,得频谱。
第二步:
取频谱幅度的平方,并除以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点FFT
Sk3=abs(Xk1).*(abs(Xk1))./N;%取频谱幅度的平方,并除以N,以此作为对xn真实功率谱的估计
Sk4=10*log(Sk3);%对估计出的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法功率谱估计
(一)算法原理简介
当我们用相关法或者周期图法对信号的功率谱进行估计时,都不是对的一致估计,主要原因是方差大。
于是就产生了周期图法的改进。
改进的主要途径是平滑和平均。
平滑是用一个适当的窗函数与计算的功率谱进行卷积,是谱线平滑。
这种方法的出的谱估计是无偏的,方差也小,但分辨率下降。
平均就是将截取的数据段再分成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*randn(size(n));%产生加有均值为0,方差为1的AWGN的信号的前半段
Xk1=fft(x1n,N);%进行N点FFT
Sk5=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点FFT
Sk6=abs(Xk2).^2./N;%取频谱幅度的平方,并除以N,以此作为对xn真实功率谱的估计
Sk7=(Sk5+Sk6)/2;%相加求平均
Sk8=10*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点FFT
Sk1=abs(Xk1).^2./N;%取频谱幅度的平方,并除以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点FFT
Sk2=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=fft(x3n,N);%进行N点FFT
Sk3=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点FFT
Sk4=abs(Xk4).*(abs(Xk4))./N;%取频谱幅度的平方,并除以N,以此作为对xn真实功率谱的估计
Sk5=(Sk1+Sk2+Sk3+Sk4)/4;%相加求平均
Sk6=10*log(Sk5);
n=0:
127;
figure(5)
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法进行信号功率谱估计的波形相比我们可以很明显的看出L大的那个波形更加平坦,这与之前在算法原理中介绍的一样,L越大,平均后的方差就越小,越能达到一致估计的目的。
五、Welch法功率谱估计
(一)算法原理简介
现在比较常用的功率谱估计改进方法是Welch法,又叫加权交叠平均法。
这种方法以加窗(加权)求取平滑,以分段重叠求得平均,因此集平均与平滑的优点于一体,同时也不可避免的带有两者的缺点。
其主要步骤如下:
第一步:
将N长的数据段分成L个小段,每小段M点,相邻小段见交叠M/2点。
第二步:
对个小段加同样的品挂窗后求傅里叶变换
第三步:
求个小段功率谱的平均,得
这里
(二)运算简要框图
矩形窗截断
X(n)分成L小段输出
(三)程序示例
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(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=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,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 古典 功率 估计