FFT algoritms
in [DSP]
|
Prev: Scilab and DSP
Next: Separating Harmonic Spectra
From: Richard Owlett on 24 Nov 2007 15:54 Martin Blume wrote: > "Martin Blume" schrieb > >>>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. >>> > > Another possible site is: > http://mathworld.wolfram.com/FastFourierTransform.html > (didn't see any mention at first glance of non-bit-reversal stuff, > but you might follow the links [*]) Links were useful. Lead me to Hartley transform which matches my problem domain better as it has only REAL inputs -- in both senses of word. It is also explained with sines and cosines which I'm more comfortable with than with e. j. theta Google then turned up _An Algorithm for the Fast Hartley Transform_ , by Ronald F. Ullmann at http://sepwww.stanford.edu/oldreports/sep38/38_29_abs.html which includes code in RATFOR. I can read FORTRAN but no C or its variants ;/ > > Regards > Martin > > [*] beware of the danger mentioned in http://xkcd.com/214/ Including fact that "Perfection" has a 58.8% relevance to fast Hartley Transform ;) see page 2 of http://en.wikipedia.org/w/index.php?title=Special:Search&search=fast+Hartley+Transform
From: Richard Owlett on 24 Nov 2007 16:18 Ron N. wrote: > On Nov 24, 7:14 am, "Martin Blume" <mbl...(a)socha.net> wrote: > >>"Richard Owlett" schrieb> 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." > > > Why do you need to avoid sorting? Cause confuseth me it doeth ;) If I actually am forced to implement my own FFT (or more likely FHT) it will likely be in Forth. The bit reversal of the address become trivial ( an XOR ). > [snip]
From: Richard Owlett on 24 Nov 2007 16:21 Ron N. wrote: >>>"Richard Owlett" schrieb >> >>I'm not interested in the computation time but in whether or not I can >>follow what's going on. Eventually I want to experiment with using >>square waves instead of sinusoids as my basis functions. I suspect it >>will be a better model for my situation. > > > For experimenting with other basis vectors, you probably > want to stick with the N^2 DFT algorithm. The N*log(N) FFT > algorithm depends on the fact that sinusoidal basis vectors > have periodicity, symmetry and factorization properties. > > IMHO. YMMV. > -- > rhn A.T nicholson d.0.t C-o-M I don't see where square waves lack periodicity and/or symmetry where sine or cosine are used in code.
From: Martin Blume on 24 Nov 2007 17:51 "Richard Owlett" schrieb > > > > I'm not interested in the computation time but in whether > or not I can follow what's going on. Eventually I want to > experiment with using square waves instead of sinusoids > as my basis functions. > Why didn't you ask for it in the first place? You might want to look at http://mathworld.wolfram.com/WalshFunction.html > > I suspect it will be a better model for my situation. > Describe your sitiuation and somebody more knowledgeable than I might help. Regards Martin
From: Ron N. on 24 Nov 2007 19:18
On Nov 24, 1:18 pm, Richard Owlett <rowl...(a)atlascomm.net> wrote: > Ron N. wrote: > > On Nov 24, 7:14 am, "Martin Blume" <mbl...(a)socha.net> wrote: > > >>"Richard Owlett" schrieb> 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." > > > Why do you need to avoid sorting? > > Cause confuseth me it doeth ;) > If I actually am forced to implement my own FFT (or more likely FHT) it > will likely be in Forth. The bit reversal of the address become trivial > ( an XOR ). An XOR does a masked bit invert, not a bit reverse. To reverse the m-bit variable k in Basic : ki = k : rem temporary shift register kr = 0 : rem initialize result variable for i=1 to m kr = kr * 2 : rem shift left if (ki mod 2) = 1 then kr = kr + 1 ki = int(ki / 2) : rem shift right next i rem kr is k bit reversed m bits print k, kr If your Basic doesn't have the mod (modulo) operator, then (x mod y) = x - y * int(x / y) Code optimizations left as an exercise for the student. Chipmunk Basic does have a mod (modulo) operator: http://www.nicholson.com/rhn/basic IMHO. YMMV. -- rhn A.T nicholson d.0.t C-o-M |