fftfilt

FFT-based FIR filtering using overlap-add method

Syntax

y = fftfilt(b,x)

y = fftfilt(b,x,nfft)

y = fftfilt(d,x)

y = fftfilt(d,x,nfft)

Description

y = fftfilt(b,x) filters the data specified in vector x. The function uses the filter described by the coefficient vector b.

y = fftfilt(b,x,nfft) uses nfft to determine the length of the FFT.

y = fftfilt(d,x) filters the data in vector x with a digitalFilter object d.

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

Examples

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. Repeat the experiment 100 times to improve the statistics.

rng("default")

N = 100;

s = 20;
l = 2000;

tfs = 0;
tls = 0;
tfl = 0;
tll = 0;

for kj = 1:N
    
    x = rand(1,1e6);

    bshrt = rand(1,s);

    tic
    sfs = fftfilt(bshrt,x);
    tfs = tfs+toc/N;

    tic
    sls = filter(bshrt,1,x);
    tls = tls+toc/N;

    blong = rand(1,l);

    tic
    sfl = fftfilt(blong,x);
    tfl = tfl+toc/N;
    
    tic
    sll = filter(blong,1,x);
    tll = tll+toc/N;

end

Compare and display the average times.

table(table(1000*[tfs;tls],1000*[tfl;tll], ...
    RowNames=["fftfilt" "filter"],VariableNames=[s;l]+"-tap"), ...
    VariableNames="Filter averages (milliseconds)")

ans=2×1 table
            Filter averages (milliseconds)        
    ______________________________________________

               20-tap    2000-tap
               ______    ________
                                                  
    fftfilt    142.65     60.027 
    filter     10.064     149.01 

Overlap-Add Filtering on the GPU

This example requires Parallel Computing Toolbox™ software. Refer to GPU Computing Requirements (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 designfilt.

d = designfilt("lowpassfir",SampleRate=Fs, ...
    PassbandFrequency=5500,StopbandFrequency=6000, ...
    PassbandRipple=0.5,StopbandAttenuation=50);
B = d.Numerator;

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)

Input Arguments

b — Filter coefficients

vector | matrix

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.

x — Input data

vector | matrix

Input data, specified as a vector. If x is a matrix, fftfilt filters its columns. 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. fftfilt works for both real and complex inputs.

nfft — FFT length

positive integer

FFT length, specified as a power-of-2 positive integer greater than or equal to the filter length.

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

• If you do not specify nfft as a power-of-2 positive integer, or if nfft is less than the filter length, fftfilt chooses an FFT length of 2^nextpow2(max(nfft,N)), where N is the filter length.

  • N = numel(b) if b is a vector.

  • N = size(b,1) if b is a matrix.

d — Digital filter

digitalFilter object

Digital filter, specified as a digitalFilter object. Use designfilt to generate d based on frequency-response specifications.

Output Arguments

y — Output data

vector | matrix

Output data, returned as a vector or matrix.

More About

Comparison to filter function

When the input signal is relatively large, fftfilt is faster than filter.

filter 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 $$\frac{1}{2}L{\log }_{2}L$$ where L is the block length. It then performs L point-wise multiplications for a total cost of $$L+L{\log }_{2}L=L(1+{\log }_{2}L)$$ multiplications. The cost ratio is therefore $$L(1+{\log }_{2}L)/(NL)=(1+{\log }_{2}L)/N$$ which is approximately log₂L / N.

Therefore, fftfilt is faster when log₂L is less than N.

Algorithms

fftfilt filters data using the efficient FFT-based method of overlap-add [1], a frequency domain filtering technique that works only for FIR filters by combining successive frequency domain filtered blocks of an input sequence.

Assume an input signal vector x with N_x elements and a filter vector b with N elements, where b = [b₁ b₂ ⋯ b_N]. The operation performed by fftfilt is described in the time domain by the difference equation:

$$y[n]=b_{\text{}1}\text{}x[n]+b_{\text{}2}\text{}x[n-1]+\cdots +b_{\text{}N}\text{}x[n-(N-1)]$$

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

$$Y(z)=\left[b_{\text{}\text{}1}+b_{\text{}2}\text{}{z}^{-1}+\cdots +b_{\text{}N}\text{}{z}^{-(N-1)}\right]X(z)$$

The fftfilt function uses fft to implement the overlap-add method. fftfilt breaks the input sequence x into k length-L data blocks, where L must be greater than the filter length N, where k = ⌈N_x/L⌉ and the ⌈⌉ symbols denote the ceiling function.

The function convolves each block with the filter b using

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

where i = 1, L+1, 2L+1, ⋯ and nfft is the FFT length. fftfilt overlaps successive output sections by N–1 points and sums them.

fftfilt chooses the key parameters L and nfft in different ways, depending on whether you specify an FFT length nfft for the filter and signal.

• If you do not specify a value for nfft (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).

    These assumptions yield y as follows:

    y = ifft(fft(b,nfft).*fft(x,nfft))
    

• If you specify a value for nfft, fftfilt chooses an FFT length of 2^nextpow2(nfft) and a data block length of nfft-length(b)+1. If nfft is less than length(b), then fftfilt chooses an FFT length of 2^nextpow(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.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

digitalFilter objects are not supported for code generation.

GPU Code Generation

Generate CUDA® code for NVIDIA® GPUs using GPU Coder™.

Usage notes and limitations:

digitalFilter objects are not supported for code generation.

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

This function fully supports GPU arrays. For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced before R2006a