7
$\begingroup$

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 ?

$\endgroup$
  • $\begingroup$ there's a lot to your question, and i don't have time to go through it all. no, you need not overlap-add if the size of your data fits nicely into the FFT. you might want to window and zero-pad it. maybe not. remember that the FFT essentially periodically extends whatever data set passed to it. so you might need to worry about the edges. $\endgroup$ – robert bristow-johnson Jan 17 '14 at 15:41
  • 6
    $\begingroup$ Unless you need to squeak out every drop of performance that you can, I would highly recommend not filtering using FFT's. Filter in the time domain. For shorter filters this is actually more computationally efficient than filtering with FFT's anyways. $\endgroup$ – Jim Clay Jan 17 '14 at 16:03
  • $\begingroup$ Thx Robert : What do you mean by "the size of the data fits in the FFT". For now, I choose let's say 2048 points of data and perform an FFT, as my data is Real I obtain an Hermitian symetric FFT of 1024 usefull complex points. What must I window or pad ? $\endgroup$ – Arnaud Jan 18 '14 at 11:19
  • $\begingroup$ Thx Jim : Do you have any hints on how to filter in the time domain ? (something I haven't wrapped my mind around at the moment) Any good links on the subject ? $\endgroup$ – Arnaud Jan 18 '14 at 11:21
  • $\begingroup$ @Arnaud, i cannot think of a simpler way to put it. if your 2048 data points fit into the 2048-point DFT, there is no issue of breaking up the data to make it fit and overlapping the results of multiple uses of the DFT. now, if you are putting 2048 non-zero samples into a 2048-point DFT, you must keep in mind the discontinuity that exists between $x[2047]$ and $x[0]$ because of the inherent periodic extension done by the DFT. FYI, the "FFT" is an efficient method to compute the DFT. all of the mathematical properties of the DFT apply. $\endgroup$ – robert bristow-johnson Jan 19 '14 at 17:24
1
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ why not include the graphic in your answer? it's not hard to do. $\endgroup$ – robert bristow-johnson Sep 14 at 4:49
  • 1
    $\begingroup$ BTW, there are more terms in the cost function to worry about than just the number of complex multiplications. this is a paper where i have generalized the cost function for fast-convolution. i didn't know this question existed. i might include this paper in an answer. $\endgroup$ – robert bristow-johnson Sep 14 at 5:03
  • $\begingroup$ I did not know about embedding the graphic here. But I'd have still used the URL, because that web page also has the Octave source code, in case anyone is interested. $\endgroup$ – Bob K Sep 14 at 10:35
  • $\begingroup$ again, there are more terms to the cost function than just the $N \log_2(N)$ term. there is another salient term, $N$, with it's own cost coefficient. may i modify your answer to embed the graphic in your answer? $\endgroup$ – robert bristow-johnson Sep 14 at 14:33
  • $\begingroup$ Thanks for asking, but no thank you. Please just post your own answer. I am not the author of the formula in question. My role was to turn it into a more useful graphic. I have no opinion about "their" formula vs yours, but it has withstood the scrutiny of other Wikipedians for years. And some of the finer points are going to be implementation-dependent... not really worth quibbling about. Just make your own case, and let the "voters" decide. $\endgroup$ – Bob K Sep 14 at 23:13
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thx hotpaw2 : I will perform only one FFT at a time. WE have a track of data, we will choose a segment of N (1024/2048/4096 pts) data and perform FFT, Filtering, IFFT. You say that the size of the FFT is >= size of data+filter impulse response. I don't understand, which FFT ? My data one, my filter one, I'm lost on this one. The resulting FFT ? $\endgroup$ – Arnaud Jan 18 '14 at 11:31
  • $\begingroup$ See the first part of my answer. Say with a N=4096 FFT you can use a spacing of the segments of L=2048, an overlap of D=1024, that is, in each FFT you take L+D=3072 samples from the signal, and a filter length of F=1024. $\endgroup$ – Lutz Lehmann Jan 18 '14 at 11:36
  • $\begingroup$ @Arnaud: For which FFT size, please read about and study the overlap-add or overlap-save fast convolution algorithms. There are too many important small details to post here. $\endgroup$ – hotpaw2 Jan 18 '14 at 17:08
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thx LutzL : I will look into wavelets transforms but considering the time already spent on FFT I have almost reached my goal, I just need a couple of last advice on how to properly apply my filter FFT on my data FFT. I have already found information about wavelets tranform on the net but wasn't sure if it's the good way to go for my "low frequency elimination problem". $\endgroup$ – Arnaud Jan 18 '14 at 11:26
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.