fftfilt

FFT-based FIR filtering using overlap-add method

Syntax

y = fftfilt(b,x)
y = fftfilt(b,x,n)
y = fftfilt(d,x)
y = fftfilt(d,x,n)
y = fftfilt(gpuArrayb,gpuArrayX,n)

Description

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 vector x with the filter described by coefficient vector b. It returns the data vector y. The operation performed by fftfilt is described in the time domain by the difference equation:

y(n)=b(1)x(n)+b(2)x(n1)++b(nb+1)x(nnb)

An equivalent representation is the Z-transform or frequency domain description:

Y(z)=(b(1)+b(2)z1++b(nb+1)znb)X(z)

By default, fftfilt chooses an FFT length and a data block length that guarantee efficient execution time.

If x is a matrix, fftfilt filters its columns. If b is a matrix, fftfilt applies the filter in each column of b to the signal vector x. If b and x are both matrices with the same number of columns, the ith column of b is used to filter the ith column of x.

y = fftfilt(b,x,n) uses n to determine the length of the FFT. See Algorithms for information.

y = fftfilt(d,x) filters the data in vector x with a digitalFilter object, d. Use designfilt to generate d based on frequency-response specifications.

y = fftfilt(d,x,n) uses n to determine the length of the FFT.

y = fftfilt(gpuArrayb,gpuArrayX,n) filters the data in the gpuArray object, gpuArrayX, with the FIR filter coefficients in the gpuArray, gpuArrayb. See Establish Arrays on a GPU for details on gpuArray objects. Using fftfilt with gpuArray objects requires Parallel Computing Toolbox™ software and a CUDA-enabled NVIDIA GPU with compute capability 1.3 or above. See http://www.mathworks.com/products/parallel-computing/requirements.html for details. The filtered data, y, is a gpuArray object. See Overlap-Add Filtering on the GPU for example of overlap-add filtering on the GPU.

fftfilt works for both real and complex inputs.

Comparison to filter function

When the input signal is relatively large, it is advantageous to use fftfilt instead of filter, which performs N multiplications for each sample in x, where N is the filter length. 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

½Llog2L

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.

Therefore, fftfilt becomes advantageous when log2Lis less than N.

Examples

expand all

fftfilt and filter for Short and Long Filters

Verify that filter is more efficient for smaller operands and fftfilt is more efficient for large operands. Filter 10^6 random numbers with two random filters, a short one, with 20 taps, and a long one, with 2000. Use tic and toc to measure the execution times.

rng default
x = rand(1,1e6);

shrt = 20;
fprintf('Filter of length %d:\n',shrt)
bshrt = rand(1,shrt);

tic
sfs = fftfilt(bshrt,x);
tsfs = toc;

tic
sls = filter(bshrt,1,x);
tsls = toc;

fprintf('fftfilt takes %f s to run, filter takes %f s to run\n',tsfs,tsls)

long = 2000;
fprintf('\nFilter of length %d:\n',long)
blong = rand(1,long);

tic
sfl = fftfilt(blong,x);
tsfl = toc;

tic
sll = filter(blong,1,x);
tsll = toc;

fprintf('fftfilt takes %f s to run, filter takes %f s to run\n',tsfl,tsll)
Filter of length 20:
fftfilt takes 1.364652 s to run, filter takes 0.013516 s to run

Filter of length 2000:
fftfilt takes 0.177251 s to run, filter takes 0.241303 s to run

Overlap-Add Filtering on the GPU

The following example requires Parallel Computing Toolbox software and a CUDA-enabled NVIDIA GPU with compute capability 1.3 or above. See http://www.mathworks.com/products/parallel-computing/requirements.html for details.

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 designfilt.

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);

More About

expand all

Algorithms

fftfilt uses fft to implement the overlap-add method [1], 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 b by

y = ifft(fft(x(i:i+L-1),nfft).*fft(b,nfft));

where nfft is the FFT length. fftfilt overlaps successive output sections by n-1 points, where n is the length of the filter, and sums them.

fftfilt chooses the key parameters L and nfft in different ways, depending on whether you supply an FFT length 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:

  • If length(x)is greater than length(b), fftfilt chooses values that minimize the number of blocks times the number of flops per FFT.

  • If length(b) is greater than or equal to length(x), fftfilt uses a single FFT of length

    2^nextpow2(length(b) + length(x) - 1)
    

    This essentially computes

    y = ifft(fft(B,nfft).*fft(X,nfft))
    

If you supply a value for n, fftfilt chooses an FFT length, nfft, of 2^nextpow2(n)and a data block length of nfft - length(b) + 1. If n is less than length(b), fftfilt sets n to length(b).

References

[1] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

Was this topic helpful?