18
$\begingroup$

I'm trying to conceptually understand what is happening when the forward and inverse Short Time Fourier Transforms (STFT) are applied to a discrete time-domain signal. I've found the classic paper by Allen and Rabiner (1977), as well as a Wikipedia article (link). I believe that there is also another good article to be found here.

I am interested in computing the Gabor transform, which is nothing more than the STFT with a Gaussian window.

This is what I understand about the forward STFT:

  1. A sub-sequence is selected from the signal, comprised of time domain elements.
  2. The sub-sequence is multiplied by a window function using point-by-point multiplication in the time domain.
  3. The multiplied sub-sequence is taken into the frequency domain using the FFT.
  4. By selecting successive overlapping sub-sequences, and repeating the procedure above, we get a matrix with m rows and n columns. Each column is the sub-sequence computed at a given time. This can be used to compute a spectrogram.

However, for the inverse STFT, the papers talk about a summation over the overlapping analysis sections. I am finding it very challenging to visualize what is really going on here. What do I have to do to be able to compute the inverse STFT (in step by step order as above)?

Forward STFT

I've created a drawing showing what I think is going on for the forward STFT. What I don't understand is how to assemble each of the sub-sequences so that I get back the original time sequence. Could someone either modify this drawing or give an equation showing how the sub-sequences are added? Forward transform

Inverse Transform

Here's what I understand about the inverse transform. Each successive window is taken back into the time domain using the IFFT. Then each window is shifted by the step size, and added to the result of the previous shift. The following diagram shows this process. The summed output is the time domain signal.

Inverse transform

Code Example

The following Matlab code generates a synthetic time domain signal, and then tests the STFT process, demonstrating that the inverse is the dual of the forward transform, within numerical round-off error. The beginning and end of the signal is zero-padded to ensure that the center of the window can be situated at the first and last elements of the time-domain signal.

NOTE that according to Allen and Rabiner (1977), if a multiplication happens in the frequency domain to change the frequency response, the length of the analysis window must be equal or greater than $N + N_0 - 1$ points, where $N_0$ is the filter response. The length is extended by zero-padding. The test code simply shows that the inverse is the dual of the forward transform. The length must be extended to prevent a circular convolution.

% The code computes the STFT (Gabor transform) with step size = 1
% This is most useful when modifications of the signal is required in
% the frequency domain

% The Gabor transform is a STFT with a Gaussian window (w_t in the code)

% written by Nicholas Kinar

% Reference:
% [1] J. B. Allen and L. R. Rabiner, 
% “A unified approach to short-time Fourier analysis and synthesis,” 
% Proceedings of the IEEE, vol. 65, no. 11, pp. 1558 – 1564, Nov. 1977.

% generate the signal
mm = 8192;                  % signal points
t = linspace(0,1,mm);       % time axis

dt = t(2) - t(1);           % timestep t
wSize = 101;                % window size


% generate time-domain test function
% See pg. 156
% J. S. Walker, A Primer on Wavelets and Their Scientific Applications, 
% 2nd ed., Updated and fully rev. Boca Raton: Chapman & Hall/CRC, 2008.
% http://www.uwec.edu/walkerjs/primer/Ch5extract.pdf
term1 = exp(-400 .* (t - 0.2).^2);
term2 = sin(1024 .* pi .* t);
term3 = exp(-400.*(t- 0.5).^2);
term4 = cos(2048 .* pi .* t);
term5 = exp(-400 .* (t-0.7).^2);
term6 = sin(512.*pi.*t) - cos(3072.*pi.*t);
u = term1.*term2  + term3.*term4 + term5.*term6; % time domain signal
u = u';

figure;
plot(u)

Nmid = (wSize - 1) / 2 + 1;    % midway point in the window
hN = Nmid - 1;                 % number on each side of center point       


% stores the output of the Gabor transform in the frequency domain
% each column is the FFT output
Umat = zeros(wSize, mm);     


% generate the Gaussian window 
% [1] Y. Wang, Seismic inverse Q filtering. Blackwell Pub., 2008.
% pg. 123.
T = dt * hN;                    % half-width
sp = linspace(dt, T, hN); 
targ = [-sp(end:-1:1) 0 sp];    % this is t - tau
term1 = -((2 .* targ) ./ T).^2;
term2 = exp(term1);
term3 = 2 / (T * sqrt(pi));
w_t = term3 .* term2;
wt_sum = sum ( w_t ); % sum of the wavelet


% sliding window code
% NOTE that the beginning and end of the sequence
% are padded with zeros 
for Ntau = 1:mm

    % case #1: pad the beginning with zeros
    if( Ntau <= Nmid )
        diff = Nmid - Ntau;
        u_sub = [zeros(diff,1); u(1:hN+Ntau)];
    end

    % case #2: simply extract the window in the middle
    if (Ntau < mm-hN+1 && Ntau > Nmid)
        u_sub = u(Ntau-hN:Ntau+hN);
    end

    % case #3: less than the end
    if(Ntau >= mm-hN+1)
        diff = mm - Ntau;
        adiff = hN - diff;
        u_sub = [ u(Ntau-hN:Ntau+diff);  zeros(adiff,1)]; 
    end   

    % windowed trace segment
    % multiplication in time domain with
    % Gaussian window  function
    u_tau_omega = u_sub .* w_t';

    % segment in Fourier domain
    % NOTE that this must be padded to prevent
    % circular convolution if some sort of multiplication
    % occurs in the frequency domain
    U = fft( u_tau_omega );

    % make an assignment to each trace
    % in the output matrix
    Umat(:,Ntau) = U;

end

% By here, Umat contains the STFT (Gabor transform)

% Notice how the Fourier transform is symmetrical 
% (we only need the first N/2+1
% points, but I've plotted the full transform here
figure;
imagesc( (abs(Umat)).^2 )


% now let's try to get back the original signal from the transformed
% signal

% use IFFT on matrix along the cols
us = zeros(wSize,mm);
for i = 1:mm 
    us(:,i) = ifft(Umat(:,i));
end

figure;
imagesc( us );

% create a vector that is the same size as the original signal,
% but allows for the zero padding at the beginning and the end of the time
% domain sequence
Nuu = hN + mm + hN;
uu = zeros(1, Nuu);

% add each one of the windows to each other, progressively shifting the
% sequence forward 
cc = 1; 
for i = 1:mm
   uu(cc:cc+wSize-1) = us(:,i) + uu(cc:cc+wSize-1)';
   cc = cc + 1;
end

% trim the beginning and end of uu 
% NOTE that this could probably be done in a more efficient manner
% but it is easiest to do here

% Divide by the sum of the window 
% see Equation 4.4 of paper by Allen and Rabiner (1977)
% We don't need to divide by L, the FFT transform size since 
% Matlab has already taken care of it 
uu2 = uu(hN+1:end-hN) ./ (wt_sum); 

figure;
plot(uu2)

% Compare the differences bewteen the original and the reconstructed
% signals.  There will be some small difference due to round-off error
% since floating point numbers are not exact
dd = u - uu2';

figure;
plot(dd);
$\endgroup$
  • 2
    $\begingroup$ Great question - but, how did you make those diagrams this quickly on the fly?... $\endgroup$ – Spacey Jul 8 '12 at 17:53
  • 2
    $\begingroup$ I used Adobe Illustrator for the diagrams and Mathtype for the Greek characters. $\endgroup$ – Nicholas Kinar Jul 8 '12 at 22:36
  • 1
    $\begingroup$ "I am interested in computing the Gabor transform, which is nothing more than the STFT with a Gaussian window." Remember that Gabor transform is a continuous integral, and that Gaussian windows extend to infinity. Typical implementations of STFT use discrete overlapped chunks and have to use finite-length windows. $\endgroup$ – endolith Jan 13 '15 at 14:43
  • $\begingroup$ Thanks for pointing that out, endolith. I tend to think in a very discrete manner when doing signal processing. $\endgroup$ – Nicholas Kinar Jan 23 '15 at 15:44
11
$\begingroup$

The STFT transform pair can be characterized by 4 different parameters:

  1. FFT size (N)
  2. Step size (M)
  3. Analysis window (size N)
  4. Synthesis window (size N)

The process is as follows:

  1. Grab N (fft size) samples from the current input location
  2. Apply analysis window
  3. Do the FFT
  4. Do whatever you want to do in the frequency domain
  5. Inverse FFT
  6. Apply synthesis window
  7. Add to the output at the current output location
  8. Advance input and output location by M (step size) samples

The overlap add algorithm is a good example for that. In this case the step size is N, the FFT size is 2*N, the analysis window is rectangular with N ones followed by N zeros and the synthesis windows is simply all ones.

There are lot of other choices for that and under certain conditions the forward/inverse transfer is fully reconstructing (i.e. you can get the original signal back).

The key point here is that each output sample typically receives additive contributions from more than one inverse FFT. The output needs to be accumulated over multiple frames. The number of contributing frames is simply given by the FFT size divided by the step size (rounded up, if necessary).

$\endgroup$
  • $\begingroup$ Thank you very much for your insightful answer. I understand the overlap-add method. What do I use for the synthesis window? Is there an equation? If I know the analysis window function (such as a Gaussian window), how do I compute the synthesis window? I understand how the overlap-add method is used for convolution, but I don't understand how it is used for the STFT. If the step size is step = 1, how do I add the frames together? Is there an equation? $\endgroup$ – Nicholas Kinar Jul 2 '12 at 4:01
  • $\begingroup$ If the analysis window function is centered at each sample with step size step = 1, do I zero-pad the beginning and end of the time-domain sequence so that the middle of the window is centered at each sample (including the first and last samples in the time-domain sequence)? $\endgroup$ – Nicholas Kinar Jul 2 '12 at 4:02
  • $\begingroup$ You can choose step size, fft size, analysis and synthesis windows depending on the specific needs of your application. One example is step size N, FFT size 2*N, analysis hanning, synthesis all ones. You can modify that to analysis sqrt(hanning) and synthesis sqrt(hanning). Either one will work. I boils down to what you do in the frequency domain and what type of artifacts such as time domain aliasing you may create. $\endgroup$ – Hilmar Jul 2 '12 at 13:37
  • $\begingroup$ @Hilmar: I need to be able to make frequency-domain modifications to the signal, and then take the IFFT to get a time domain signal. I would like to minimize time domain aliasing. I still don't understand how to take each sub-sequence back into the time domain and then add them together. $\endgroup$ – Nicholas Kinar Jul 2 '12 at 18:18
  • $\begingroup$ I've written some test code and then updated my original question. $\endgroup$ – Nicholas Kinar Jul 2 '12 at 23:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.