古典法功率谱估计Word格式.docx
《古典法功率谱估计Word格式.docx》由会员分享，可在线阅读，更多相关《古典法功率谱估计Word格式.docx（16页珍藏版）》请在冰豆网上搜索。
相关法谱估计运算简要框图
图中快速相关的输出时从-（N-1）到（N-1）的2N-1点，加窗后截取的是-（M-1）到（M-1）的，最后做（2M-1）点FFT，即可得到结果。
（三）程序示例
程序的主要思路就是按照运算框图一步一步进行计算，下面附程序并进行简要解释：
N=512,n=0:
N-1;
%N是FFT的变换区间
%产生加有均值为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'
（b）相关法估计功率谱密度'
Sk2=10*log（Sk1）;
%对估计出的Sk取对数，使画出的图更加突出特点
subplot（2,1,2）;
plot（n/N,Sk2）;
10log（PSD）'
（四）结果分析
下面是程序运行后的结果
从上图中我们可以较为明显的看到信号中有两个频率分量，一个在0.2处，一个在0.4处，与产生的信号相一致。
但是我们不难看出，估计出的功率谱谱线非常不平坦，有很多起伏。
三、周期图法谱估计
（一）算法原理简介
周期图法又称直接法。
它是从随机信号x（n）中截取N长的一段，把它视为能量有限x（n）真实功率谱的估计的抽样。
其具体步骤如下：
第一步：
由获得的N点数据构成有限长序列直接求傅里叶变换，得频谱。
第二步：
取频谱幅度的平方，并除以N，以此作为对x（n）真实功率谱的估计。
事实上，周期图法谱估计与自相关法谱估计的差异只是估计自相关函数的方法不同。
矩形窗（长度N）截断
X（n）
图中用FFT来代替傅里叶变换
Xk1=fft（xn,512）;
%进行N点FFT
Sk3=abs（Xk1）.*（abs（Xk1））./N;
%取频谱幅度的平方，并除以N，以此作为对xn真实功率谱的估计
Sk4=10*log（Sk3）;
figure（3）
plot（n/N,Sk3）;
（c）周期图法估计功率谱密度'
plot（n/N,Sk4）;
下面是程序运行后的结果
从上图中我们同样可以看到信号中有两个频率分量，一个在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）;
Sk5=abs（Xk1）.^2./N;
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）;
Sk6=abs（Xk2）.^2./N;
Sk7=（Sk5+Sk6）/2;
%相加求平均
Sk8=10*log（Sk7）;
n=0:
figure（4）
plot（n/N,Sk7）;
（d）Bartlett法估计功率谱密度L=2'
plot（n/N,Sk8）;
%L=4时bartlett法
N=128,n=0:
127;
%产生加有均值为0，方差为1的AWGN的信号的1/4段
%进行N点FFT
Sk1=abs（Xk1）.^2./N;
n=128:
Sk2=abs（Xk2）.*（abs（Xk2））./N;
N=128,n=256:
383;
x3n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;
Xk3=fft（x3n,N）;
Sk3=abs（Xk3）.^2./N;
n=384:
x4n=2*cos（2*pi*0.2*n）+cos（2*pi*0.4*n）+2*randn（size（n））;
Xk4=fft（x4n,N）;
Sk4=abs（Xk4）.*（abs（Xk4））./N;
Sk5=（Sk1+Sk2+Sk3+Sk4）/4;
Sk6=10*log（Sk5）;
figure（5）
plot（n/N,Sk5）;
（e）Bartlett法估计功率谱密度L=4'
plot（n/N,Sk6）;
上图分别是L=2和L=4时用bartlett法进行信号功率谱估计的波形。
但是我们不难看出，估计出的功率谱谱线与之前相关法功率谱估计和周期图法功率谱估计相比，波形相对平坦了一些。
L=2和L=4时用bartlett法进行信号功率谱估计的波形相比我们可以很明显的看出L大的那个波形更加平坦，这与之前在算法原理中介绍的一样，L越大，平均后的方差就越小，越能达到一致估计的目的。
五、Welch法功率谱估计
现在比较常用的功率谱估计改进方法是Welch法，又叫加权交叠平均法。
这种方法以加窗（加权）求取平滑，以分段重叠求得平均，因此集平均与平滑的优点于一体，同时也不可避免的带有两者的缺点。
其主要步骤如下：
将N长的数据段分成L个小段，每小段M点，相邻小段见交叠M/2点。
对个小段加同样的品挂窗后求傅里叶变换
第三步：
求个小段功率谱的平均，得
这里
矩形窗截断
X（n）分成L小段输出
N=128;
%产生加有均值为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=64:
191;
x22n=x2n.*wn'
Xk2=fft（x22n,128）;
Sk2=abs（Xk2）.^2.*（1/U）;
N=128,n=128:
x33n=x3n.*wn'
Xk3=fft（x33n,128）;
Sk3=abs（Xk3）.^2.*（1/U）;
N=128,n=192:
319;
x44n=x4n.*wn'
Xk4=fft（x44n,128）;
Sk4=abs（Xk4）.^2.*（1/U）;
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,
升级会员 