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

一、直接法

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;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';%频率间隔为[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;

title('Welch法海明窗');

ylabel('Am/dB');

set(gca,'xlim',[1 100]);

hold on;

plot(f,plot_Pxx2,'r');

title('Welch法海明窗');

ylabel('Am/dB');

xlabel('Frequency/Hz');