经典功率谱设计.docx
- 文档编号:8891501
- 上传时间:2023-05-15
- 格式:DOCX
- 页数:10
- 大小:216.94KB
经典功率谱设计.docx
《经典功率谱设计.docx》由会员分享,可在线阅读,更多相关《经典功率谱设计.docx(10页珍藏版)》请在冰点文库上搜索。
经典功率谱设计
DSP实验报告
实验题目:
实验功率谱估计
一、实验要求:
(1)理解功率谱估计的基本概念;
(2)掌握经典功率谱估计方法——直接法和间接法;
(3)掌握改进的经典功率谱估计方法,例如Welch法。
二、实验内容与原理:
功率谱估计就是基于有限的数据寻找信号、随机过程或系统的频率成分。
它表示随机信号频率域的统计特性。
随机信号是无始无终具有无限能量的,所以其傅立叶变换并不存在,因为它不满足绝对可积的条件。
因此需要研究其在频率域上的功率分布情况,即功率谱密度或功率谱。
根据实验要求,完成该实验首先要正确的生成被估计信号
。
数据的长度和FFT所用的数据长度都设为1024。
1.周期图法:
直接法,即周期图法,是由傅立叶变换得到的:
将随机信号
的N点样本值
看作能量有限信号,取其傅立叶变换,得到
;然后再取其幅值的平方,并除以N作为
的真实功率谱
的估计,即
实验中,将随机信号x(n)的N点样本值看作xN(n)能量有限信号,取其傅立叶变换,得到X,然后再取其幅值的平方,并除以N作为x(n)的真实功率谱P的估计。
2.间接法:
间接法,又称为自相关法或BT法,是由随机信号N个观察值
,估计出自相关函数
,然后再求
的傅立叶变换作为功率谱的估计:
即如下计算:
实验中由随机信号N个观察值估计出自相关函数R(m),然后再求R(m)傅立叶变换作为功率谱的估计PBT。
直接法和间接法的方差性能很差,而且当数据长度太大时,谱曲线起伏加剧;若数据长度太小,则谱的分辨率又不好,所以需要改进。
改进的直接谱估计方法由Bartlett法和Welch法。
3.BARTLETT算法
Bartlett法将采样数据
分成L段,每段的长度都是M,即N=LM,对每段数据加矩形窗,再计算其各自的功率谱
,把
对应相加,再取平均,得到平均周期图
。
即如下过程:
首先将观测数据分为L段,每段长M,分段后的数据可以用式
(1)表示:
;
;
(1)
其中
平均的周期图为:
BARTLETT法是周期图算法的一种改进。
由概率论的知识知道,如果
是N个不相关的随机变量,每个随机变量的期望值为
,方差为
,那么将这N个随机变量求平均,它的期望仍为
,方差变为
。
BARTLETT法即是受此启发,将观测数据分段,先求每段数据的周期图,再求平均的周期图,当分段较多时,估计出的功率谱较平滑,频率分辨率较差;当分段较少时,估计出的功率谱起伏较大,频率分辨率较好。
4.WELCH算法
Welch法是对Bartlett法的改进:
一,在对
分段时,可允许每段数据有部分重叠;二,每段数据窗口可以不是矩形窗口,例如使用汉宁窗或哈明窗,记为
。
然后按Bartlett法求每一段的功率谱,记为
=
,其中
。
平均后的功率谱为:
(段数
)
加窗的优点是使得无论对于什么样的窗函数均可以谱估计为非负值;二是在分段时,各段之间有重叠,这样会使方差减小。
三、实验目的:
分别用四种不同的发放进行功率谱估计,并对比结果。
四、实验程序:
clear;
%数据的长度和FFT所用的数据长度
nfft=1024;N=1024;
%每段长度
Ns=256;
%产生含有噪声的序列xn
n=[0:
N-1];
w1=2*pi*0.02;
w2=2*pi*0.28;
wn=randn(1,N);
xn=sin(w1*n)+2*cos(w2*n)+wn;
%直接法求功率谱
%将随机信号x(n)的N点样本值看作xN(n)能量有限信号,取其傅立叶变换,得到X;然后再取其幅值的平方,并除以N作为x(n)
%的真实功率谱P的估计。
%计算序列的DFT
XN=fft(xn,nfft);
%对序列取绝对值后平方
PER=abs(XN).^2/N;
%并转化为dB
PERdb=10*log10(PER);
%给出频率序列
f=(0:
length(PERdb)-1)/length(PERdb);
%绘制功率谱图形
figure
(1);
plot(f,PERdb);
xlabel('频率/Hz');
ylabel('功率谱/dB');
title('直接法N=1024');
grid;
%间接法求功率谱
%又称为自相关法或BT法,是由随机信号N个观察值估计出自相关函数R(m),然后再求R(m)傅立叶变换作为功率谱的估计PBT
%计算序列的自相关函数Rm
Rm=xcorr(xn,'unbiased');
%计算自相关函数Rm的DTFT
PBT=fft(Rm,nfft);
%把PBT转化为dB
PBTdb=10*log10(abs(PBT));
Fbt=(0:
length(PBTdb)-1)/length(PBTdb);
%绘制功率谱图形
figure
(2);
plot(Fbt,PBTdb);
xlabel('频率/Hz');
ylabel('功率谱/dB');
title('间接法N=1024');
grid;
%bartlett法求功率谱
%Bartlett平均周期图的方法是将N点的有限长序列x(n)分段,对各段用周期图法求解功率后再平均。
%加矩形窗
window=ones(1,Ns);
normvalue=norm(window);
PBAR1=abs(fft(window.*xn(1:
256),Ns).^2)/normvalue^2;%第一段功率谱
PBAR2=abs(fft(window.*xn(257:
512),Ns).^2)/normvalue^2;%第二段功率谱
PBAR3=abs(fft(window.*xn(513:
768),Ns).^2)/normvalue^2;%第三段功率谱
PBAR4=abs(fft(window.*xn(769:
1024),Ns).^2)/normvalue^2;%第四段功率谱
%求Fourier振幅谱的平均值,并转化为dB
PBAR=10*log10((PBAR1+PBAR2+PBAR3+PBAR4)/4);
%给出频率序列
Fpbar=(0:
length(PBAR)-1)/length(PBAR);
%绘制功率谱曲线
figure(3);
plot(Fpbar,PBAR);
xlabel('频率/Hz');ylabel('功率谱/dB');
title('bartlett法4*256');
grid;
%welch法求功率谱
%Welch法对Bartlett法进行了两方面的修正,一是选择适当的窗函数w(n),并再周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。
%二是在分段时,可使各段之间有重叠,这样会使方差减小。
%加汉宁窗
hanningwindow=hanning(Ns)';
normvalue=norm(window);
PWel1=abs(fft(hanningwindow.*xn(1:
256),Ns).^2)/normvalue^2;%第一段功率谱
PWel2=abs(fft(hanningwindow.*xn(129:
384),Ns).^2)/normvalue^2;%第二段功率谱
PWel3=abs(fft(hanningwindow.*xn(257:
512),Ns).^2)/normvalue^2;%第三段功率谱
PWel4=abs(fft(hanningwindow.*xn(385:
640),Ns).^2)/normvalue^2;%第四段功率谱
PWel5=abs(fft(hanningwindow.*xn(513:
768),Ns).^2)/normvalue^2;%第五段功率谱
PWel6=abs(fft(hanningwindow.*xn(641:
896),Ns).^2)/normvalue^2;%第六段功率谱
PWel7=abs(fft(hanningwindow.*xn(769:
1024),Ns).^2)/normvalue^2;%第七段功率谱
%求Fourier振幅谱的平均值,并转化为dB
PWel=10*log10((PWel1+PWel2+PWel3+PWel4+PWel5+PWel6+PWel7)/7);
%给出频率序列
FWel=(0:
length(PWel)-1)/length(PWel);
%绘制功率谱曲线
figure(4);
plot(FWel,PWel);
xlabel('频率/Hz');ylabel('功率谱/dB');
title('welch法7*256');
grid;
五、实验结果:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 经典 功率 设计