f8_6.m

matlab时频分析配套源码第8章的程序

%实验信号
N=1024;
t=1:N;
%齿轮裂纹故障信号
fid=fopen('gearch1.dat','r'); 
sig=fread(fid,N,'int16');
fclose(fid);
%归一化
sig=(sig-mean(sig))/std(sig,1);
%采样频率
fs=10000;
%时域波形
figure(1)
plot(t/fs,sig);
xlabel('时间 t/s');
ylabel('幅值 A/V');
%功率谱
nfft=10240;
y=abs(psd(sig,nfft));
figure(2)
plot((0:nfft/2-1)/nfft*fs,y(1:nfft/2));
xlabel('频率 f/Hz');
ylabel('功率谱 S/V');
%计算伪Wigner-Ville分布
sig1=hilbert(sig);
[tfr,t,f]=tfrpwv(sig1);
figure(3)
contour(t/fs,f(1:N)*fs,abs(tfr));
xlabel('时间 t/s');
ylabel('频率 f/Hz');

%带通滤波
wp=2*[1500 2400]/fs;
ws=2*[1400 2500]/fs;
rp=3;
rs=30;
[n wn]=buttord(wp,ws,rp,rs);
[b,a]=butter(n,wn);
sig2=filter(b,a,sig);
%计算伪Wigner-Ville分布
sig=hilbert(sig2);
[tfr,t,f]=tfrpwv(sig2);
figure(4)
contour(t/fs,f(1:N)*fs,abs(tfr));
xlabel('时间 t/s');
ylabel('频率 f/Hz');