selective fft

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 &#4294967295;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 &#4294967295;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