经典功率谱分析Matlab程序.docx

一、直接法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,1 120); ;ylabel(Am/dB);xlabel(Frequency/Hz);四、Welch法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=128;range=onesided; %频率间隔为0 Fs/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,1 120); 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