15

I've got an algorithm that zero pads a sequence to 4N, does an FFT, and only uses the lowest frequency N points out of the generated 4N.

This seems like a lot of wasted work, any ideas how this can be done faster?

Mark Borgerding
  • 3,020
  • 19
  • 26
  • @Dilip. I'll be using FFTW or IMKL libraries. I could of course use my kissfft library, but it's starting out at a speed disadvantage vs the others – Mark Borgerding Nov 09 '11 at 17:30
  • 2
    I deleted the comment to which you responded since I meant to say decimation-in-frequency but wrote decimation-in-time instead. But look at the butterfly diagram here. If you write some code for the first two stages for the $4N$-FFT to take into account the large number of zeroes and skip the corresponding multiplications, you can then call the FFT library subroutine $4$ times for $N$-FFTs in which the input vectors are "full". Of course, you need only $N/4$ of the outputs from each subroutine call. – Dilip Sarwate Nov 09 '11 at 17:58

3 Answers3

2

If you only a few bins then the following may be very efficient for you:
1. Simply do the DFT of at each frequency you need.
2. Use the Goertzel algorithm for each frequency in question.

Jacob
  • 361
  • 2
  • 12
  • Mark said he needed $N$ bins out of $4N$, so 1) seems not to be a reasonable option. The Goertzel algorithm has advantages such as on-line computation as the data are received, small storage, etc, but needs $2N+4$ multiplications per bin, while each bin computed as a polynomial evaluation via Horner's rule needs only $N$ multiplications. Thus, 2) does not appear to be a particularly reasonable option either. – Dilip Sarwate Nov 10 '11 at 23:01
  • You're right, while reading the question I somehow missed the details. While I was answering I was thinking, "Gee, it would be nice to know how many bins he wants..." Guess I should re-read the question before answering. – Jacob Nov 11 '11 at 15:29
2

Zero padding to 4X length, computing the longer FFT, and then using only the bottom 1/4th bins produces almost identical results to windowed Sinc interpolation of the original length FFT.

So just use the original FFT length and interpolate using a 3 phase Sinc interpolation kernel with a suitable window width.

hotpaw2
  • 35,346
  • 9
  • 47
  • 90
0

Zero padding in the time domain gives you higher frequency solution but no new information, so it provides essentially interpolating in the frequency domain. Depending on the nature of your signals and the precision required you may be able to get the additional frequency points with a regular FFT of N points and doing a suitable interpolation (linear, spline, pchip, sinc etc).

Hilmar
  • 9,472
  • 28
  • 35
  • Let $x(z)=\sum_{i=0}^{N-1}x_iz^{i}$ be a polynomial (possibly with complex coefficients $x_i$) of degree $N-1$. We evaluate it at the $N$ points $\alpha^n, 0\leq n\leq N-1$ where $\alpha=\exp(-j2\pi/N)$ is a $N$-th root of unity to get $N$ numbers $X_n=x(\alpha^n)$. These are values of $x(z)$ at $N$ equally spaced points on the unit circle. What we really want is the values of $x(z)$ at $\beta^n, 0\leq n\leq N-1$ where $\beta=\exp(-j2\pi/4N)$, which are $N$ points on the first quadrant of the unit circle. I don't see how linear, spline etc interpolation is going to work. Please explain. – Dilip Sarwate Nov 10 '11 at 19:31
  • Sorry, that penultimate sentence in my previous comment should have said fourth quadrant of the unit circle. Since $\beta^4 = \alpha$, every fourth desired value $x(\beta^{4k})$ has already been computed by the FFT: $x(\beta^{4k})=x(\alpha^k)$. – Dilip Sarwate Nov 10 '11 at 20:07
  • I suspect it would be difficult to do a decent interpolation faster than doing the larger FFT. – Mark Borgerding Nov 10 '11 at 20:30
  • Let's say you have a 128 point FFT and 12800Hz sample rate. A 128 point FFTs gives values at 0Hz, 100Hz, 200 Hz, 300Hz, etc. What the zero padding does is to increase the frequency resolution to 0 Hz, 25Hz,50 Hz, 100Hz etc. This can be viewed as an interpolation problem. To me mathematically precised you need to do circular sinc intperpolation of 128th order. That certainly isn't worth the bother but depending on application and precision required a much lower order interpolation would be good enough – Hilmar Nov 11 '11 at 14:32