FFT-based FIR filtering using overlap-add method
filter is more efficient for smaller operands and
fftfilt is more efficient for large operands. Filter random numbers with two random filters: a short one, with 20 taps, and a long one, with 2000. Use
toc to measure the execution times. Repeat the experiment 100 times to improve the statistics.
rng default N = 100; shrt = 20; long = 2000; tfs = 0; tls = 0; tfl = 0; tll = 0; for kj = 1:N x = rand(1,1e6); bshrt = rand(1,shrt); tic sfs = fftfilt(bshrt,x); tfs = tfs+toc/N; tic sls = filter(bshrt,1,x); tls = tls+toc/N; blong = rand(1,long); tic sfl = fftfilt(blong,x); tfl = tfl+toc/N; tic sll = filter(blong,1,x); tll = tll+toc/N; end
Compare the average times.
fprintf('%4d-tap filter averages: fftfilt: %f s; filter: %f s\n',shrt,tfs,tls)
20-tap filter averages: fftfilt: 0.178331 s; filter: 0.012774 s
fprintf('%4d-tap filter averages: fftfilt: %f s; filter: %f s\n',long,tfl,tll)
2000-tap filter averages: fftfilt: 0.063776 s; filter: 0.150068 s
This example requires Parallel Computing Toolbox™ software. Refer to GPU Support by Release (Parallel Computing Toolbox) for a list of supported GPUs.
Create a signal consisting of a sum of sine waves in white Gaussian additive noise. The sine wave frequencies are 2.5, 5, 10, and 15 kHz. The sampling frequency is 50 kHz.
Fs = 50e3; t = 0:1/Fs:10-(1/Fs); x = cos(2*pi*2500*t) + 0.5*sin(2*pi*5000*t) + 0.25*cos(2*pi*10000*t)+ ... 0.125*sin(2*pi*15000*t) + randn(size(t));
Design a lowpass FIR equiripple filter using
d = designfilt('lowpassfir','SampleRate',Fs, ... 'PassbandFrequency',5500,'StopbandFrequency',6000, ... 'PassbandRipple',0.5,'StopbandAttenuation',50); B = d.Coefficients;
Filter the data on the GPU using the overlap-add method. Put the data on the GPU using
gpuArray. Return the output to the MATLAB® workspace using
gather and plot the power spectral density estimate of the filtered data.
y = fftfilt(gpuArray(B),gpuArray(x)); periodogram(gather(y),rectwin(length(y)),length(y),50e3)
b— Filter coefficients
Filter coefficients, specified as a vector. If
b is a matrix,
fftfilt applies the filter in each column of
b to the signal vector
x— Input data
Input data, specified as a vector. If
x is a matrix,
fftfilt filters its columns. If
x are both matrices with the same number of columns, the
ith column of
b is used to filter the
ith column of
fftfilt works for both real and complex inputs.
n— FFT length
FFT length, specified as a positive integer. By default,
fftfilt chooses an FFT length and a data block length that
guarantee efficient execution time.
gpuArrayX— GPU arrays
GPU arrays, specified as a
gpuArrayb contains the filter coefficients, and
gpuArrayX is the input data. See Run MATLAB Functions on a GPU (Parallel Computing Toolbox) for details on
Computing Toolbox™ software. Refer to GPU Support by Release (Parallel Computing Toolbox) for a list of
supported GPUs. The filtered data,
y, is a
gpuArray object. See Overlap-Add Filtering on the GPU for an example of
overlap-add filtering on the GPU.
y— Output data
Output data, returned as a vector, matrix or
When the input signal is relatively large,
is faster than
filter performs N
multiplications for each sample in
x, where N is the
fftfilt performs 2 FFT operations — the FFT of
the signal block of length L plus the inverse FT of the product of the
FFTs — at the cost of where L is the block length. It then performs
L point-wise multiplications for a total cost of multiplications. The cost ratio is therefore which is approximately
log2L / N.
fftfilt is faster when
log2L is less than N.
fftfilt filters data using the efficient FFT-based method of
, a frequency domain
filtering technique that works only for FIR filters by combining successive frequency domain
filtered blocks of an input sequence. The operation performed by
is described in the time domain by the difference equation:
An equivalent representation is the Z-transform or frequency domain description:
fft to implement the overlap-add method.
fftfilt breaks an
x into length L data blocks, where
L must be greater than the filter length N.
and convolves each block with the filter
y = ifft(fft(x(i:i+L-1),nfft).*fft(b,nfft));
nfft is the FFT length.
successive output sections by
n-1 points, where
n is the
length of the filter, and sums them.
fftfilt chooses the key parameters
nfft in different ways, depending on whether you supply an FFT length
n for the filter and signal. If you do not specify a value for
n (which determines FFT length),
these key parameters automatically:
length(x) is greater than
fftfilt chooses values that minimize the number of blocks times the
number of flops per FFT.
length(b) is greater than or equal to
fftfilt uses a single FFT of
2^nextpow2(length(b) + length(x) - 1)
y = ifft(fft(B,nfft).*fft(X,nfft))
If you supply a value for
fftfilt chooses an
2^nextpow2(n) and a data block
n is less than
 Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.
Usage notes and limitations:
Digital filter objects are not supported for code generation.
This function fully supports GPU arrays. For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).