短时傅里叶变换（Short-time Fourier Transform, STFT）Matlab代码

function [S, f, t] = stft(x, fs, window, nfft, noverlap)
% STFT - Short-time Fourier Transform
% [S, f, t] = stft(x, fs, window, nfft, noverlap)
% x: input signal
% fs: sampling frequency
% window: window function (default: hamming)
% nfft: number of FFT points (default: length of window)
% noverlap: number of samples overlapped between adjacent frames (default: 0)
% S: STFT matrix (nfft/2+1 x nframes)
% f: frequency vector
% t: time vector
% Example:
% [x, fs] = audioread('speech.wav');
% [S, f, t] = stft(x, fs, hamming(512), 512, 256);
% imagesc(t, f, 20*log10(abs(S)));
% axis xy; colormap(jet); colorbar; xlabel('Time (s)'); ylabel('Frequency (Hz)');
% Written by Yuancheng Zhu (yzhu@nd.edu)
% Last update: 2021/9/9

% Check inputs
narginchk(2, 5);
if nargin < 3 || isempty(window)
window = hamming(256);
end
if nargin < 4 || isempty(nfft)
nfft = length(window);
end
if nargin < 5 || isempty(noverlap)
noverlap = 0;
end
if ~isvector(window) || ~isnumeric(window)
error('Window function must be a numeric vector.');
end
if ~isscalar(fs) || ~isnumeric(fs) || fs <= 0
error('Sampling frequency must be a positive scalar.');
end
if ~isscalar(nfft) || ~isnumeric(nfft) || nfft <= 0
error('Number of FFT points must be a positive scalar.');
end
if ~isscalar(noverlap) || ~isnumeric(noverlap) || noverlap < 0
error('Number of overlapped samples must be a non-negative scalar.');
end
if noverlap >= length(window)
error('Number of overlapped samples must be less than the window length.');
end

% Calculate STFT parameters
nframes = fix((length(x)-noverlap)/(length(window)-noverlap));
if nframes < 1
error('Signal is too short for STFT with the given parameters.');
end
x = x(:);
window = window(:);
S = zeros(nfft/2+1, nframes);
t = (nfft/2:nfft/2+nframes-1) / fs;

% Compute STFT
for i = 1:nframes
idx = (1:length(window)) + (i-1)*(length(window)-noverlap);
if idx(end) > length(x)
idx = idx(idx <= length(x));
xw = [x(idx); zeros(length(window)-length(idx), 1)];
else
xw = x(idx) .* window;
end
X = fft(xw, nfft);
S(:, i) = X(1:nfft/2+1);
end
f = (0:nfft/2)' / nfft * fs;

end