I'm not sure what this is. I asked for a STFT and ISTFT in Python and this is some collected implementations? I guess I was trying to write a fully-reversible STFT/ISTFT in Python?
First Matlab version from from Time-frequency Automatic Gain Control (AGC) by Dan Ellis Matlab updates from physio-acoustic-simulator's stft.m
pytfd (Time-Frequency distributions in Python) has a stft.py
| function x = istft(d, ftsize, w, h) | |
| % X = istft(D, F, W, H) Inverse short-time Fourier transform. | |
| % Performs overlap-add resynthesis from the short-time Fourier transform | |
| % data in D. Each column of D is taken as the result of an F-point | |
| % fft; each successive frame was offset by H points (default | |
| % W/2, or F/2 if W==0). Data is hann-windowed at W pts, or | |
| % W = 0 gives a rectangular window (default); | |
| % W as a vector uses that as window. | |
| % This version scales the output so the loop gain is 1.0 for | |
| % either hann-win an-syn with 25% overlap, or hann-win on | |
| % analysis and rect-win (W=0) on synthesis with 50% overlap. | |
| % dpwe 1994may24. Uses built-in 'ifft' etc. | |
| % $Header: /home/empire6/dpwe/public_html/resources/matlab/pvoc/RCS/istft.m,v 1.5 2010/08/12 20:39:42 dpwe Exp $ | |
| if nargin < 2; ftsize = 2*(size(d,1)-1); end | |
| if nargin < 3; w = 0; end | |
| if nargin < 4; h = 0; end % will become winlen/2 later | |
| s = size(d); | |
| if s(1) ~= (ftsize/2)+1 | |
| error('number of rows should be fftsize/2+1') | |
| end | |
| cols = s(2); | |
| if length(w) == 1 | |
| if w == 0 | |
| % special case: rectangular window | |
| win = ones(1,ftsize); | |
| else | |
| if rem(w, 2) == 0 % force window to be odd-len | |
| w = w + 1; | |
| end | |
| halflen = (w-1)/2; | |
| halff = ftsize/2; | |
| halfwin = 0.5 * ( 1 + cos( pi * (0:halflen)/halflen)); | |
| win = zeros(1, ftsize); | |
| acthalflen = min(halff, halflen); | |
| win((halff+1):(halff+acthalflen)) = halfwin(1:acthalflen); | |
| win((halff+1):-1:(halff-acthalflen+2)) = halfwin(1:acthalflen); | |
| % 2009-01-06: Make stft-istft loop be identity for 25% hop | |
| win = 2/3*win; | |
| end | |
| else | |
| win = w; | |
| end | |
| w = length(win); | |
| % now can set default hop | |
| if h == 0 | |
| h = floor(w/2); | |
| end | |
| xlen = ftsize + (cols-1)*h; | |
| x = zeros(1,xlen); | |
| for b = 0:h:(h*(cols-1)) | |
| ft = d(:,1+b/h)'; | |
| ft = [ft, conj(ft([((ftsize/2)):-1:2]))]; | |
| px = real(ifft(ft)); | |
| x((b+1):(b+ftsize)) = x((b+1):(b+ftsize))+px.*win; | |
| end; |
| """This file defines the Short-Time Fourier Transfor (STFT) | |
| or more precisely the Discrete Time and Frequency STFT. | |
| """ | |
| from __future__ import division | |
| from numpy import * | |
| from numpy.fft import fft | |
| from pytfd import helpers as h | |
| def stft(x, w, L=None): | |
| # L is the overlap, see http://cnx.org/content/m10570/latest/ | |
| N = len(x) | |
| #T = len(w) | |
| if L is None: | |
| L = N | |
| # Zerro pad the window | |
| w = h.zeropad(w, N) | |
| X_stft = [] | |
| points = range(0, N, N//L) | |
| for i in points: | |
| x_subset = h.subset(x, i, N) | |
| fft_subset = fft(x_subset * w) | |
| X_stft.append(fft_subset) | |
| X_stft = array(X_stft).transpose() | |
| return X_stft | |
| def spec(x, w): | |
| return abs(stft(x, w)) | |
| spectogram = spec | |
| __all__ = ['stft', 'spec', 'spectogram'] |
| function D = stft(x, f, w, h, sr) | |
| % D = stft(X, F, W, H, SR) Short-time Fourier transform. | |
| % Returns some frames of short-term Fourier transform of x. Each | |
| % column of the result is one F-point fft (default 256); each | |
| % successive frame is offset by H points (W/2) until X is exhausted. | |
| % Data is hann-windowed at W pts (F), or rectangular if W=0, or | |
| % with W if it is a vector. | |
| % Without output arguments, will plot like sgram (SR will get | |
| % axes right, defaults to 8000). | |
| % See also 'istft.m'. | |
| % dpwe 1994may05. Uses built-in 'fft' | |
| % $Header: /home/empire6/dpwe/public_html/resources/matlab/pvoc/RCS/stft.m,v 1.4 2010/08/13 16:03:14 dpwe Exp $ | |
| if nargin < 2; f = 256; end | |
| if nargin < 3; w = f; end | |
| if nargin < 4; h = 0; end | |
| if nargin < 5; sr = 8000; end | |
| % expect x as a row | |
| if size(x,1) > 1 | |
| x = x'; | |
| end | |
| s = length(x); | |
| if length(w) == 1 | |
| if w == 0 | |
| % special case: rectangular window | |
| win = ones(1,f); | |
| else | |
| if rem(w, 2) == 0 % force window to be odd-len | |
| w = w + 1; | |
| end | |
| halflen = (w-1)/2; | |
| halff = f/2; % midpoint of win | |
| halfwin = 0.5 * ( 1 + cos( pi * (0:halflen)/halflen)); | |
| win = zeros(1, f); | |
| acthalflen = min(halff, halflen); | |
| win((halff+1):(halff+acthalflen)) = halfwin(1:acthalflen); | |
| win((halff+1):-1:(halff-acthalflen+2)) = halfwin(1:acthalflen); | |
| end | |
| else | |
| win = w; | |
| end | |
| w = length(win); | |
| % now can set default hop | |
| if h == 0 | |
| h = floor(w/2); | |
| end | |
| c = 1; | |
| % pre-allocate output array | |
| d = zeros((1+f/2),1+fix((s-f)/h)); | |
| for b = 0:h:(s-f) | |
| u = win.*x((b+1):(b+f)); | |
| t = fft(u); | |
| % FILTERS APPLIED HERE | |
| lo = floor(f / 2 * 0 / (sr/2)) + 1; | |
| hi = floor(f / 2 * 1000 / (sr/2)) + 1; | |
| t(lo:hi) = t(lo:hi) * 0; | |
| d(:,c) = t(1:(1+f/2))'; | |
| c = c+1; | |
| end; | |
| % If no output arguments, plot a spectrogram | |
| if nargout == 0 | |
| tt = [0:size(d,2)]*h/sr; | |
| ff = [0:size(d,1)]*sr/f; | |
| imagesc(tt,ff,20*log10(abs(d))); | |
| axis('xy'); | |
| xlabel('time / sec'); | |
| ylabel('freq / Hz') | |
| % leave output variable D undefined | |
| else | |
| % Otherwise, no plot, but return STFT | |
| D = d; | |
| end |
| # Originally from http://stackoverflow.com/a/6891772/125507 | |
| import scipy, pylab | |
| def stft(x, fs, framesz, hop): | |
| """x is the time-domain signal | |
| fs is the sampling frequency | |
| framesz is the frame size, in seconds | |
| hop is the the time between the start of consecutive frames, in seconds | |
| """ | |
| framesamp = int(framesz*fs) | |
| hopsamp = int(hop*fs) | |
| w = scipy.hamming(framesamp) | |
| X = scipy.array([scipy.fft(w*x[i:i+framesamp]) | |
| for i in range(0, len(x)-framesamp, hopsamp)]) | |
| return X | |
| def istft(X, fs, T, hop): | |
| """X is the short-time Fourier transform | |
| fs is the sampling frequency | |
| T is the total length of the time-domain output in seconds | |
| hop is the the time between the start of consecutive frames, in seconds | |
| """ | |
| x = scipy.zeros(T*fs) | |
| framesamp = X.shape[1] | |
| hopsamp = int(hop*fs) | |
| for n,i in enumerate(range(0, len(x)-framesamp, hopsamp)): | |
| x[i:i+framesamp] += scipy.real(scipy.ifft(X[n])) | |
| return x | |
| if __name__ == '__main__': | |
| f0 = 440 # Compute the STFT of a 440 Hz sinusoid | |
| fs = 8000 # sampled at 8 kHz | |
| T = 5 # lasting 5 seconds | |
| framesz = 0.050 # with a frame size of 50 milliseconds | |
| hop = 0.020 # and hop size of 20 milliseconds. | |
| # Create test signal and STFT. | |
| t = scipy.linspace(0, T, T*fs, endpoint=False) | |
| x = scipy.sin(2*scipy.pi*f0*t) | |
| X = stft(x, fs, framesz, hop) | |
| # Plot the magnitude spectrogram. | |
| pylab.figure() | |
| pylab.imshow(scipy.absolute(X.T), origin='lower', aspect='auto', | |
| interpolation='nearest') | |
| pylab.xlabel('Time') | |
| pylab.ylabel('Frequency') | |
| pylab.show() | |
| # Compute the ISTFT. | |
| xhat = istft(X, fs, T, hop) | |
| # Plot the input and output signals over 0.1 seconds. | |
| T1 = int(0.1*fs) | |
| pylab.figure() | |
| pylab.plot(t[:T1], x[:T1], t[:T1], xhat[:T1]) | |
| pylab.xlabel('Time (seconds)') | |
| pylab.figure() | |
| pylab.plot(t[-T1:], x[-T1:], t[-T1:], xhat[-T1:]) | |
| pylab.xlabel('Time (seconds)') |