Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
STFT ISTFT Matlab Python
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)')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.