%產生信號
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');