%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Short-Time Fourier Transform %
% with MATLAB Implementation %
% %
% Author: M.Sc. Eng. Hristo Zhivomirov 12/21/13 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [stft, f, t] = stft(x, wlen, h, nfft, fs)
% function: [stft, t, f] = mystftfun(x, fs, wlen, h, nfft)
% x - signal in the time domain
% wlen - length of the hamming window
% h - hop size
% nfft - number of FFT points
% fs - sampling frequency, Hz
% f - frequency vector, Hz
% t - time vector, s
% stft - STFT matrix (time across columns, freq across rows)
% represent x as column-vector if it is not
if size(x,2) > 1
x = x';
end
% length of the signal
xlen = length(x);
% form a periodic hamming window
%win = hamming(wlen, 'periodic');
win = hamming(nfft);
% form the stft matrix
rown = ceil((1+nfft)/2); % calculate the total number of rows
coln = 1+fix((xlen-wlen)/h); % calculate the total number of columns
stft = zeros(rown, coln); % form the stft matrix
% initialize the indexes
indx = 0;
col = 1;
% perform STFT
while indx + wlen <= xlen
% windowing
xw = x(indx+1:indx+wlen).*win;
% FFT
X = fft(xw, nfft);
% update the stft matrix
stft(:,col) = X(1:(rown));
% update the indexes
indx = indx + h;
col = col + 1;
end
% calculate the time and frequency vectors
t = (wlen/2:h:xlen-wlen/2-1)/fs;
f = (0:rown-1)*fs/nfft;
end