stft.m

对信号进行短时傅立叶变换以及Wigner-Ville分布分析的matlab程序

%产生信号
s=data;
%采集点数
N=8192;
n=0:N-1;
%采样频率
Fs=15360;
t=n/Fs;
F=n*Fs/N;
%滤波器阶数;
n=10;
%低通滤波器的截止频率设为500Hz(直接给出截止频率即可);
Wn=500/(Fs/2);
%计算滤波器系数;
[a,b]=butter(n,Wn);
%对信号进行低通滤波;
a=filter(a,b,s);

%时域波形
figure(1);plot(t,a);title('时域图');
xlabel('时间t/s');
ylabel('幅值A/m/s^2');
%时域指标
X=max(a);
y=var(a);
u=mean(a);
s=std(a);
Xrms=(1/N*sum(a.^2))^(1/2);
I=X/abs(u);
Kv=kurtosis(a);
CLf=skewness(a);
C=max(a)/(N*(u^2));
K=Xrms/abs(u);
%幅值图
b=fft(a);
Fp=abs(b);
figure(2);plot(F(1:N/2),Fp(1:N/2));title('幅值图');
xlabel('频率Hz');
ylabel('振幅');
%功率谱
Sn=1/N*abs(Fp.^2);
figure(3);plot(F(1:N/2),Sn(1:N/2));title('功率谱');
%设置窗函数1
h1=window(@guasswin,52);
%计算短时傅立叶变化1
a1=a';
sig=hilbert(a1);
[tfr1,t,f]=tfrstft(sig,1:1024,1024,h1);
figure(4);
contour(t,f,abs(tfr1));title('短时傅立叶加guasswin窗函数');
xlabel('时间t');
ylabel('频率f');
%设置窗函数2
h2=window(@hamming,65);
%计算短时傅立叶变化2
a1=a';
sig=hilbert(a1);
[tfr2,t,f]=tfrstft(sig,1:1024,1024,h2);
figure(5);
contour(t,f,abs(tfr2));title('短时傅立叶加hamming窗函数');
xlabel('时间t');
ylabel('频率f');
%设置窗函数3
h3=window(@hann,65);
%计算短时傅立叶变化3
a1=a';
sig=hilbert(a1);
[tfr3,t,f]=tfrstft(sig,1:1024,1024,h3);
figure(6);
contour(t,f,abs(tfr3));title('短时傅立叶加hann窗函数');
xlabel('时间t');
ylabel('频率f');
%设置窗函数4
h4=window(@flattopwin,65);
%计算短时傅立叶变化4
a1=a';
sig=hilbert(a1);
[tfr4,t,f]=tfrstft(sig,1:1024,1024,h4);
figure(7);
contour(t,f,abs(tfr4));title('短时傅立叶加flattopwin窗函数');
xlabel('时间t');
ylabel('频率f');
%设置窗函数5
h5=window(@rectwin,65);
%计算短时傅立叶变化5
a1=a';
sig=hilbert(a1);
[tfr5,t,f]=tfrstft(sig,1:1024,1024,h5);
figure(8);
contour(t,f,abs(tfr5));title('短时傅立叶加rectwin窗函数');
xlabel('时间t');
ylabel('频率f');
%计算Wigner-Ville分布1
a1=a';
sig=hilbert(a1);
[tfr5,t,f]=tfrwv(sig,1:256,256);
figure(9);
contour(t,f,tfr5);title('Wigner-Ville分布(等高线)');
xlabel('时间t');
ylabel('频率f');
figure(10);
mesh(t,f,tfr5);title('Wigner-Ville分布(三维图)');
xlabel('时间t');
ylabel('频率f');
zlabel('幅值A');
%计算伪Wigner-Ville分布
a1=a';
sig=hilbert(a1);
[tfr6,t,f]=tfrpwv(sig,1:256,256);
figure(11);
contour(t,f,abs(tfr6));title('伪Wigner-Ville分布(等高线)');
xlabel('时间t');
ylabel('频率f');
%计算平滑伪Wigner-Ville分布
a1=a';
sig=hilbert(a1);
[tfr8,t,f]=tfrspwv(sig,1:256,256);
figure(12);
contour(t,f,abs(tfr8));title('平滑伪Wigner-Ville分布(等高线)');
xlabel('时间t');
ylabel('频率f');