【20211206】【信号处理】时频分析 —— 短时傅里叶变换（STFT）

一、时频分析（JTFA）

        分析时域可以得到信号随时间变化的信息，分析频域可以得到信号随频率变化的信息，这两者都只能分析时域或频域，而不能同时观察时域和频域。

        时频分析是时频联合域分析的简称，是分析时变非平稳信号的有力工具，是一个同时观察时域、频域信息的工具。时频分析法提供了时间域和频率域的联合分布信息，清楚地描述了信号频率随时间变化的关系。

        时频分析的基本思想是：设计时间和频率的联合函数，用它同时描述信号在不同时间和频率的能量密度或强度。

        时间和频率的这种联合函数简称为时频分布。利用时频分布来分析信号，能给出各个时刻的瞬时频率和幅值。

        常见的时频分布有：短时傅里叶变换、小波变换、希尔伯特黄变换等。

        可以这么理解时频分析：图（a）是时域的一个时序信号，可以看出在 0~1s 信号频率低，1s 之后的信号频率增大；图（b）是该时序信号的频谱，能够看出该信号有两个主要频率，分别是 10Hz 和 20Hz；图（c）是该信号的时频图，横坐标是时间、纵坐标是频率，在 0~1s 的频率为 10Hz，1s 之后信号频率是 20Hz。

        （参考：时频分析）

        （参考：短时傅里叶变换）

二、短时傅里叶变换（STFT）

        短时傅里叶变换是线性时频分析中的一种。

        实现方式是：在时域用窗函数截断，对窗信号做傅里叶变换，即得到该时刻的傅里叶变换，不断地移动窗口的中心位置，即可得到不同时刻的傅里叶变换，这些傅里叶变换的集合就是短时傅里叶变换，得到的图形叫时频图，横坐标是时间，纵坐标是频率，z 轴是信号能量/功率。

        短时傅里叶变换选用短窗口，会有较高的时间分辨率，但频率分辨率较差；选用长窗口，会有较高的频率分辨率，但时间分辨率就差。所以，在应用中要注意时窗和频窗宽度的折衷。

        处理步骤：

        （1）选择一个窗函数；

        （2）截取窗长度的信号，并和窗函数相乘；

        （3）进行 FFT；

        （4）滑窗得到每个时刻的 STFT。

        简单来说，短时傅里叶变换就是通过窗内的一段信号来表示某一时刻的信号特征，傅里叶变换就是先把一个函数和窗函数进行相乘，然后再进行一维的傅里叶变换，并通过窗函数的滑动得到一系列的频谱函数，所以最终得到一个二维的时频图。

       （时频联合分析参考：时频联合分析）

        （参考：信号处理 | 傅里叶变换、短时傅里叶变换、小波变换、希尔伯特变换、希尔伯特黄变换）

        （参考：时频分析：短时傅里叶变换（1））

三、举个栗子

%% 参数设置
fs = 50;  % 采样频率
L = 10;  % 时间长度
t = 0 : 1/fs : L;  % 时间坐标
%% 生成信号
fsignal = 10;  % 信号频率
A = 10;  % 信号幅值
datalength = length(t);
s = A * sin(2 * pi * fsignal * t) + randn(1, datalength);
numfft = 256;   % fft 点数
f = (0 : numfft/2-1) / numfft * fs;  % 频率（STFT 纵轴坐标） 
%% 滑窗进行 stft
winLen = 20;
stepLen = 2;
winFun = chebwin(winLen);
win_ind_packet = 1 : winLen;
s_STFT_dB = zeros(numfft/2, datalength);
cnt = 1;
while win_ind_packet(end) < datalength
    subData = s(win_ind_packet)';
    subDataWin = subData .* winFun;
    subDataWinFFT = fft(subDataWin, numfft, 1);
    subDataWinFFTPow = abs(subDataWinFFT(1:numfft/2)) .^ 2;
    subDataWinFFTPowdB = pow2db(subDataWinFFTPow);   
    s_STFT_dB(:, cnt) = subDataWinFFTPowdB;
    
    win_ind_packet = win_ind_packet + stepLen;
    cnt = cnt + 1;
end
s_STFT_dB(:, cnt+1:datalength) = nan(numfft/2, datalength-cnt);  % 没有进行傅里叶变换的窗数据补nan
 
%% 作图
packet_ind = 1 : datalength;
figure(2);
surf(packet_ind, f, s_STFT_dB, 'edgecolor', 'none'); colormap('jet'); axis('tight');
xlabel('时间'); ylabel('频率（Hz）'); title('信号s的时频图');