%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inverse Short-Time Fourier Transform %
% with MATLAB Implementation %
% %
% Author: M.Sc. Eng. Hristo Zhivomirov 12/26/13 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, t] = istft(stft, wlen, hop, nfft, fs)
% function: [x, t] = istft(stft, wlen, hop, nfft, fs)
% stft - STFT matrix (only unique points, time across columns, freq across rows)
% wlen - length of the sinthesis Hamming window
% hop - hop size
% nfft - number of FFT points
% fs - sampling frequency, Hz
% x - signal in the time domain
% t - time vector, s
% signal length estimation and preallocation
coln = size(stft, 2);
xlen = wlen + (coln-1)*hop;
x = zeros(1, xlen);
% form a periodic hamming window
win = hamming(wlen, 'periodic');
% initialize the signal time segment index
indx = 0;
% perform ISTFT (via IFFT and Weighted-OLA)
if rem(nfft, 2) % odd nfft excludes Nyquist point
for col = 1:coln
% extract FFT points
X = stft(:, col);
X = [X; conj(X(end:-1:2))];
% IFFT
xprim = real(ifft(X));
xprim = xprim(1:wlen); %?????????????????????????????
% weighted-OLA
x((indx+1):(indx+wlen)) = x((indx+1):(indx+wlen)) + (xprim.*win)';
% update the index
indx = indx + hop;
end
else % even nfft includes Nyquist point
for col = 1:coln
% extract FFT points
X = stft(:, col);
X = [X; conj(X(end-1:-1:2))];
% IFFT
xprim = real(ifft(X));
xprim = xprim(1:wlen);
% weighted-OLA
x((indx+1):(indx+wlen)) = x((indx+1):(indx+wlen)) + (xprim.*win)';
% update the index
indx = indx + hop;
end
end
% scale the signal
W0 = sum(win.^2);
x = x.*hop/W0;
% generate time vector
t = (0:xlen-1)/fs;
end