On Nov 27, 9:24 am, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> On Nov 24, 9:45 am, Richard Owlett <rowl...@atlascomm.net> wrote:
>
> > In "Chapter 12 - The Fast Fourier Transform / FFT Programs" of _The
> > Scientist and Engineer's Guide to Digital Signal Processing_
> > Smith says "There are also FFT routines that completely eliminate the
> > bit reversal sorting."
>
> > Could someone point me to a description of one, preferably with BASIC or
> > FORTRAN sample code.
>
> A number of references are given in the Wikipedia article on the
> Cooley-Tukey FFT algorithm.
>
> http://en.wikipedia.org/wiki/Cooley-Tukey_FFT_algorithm
>
> The Cooley-Tukey FFT (the most common algorithm), along with other FFT
> algorithms for that matter, intrinsically requires some sort of re-
> ordering, in the sense that its operations map consecutive inputs to
> non-consecutive outputs and vice-versa.
>
> The most well-known approach is to do the re-ordering all at once,
> with a single bit-reversal pass (or, more generally, a mixed-base
> digit reversal). Another approach that still works in-place is to
> arrange the sequence of radices so that one can do a sequence of small
> transpositions interleaved with the transform (alternating DIT/DIF
> stages work well here). The simplest approaches work out-of-place.
> e.g. the Stockham algorithm goes back and forth between two arrays
> with each pass, transposing one digit at a time.
>
> The easiest approach, in my mind, is simply to do the algorithm
> recursively out-of-place; if you directly transcribe the Cooley-Tukey
> algorithm from the mathematical description, working recursively, all
> of the outputs are "automatically" put in the right place. (This was
> the approach we used in the original FFTW; the current version of FFTW
> implements a variety of algorithms as described in our Proc. IEEE
> paper.)
>
> For example, here is a very simple 15-line radix-2 FFT, working
> recursively and out-of-place, that does not require any separate bit-
> reversal pass. (To make it reasonably fast, one would want to use a
> larger base case than the size-1 base-case used here. It also uses a
> relatively inaccurate trigonometric recurrence. But of course, if you
> were actually interested in such practical issues you should be using
> a pre-existing FFT library.)
>
> #include <complex.h>
> #define TWOPI 6.2831853071795864769252867665590057683943388
> /* call with rec_fft(+1 or -1, N, input, output, 1, 1), where input
> and output are arrays of length N ...
> requires C99 complex-number support */
> void rec_fft(int sign, int N, complex double *x, complex double *y,
> int is, int os)
> {
> if (N == 1) *y = *x;
> else {
> int N2 = N/2;
> rec_fft(sign, N2, x, y, is*2, os);
> rec_fft(sign, N2, x+is, y+os*N2, is*2, os);
> complex double omega = 1.0, omega1 = cexp(sign * I * TWOPI /
> N);
> for (int i = 0; i < N2; ++i) {
> complex double a = y[os*i], b = omega * y[os*(i + N2)];
> y[os*i] = a + b;
> y[os*(i + N2)] = a - b;
> omega *= omega1;
> }
> }
>
> }
>
> Regards,
> Steven G. Johnson
Richard wanted the FFT routines in BASIC, so I converted
the above routine into Basic for Chipmunk Basic, changing
pointer arithmetic into array access, converting complex
data types back to standard float variables, and removing
the trig recursion. A Basic implementation that allow
recursive subroutines with locally scoped variables is
required.
I found that the line in the above C code:
> if (N == 1) *y = *x;
does pretty much the same thing as does the final sort of
a generic in-place decimation in frequency FFT algorithm,
except distributing the sort with the butterflys, and
calculating the bit-reversed indexes recursively.
rem *** recursive out-of-place radix-2 FFT in Basic ***
rem
rem *** call as: rec_fft(flag, n, xr,xi,0, yr,yi,0, 1,1)
rem flag = 1 for FFT, -1 for IFFT
rem n is FFT length
rem xr,xi are input arrays (must be dimensioned to length >= n)
rem yr,yi are result arrays (must be dimensiond to length >= n)
rem (kx,ky,ks,os) are sub-vector parameter state variables
rem (n2,i,c,s,k1,k2,ar,ai,br,bi) are local variables
rem
sub rec_fft(flag,n,xr,xi,kx,yr,yi,ky,ks,os,
n2,i,c,s,k1,k2,ar,ai,br,bi)
if (n = 1)
rem ** this does a bit-reversed-index copy and scale **
rem print ky,kx to see post (DiF) sorting
s = 1 : if (flag = -1) then s = 1/ks
yr(ky) = xr(kx) * s
yi(ky) = xi(kx) * s
else
n2 = n/2
rec_fft(flag,n2, xr,xi,kx , yr,yi,ky, ks*2,os)
rec_fft(flag,n2, xr,xi,kx+ks, yr,yi,ky+os*n2, ks*2,os)
for i = 0 to n2-1
c = cos(i * 2*pi/n)
s = sin(i * flag * 2*pi/n)
k1 = ky+os*(i)
k2 = ky+os*(i+n2)
ar = yr(k1)
ai = yi(k1)
br = c * yr(k2) - s * yi(k2)
bi = c * yi(k2) + s * yr(k2)
yr(k1) = ar + br
yi(k1) = ai + bi
yr(k2) = ar - br
yi(k2) = ai - bi
next i
endif
end sub : rem end rec_fft
If the above code formats wrong (old-fashioned Basic is
line-break sensitive), a copy is here:
http://www.nicholson.com/dsp.fft1.html
IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M