基于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点点，得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段平均后求的的就是谱估计。(5)韦尔奇（Welch）谱估计法韦尔奇谱估计法将加窗平滑法和平均周期法二者相结合，其原理在于：a)把N个数据分成k段，每段可以互相独立(如平均周期图法)b)再把每段数据乘上窗函数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: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;endxpsd=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:49x=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;endxpsd=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:49x=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;endxpsd=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:49x=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;endxpsd=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: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;endxpsd=10*log10(sum);plot(xpsd);axis([0,256,-15,10]);title('welch法法求功率谱密度函数');xlabel('frequency(Hz)');ylabel('power/frequency(dB/Hz)')Bartlett法：noverlap=0时的welch法