FFT-based FIR filtering
FFT-based FIR filtering using the overlap-add method.
fftfilt(b, x, n = NULL) ## Default S3 method: fftfilt(b, x, n = NULL) ## S3 method for class 'Ma' fftfilt(b, x, n = NULL)
b: moving average (Ma) coefficients of a FIR filter, specified as a vector.
x: the input signal to be filtered. If x is a matrix, its columns are filtered.
n: FFT length, specified as a positive integer. The FFT size must be an even power of 2 and must be greater than or equal to the length of filt. If the specified n does not meet these criteria, it is automatically adjusted to the nearest value that does. If n = NULL
(default), then the overlap-add method is not used.
The filtered signal, returned as a vector or matrix with the same dimensions as x.
This function combines two important techniques to speed up filtering of long signals, the overlap-add method, and FFT convolution. The overlap-add method is used to break long signals into smaller segments for easier processing or preventing memory problems. FFT convolution uses the overlap-add method together with the Fast Fourier Transform, allowing signals to be convolved by multiplying their frequency spectra. For filter kernels longer than about 64 points, FFT convolution is faster than standard convolution, while producing exactly the same result.
The overlap-add technique works as follows. When an N length signal is convolved with a filter kernel of length M, the output signal is N + M - 1 samples long, i.e., the signal is expanded 'to the right'. The signal is then broken into k smaller segments, and the convolution of each segment with the f kernel will have a result of length N / k + M -1. The individual segments are then added together. The rightmost M - 1 samples overlap with the leftmost M - 1 samples of the next segment. The overlap-add method produces exactly the same output signal as direct convolution.
FFT convolution uses the principle that multiplication in the frequency domain corresponds to convolution in the time domain. The input signal is transformed into the frequency domain using the FFT, multiplied by the frequency response of the filter, and then transformed back into the time domain using the inverse FFT. With FFT convolution, the filter kernel can be made very long, with very little penalty in execution time.
t <- seq(0, 1, len = 10000) # 1 second sample x <- sin(2* pi * t * 2.3) + 0.25 * rnorm(length(t)) # 2.3 Hz sinusoid+noise filt <- rep(0.1, 10) # filter kernel y1 <- filter(filt, 1, x) # use normal convolution y2 <- fftfilt(filt, x) # FFT convolution plot(t, x, type = "l") lines(t, y1, col = "red") lines(t, y2, col = "blue") ## use 'filter' with different classes t <- seq(0, 1, len = 10000) # 1 second sample x <- sin(2* pi * t * 2.3) + 0.25 * rnorm(length(t)) # 2.3 Hz sinusoid+noise ma <- Ma(rep(0.1, 10)) # filter kernel y1 <- filter(ma, x) # convulution filter y2 <- fftfilt(ma, x) # FFT filter all.equal(y1, y2) # same result
https://en.wikipedia.org/wiki/Overlap-add_method.
filter
Kurt Hornik, Kurt.Hornik@wu-wien.ac.at ,
adapted by John W. Eaton.
Conversion to R by Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com .