FFT-based FIR filtering using overlap-add method
y = fftfilt(b,x)
y = fftfilt(b,x,n)
y = fftfilt(d,x)
y = fftfilt(d,x,n)
y = fftfilt(gpuArrayb,gpuArrayX,n)
fftfilt filters data using the efficient FFT-based method of
overlap-add, a frequency domain filtering technique that works
only for FIR filters.
y = fftfilt(b,x) filters the data in
x with the filter described by coefficient vector
b. It returns the data vector
y. The operation
fftfilt is described in the time
domain by the difference equation:
An equivalent representation is the Z-transform or frequency domain description:
fftfilt chooses an FFT length and a data block length
that guarantee efficient execution time.
x is a matrix,
fftfilt filters its columns.
b is a matrix,
fftfilt applies the filter in
each column of
b to the signal vector
x are both matrices with the same number
of columns, the ith column of
b is used to filter
the ith column of
y = fftfilt(b,x,n) uses
n to determine the length of the FFT. See Algorithms for information.
y = fftfilt(d,x,n) uses
n to determine the
length of the FFT.
y = fftfilt(gpuArrayb,gpuArrayX,n) filters the data in the
gpuArrayX, with the FIR filter coefficients in the
gpuArrayb. 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) to see what GPUs
are supported. The filtered data,
y, is a
gpuArray object. See Overlap-Add Filtering on the GPU for example of overlap-add filtering on
fftfilt works for both real and complex inputs.
When the input signal is relatively large, it is advantageous to use
fftfilt instead of
filter, which performs
N multiplications for each sample in
where N is the filter length.
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
L + Llog2L = L(1 + log2L)
multiplications. The cost ratio is therefore
L(1+log2L) / (NL) = (1 + log2L)/N
which is approximately log2L / N.
fftfilt becomes advantageous when
log2Lis less than
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.289484 s; filter: 0.024282 s
fprintf('%4d-tap filter averages: fftfilt: %f s; filter: %f s\n',long,tfl,tll)
2000-tap filter averages: fftfilt: 0.096268 s; filter: 0.354657 s
This example requires Parallel Computing Toolbox™ software. Refer to GPU Support by Release (Parallel Computing Toolbox) to see what GPUs are supported.
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)
fft to implement the overlap-add method
, a technique that combines successive frequency domain filtered blocks
of an input sequence.
fftfilt breaks an input sequence
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
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
n and on the lengths of the filter and signal. If you do not
specify a value for
n (which determines FFT length),
fftfilt chooses 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
2^nextpow2(length(b) + length(x) - 1)
This essentially computes
y = ifft(fft(B,nfft).*fft(X,nfft))
If you supply a value for
an FFT length,
2^nextpow2(n)and a data
block length of
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.
This function fully supports GPU arrays. For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).