Forums

selective fft

Started by amara vati September 21, 2004
my application is to select and decode only a few sub-carriers in a
OFDM sytem. so it is neither continuous nor has any regularity.

amar

"Robert Lacoste" <see-www-alciom-com-for-email-adress> wrote in message news:<415284b9$0$15748$7a628cd7@news.club-internet.fr>...
> Dear Amar, > you can check what could give a "downconverter" approach : > > - Do a digital downconvertion of your band of interest to DC or to a low IF > (mixing your input signal with a DDS-generated digital sine wave with a > fixed frequency close to your band of interest, then low pass filter it and > decimate it with a CIC and/or FIR filter. This way you get a signal that > includes only the information you are looking for... > - Do a FFT on the downconverted signal, with a sample rate 10 times lower > than your input signal... > > For 10% I'm quite sure it will not be cost effective in terms of Mips with a > full FFT, however for very narrow analysis bandwidth it could be efficient, > no ? > > Cheers, > > Robert Lacoste > ALCIOM - The mixed signal experts > www.alciom.com > > > "amara vati" <amaraavati@yahoo.com> a &#2013265929;crit dans le message de > news:f89b870.0409202130.deef21f@posting.google.com... > > which is the best fft method if want the spectrum only for a selective > > set of points, say only 10-20% of the points in the entire spectrum. > > > > amar
Thanks Andor! - I didn't know that and now I'll look it up.
Best  of  Luck - Mike

"Andor Bariska" <an2or@nospam.net> wrote in message
news:4152b3bc@pfaff2.ethz.ch...
> Mike Yarwood wrote: > > "Andor Bariska" <an2or@nospam.net> wrote in message > > news:41517f90@pfaff2.ethz.ch... > > > >>amara vati wrote: > >> > >>>which is the best fft method if want the spectrum only for a selective > >>>set of points, say only 10-20% of the points in the entire spectrum. > >> > >>If the points are equally spaced on the unit circle, you can use the > >>chirp-z transform. > >> > >>Regards, > >>Andor > >> > > > > > > Don't you get the same thing just doing a shorter FFT? > > I don't understand your question: shorter than what? > > The N point DFT evaluates the spectrum at N evenly spaced points on the > whole unit circle. The chirp z transform can evaluate the spectrum at N > points evenly spaced on any subinterval of the unit circle (perhaps you > have heard of zoom FFT?). > > The N point chirp z transform can be implemented using two N point FFTs > plus an N point convolution (I think). > > Regards, > Andor >
Let me try to distill some of the wisdom in this thread and mix in a 
little foolishness of my own.


Q: How can I get just a few bins of a DFT?


A0: For a continous subset of frequency bins ...
  use the "zoom fft" technique
	0) mix (freq shift) the center of the band of interest to DC
	1) anti-alias filter and decimate
	2) optionally repeat step 1 multiple times
	3) perform an FFT at the lower sample rate

A1:  For a few erratically spaced bins ...
  it may be fastest to use a Goertzel filter for each bin.  This 
technique has O(n^2) efficiency so it quickly  becomes more expensive 
than just performing the entire FFT.


A2:  For a few equally spaced bins ...
  it sounds like chirp-Z is recommended (anyone wish to elaborate?).

A3: Also for a few equally spaced bins ... (here's my own bit of 
foolishness)
	0) mix the input signal to align the first bin of interest to DC
	1) block-accumulate the input signal into a smaller sequence
	2) perform the fft of the smaller sequence

	I believe the cost is less than Andor mentioned for a chirp-Z transform.


-- Mark


amara vati wrote:
> my application is to select and decode only a few sub-carriers in a > OFDM sytem. so it is neither continuous nor has any regularity. > > amar > > "Robert Lacoste" <see-www-alciom-com-for-email-adress> wrote in message news:<415284b9$0$15748$7a628cd7@news.club-internet.fr>... > >>Dear Amar, >>you can check what could give a "downconverter" approach : >> >>- Do a digital downconvertion of your band of interest to DC or to a low IF >>(mixing your input signal with a DDS-generated digital sine wave with a >>fixed frequency close to your band of interest, then low pass filter it and >>decimate it with a CIC and/or FIR filter. This way you get a signal that >>includes only the information you are looking for... >>- Do a FFT on the downconverted signal, with a sample rate 10 times lower >>than your input signal... >> >>For 10% I'm quite sure it will not be cost effective in terms of Mips with a >>full FFT, however for very narrow analysis bandwidth it could be efficient, >>no ? >> >>Cheers, >> >>Robert Lacoste >>ALCIOM - The mixed signal experts >>www.alciom.com >> >> >>"amara vati" <amaraavati@yahoo.com> a &#2013265929;crit dans le message de >>news:f89b870.0409202130.deef21f@posting.google.com... >> >>>which is the best fft method if want the spectrum only for a selective >>>set of points, say only 10-20% of the points in the entire spectrum. >>> >>>amar
On Sun, 26 Sep 2004 19:13:15 GMT, Mark Borgerding
<mark@borgerding.net> wrote:

>Let me try to distill some of the wisdom in this thread and mix in a >little foolishness of my own. > > >Q: How can I get just a few bins of a DFT? > > >A0: For a continous subset of frequency bins ... > use the "zoom fft" technique > 0) mix (freq shift) the center of the band of interest to DC > 1) anti-alias filter and decimate > 2) optionally repeat step 1 multiple times > 3) perform an FFT at the lower sample rate > >A1: For a few erratically spaced bins ... > it may be fastest to use a Goertzel filter for each bin. This >technique has O(n^2) efficiency so it quickly becomes more expensive >than just performing the entire FFT. > > >A2: For a few equally spaced bins ... > it sounds like chirp-Z is recommended (anyone wish to elaborate?). > >A3: Also for a few equally spaced bins ... (here's my own bit of >foolishness) > 0) mix the input signal to align the first bin of interest to DC > 1) block-accumulate the input signal into a smaller sequence > 2) perform the fft of the smaller sequence > > I believe the cost is less than Andor mentioned for a chirp-Z transform. > > >-- Mark >
Hi Mark, Perhaps the "Sliding DFT" might be useful to Amara. That DFT allows you to compute single-bin DFT samples. The nice part is that upon arrival of a new time sample, the new single-bin DFT sample can be computed with (roughly) a half dozen arithmetic operations. The best places for Amara to learn about the "Sliding DFT" are: "The Sliding DFT, an Update", IEEE Signal Processing Magazine, DSP Tips & Tricks column, Vol. 21, No. 1, Jan. 2004, by E.Jacobsen & R. Lyons "The Sliding DFT", IEEE Signal Processing Magazine, DSP Tips & Tricks column, Vol. 20, No. 2, March 2003, by E. Jacobsen & R. Lyons See Ya', [-Rick-]
By the way, an FFT where you only want a few outputs (or, conversely,
where only a few inputs are nonzero) is called a "pruned FFT" in the
literature.  See:

                  http://www.fftw.org/pruned.html

Mark Borgerding <mark@borgerding.net> wrote...>
>A0: For a continous subset of frequency bins ... > use the "zoom fft" technique
Instead, for an exact method (in exact arithmetic), see the URL above. Essentially, for a block of K outputs, where K divides N, you perform a single step of a radix-K decimation-in-frequency Cooley-Tukey FFT. I should also mention that there are also several approximate methods in the literature based on wavelets etc. to approximately compute a subset of the DFT outputs (Guo and Burrus, 1996; Shentov, 1995).
> A1: For a few erratically spaced bins ... > it may be fastest to use a Goertzel filter for each bin. This > technique has O(n^2) efficiency so it quickly becomes more expensive > than just performing the entire FFT.
Note that Goertzel's algorithm has *much* worse numerical errors than the FFT, or even compared to the naive DFT formula, so it is probably not appropriate for large N. (I think Mark meant to say that it has O(N K) complexity, where K is the number of outputs desired and N is the number of inputs, so that it approaches O(N^2) as K goes to N. Similarly for the naive DFT formula, of course.)
> A2: For a few equally spaced bins ... > it sounds like chirp-Z is recommended (anyone wish to elaborate?).
No, I believe the chirp-z algorithm (http://en.wikipedia.org/wiki/Bluestein's_FFT_algorithm) is generally more costly than computing the ordinary FFT of the whole data set, even for K < N outputs, assuming the size N is highly composite. chirp-z is mostly useful if you want a different set of frequency points in the complex plane than what the DFT gives you, or if your data size is prime. Cordially, Steven G. Johnson
On 27 Sep 2004 08:16:42 -0700, stevenj@alum.mit.edu (Steven G.
Johnson) wrote:

  (snipped)
> >Note that Goertzel's algorithm has *much* worse numerical errors than >the FFT, or even compared to the naive DFT formula, so it is probably >not appropriate for large N. >
(snipped)
> >Cordially, >Steven G. Johnson
Hi Steven, your post is interesting. I've modeled the Goertzel algorithm (using the floating-point precision of MATLAB) and have not seen the numerical errors that you mention. I'm wondering under what condition(s), or implementation(s), the Goertzel algorithm exhibits *much* worse numerical errors than FFT computations. Can I bug you tell us a little more? Thanks Steven. Regards, [-Rick-]
Rick Lyons wrote:

> ... DFT allows you to compute > single-bin DFT samples. The nice part is that > upon arrival of a new time sample, the new single-bin > DFT sample can be computed with (roughly) > a half dozen arithmetic operations. >
The question was posed for case of only a few frequencies being of interest. The answer got me musing. Having already calculated an m point DFT/FFT for points n thru n+m when point n+m+1 arrives is there a simple algorithm for calculating DFT/FFT of points n+1 thru n+m+1?
r.lyons@_BOGUS_ieee.org (Rick Lyons) wrote...
> your post is interesting. I've modeled the > Goertzel algorithm (using the floating-point > precision of MATLAB) and have not seen the > numerical errors that you mention. > > I'm wondering under what condition(s), or > implementation(s), the Goertzel algorithm exhibits > *much* worse numerical errors than FFT computations. > Can I bug you tell us a little more?
See: W. M. Gentleman, "An error analysis of Goertzel's (Watt's) method for computing Fourier coefficients," Computer Journal 12 (2), 160-165 (1969). Abstract: Goertzel's method is commonly used when only a small number of coefficients is desired for a given sequence. The paper gives a floating-point error analysis of the technique, and shows why it should be avoided, particularly for low frequencies. His error analysis indicates that Goertzel's method leads to floating point errors that have an upper bound O(N^(5/2)), with the worst case being low frequencies. Compare this to an error bound of O(log N) for the FFT, and mean errors of O(sqrt(log N)). It's perhaps not that surprising -- Goertzel's algorithm does not explicitly compute N twiddle factors exp(ikn/N), but instead contains something very similar to a trigonometric recurrence formula...and the recurrence that it most resembles has horrible O(N^2) relative error growth (both worst case *and* mean). Here are some sample numbers I got from a C implementation of Goertzel's algorithm, compiled with gcc on a Pentium 4 (double precision; code is structured to *not* stay within x86 extended-precision registers). Computing the first 10 bins of a size-N DFT with Goertzel and calculating the maximum relative error w.r.t. an FFT (defined as max |difference| / max |value|) for uniform random inputs, I get: N, relative error 100, 1.2e-14 1000, 9.3e-13 10000, 1.4e-10 100000, 3.3e-8 1000000, 7.7e-7 10000000, 7.1e-5 which is almost exactly N^2 growth in the error. In contrast, the errors accumulated by the naive DFT formula in a similar experiment grow only as ~sqrt(N) as expected, and are still only 1.2e-13 for N=10^7 (but are still worse than the FFT's O(sqrt(log N))). I'd be interested to hear if there were some trick to make Goertzel's algorithm more accurate without sacrificing performance (or if your Matlab code achieves better errors for a similar experiment to the above). Of course, if you use e.g. extended-precision registers on the Pentium, matters will improve somewhat (by a constant factor). Cordially, Steven G. Johnson
Richard Owlett wrote:
> Rick Lyons wrote: > >> ... DFT allows you to compute >> single-bin DFT samples. The nice part is that upon arrival of a new >> time sample, the new single-bin DFT sample can be computed with >> (roughly) a half dozen arithmetic operations. >> > > The question was posed for case of only a few frequencies being of > interest. > > The answer got me musing. > > Having already calculated an m point DFT/FFT for points n thru n+m > > when point n+m+1 arrives > > is there a simple algorithm for calculating DFT/FFT of points n+1 thru > n+m+1? > >
Instead of musing I should have looked in Rick's book (sect 13.18) ;]
Steven,
see at the bottom


"Steven G. Johnson" <stevenj@alum.mit.edu> wrote in message
news:27cfb406.0409281803.69328485@posting.google.com...
> r.lyons@_BOGUS_ieee.org (Rick Lyons) wrote... > > your post is interesting. I've modeled the > > Goertzel algorithm (using the floating-point > > precision of MATLAB) and have not seen the > > numerical errors that you mention. > > > > I'm wondering under what condition(s), or > > implementation(s), the Goertzel algorithm exhibits > > *much* worse numerical errors than FFT computations. > > Can I bug you tell us a little more? > > See: > > W. M. Gentleman, "An error analysis of Goertzel's (Watt's) method > for computing Fourier coefficients," Computer Journal 12 (2), > 160-165 (1969). > > Abstract: Goertzel's method is commonly used when only a small > number of coefficients is desired for a given sequence. The > paper gives a floating-point error analysis of the technique, > and shows why it should be avoided, > particularly for low frequencies. > > His error analysis indicates that Goertzel's method leads to floating > point errors that have an upper bound O(N^(5/2)), with the worst case > being low frequencies. Compare this to an error bound of O(log N) for > the FFT, and mean errors of O(sqrt(log N)). > > It's perhaps not that surprising -- Goertzel's algorithm does not > explicitly compute N twiddle factors exp(ikn/N), but instead contains > something very similar to a trigonometric recurrence formula...and the > recurrence that it most resembles has horrible O(N^2) relative error > growth (both worst case *and* mean). > > Here are some sample numbers I got from a C implementation of > Goertzel's algorithm, compiled with gcc on a Pentium 4 (double > precision; code is structured to *not* stay within x86 > extended-precision registers). Computing the first 10 bins of a > size-N DFT with Goertzel and calculating the maximum relative error > w.r.t. an FFT (defined as max |difference| / max |value|) for uniform > random inputs, I get: > > N, relative error > 100, 1.2e-14 > 1000, 9.3e-13 > 10000, 1.4e-10 > 100000, 3.3e-8 > 1000000, 7.7e-7 > 10000000, 7.1e-5 > > which is almost exactly N^2 growth in the error. In contrast, the > errors accumulated by the naive DFT formula in a similar experiment > grow only as ~sqrt(N) as expected, and are still only 1.2e-13 for > N=10^7 (but are still worse than the FFT's O(sqrt(log N))). > > I'd be interested to hear if there were some trick to make Goertzel's > algorithm more accurate without sacrificing performance (or if your > Matlab code achieves better errors for a similar experiment to the > above). Of course, if you use e.g. extended-precision registers on > the Pentium, matters will improve somewhat (by a constant factor). > > Cordially, > Steven G. Johnson
Hello Steven, A trick I use for low frequencies that is Goertzel like uses a near quadrature oscillator. The "c" like code is below: It does cost two mults per iter, but like the Goertzel algo, the twiddle factors are also computed on the fly. The advantage of this method is when viewed in matrix form, the condition number for the matrix approaches 1 for low frequencies. Goerzel uses a biquad oscillator which needs to step 90 degrees per sample to achieve the low condition number. The input data is in x[], the data has N samples. And the bin number is w. y1=0; y2=0; k=2*sin(pi*w/N); // not 2 pi !! for (j=0;j<N;j++) { y2=y2-k*y1+x[j]; y1=y1+k*y2; } And the energy is simply E = y1*y1 + y2*y2 - k*y1*y2; I would be interested in how this algo fairs in your comparison tests. Clay S. Turner