Filter size vs. FFT size and Overlap add

Question

First of all I would like to say that :

1. I'm new to DSP and in a "learning curve" hence the "basicness" of my questions.
2. I've read as many posts/articles as possible on the internet but still need "guidance".
3. I'm french hence the bad formulation sometimes.

My mission/data context :

Performing data analysis on real data sets (oil-rig data). Sometimes data collected on the rig present low frequency component due to rig movements. Data are collected every 5 seconds so the sampling rate is 0.2 Hz. The size of the data to be analysed will be 1024/2048/4096 points let's call it N is this post. My mission is to perform a filtering on these data to eliminate the low frequency component. This is not a real time application. The range of data to be analysed will be chosen and then the analysis performed.

Chosen approach :

Process these data with a FFT, filtering the FFT and perform an IFFT to obtain and clean(er) set of data. Basically the lowest frequency found in the FFT is to be eliminated. Many users suggested me to use the overlap-add method to obtain the desired result. First I was thinking of just zeroing low frequency related bins in my FFT but my readings on the subject suggested to rather use a high pass filtering (with spectral inverted windowed-sinc).

Designing this filter and applying it is what I post here for because it turned out to be way more difficult than expected.

Questions :

Do I have to use the overlap-add method considering that my application is not a real time one ?

For what I have understood I need to design a high pass filter kernel with a spectral inverted sinc function. Then I have to perform a FFT on my filter kernel and multiply points by points my filter FFT(1) coefficients with my data FFT(2) coefficient before performing a IFFT on the resulting FFT. Am I correct on this one ?

1. According to my readings I can obtain a spectral inverted sinc function by inverting (-sin(x)/x) it and adding 1 to the center of symmetry. I am right, do you guys confirm ? Do I add 1 on the one center point only ?

2. What is the right size to choose for my filter kernel ? Is it the same size N as my data ? As I have to multiply my FFT(1) with my FFT(2) I would think that I have to choose the same size but I'm probably wrong on this one and the overlap add method may be the answer.

3. The kernel filter designing is tricky, I don't know where to place the central lobe of my spectral inverted sinc, in 0, on a specific offset in the time domain ?

Answer 1

For your non-real-time application, you do not need OS or OA. Direct form convolution (such as Octave function conv()) is easy and gives the same result.

But since the title of your post is "Filter size vs. FFT size and Overlap add", this graphic directly addresses that issue.

Answer 2

Choosing the FFT size for optimally efficient fast convolution

The result of linear convolution

$$ y[n] = \sum\limits_{i=0}^{L-1} h[i] x[n-i] $$

of a block of input, $x[n]$, of length $B$

$$ x[n] = 0 \qquad \text{for } n < 0 \text{ or } n \ge B $$

with a finite impulse response (FIR), $h[n]$, of length $L$

$$ h[n] = 0 \qquad \text{for } n < 0 \text{ or } n \ge L $$

has a non-zero length $B+L−1$.

$$ y[n] = 0 \qquad \text{for } n < 0 \text{ or } n \ge B+L-1 $$

Whether done by overlap-add or overlap-save (sometimes referred to as “overlap scrap”) when using the FFT and circular convolution to implement an FIR on a block of samples, the FFT length, $N$, is the period of the periodic result.

$$ y[n+N] = y[n] \qquad \forall n \in \mathbb{Z} $$

To prevent overlapping and time-aliasing in $y[n]$, this period, $N$, must be at least as long as the non-zero length of the FIR result. The block length, $B$, FIR length, $L$, and FFT size, $N$, are constrained by

$$ N \ge B+L−1 \ .$$

Given an FIR length, $L$, and an FFT size, $N$, the most samples we can hope to process per block is $B=N−L+1$. Performing fast convolution using a radix-2 FFT, the computational cost, in instruction execution time per block, is

$$ C = C_1 N\log_2(N) + C_2 N + C_3 \log_2(N) + C_4 $$

for constants $C_1$, $C_2$, $C_3$, and $C_4$ to be determined from analysis or profiling of the code performing the FFT and the other tasks of the circular convolution.

$C_1$ is the cost of performing an FFT butterfly. $C_2$ is the cost of all operations performed on each FFT bin such as multiplication of $H[k]$ to $X[k]$. $C_3$ is the overhead cost of setting up each FFT pass and $C_4$ is the total constant overhead of setting up each block or frame. It is likely that $C_3$ and $C_4$ can be ignored since they multiply much smaller numbers than $C_1$ and $C_2$.

This circular convolution operation will have to be applied every $B$ samples and the optimal choice of $N$ is such that minimizes that cost per sample.

So that a single simple implementation of the FFT can be used for a variety of FFT sizes, only radix-2 FFTs will be considered and $N$ must be a power of $2$.

$$ N = 2^p $$ $$ p = \log_2(N) $$

The cost per sample is

$$ \frac{C}{B} = \frac{C_1 N\log_2(N) + C_2 N + C_3 \log_2(N) + C_4}{N-L+1} $$

or

$$ \frac{C}{B} = \frac{C_1 2^p p + C_2 2^p + C_3 p + C_4}{2^p-L+1} $$

Given an FIR length, $L$, to optimize the FFT size, $N=2^p$ , we choose $N$ to be such a power of $2$ so to minimize the per-sample cost $\frac{C}{B}$. To do that we find the minimum value of $N=2^p$ or minimum integer $p$ so that

$$ \frac{C_1 2^{p+1} (p+1) + C_2 2^{p+1} + C_3 (p+1) + C_4}{2^{p+1}-L+1} \ > \ \frac{C_1 2^p p + C_2 2^p + C_3 p + C_4}{2^p-L+1} $$

or

$$ \big( C_1 2^{p+1} (p+1) + C_2 2^{p+1} + C_3 (p+1) + C_4 \big) \big(2^p-L+1\big) \ > \ \big( C_1 2^p p + C_2 2^p + C_3 p + C_4 \big) \big( 2^{p+1}-L+1 \big) $$

or

$$ \big( 2 C_1 2^p (p+1) + 2 C_2 2^p + C_3 (p+1) + C_4 \big) \big(2^p-(L-1)\big) \\ > \ \big( C_1 2^p p + C_2 2^p + C_3 p + C_4 \big) \big( 2 \cdot 2^p-(L-1) \big) $$

or

$$ 2^p\big( 2 C_1 2^p (p+1) + 2 C_2 2^p + C_3 (p+1) + C_4 \big) - 2 \cdot 2^p\big( C_1 2^p p + C_2 2^p + C_3 p + C_4 \big) \\ > \ (L-1)\big( 2 C_1 2^p (p+1) + 2 C_2 2^p + C_3 (p+1) + C_4 \big) - (L-1)\big( C_1 2^p p + C_2 2^p + C_3 p + C_4 \big) $$

or

$$ 2^p\big( 2 C_1 2^p - C_3 p + C_3 - C_4 \big) \ > \ (L-1)\big( C_1 2^p p + 2 C_1 2^p + C_2 2^p + C_3 \big) $$

or

$$ 2 C_1 2^p - C_3 p + C_3 - C_4 \ > \ (L-1)\big( C_1 p + 2 C_1 + C_2 + C_3 2^{-p} \big) $$

Simplification results in

$$ L-1 \ < \ \frac{2 C_1 2^p - C_3 (p-1) - C_4}{C_1 (p+2) + C_2 + C_3 2^{-p}} $$

or

$$ L-1 \ < \ \frac{2 C_1 N - C_3 (\log_2(N)-1) - C_4}{C_1 (\log_2(N)+2) + C_2 + \frac{C_3}{N}} $$

or

$$ L \ < \ 1 + \frac{2 C_1 N}{C_1 (\log_2(N)+2) + C_2 + \frac{C_3}{N}} - \frac{C_3 (\log_2(N)-1) + C_4}{C_1 (\log_2(N)+2) + C_2 + \frac{C_3}{N}} $$

Then, given an FFT size that is a power of two, $N=2^p$ , we can see what is the largest value of FIR length, $L$, for which that FFT size is optimally efficient. If the terms with $C_3$ and $C_4$ are ignored (because they multiply numbers that are much smaller than what $C_1$ and $C_2$ multiply), then this comes out approximately as

$$ L \ < \ \frac{2 N}{ \log_2(N)+2 + \frac{C_2}{C_1} } . $$

I believe that the cost factor $C_2$ (essentially the cost of a complex multiplication) will always be less than $C_1$ (the cost of an FFT butterfly). Examining this result indicates that, say for $C_1=C_2$ and $N=256$, the largest optimal $L$ is $46$ and the corresponding block length $B$ is $211$. Choosing a smaller block or buffer length for that size of FFT is not as computationally efficient.

Answer 3

If you have a data sequence that is larger than the size of data window used in each FFT, then you will need to use an algorithm such as overlap-add or save. The size of the FFT should be equal to or larger than the sum of the data window plus the impulse response of your filter (your windowed Sinc). The size of your windowed Sinc will be related to your desired filter transition width and frequency.

Answer 4

To do the computation of overlap-add: You divide the signal into segments starting at indices $kL$ using an overlap of size $D$, that is, between kL and kL+D there is a cross-fade adding to 1 of the segment $(k-1)$ and segment $k$. Applying a filter of finite length is equivalent to polynomial multiplication. The windowed $k$th segment of length $L+D$ and the filter of length $F$ result in a filtered segment of length $L+D+F-1$. This is the size that has to fit in the dyadic FFT-length $N$. After zeropadding, FFT, multiplication with the FFT of the filter and iFFT the filtered segment is added in its full length $N$ to the output signal starting at positon $kL$.

---

The combination of filtering and downsampling can be factored into steps in the so-called lifting scheme. For short filters like those of the discrete wavelet transform, this is faster than overlap-add. Using the idea of the Remez-algorithm, relatively short wavelet transforms with good frequency separation can be constructed.