FFT algoritms

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

Ron N. wrote:

> On Nov 24, 7:14 am, "Martin Blume" <mbl...@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]

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.

"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

On Nov 24, 1:18 pm, Richard Owlett <rowl...@atlascomm.net> wrote:
> Ron N. wrote:
> > On Nov 24, 7:14 am, "Martin Blume" <mbl...@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

On Nov 24, 2:13 pm, Richard Owlett <rowl...@atlascomm.net> wrote:
>
> > PS: My ancient dead-tree version (from 1987) mentions other
> > algorithms (other than bit-reversal) briefly (2 pages), but
> > says that these are used for convolution using the FFT and
> > that not doing the bit-reversal step doesn't save much time.
>
> 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.

there are drawings of it in Oppenheim and Schafer (any era of O&S) for
the case of N=8.  i think i can see a pattern in them but, they do not
appear to me to be as efficient, for most instruction sets, as the
simple FFT where you have to bit reverse before or after.  you need to
maintain many more pointers and i'm not completely sure how you
manipulate them all for the different FFT passes.

for some reason, way back when, when i wrote FFT routines in different
contexts, i remember coming to the conclusion that the FFT form that
has normal order in and out, was more of a mess than the regular FFT
where bit reversing is sometimes needed.  it was not needed in the
case of fast-convolution, since i had a different iFFT from the FFT
where one was decimation-in-time form and the other was decimation-in-
frequency.  the FFT had normal order in, bit-reversed out, and the
iFFT had bit-reversed in, normal out.  the transfer function could be
pre-bit-reversed in order, so multiplying by it can be done with
normal order indexing in memory (which is cheap in native C, some DSPs
have built in bit reversing that is as cheap as incrementing in normal
order).

anyway, i'm still not sure that the normal-in, normal-out FFT is
useful for anything.  might be.

r b-j

Richard Owlett wrote:

(snip)

> 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.

Square waves are not orthogonal, which makes the transform
harder, but not impossible.

Still, you might look at the ICT.  It is a transform specifically
designed to make the forward transform fast on a machine with no
multiply instructions.   You might find from there some of the
properties of transforms on non-orthogonal basis functions.

The ICT assumes that the inverse transform is done on a much
faster processor (on earth).  The forward case can be done
fast on a CDP1802.

-- glen

Martin Blume 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. 
>>
> 
> Why didn't you ask for it in the first place?

My first goal was to "grok" FFT.
*THEN* I would go on to investigating if square waves were appropriate 
to a particular problem. My prime goal is understanding FFT. A secondary 
goal is gaining skill in programming a solution. The way my mind works, 
I have to pose a "problem" to work on. In this case it might be stated 
as "Why does everybody use sinusoids to analyze the output of a system 
excited by a two state source (larynx)?"


> 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.

Well the more they know about the "problem", the more likely I'll be 
told "don't do that" ;)
As I said above, the "problem" is somewhat a straw man. It's purpose is 
to serve as a framework to explore something I want to learn. Who knows, 
there may be useful fallout. Consider vulcanization and peneccillin.

glen herrmannsfeldt wrote:

> Richard Owlett wrote:
> 
> (snip)
> 
>> 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.
> 
> 
> Square waves are not orthogonal, which makes the transform
> harder, but not impossible.

Why would they be any less orthogonal than sets of sines and cosines?
They should also span the space as well as sines and cosines.

> 
> Still, you might look at the ICT.  It is a transform specifically
> designed to make the forward transform fast on a machine with no
> multiply instructions.   You might find from there some of the
> properties of transforms on non-orthogonal basis functions.
> 
> The ICT assumes that the inverse transform is done on a much
> faster processor (on earth).  The forward case can be done
> fast on a CDP1802.
> 

By ICT do you mean Inverse Cosine Transform? That was the only relevant 
Google hit I found. My signal is real but not necessarily even (or odd).

I'd be interested in anything that eliminated multiplications.

On Nov 25, 4:30 am, Richard Owlett <rowl...@atlascomm.net> wrote:
> glen herrmannsfeldt wrote:
> > Richard Owlett wrote:
>
> > (snip)
>
> >> 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.
>
> > Square waves are not orthogonal, which makes the transform
> > harder, but not impossible.
>
> Why would they be any less orthogonal than sets of sines and cosines?
> They should also span the space as well as sines and cosines.

It depends on which waves.  Even waves are orthogonal to odd
waves (sinusoidal or square), periodic waves are orthogonal
to waves at 2^N times the frequency (sinusoidal or square).
But square waves are not orthogonal to square waves at three
times the frequency (or other odd multiples), whereas
sinusoids are.

There are Walsh functions which do span a space and are
orthogonal, but they are not purely "square" waves, if square
means periodic with a constant 50% duty cycle.  Since they
are not completely constant in periodicity or duty cycle,
using them to determine "frequency content", for more common
usages of that term, might be problematic.




IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M