function TF_grayimage(file_name,var_name)
% TF_grayimage() 繪制出指定時域信號的時頻灰度圖；
% 輸入?yún)?shù)為保存時域信號的mat文件名,eg. TF_grayimage('mydata.mat')；
% 處理過程如下：
% 給指定信號加窗；
% 再進(jìn)行短時傅立葉變換；
% 繪制時頻灰度圖；
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；

%% 設(shè)計窗函數(shù)；
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;            %獲取時域信號全部加窗段數(shù)目 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