matlab程序中功率谱分析的经典常用方法.doc

一、直接法 clear;clc;close all; %清除变量；清屏；关闭当前图形窗口 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,1 120); 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,1 120); xlabel(Frequency/Hz); 三、Bartlett法 clear;clc;close all; %清除变量；清屏；关闭当前图形窗口 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;close all;%清除变量；清屏；关闭当前图形窗口 Fs=2000; %采样频率 load Chanel8Xia2Data.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);grid on; title(直接法 ); set(gca,xlim,1 100); ylabel(Am/dB); hold on; plot(f,plot_Pxx2,r); axis(0,100,-50,40);grid on; title(直接法 ); xlabel(Frequency/Hz);ylabel(Am/dB); 二、Welch法 clear;clc;close all; %清除变量；清屏；关闭当前图形窗口 Fs=2000; load Chanel8Xia2Data.mat; x1=Chanel8Xia2Gray; x2=Chanel8Xia2Sti; window1=hamming(1024); %海明窗 noverlap=256; %数据无重叠 range=onesided; %频率间隔为0 Fs/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);grid on;