Project: Short-Time Fourier Transform (Small)

Table of Contents

Abstract
Fourier Transform
Short Time Fourier Transform
References

Developer	sonots
First Edition	04/2006
Last Modified	06/2006
Language	Matlab

Abstract

Short-Time Fourier Transform is a well studied filter bank. It can be seen in various ways, simply taking fourier transform in short time, low-pass filter applied for modulated signal, filter bank.

Fourier Transform

First, examine fourier transform.

Matlab Code

% NAME
%   spFft - Fourier Transform of a signal
% SYNOPSIS
%   [Y, YMag, YPhase, F] = spFft(X, window, show)
% DESCRIPTION
%   Get Frequency Response (and its magnitude and phase) of a signal and
%   plot
% INPUTS
%   X      (vector) of size Nx1 which contains signal
%   [window]
%          (string) the window function such as 'rectwin', 'hamming'. 
%           The default is 'rectwin' (or no window). 
%   [show] (boolean)   plot or not. The default is 0.
% OUTPUTS
%   Y      (vector) of size Nx1 which contains Frequency Response. 
%   [YMag] (vector) of size Nx1 which contains Magnitude Response. 
%   [YPhase] 
%          (vector) of size Nx1 which contains Phase Response. 
%   [F]    (vector) of size Nx1 which contains the x-axis (normalized 
%           frequency). 
% AUTHOR
%   Naotoshi Seo, April 2006
function [Y, YMag, YPhase, F] = spFft(X, window, show)
 %% Initialization
 N = length(X);
 X = X(:); % assure column vector
 if ~exist('show', 'var')
     show = 0;
 end
 if ~exist('window', 'var')
     window = 'rectwin';
 end

 %% do fourier transform of a windowed signal
 if ~strcmp(window, 'rectwin')
     X = X .* eval(sprintf('%s(N)', window)); % X .* hamming(N)
 end
 Y = fft(X, N);

 % Magnitude Response in dB (20 * log10)
 YMag = 20 * log10(abs(Y));
 % Phase Response
 YPhase = unwrap(angle(Y));
 % Normalized Frequency
 F = (0:N-1)/N; % * fs;

 %% Plot
 if show
     % plot magnitude
     figure;
     subplot(2,1,1);
     plot(F, YMag); grid on;
     xlabel('Normalized Frequency (\times\pi rad/sample)');
     ylabel('Magnitude (dB)');
     title('Magnitude Response');
     % plot phase
     subplot(2,1,2);
     plot(F, YPhase); grid on;
     xlabel('Normalized Frequency (\times\pi rad/sample)');
     ylabel('Phase (radian)');
     title('Phase Response');
 end
end

Experiments

function spFftDemo
 [x, fs] = wavread('bee.wav');
 %x = wgn(1, 1000, 2);
 %t = linspace(-10,10);
 %x = sinc(t);
 y = spFft(x, 'rectwin', 'plot');

Input data

bee.wav

Result

FYI: fvtool(x); allows you to do DFT plot, too. The Fourier Analysis of LTI system can be done by freqz.

Short Time Fourier Transform

Math Equations

Refer [1] or [2]. Refer [2] especially for inverse STFT.

Matlab Code

% NAME
%   spStft: Short-Time Fourier Transform. 
% SYNOPSIS
%   [S, F, T] = spStft(x, [window], [frame_overlap], [frame_length],
%               [fs], [show])
% DESCRIPTION
%   Short-Time Fourier Transform
% INPUTS
%   x      (vector) of size Nx1.
%   window (string) the window function such as rectwin, hamming.  
%           if not specified, equivalent to hamming
%   frame_overlap 
%          (scalar) the length of each frame overlaps in micro second. 
%           The default is frame_length / 2. 
%   frame_length
%          (scalar) the length of each frame in micro second. 
%           The default is 20ms.  
%   fs     (scalar) the sampling rate in Hz. The default assumes 16kHz
%   [show] (boolean) plot or not. The default is 0.
% OUTPUTS
%   S      (matrix) of size MxK where M is the size of each frame(fft))
%           and K is the number of frames. S is spectrogram (complex number)
%   F      (vector) of size Mx1 that is a vector of frequencies (y-axis)
%           in Hz. 
%   T      (vector) of size 1xK whose value corresponds to the center of 
%           each frame (x-axis) in sec.
% AUTHOR
%   Naotoshi Seo, April 2006
% SEE ALSO
%   spectrogram.m (digital signal processing toolbox) 
%    - arguments order was arranged similar to this (but in micro second). 
function [S, F, T] = spStft(x, window, frame_overlap, frame_length, fs, show)
 %% Initialization
 N = length(x);
 if ~exist('window', 'var') || isempty(window)
     window = 'hamming';
 end
 if ~exist('fs', 'var') || isempty(fs)
     fs = 16000;
 end
 if ~exist('frame_length', 'var') || isempty(frame_length)
     frame_length = 20;
 end
 if ~exist('frame_overlap', 'var') || isempty(frame_overlap)
     frame_overlap = frame_length / 2;
 end
 if ~exist('show', 'var') || isempty(show)
     show = 0;
 end
 nsample = round(frame_length  * fs / 1000); % convert ms to points
 noverlap = round(frame_overlap * fs / 1000); % convert ms to points
 window   = eval(sprintf('%s(nsample)', window)); % e.g., hamming(nsample)

 %% spectrogram
 [S, F, T] = gen_spectrogram(x, window, noverlap, nsample, fs); % below

 %% plot
 if show
     plot_spectrogram(S, F, T); % below
 end

function [S, F, T] = gen_spectrogram(x, window, noverlap, nsample, fs)
 N = length(x);
 S = [];
 pos = 1;
 while (pos+nsample <= N)
     frame = x(pos:pos+nsample-1);
     pos = pos + (nsample - noverlap);
     Y = fft(frame .* window, nsample);
     % see also goertzel. signal/spectrogram.m is using it. 
     S = [S Y(1:round(nsample/2), 1)]; % half is enough, another half is just mirror
 end
 [M, K] = size(S);
 F = (0:round(nsample/2)-1)' / nsample * fs; % [0, fs/2) Hz 
 % F = psdfreqvec(nsample, fs, 'half');
 T = (round(nsample/2):(nsample-noverlap):N-1-round(nsample/2))/fs;

% equals to spectrogram(..., 'yaxis'), no argout
function plot_spectrogram(S, F, T)
 SdB = 20 * log10(abs(S)); % dB
 ax = imagesc(T, F, SdB);
 set(get(ax, 'Parent'), 'YDir', 'normal');
 xlabel('Time')
 ylabel('Frequency (Hz)')
 % colorbar

Experiments

function spStftDemo
 [x, fs] = wavread('bee.wav');
 %x = wgn(1, 1000, 2); fs = 16000;
 frame_overlap = 10; % ms
 frame_length  = 20;
 spStft(x(:), 'hamming', frame_overlap, frame_length, fs, 'plot');

spectrogram.m (signal processing toolbox) version

function spStftDemo2
% use spectrogram (digital signal processing toolbox)
 [x, fs] = wavread('bee.wav');
 frame_overlap = 10; % ms
 frame_length  = 20;
 window        = 'rectwin';

 nfft = round(frame_length  * fs / 1000); % convert ms to points
 noverlap = round(frame_overlap * fs / 1000); % convert ms to points
 window   = eval(sprintf('%s(nfft)', window)); % e.g., hamming(nfft)

 [S, F, T, P] = spectrogram(x, window, noverlap, nfft, fs);
 spectrogram(x, window, noverlap, nfft, fs, 'yaxis'); % plot

Input data

bee.wav

Results

spStft.m

spectrogram.m (signal processing toolbox)

They are a bit different because function spectrogram uses goertzel function which computes the discrete Fourier transform (DFT) using second-order goertzel algorithm although my spectrogram uses fft function.

Discussion

As the frame (window, segment) length increases, frequency resolutions are increased, however, time resolutions are decreased. As the frame (window, segment) length decreases, time resolutions are increased, however, frequency resolutions are decreased.

In the case of wavelet, we do not need to care of this.

Short report

STFT.pdf

References

[1] Xuedong Huang, Alex Acero, Hsiao-wuen hon, "Spoken Language Processing: A guide to Theory, Algorithm, and System Development," Prentice Hall, January 2001, p276-83
[2] J. B. Allen, \Applications of the Short Time Fourier Transform to Speech Processing and Spectral Analysis," Acoustics, Speech, and Signal Processing, IEEE International Conference on ICASSP '82. , vol.7, no., pp. 1012-1015, May 1982

Naotoshi Seo

Announce