您好,欢迎访问三七文档
当前位置:首页 > 行业资料 > 其它行业文档 > 经典功率谱分析Matlab程序
一、直接法clear;clc;closeall;%清除变量;清屏;关闭当前图形窗口Fs=1000;t=0:1/Fs:1;nfft=2048;%改变nfft的值可对比不同采样值时的谱估计效果%****************生成信号、噪声**************%x1=cos(2*pi*40*t)+3*cos(2*pi*45*t);%信号x2=randn(size(t));%噪声x3=x1+x2;%信号+噪声[Pxx,f]=periodogram(x3,window,nfft,Fs);%直接法plot(f,10*log10(Pxx));title('直接法nfft=2048');;set(gca,'xlim',[1120]);;ylabel('Am/dB');xlabel('Frequency/Hz');二、间接法Fs=1000;%采样频率n=0:1/Fs:1;%产生含有噪声的序列x1=cos(2*pi*40*n)+3*cos(2*pi*45*n);%信号x2=randn(size(n));%噪声x3=x1+x2;%信号+噪声nfft=1024;cxn=xcorr(x3);%计算序列的自相关函数CXk=fft(cxn);Pxx=abs(CXk);index=0:round(nfft/2-1);f=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));figure(1)plot(f,plot_Pxx);title('间接法nfft=1024');ylabel('Am/dB');set(gca,'xlim',[1120]);xlabel('Frequency/Hz');三、Bartlett法clear;clc;closeall;%清除变量;清屏;关闭当前图形窗口Fs=1000;t=0:1/Fs:1;nfft=1024;%****************生成信号、噪声**************%x1=cos(2*pi*40*t)+3*cos(2*pi*45*t);%信号x2=randn(size(t));%噪声x3=x1+x2;%信号+噪声window=hamming(512);%海明窗noverlap=0;%数据无重叠p=0.9;%置信概率[Pxx,Pxxc]=psd(x3,nfft,Fs,window,noverlap,p);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));plot_Pxxc=10*log10(Pxxc(index+1));figure(1)plot(k,plot_Pxx);title('Bartlett法海明窗');;set(gca,'xlim',[1120]);;ylabel('Am/dB');xlabel('Frequency/Hz');四、Welch法clear;clc;closeall;%清除变量;清屏;关闭当前图形窗口Fs=1000;t=0:1/Fs:1;nfft=1024;%****************生成信号、噪声**************%x1=cos(2*pi*40*t)+3*cos(2*pi*45*t);%信号x2=randn(size(t));%噪声x3=x1+x2;%信号+噪声window=hamming(512);%海明窗noverlap=128;range='onesided';%频率间隔为[0Fs/2],只计算一半的频[Pxx1,f]=pwelch(x3,window,noverlap,nfft,Fs,range);plot_Pxx1=10*log10(Pxx1);figure(1);plot(f,plot_Pxx1);title('Welch法海明窗');ylabel('Am/dB');set(gca,'xlim',[1120]);xlabel('Frequency/Hz');对所给的实测信号进行谱估计,本文采用了周期图法和Welch法。其程序如下所示:一、周期图法clear;clc;closeall;%清除变量;清屏;关闭当前图形窗口Fs=2000;%采样频率loadChanel8Xia2Data.mat;x1=Chanel8Xia2Gray;x2=Chanel8Xia2Sti;Mlag=length(x1);%****************求功率谱密度**************%[Pxx1,f]=periodogram(x1,window,length(x1),Fs);%直接法[Pxx2,f]=periodogram(x2,window,length(x1),Fs);plot_Pxx1=10*log10(Pxx1);plot_Pxx2=10*log10(Pxx2);%****************显示功率谱密度曲线**************%figure(1);plot(f,plot_Pxx1,'b');axis([0,100,-50,40]);gridon;title('直接法');;set(gca,'xlim',[1100]);ylabel('Am/dB');holdon;plot(f,plot_Pxx2,'r');axis([0,100,-50,40]);gridon;title('直接法');xlabel('Frequency/Hz');ylabel('Am/dB');二、Welch法clear;clc;closeall;%清除变量;清屏;关闭当前图形窗口Fs=2000;loadChanel8Xia2Data.mat;x1=Chanel8Xia2Gray;x2=Chanel8Xia2Sti;window1=hamming(1024);%海明窗noverlap=256;%数据无重叠range='onesided';%频率间隔为[0Fs/2],只计算一半的频、%****************求功率谱密度**************%[Pxx1,f]=pwelch(x1,window1,noverlap,length(x1),Fs,range);[Pxx2,f]=pwelch(x2,window1,noverlap,length(x1),Fs,range);plot_Pxx1=10*log10(Pxx1);plot_Pxx2=10*log10(Pxx2);%****************显示功率谱密度曲线**************%figure(1);plot(f,plot_Pxx1,'b');axis([0,100,-30,30]);gridon;title('Welch法海明窗');ylabel('Am/dB');set(gca,'xlim',[1100]);holdon;plot(f,plot_Pxx2,'r');title('Welch法海明窗');ylabel('Am/dB');xlabel('Frequency/Hz');
本文标题:经典功率谱分析Matlab程序
链接地址:https://www.777doc.com/doc-6134819 .html