 ## Announce

Puki contents have been moved into SONOTS Plugin (20070703)

## Project: Short-Time Fourier Transform (Small)

Developer sonots 04/2006 06/2006 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;
ylabel('Magnitude (dB)');
title('Magnitude Response');
% plot phase
subplot(2,1,2);
plot(F, YPhase); grid on;
title('Phase Response');
end
end

### Experiments

function spFftDemo
%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  or . Refer  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
%   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);
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 = 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)
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

•  Xuedong Huang, Alex Acero, Hsiao-wuen hon, "Spoken Language Processing: A guide to Theory, Algorithm, and System Development," Prentice Hall, January 2001, p276-83
•  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