Next  |  Prev  |  Up  |  Top  |  Index  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search


Example 1: Low-Pass Filtering by FFT Convolution

In this example, we design and implement a length $ L=257$ FIR lowpass filter having a cut-off frequency at $ f_c = 600$ Hz. The filter is tested on an input signal $ x(n)$ consisting of a sum of sinusoidal components at frequencies $ (440, 880, 1000, 2000)$ Hz. We'll filter a single input frame of length $ M=256$ , which allows the FFT to be $ N=512$ samples (no wasted zero-padding).

% Signal parameters:
f = [ 440 880 1000 2000 ];      % frequencies   
M = 256;                        % signal length 
Fs = 5000;                      % sampling rate

% Generate a signal by adding up sinusoids:
x = zeros(1,M); % pre-allocate 'accumulator'
n = 0:(M-1);    % discrete-time grid 
for fk = f; 
    x = x + sin(2*pi*n*fk/Fs); 
end

Next we design the lowpass filter using the window method:

% Filter parameters:
L = 257;    % filter length 
fc = 600;   % cutoff frequency 

% Design the filter using the window method:
hsupp = (-(L-1)/2:(L-1)/2);
hideal = (2*fc/Fs)*sinc(2*fc*hsupp/Fs);
h = hamming(L)' .* hideal; % h is our filter

Figure 8.3: FIR filter impulse response (top) and amplitude response (bottom).
\includegraphics[width=\twidth]{eps/filter}

Figure 8.3 plots the impulse response and amplitude response of our FIR filter designed by the window method. Next, the signal frame and filter impulse response are zero-padded out to the FFT size and transformed:

% Choose the next power of 2 greater than L+M-1 
Nfft = 2^(ceil(log2(L+M-1))); % or 2^nextpow2(L+M-1)

% Zero pad the signal and impulse response:
xzp = [ x zeros(1,Nfft-M) ];
hzp = [ h zeros(1,Nfft-L) ];

X = fft(xzp); % signal
H = fft(hzp); % filter

Figure 8.4 shows the input signal spectrum and the filter amplitude response overlaid. We see that only one sinusoidal component falls within the pass-band.

Figure 8.4: Overlay of input signal spectrum and desired lowpass filter pass-band.
\includegraphics[width=\twidth,height=1.8in]{eps/signal_transform}

Figure 8.5: Output signal magnitude spectrum = magnitude of input spectrum times filter frequency response.
\includegraphics[width=\twidth,height=1.8in]{eps/filtered_transform}

Now we perform cyclic convolution in the time domain using pointwise multiplication in the frequency domain:

Y = X .* H;
The modified spectrum is shown in Fig.8.5.

The final acyclic convolution is the inverse transform of the pointwise product in the frequency domain. The imaginary part is not quite zero as it should be due to finite numerical precision:

y = ifft(Y);
relrmserr = norm(imag(y))/norm(y) % check... should be zero
y = real(y);

Figure 8.6: Filtered output signal, with close-up showing the filter start-up transient (``pre-ring'').
\includegraphics[width=\twidth]{eps/filteredSignalAnn}

Figure 8.6 shows the filter output signal in the time domain. As expected, it looks like a pure tone in steady state. Note the equal amounts of ``pre-ringing'' and ``post-ringing'' due to the use of a linear-phase FIR filter.9.5

For an input signal approximately $ 4000$ samples long, this example is 2-3 times faster than the conv function in Matlab (which is precompiled C code implementing time-domain convolution).


Next  |  Prev  |  Up  |  Top  |  Index  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search

[How to cite this work]  [Order a printed hardcopy]  [Comment on this page via email]
[Lecture Video]  [Exercises]  [Examination]  
``Spectral Audio Signal Processing'', by Julius O. Smith III, W3K Publishing, 2011, ISBN 978-0-9745607-3-1.
Copyright © 2014-03-23 by Julius O. Smith III
Center for Computer Research in Music and Acoustics (CCRMA),   Stanford University
CCRMA