基于matlab的功率谱估计.docx
- 文档编号:7386213
- 上传时间:2023-01-23
- 格式:DOCX
- 页数:23
- 大小:696.66KB
基于matlab的功率谱估计.docx
《基于matlab的功率谱估计.docx》由会员分享,可在线阅读,更多相关《基于matlab的功率谱估计.docx(23页珍藏版)》请在冰豆网上搜索。
基于matlab的功率谱估计
题目基于matlab的功率谱估计
学院通信工程学院
专业通信与信息系统
第一章功率谱估计分析及比较
1.实验目的
(1)掌握Welch算法的概念、应用及特点;
(2)了解谱估计在信号分析中的作用;
(3)能够利用Welch法对信号作谱估计,对信号的特点加以分析。
2.实验内容
(1)读入实验数据。
(2)编写一利用Welch法作估计的算法程序。
(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,加汉明窗
取平均后的图形为:
第三章结果分析
(1)周期图法
周期图法所求得的功率谱振荡剧烈,信号方差较大,不利于对功率信号的分析。
缺少了统计平均时,记录的信号序列长度一定的条件下,要保证足够高的谱分辨率,谱估计的方差就会很大,谱的正确性会很差。
当数据长度N太大时,谱曲线呈现较大的起伏;当数据长度N太小时,谱的分辨率又不好。
(2)自相关法
先根据实验所给数据求解出自相关函数,然后对自相关函数进行傅里叶变换,从而得到功率谱估计。
从图中可以看出,自相关法得到的结果与周期图法相似,同样是方差较大,信号振荡大。
(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)');
加窗周期图法:
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));
a=length(x);
window=rectwin(a);
Rx=periodogram(x,window,nfft,Fs);
t=length(Rx);
sum=zeros(t,1);
fori=0:
49
x=cos(100*pi*n)+cos(200*pi*n)+0.25*cos(400*pi*n)+randn(size(n));
Rx=periodogram(x,window,nfft,Fs);
sum=sum+Rx;
end
xpsd=10*log10(sum);
plot(xpsd)
axis([0,250,-15,15]);
title('周期图法求功率谱密度');
xlabel('frequency(Hz)');
ylabel('power/frequency(dB/Hz)');
未加窗时自相关法:
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));
rx=xcorr(x,'biased');
y=fft(rx);
t=length(y);
sum=zeros(1,t);
fori=0:
49
x=cos(100*pi*n)+cos(200*pi*n)+0.25*cos(400*pi*n)+randn(size(n));
rx=xcorr(x,'biased');
y=fft(rx);
Rx=abs(y);
sum=sum+Rx;
end
xpsd=10*log10(sum);
plot(xpsd)
axis([0,500,10,45])
title('自相关法求功率谱密度');
xlabel('frequency(Hz)');
ylabel('power/frequency(dB/Hz)');
B-T法:
clear;
Fs=1000;
n=0:
1/Fs:
1;
nfft=1024;
L=256;
x=cos(100*pi*n)+cos(200*pi*n)+0.25*cos(400*pi*n)+randn(size(n));
window=hamming(L);
rxx=xcorr(x,'biased');
rxxx=rxx(1:
L);
rx=rxxx*window;
rxx(1:
L)=rx;
y=fft(rxx,nfft);
t=length(y);
sum=zeros(1,t);
fori=0:
49
x=cos(100*pi*n)+cos(200*pi*n)+0.25*cos(400*pi*n)+randn(size(n));
rxx=xcorr(x,'biased');
rxxx=rxx(1:
L);
rx=rxxx*window;
rxx(1:
L)=rx;
y=fft(rxx,nfft);
Rx=abs(y);
sum=sum+Rx;
end
xpsd=10*log10(sum);
plot(xpsd);
title('B-T法求功率谱密度函数');
xlabel('frequency(Hz)');
ylabel('power/frequency(dB/Hz)')
axis([0,256,15,40]);
welch法:
clear;
Fs=1000;
n=0:
1/Fs:
1;
nfft=1024;
L=256;
noverlap=L/2;
window=hamming(L);
x=cos(100*pi*n)+cos(200*pi*n)+0.25*cos(400*pi*n)+randn(size(n));
Rx=pwelch(x,window,noverlap,nfft,Fs);
t=length(Rx);
sum=zeros(t,1);
fori=0:
49
x=cos(100*pi*n)+cos(200*pi*n)+0.25*cos(400*pi*n)+randn(size(n));
Rx=pwelch(x,window,noverlap,nfft,Fs);
sum=sum+Rx;
end
xpsd=10*log10(sum);
plot(xpsd);
axis([0,256,-15,10]);
title('welch法法求功率谱密度函数');
xlabel('frequency(Hz)');
ylabel('power/frequency(dB/Hz)')
Bartlett法:
noverlap=0时的welch法
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 基于 matlab 功率 估计