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

Table of Contents |

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.

Tag: Scientific SoundProcessing SignalProcessing STFT Matlab FFT

## 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

### 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

### 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

## 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