tf_grayimage.m

小波方法应用到故障诊断

function TF_grayimage(file_name,var_name)
% TF_grayimage() 绘制出指定时域信号的时频灰度图；
% 输入参数为保存时域信号的mat文件名,eg. TF_grayimage('mydata.mat')；
% 处理过程如下：
% 给指定信号加窗；
% 再进行短时傅立叶变换；
% 绘制时频灰度图；
echo all on 
fs=10000;                               % 采样频率；

%% 载入时域信号
%    Load time-domain signal in x_t
isexist=exist(file_name,'file');
if isexist~=2
    error('Data-file not existed!');
    return;
end
x_t=load(file_name);                 % here x_t is a struct containing all varables in file_name, 
                                                 % eg. x_t.var1,  x_t.var2……；
eval(['x_t=x_t.' var_name ';']);   % get file_name.var_name in x_t；

%% 设计窗函数；
length_win=102;                        % round(fs/fi), fs--sampling frequency, fi--fault frequency；
window=hamming(length_win)';
%sum_win=sum(window);
%window=2*window/sum_win;

%% 短时傅立叶变换算法；
length_fft=128;
frame=2;
sections=(length(x_t)-length_win)/frame;            %获取时域信号全部加窗段数目 numbers : sections of stft；
for i=1:sections
    win_start=frame*(i-1)+1;
    win_end=frame*(i-1)+length_win;
    stft(i,:)=abs(fft(x_t(win_start:win_end).*window,length_fft)/length_fft);
end
stft=2*stft(1:sections,1:length_fft/2)';
axis_freq=fs*(1:length_fft/2)/length_fft;
axis_time=1000*(0:frame:(sections-1)*frame)/fs;
save STFT_stft stft axis_freq axis_time
%% 灰度图归一化；
[row column]=size(stft);
max_element=max(max(stft));
for i=1:row
    for j=1:column
        stft(i,j)=255-abs(stft(i,j))*255/max_element;
    end
end

%% 绘制灰度图；
figure();
image(axis_time,axis_freq,stft);
colormap(gray);
set(gca,'YDir','Normal');
echo all off