基于matlab的功率谱估计Word下载.docx
《基于matlab的功率谱估计Word下载.docx》由会员分享，可在线阅读，更多相关《基于matlab的功率谱估计Word下载.docx（23页珍藏版）》请在冰豆网上搜索。
（3）将计算结果表示成图形的形式，给出信号谱的分布情况图。
3.谱估计方法简介
（1）周期图法
周期图法是直接建立在功率谱的定义式上的，也称之为直接法。
原理计算如下：
a）取N点数据的DTFT（DFT）；
b）求模之平方并除以N；
（2）自相关法
自相关法的原理是由维纳-辛钦公式，经自相关函数间接获得的。
a）x（n），N点Þ
2N-1点，得Rxx（m）；
b）按2N-1点对Rxx（m）；
作DFT，Rxx（m）能推出SBT（k）；
（3）加窗平滑法（BT法）
加窗平滑法的原理是先做自相关估计，在选择合适的窗函数相乘，也即截断，然后作DFT。
（4）平均周期图法（Bartlett法）
平均周期图法的原理是个独立同分布的随机变量的均值之方差，等于单个变量方差的1/k。
具体的方法步骤是长数据N分成k段，每段M=N/k，针对每段分别用周期图法求谱，然后k段平均后求的的就是谱估计。
ww
（5）韦尔奇（Welch）谱估计法
韦尔奇谱估计法将加窗平滑法和平均周期法二者相结合，其原理在于：
a）把N个数据分成k段，每段可以互相独立（如平均周期图法）
再把每段数据乘上窗函数w（n）（如加窗平滑法）后作DFT。
第二部分仿真结果图
未加窗时周期图：
取平均后的图形为：
加窗周期图法：
加矩形窗：
加汉明窗后的图形为：
下图为采样频率为100，窗函数分别取矩形窗、汉宁窗和汉明窗时的图形：
未加窗时自相关法：
B-T法：
加汉明窗时：
welch法：
L=512,加汉明窗，50%重叠
L=512,加矩形窗，50%重叠
L=256,加汉明窗，50%重叠
L=256,加矩形窗，50%重叠
L=128,加汉明窗，50%重叠
Bartlett法：
L=512，加汉明窗
L=256，加汉明窗
L=128,加汉明窗
第三章结果分析
周期图法所求得的功率谱振荡剧烈，信号方差较大，不利于对功率信号的分析。
缺少了统计平均时，记录的信号序列长度一定的条件下，要保证足够高的谱分辨率，谱估计的方差就会很大，谱的正确性会很差。
当数据长度N太大时，谱曲线呈现较大的起伏；
当数据长度N太小时，谱的分辨率又不好。
先根据实验所给数据求解出自相关函数，然后对自相关函数进行傅里叶变换，从而得到功率谱估计。
从图中可以看出，自相关法得到的结果与周期图法相似，同样是方差较大，信号振荡大。
（3）加窗平滑法（BT）
根据加窗平滑法的求解步骤进行编程，取窗函数为矩形窗。
但是可以看出在L取较大时性能不是很好。
（4）平均周期图法（Bartlett）
将数据平均分为K段，Bartlett法很好地改善了直接法的方差特性，但是它是以牺牲偏差和分辨率为代价的。
（5）Welch法
Welch法各段允许交叠，从而增大了段数L，这样可以更好地改善方差特性。
但是，数据的交叠又减小了每一段的不相关性，使方差的减小不会达到理论计算的程度。
另外，选择合适的窗函数可以减小频谱的泄漏，改善分辨率。
另外，由矩形窗处理的谱估计的主瓣宽度最窄，分辨率最好，但是其旁瓣比其他窗函数的旁瓣要高，因此其正弦谱线附近的旁瓣泄漏比较严重，而且其起伏性较大，所以其方差特性最差。
由海明窗处理谱估计的主瓣宽度最宽，因此其分辨率相对较差，但其旁瓣较小，大大改善了由矩形窗处理的谱估计旁瓣较大所产生的谱失真。
究其原因，选择不同的窗函数其主瓣宽度不一样，造成谱估计的分辨率也不相同；
另外，选择不同的窗函数旁瓣的衰减速度也不相同，因而谱估计旁瓣的泄漏程度也不一样。
附录：
clear;
Fs=1000;
n=0:
1/Fs:
1;
nfft=1024;
x=cos（100*pi*n）+cos（200*pi*n）+0.25*cos（400*pi*n）+randn（size（n））;
y=fft（x）;
Rx=abs（y）.^2/length（n）;
a=length（Rx）;
sum=zeros（1,a）;
fori=0:
49
x=cos（100*pi*n）+cos（200*pi*n）+0.25*cos（400*pi*n）+randn（size（n））;
y=fft（x）;
Rx=abs（y）.^2/length（n）;
sum=sum+Rx;
end
xpsd=10*log10（sum）;
plot（xpsd）
axis（[0,280,15,45]）;
title（'
周期图法求功率谱密度函数'
）;
xlabel（'
frequency（Hz）'
ylabel（'
power/frequency（dB/Hz）'
a=length（x）;
window=rectwin（a）;
Rx=periodogram（x,window,nfft,Fs）;
t=length（Rx）;
sum=zeros（t,1）;
Rx=periodogram（x,window,nfft,Fs）;
axis（[0,250,-15,15]）;
周期图法求功率谱密度'
rx=xcorr（x,'
biased'
y=fft（rx）;
t=length（y）;
sum=zeros（1,t）;
rx=xcorr（x,'
y=fft（rx）;
Rx=abs（y）;
axis（[0,500,10,45]）
自相关法求功率谱密度'
B-T法:
L=256;
window=hamming（L）;
rxx=xcorr（x,'
rxxx=rxx（1:
L）;
rx=rxxx*window;
rxx（1:
L）=rx;
y=fft（rxx,nfft）;
rxx=xcorr（x,'
rxxx=rxx（1:
rx=rxxx*window;
rxx（1:
y=fft（rxx,nfft）;
plot（xpsd）;
B-T法求功率谱密度函数'
）
axis（[0,256,15,40]）;
noverlap=L/2;
Rx=pwelch（x,window,noverlap,nfft,Fs）;
Rx=pwelch（x,window,noverlap,nfft,Fs）;
axis（[0,256,-15,10]）;
welch法法求功率谱密度函数'
noverlap=0时的welch法
升级会员 