《基于matlab的功率谱估计》由会员分享，可在线阅读，更多相关《基于matlab的功率谱估计（24页珍藏版）》请在金锄头文库上搜索。
1、2013 毕业设计论文2013 毕业设计论文题题 目目 基于 matlab 的功率谱估计 学学 院院 通信工程学院 专专 业业 通信与信息系统 2013 毕业设计论文2013 毕业设计论文第一章第一章 功率谱估计分析及比较功率谱估计分析及比较1. 实验目的实验目的(1) 掌握 Welch 算法的概念、应用及特点；(2) 了解谱估计在信号分析中的作用；(3) 能够利用 Welch 法对信号作谱估计，对信号的特点加以分析。2. 实验内容实验内容(1) 读入实验数据。(2) 编写一利用 Welch 法作估计的算法程序。(3) 将计算结果表示成图形的形式，给出信号谱的分布情况图。3. 谱估计方法简介谱
2、估计方法简介(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 法） 平均周期图法的原理是个独立
3、同分布的随机变量的均值之方差，等于单个变量方差的 1/ k。 具体的方法步骤是长数据 N 分成 k 段，每段 M =N/k，针对每段分别用周期图法求谱，然后 k 段平均后求的的就是谱估计。 (5) 韦尔奇（Welch）谱估计法韦尔奇谱估计法将加窗平滑法和平均周期法二者相结合，其原理在于：2013 毕业设计论文2013 毕业设计论文a) 把 N 个数据分成 k 段，每段可以互相独立(如平均周期图法)b) 再把每段数据乘上窗函数 w(n) (如加窗平滑法)后作 DFT。2013 毕业设计论文2013 毕业设计论文第二部分第二部分 仿真结果图仿真结果图未加窗时周期图：未加窗时周期图：取平均后的图形为
4、：2013 毕业设计论文2013 毕业设计论文加窗周期图法：加窗周期图法：加矩形窗：取平均后的图形为：2013 毕业设计论文2013 毕业设计论文加汉明窗后的图形为：取平均后的图形为：2013 毕业设计论文2013 毕业设计论文下图为采样频率为 100，窗函数分别取矩形窗、汉宁窗和汉明窗时的图形：2013 毕业设计论文2013 毕业设计论文未加窗时自相关法：未加窗时自相关法：取平均后的图形为：2013 毕业设计论文2013 毕业设计论文B-T 法：法：加汉明窗时：取平均后的图形为：2013 毕业设计论文2013 毕业设计论文welch 法：法：L=512,加汉明窗，50%重叠取平均后的图形为：
5、2013 毕业设计论文2013 毕业设计论文L=512,加矩形窗，50%重叠取平均后的图形为：2013 毕业设计论文2013 毕业设计论文L=256,加汉明窗，50%重叠取平均后的图形为：2013 毕业设计论文2013 毕业设计论文L=256,加矩形窗，50%重叠取平均后的图形为：2013 毕业设计论文2013 毕业设计论文L=128,加汉明窗，50%重叠取平均后的图形为：2013 毕业设计论文2013 毕业设计论文Bartlett 法：法：L=512，加汉明窗取平均后的图形为：2013 毕业设计论文2013 毕业设计论文L=256，加汉明窗取平均后的图形为：2013 毕业设计论文2013 毕
6、业设计论文L=128,加汉明窗取平均后的图形为：2013 毕业设计论文2013 毕业设计论文2013 毕业设计论文2013 毕业设计论文第三章第三章 结果分析结果分析（1）周期图法 周期图法所求得的功率谱振荡剧烈，信号方差较大，不利于对功率信号的分析。缺少了统计平均时，记录的信号序列长度一定的条件下，要保证足够高的谱分辨率，谱估计的方差就会很大，谱的正确性会很差。当数据长度 N 太大时，谱曲线呈现较大的起伏；当数据长度 N 太小时，谱的分辨率又不好。(2) 自相关法 先根据实验所给数据求解出自相关函数，然后对自相关函数进行傅里叶变换， 从而得到功率谱估计。从图中可以看出，自相关法得到的结果与周
7、期图法相似，同样是方差较大，信号振荡大。(3) 加窗平滑法（BT） 根据加窗平滑法的求解步骤进行编程，取窗函数为矩形窗。但是可以看出在 L 取较大时性能不是很好。(4) 平均周期图法（Bartlett） 将数据平均分为 K 段，Bartlet t 法很好地改善了直接法的方差特性，但是它是以牺牲偏差和分辨率为代价的。(5) Welch 法 Welch 法各段允许交叠，从而增大了段数 L，这样可以更好地改善方差特性。但是，数据的交叠又减小了每一段的不相关性，使方差的减小不会达到理论计算的程度。另外，选择合适的窗函数可以减小频谱的泄漏，改善分辨率。另外，由矩形窗处理的谱估计的主瓣宽度最窄，分辨率最好
8、，但是其旁瓣比其他窗函数的旁瓣要高，因此其正弦谱线附近的旁瓣泄漏比较严重，而且其起伏性较大，所以其方差特性最差。由海明窗处理谱估计的主瓣宽度最宽，因此其分辨率相对较差，但其旁瓣较小，大大改善了由矩形窗处理的谱估计旁瓣较大所产生的谱失真。究其原因，选择不同的窗函数其主瓣宽度不一样，造成谱估计的分辨率也不相同；另外，选择不同的窗函数旁瓣的衰减速度也不相同，因而谱估计旁瓣的泄漏程度也不一样。2013 毕业设计论文2013 毕业设计论文附录：附录：未加窗时周期图：未加窗时周期图：clear;Fs=1000;n=0:1/Fs:1;nfft=1024;x=cos(100*pi*n)+cos(200*pi*
9、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);for i=0:49x=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;ensd=10*log10(sum);plot(xpsd)axis(0,280,15,45);title(周期图法求功率谱密度函数);xlabel(frequency(Hz);yl
10、abel(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);2013 毕业设计论文2013 毕业设计论文a=length(x);window=rectwin(a);Rx=periodogram(x,window,nfft,Fs);t=length(Rx);sum=zeros(t,1);for i=0:49x=cos(100*pi*n)+cos(200*pi*n)+0.25*
11、cos(400*pi*n)+randn(size(n);Rx=periodogram(x,window,nfft,Fs);sum=sum+Rx;ensd=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(si
12、ze(n);rx=xcorr(x,biased);y=fft(rx);t=length(y);sum=zeros(1,t);for i=0:49x=cos(100*pi*n)+cos(200*pi*n)+0.25*cos(400*pi*n)+randn(size(n);rx=xcorr(x,biased);2013 毕业设计论文2013 毕业设计论文y=fft(rx);Rx=abs(y);sum=sum+Rx;ensd=10*log10(sum);plot(xpsd)axis(0,500,10,45)title(自相关法求功率谱密度);xlabel(frequency(Hz);ylabel(p
13、ower/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);for i=0:49x=cos(100*pi*n)+cos(200*pi*n)+0.2
14、5*cos(400*pi*n)+randn(size(n);rxx=xcorr(x,biased);rxxx=rxx(1:L);2013 毕业设计论文2013 毕业设计论文rx=rxxx* window;rxx(1:L)=rx;y=fft(rxx,nfft);Rx=abs(y);sum=sum+Rx;ensd=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/
15、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);for i=0:49x=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;end2013 毕业设计论文2013 毕业设计论文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法
