Non Power-of-2 length FFT
Gerald T Beauregard
27/3/2000
"Steven G. Johnson" <ste...@gil-galad.mit.edu> wrote in message
news:slrn8e0gd0....@gil-galad.mit.edu...
> Gerald T Beauregard <gb...@mindmaker.com> wrote:
> >For efficiency, normally people use an FFT, so the frames must be a power
of
> >2 samples long.
>
> An aside: FFT algorithms in general do not require power-of-two lengths,
> although some implementations do. (Sorry to interject; this misconception
> is a pet peeve of mine. :-)
Steven: I started a new thread on this because I think this topic by itself
should generate quite a bit of interest. I've also posted it to comp.dsp.
Many (most?) FFT implementations do assume power of 2 lengths. With many
FFTs you don't even pass in the length, but instead log2 of the length. The
Intel Signal Processing Library works that way (to cite just one example).
However, you are *absolutely* correct that not all implementations require
power-of-two lengths! (Though I must admit, I didn't realize this until
reading your response to my post). I just took a peak at some of the
documentation at www.fftw.org, and came across this:
"The first argument, n, is the size of the transform you are trying to
compute. The size n can be any positive integer, but sizes that are products
of small factors are transformed most efficiently. "
and also...
"FFTW is best at handling sizes of the form 2^a 3^b 5^c 7^d 11^e 13^f, "
The most common audio sample rates are highly composite numbers (is that the
right term?). Specifically:
44100=(2*3*5*7)^2, and
48000=2^7*3*5^3
This means that it's possible to do *Fast* Fourier Transforms on nice, exact
sub-second audio segments.
For example, suppose you're operating at 44.1kHz, and you want to use
analysis frames that are *exactly* 20ms long. Each frame will be 882
samples long. Before seeing your post, I would have assumed that to do a
DFT on 882 samples, I'd need to use a slow algorithm; and that if I were
concerned about efficiency, I'd just bump up the frame length to 1024
samples. However, since 882 is a product of small factors (882=2*3*3*7*7),
it looks like it can be done efficiently with FFTW's more flexible
algorithm. Very cool!
Does any know offhand what the tradeoffs are for using non power-of-two
lengths? For example, would an 882 point FFT be just as quick as a 1024
point one?
-Gerry
PS: Are non power-of-2 length FFTs common knowledge, or did I just fall
asleep during too many lectures in my DSP courses??
---------------------------------------------------------
Gerald T Beauregard Mindmaker, Pte Ltd
Tel: (65) 425-2280 #03-04 Creative Resource
Fax: (65) 425-2278 31 International Business Park
gb...@mindmaker.com Singapore 609921
---------------------------------------------------------
Rick Lyons
<gb...@mindmaker.com> wrote:
>In a thread on "Spectral Analysis", in comp.speech.research, started on
>27/3/2000
>"Steven G. Johnson" <ste...@gil-galad.mit.edu> wrote in message
>news:slrn8e0gd0....@gil-galad.mit.edu...
>> Gerald T Beauregard <gb...@mindmaker.com> wrote:
>> >For efficiency, normally people use an FFT, so the frames must be a power
>of
>> >2 samples long.
>>
>> An aside: FFT algorithms in general do not require power-of-two lengths,
>> although some implementations do. (Sorry to interject; this misconception
>> is a pet peeve of mine. :-)
>
>Steven: I started a new thread on this because I think this topic by itself
>should generate quite a bit of interest. I've also posted it to comp.dsp.
>
> (good stuff snipped)
>
>Does any know offhand what the tradeoffs are for using non power-of-two
>lengths? For example, would an 882 point FFT be just as quick as a 1024
>point one?
>
>-Gerry
>
>PS: Are non power-of-2 length FFTs common knowledge, or did I just fall
>asleep during too many lectures in my DSP courses??
>
Hi Gerry,
I think the existence of "non-power of two" FFTs
are fairly well-known, but their implementation details
are known to but a few (me included).
I have a fairly long bibliography of FFT papers
(mostly those darned IEEE journals). If you have
access to a technical library, and you're feeling
energetic, I'd be happy to e-mail the
bibliography to you. I'll bet some of those papers
discuss the tradeoffs for using non power-of-two
lengths.
[-Rick-]
John B. Sampson
>Does any know offhand what the tradeoffs are for using non power-of-two
>lengths? For example, would an 882 point FFT be just as quick as a 1024
>point one?
According to the following Matlab program:
n = 1024;
x = round(32767*(2*rand(1,n)-ones(1,n)));
tic;
for i = 1:100
y = fft(x);
end
toc/100
It takes twice as long on average for Matlab to do N=882 as N=1024 (500 MHz
P3, Windows NT).
>
>-Gerry
>
>PS: Are non power-of-2 length FFTs common knowledge, or did I just fall
>asleep during too many lectures in my DSP courses??
>
Powers of 2 are just a special case. The basic idea of the FFT is to factor
N = r1 x r2 x r3 x ... x rk were the r's are small numbers. If all the r's
are the same then you have a radix r FFT. r=2 and 4 are common. If the r's
are not the same then you have a mixed radix FFT.
It can take longer to do smaller FFTs because the amount of work required to
do the 3s, 5s, 7s is on many processors more than doing just 2s but more of
them.
John
Rick Lyons
wrote:
> (snipped)
>
>Hi Gerry,
> I think the existence of "non-power of two" FFTs
>are fairly well-known, but their implementation details
>are known to but a few (me included).
>
Arrghh! I typed that wrong! Meant that I'm
*NOT* one of the people who understand those
algorithms!!!
[-Rick-]
E. Robert Tisdale
> Meant that I'm *NOT* one of the people
> who understand those algorithms!
Take a look at
The C++ Scalar, Vector, Matrix and Tensor class library
http://www.netwood.net/~edwin/svmt/
The portable document file in .../svmt/doc/fft.pdf
derives a simple recursive formula
for fast discrete Fourier transforms of any length n
using a [least] prime factor algorithm --
what some people call a Cooley-Tukey FFT algorithm.
It shows how to estimate the computational complexity
of the algorithm by counting multiplications.
It also describes a digit reverse algorithm
which can be used to transpose the elements in place
either before or after the transform.
Altogether, it's just six easy pages of reading.
Hope this helps,
E. Robert Tisdale <ed...@netwood.net>
Mike Rosing
> PS: Are non power-of-2 length FFTs common knowledge, or did I just fall
> asleep during too many lectures in my DSP courses??
Yes and no. It's common knowledge if you have done it before :-)
Check out Nussbaum's book. I've got it at home, so I can't give the
full reference. He goes thru all the number theory and shows lots of
good quick algorithms for computing 3, 5, 7, as well as 4, and 8 point
FFT's, and then how to link them to make any power you want.
It's actually pretty easy to implement. I did it over 10 years ago on
a VAX in FORTRAN. Ought to be trivial in assembler :-)
Patience, persistence, truth,
Dr. mike
Mike Rosing
> Check out Nussbaum's book. I've got it at home, so I can't give the
> full reference. He goes thru all the number theory and shows lots of
> good quick algorithms for computing 3, 5, 7, as well as 4, and 8 point
> FFT's, and then how to link them to make any power you want.
found it on the library pages. Spelled his name wrong too, sorry.
Author:
Nussbaumer, Henri J., 1931-
Title:
Fast Fourier transform and convolution algorithms / Henri
J. Nussbaumer.
Edition:
2nd corr. and updated ed.
Publisher:
Berlin ; New York : Springer-Verlag, 1982.
Description:
xii, 276 p. : ill. ; 24 cm.
Series:
Springer series in information sciences ; 2
Location:
Physics Library
Call Number:
QA403.5 N87 1982
Copy Number:
1
Status:
Renewed - Due on 05/14/2000
So neither of us can check it out :-\
Steven G. Johnson
>>Does any know offhand what the tradeoffs are for using non power-of-two
>>lengths? For example, would an 882 point FFT be just as quick as a 1024
>>point one?
>It takes twice as long on average for Matlab to do N=882 as N=1024 (500 MHz
>P3, Windows NT).
Matlab's FFT implementation is (currently) not really that optimized,
especially for non-power-of-two cases.
Power-of-two sizes tend to be slightly quicker than sizes with other
factors for very similar sizes--or at least, they require fewer operations
(for the Cooley-Tukey algorithm, anyway...things become more subtle if you
consider Prime-Factor algorithms). (When you start doing large transforms
and taking the cache into account, things become hairier, because powers
of two are especially hard on direct-mapped caches.) When the sizes start
becoming significantly different, then highly-composite sizes may be
faster than a larger power-of-two, but it depends a lot on the
factorization, the implementation, and the machine.
(882 has two factors of 7; a size-7 DFT in FFTW takes 70% more arithmetic
operations than a size-8 DFT, although implementations with somewhat fewer
operation counts are possible. FFTW also comes pre-packaged with more
hard-coded transforms for power-of-two sizes, although you can change
this.)
For example, on a 400MHz Pentium-II (Linux/gcc) with FFTW and
single-precision complex transforms, a size-1024 FFT takes 172 us, a
size-882 FFT takes 215 us, and a size-864 (2^5*3^3) FFT takes 167 us.
(For purely real inputs, the numbers are 87/108/86 us, respectively.)
YMMV.
>Powers of 2 are just a special case. The basic idea of the FFT is to
>factor N = r1 x r2 x r3 x ... x rk were the r's are small numbers. If all
>the r's are the same then you have a radix r FFT. r=2 and 4 are common.
>If the r's are not the same then you have a mixed radix FFT.
This is the basic idea of the Cooley-Tukey FFT algorithm, anyway. There
are many other O(N lg N) FFT algorithms, some of which factorize N, and
some of which don't (there are FFTs even for prime sizes).
Cordially,
Steven G. Johnson
Steven G. Johnson
>for fast discrete Fourier transforms of any length n
>using a [least] prime factor algorithm --
>what some people call a Cooley-Tukey FFT algorithm.
...what everyone except for you calls a Cooley-Tukey FFT algorithm. =)
The literature exclusively (as far as I have seen) reserves the term
"Prime Factor Algorithm" for decompositions that only work for
relatively-prime factors but which avoid the twiddle factors of
Cooley-Tookey (e.g. the Good-Thomas algorithm). Do you have any
references to the contrary? (Internet posts don't count...people confuse
split-radix and mixed-radix online all the time too.)
Steven
Tom
Gerald T Beauregard wrote -
>Does any know offhand what the tradeoffs are for using non power-of-two
>lengths? For example, would an 882 point FFT be just as quick as a 1024
>point one?
A lot depends on how you want to implement them - in general terms, the
number of cycles for the complex arithmetic part of the calculation is
minimised compared with the Radix-2 FFT, whilst the control part of the
process is more complicated than for the Radix-2. So they can be more
difficult to implement on a PC, where the control overhead can be a problem
(generally have to write in-line code rather than loops), but because of
their reduced arithmetic complexity, they can be very useful for hardware
implementations using for example FPGAs.
There are a several points worth bearing in mind with mixed radix/low
arithmetic complexity FFTs -
1. If the block lengths chosen for the small transforms are relatively
prime, then Goods decomposition can be used to re-arrange the order in
which the complex samples are used for each pass of the transform in such a
way that the twiddle factor complex multiplies are eliminated. This reduces
computational complexity considerably.
2. There are a number of algorithms to calculate the small block length
transforms very efficiently: for example Winograd (amongst others) gives
minimum multiplicative complexity 'recipes' for a number of small length
transforms.
3. If you are really heavily into maths, then it might be worth considering
Number Theoretic transforms - e.g.. Fermat Number Transforms. These reduce
the arithmetic complexity to an absolute minimum (but need a Ph.D. in
Number Theory to figure out how they work!!).
The downside of these reduced arithmetic complexity algorithms is that the
control structure is no longer as regular as it is for the radix2 FFT.
However, they can be implemented in hardware moderately quickly. Back in
the early 1980s, I was using 221-point complex transforms (13 x 17)
implemented in TTL MSI that took less than 2 microseconds per FFT and
complete 899-point complex transforms (29 x 31) calculated in around 1
millisecond on a custom gate array (5k gates, 3 micron geometry): both of
these designs used Prime Radix transforms that can be implemented very
efficiently in hardware. (These were pretty good speeds twenty years
ago!). Currently, I use 1023-point complex transforms (31 x 33) on a small
Xilinx Virtex FPGA that run in under 10 microseconds (pretty much
comparable in performance with current 'state of the art' FFT chips). All
of these designs used 24 or 32 bit arithmetic, with output dynamic ranges
in excess of 140 dB.
Tom C.
**********************************************
Tom Curtis
Curtis Technology (UK) Ltd
to...@curtistech.co.uk
**********************************************
Gordon Sande
<ed...@netwood.net> wrote:
>Rick Lyons wrote:
>
>> Meant that I'm *NOT* one of the people
>> who understand those algorithms!
>
>Take a look at
>The C++ Scalar, Vector, Matrix and Tensor class library
>
> http://www.netwood.net/~edwin/svmt/
>
>The portable document file in .../svmt/doc/fft.pdf
>derives a simple recursive formula
>for fast discrete Fourier transforms of any length n
>using a [least] prime factor algorithm --
>what some people call a Cooley-Tukey FFT algorithm.
In the usual terminology the Prime Factor Algorithm
is attributed to Good. PFA uses the Chinese Remainder
Theorem in its quite nontrivial indexing which has
the additional benefit of not having the so-called
twiddle factors. PFA use factors which are mutually
prime so 24=3*8 and it can not itself use the fact
that 8=2*2*2 although whatever is used to do the
8 factor may internally use that knowledge. Radix
factoring algorithms try to factor things into small
primes and would use 24=2*2*2*3 for four factors
which need not be relatively prime. Cooley-Tukey is
the most commonly known example of a radix factoring
algorithm.
When the factors are no longer small primes the adjective
fast is misplaced. This can be cured by either the
Rader or Bluestein algorithms which recover the n*log(n)
costing (with larger coeficient) even for large primes
by the use of additional storeage and logic.
>It shows how to estimate the computational complexity
>of the algorithm by counting multiplications.
>It also describes a digit reverse algorithm
>which can be used to transpose the elements in place
>either before or after the transform.
>Altogether, it's just six easy pages of reading.
Or if it is not in place one can forgo all the
digit reversed indexing fuss.
E. Robert Tisdale
> E. Robert Tisdale <ed...@netwood.net> wrote:
>
> > for fast discrete Fourier transforms of any length n
> > using a [least] prime factor algorithm --
> > what some people call a Cooley-Tukey FFT algorithm.
>
> ...what everyone except for you calls a Cooley-Tukey FFT algorithm. =)
>
> The literature exclusively (as far as I have seen)
> reserves the term "Prime Factor Algorithm"
> for decompositions that only work for relatively-prime factors
> but which avoid the twiddle factors of Cooley-Tookey
> (e.g. the Good-Thomas algorithm).
> Do you have any references to the contrary?
> (Internet posts don't count...
> people confuse split-radix and mixed-radix online all the time too.)
Dear subscriber,
Steven and I have this discussion from time to time.
I always know exactly what Steven is going to say
and he always knows exactly what I'm going to say.
The discussion is meant to be helpful and not confuse
other subscribers.
Authors don't always have the presence of mind
to give their algorithms a descriptive name before they publish them
so other people use the authors' names to refer to the algorithm
and use the same name to identify other similar algorithms.
I think most DSP engineers associate the Cooley-Tukey algorithm
with the radix 2 algorithm but Cooley-Tukey is more generally
associated with any least prime factorization of an arbitrary length n.
An algorithm based upon a _mutually_ prime factorization of length n
is simply called a Prime Factor Algorithm (PFA)
which unfortunately leaves out the essential distinction
that two or more factors are simply mutually prime
and not necessarily prime numbers themselves.
Mutually prime factor algorithms can save a lot of multiplications
but require a significant amount of overhead up front
to determine exactly how to factor the length n
in order to take maximum advantage of the algorithm.
They work best when the factorization is done off-line for a fixed n
or when the cost of the factorization can be amortized
over several subsequent invocations of an FFT of length n.
Steven G. Johnson
> "Steven G. Johnson" wrote:
> > E. Robert Tisdale <ed...@netwood.net> wrote:
> > > for fast discrete Fourier transforms of any length n
> > > using a [least] prime factor algorithm --
> > > what some people call a Cooley-Tukey FFT algorithm.
> >
> > ...what everyone except for you calls a Cooley-Tukey FFT algorithm. =)
>
I think I commented on it once before; it's not really much of a
discussion, though, because you never really respond to my basic point.
(This is why many arguments on Usenet seem to go in circles...sigh.)
(See also http://www.deja.com/=dnc/[ST_rn=ps]/getdoc.xp?AN=488913054)
> Authors don't always have the presence of mind
> to give their algorithms a descriptive name before they publish them
> [...various comments on the workings of various FFT algorithms...]
The point, Robert, is that everyone in the literature has agreed to use
"Prime Factor Algorithm" to describe a certain kind of algorithm for
mutually-prime factors, and "mixed-radix" Cooley-Tukey to describe the
(different) general factorization that C-T described when applied to
arbitrary factors (and which you describe on your web page).
While people *could* have chosen "prime factor algorithm" to describe C-T
(which is the point you try to make, as far as I can decipher), the fact
is that they *didn't*, and moreover they used that phrase for something
else. Why you insist on gratuitously turning the conventional terminology
on its head is beyond me.
What do you have against poor Cooley and Tukey? =)
Steven
Gordon Sande
<ed...@netwood.net> wrote:
>"Steven G. Johnson" wrote:
>
>> E. Robert Tisdale <ed...@netwood.net> wrote:
>>
>> > for fast discrete Fourier transforms of any length n
>> > using a [least] prime factor algorithm --
>> > what some people call a Cooley-Tukey FFT algorithm.
>>
>> ...what everyone except for you calls a Cooley-Tukey FFT algorithm. =)
>>
>> The literature exclusively (as far as I have seen)
>> reserves the term "Prime Factor Algorithm"
>> for decompositions that only work for relatively-prime factors
>> but which avoid the twiddle factors of Cooley-Tookey
>> (e.g. the Good-Thomas algorithm).
>> Do you have any references to the contrary?
>> (Internet posts don't count...
>> people confuse split-radix and mixed-radix online all the time too.)
>
>Dear subscriber,
>
>Steven and I have this discussion from time to time.
>I always know exactly what Steven is going to say
>and he always knows exactly what I'm going to say.
>The discussion is meant to be helpful and not confuse
>other subscribers.
>
>Authors don't always have the presence of mind
>to give their algorithms a descriptive name before they publish them
>so other people use the authors' names to refer to the algorithm
>and use the same name to identify other similar algorithms.
>I think most DSP engineers associate the Cooley-Tukey algorithm
>with the radix 2 algorithm but Cooley-Tukey is more generally
>associated with any least prime factorization of an arbitrary length n.
Breaking a 1024 point Fourier transform into two factors of 32 each
(1024=32*32) would certainly fall within the standard notion of
what is meant by Cooley-Tukey but does not seem to have any least
primes factorization involved. Doing a 10 point Fourier transform
with 2*5 PFA would seem to fit the requirement of being "any least
prime factorization" but seems outside the reach of Cooley-Tukey.
The usual terminology seems to fall out with just plain factors
describing a radix based indexing form and prime factors describing
a Chinese Remainder Theorem (modular arithmetic) form of indexing.
The radix indexing ends most commonly with the name of "Cooley-Tukey"
and the CRT indexing with the name "Prime Factor Algorithm". The
primeness (or not) of the factors is irrelevant for radix indexing
and so is never mentioned. It is important for the CRT indexing and
so is mentioned. It is unfortunately common that some technical
distinctions/conditions get dropped so that primeness should really
be relative primeness which would then include a PFA for 210=10*21.
You are of course free to give any algorithm any name you choose for
your purposes but that is not prone to help communication.
E. Robert Tisdale
> The point, Robert, is that
> everyone in the literature has agreed to use "Prime Factor Algorithm"
> to describe a certain kind of algorithm for mutually-prime factors,
> and "mixed-radix" Cooley-Tukey (C-T) to describe the (different)
> general factorization that C-T described
> when applied to arbitrary factors (and which you describe on your web page).
>
> While people *could* have chosen "prime factor algorithm" to describe C-T
> (which is the point you try to make, as far as I can decipher),
> the fact is that they *didn't*,
> and moreover they used that phrase for something else.
> Why you insist on gratuitously turning
> the conventional terminology on its head is beyond me.
No. I think that you have just explained very well why I
"insist on turning the conventional terminology on its head."
> What do you have against poor Cooley and Tukey? =)
They neglected to give their algorithm a short descriptive name
so that people would be able to identify it without confusion?
The naming conventions are important
when searching through the literature
but I don't think that
very many subscribers to the comp.dsp newsgroup
besides you and I do that.
It just strikes me as a bit pedantic
when you correct someone's usage
of the names of FFT algorithms
especially when you don't bother with any explanation
of how they got the names that they have.
Steven G. Johnson
> "insist on turning the conventional terminology on its head."
You wrote a paragraph on the differences between the Cooley-Tukey and PFA
FFT algorithms; I don't see how this constitutes a defense for rejecting
conventional nomenclature.
> > What do you have against poor Cooley and Tukey? =)
>
> They neglected to give their algorithm a short descriptive name
> so that people would be able to identify it without confusion?
Instead, you use a short descriptive name which (a) is already used and
well-accepted for something else and (b) is not very accurate (mixed-radix
Cooley-Tukey does not require that its factorization be into prime
numbers).
Do you rename all equations that are named after people? Newton's
Equations? Kepler's Laws? Schrodinger's Equation? The Pythagorean
Theorem? It's one thing to have personal names for your own use; it's a
disservice to teach your names to others without explanation--especially
if your personal names are normally used for something else!
> The naming conventions are important
> when searching through the literature
Or when talking with anyone else who is familiar with these algorithms.
Or when reading any textbook which describes these algorithms. Or when
communicating with anyone who has done those things.
Common language is important in mathematics just as in other fields.
> It just strikes me as a bit pedantic
> when you correct someone's usage
> of the names of FFT algorithms
When you respond to a person's questions by promoting your own
terminology, at odds with what everyone else uses and what is used in all
the textbooks and literature, I think you are inviting confusion. When
this occurs in a thread that I am participating in, I try to correct the
confusion.
> especially when you don't bother with any explanation
> of how they got the names that they have.
Etymology is interesting, but it is not always relevant. If I am
proofreading someone's paper and I see a spelling error, I don't respond
with a discourse on the history of the word in English--I just correct the
mistake. In this case, the only relevant point is that "Prime Factor
Algorithm" is accepted to refer to a completely different FFT algorithm
than mixed-radix Cooley-Tukey.
Steven
Steven G. Johnson
> algorithm. The prime factor algorithm is a special case of the mixed
> radix algorithm which requires each of the DFT kernels to be relatively
> prime with respect to the all other kernels in the FFT. By meeting that
Hmm, we're talking about the same thing, but I'm not sure I would describe
it in those terms. The Prime Factor Algorithm isn't really a "special
case" of the mixed-radix algorithm (which usually refers to Cooley-Tukey),
it is a different algorithm entirely. e.g. you can use both PFA and
mixed-radix (C-T) to decompose 10=5*2, but the decompositions are
different, and mixed-radix (i.e. Cooley-Tukey) has twiddle factors while
PFA doesn't.
Steven
Steven G. Johnson
> couple of rules of thumb for FFTW (and other implementations on general
> purpose processors?):
>
> - Given two similar sizes, you'll get better performance using the size that
> has smaller factors.
Well...
That rule seems to work fairly well most of the time, although as I
remarked there can be some complications if you use non-Cooley-Tukey
algorithms (e.g. Prime Factor) or for large transforms on machines with
direct-mapped caches. For example, on a 192MHz MIPS R10000 (Origin2000),
a size 65536=2^16 transform takes 49ms with FFTW (double-precision
complex, SGI cc -O3), while a (larger!) size 69984=2^5*3^7 transform takes
only 34ms, 44% faster! (That machine seems inordinately sensitive to
cache voodoo.)
(I've also noticed that cache conflicts for powers of 2 can be especially
hard on shared-memory parallel transforms, but that's another story.)
> - It's hardly worth using a NonPof2 size for performance reasons alone,
> since at best the NonPof2 FFT will be only slighly faster than the next
> bigger Pof2 size.
Well...
It depends on what you consider "slightly." e.g. a size 34992=2^4*3^7
transform on my machine finishes in 17 ms, about 140% faster than the next
power of 2 (65536, 41 ms). For multidimensional transforms the difference
is even more dramatic, and non-power-of-2 transforms are essential. e.g.
a 64x64x64 FFT on my machine takes 118 ms, while a 36x36x36 FFT takes 18
ms.
The moral of this story is that computers these days are a mess, and the
only reliable way to optimize performance is by experimentation. My
generic advice, when speed is important, is to restrict yourself to
products of 2, 3, and 5, as this gives good results under most
circumstances (at least with FFTW).
Steven
E. Robert Tisdale
> "E. Robert Tisdale" <ed...@netwood.net> wrote:
>
> > No. I think that you have just explained very well why I
> > "insist on turning the conventional terminology on its head."
>
> You wrote a paragraph on the differences
> between the Cooley-Tukey and PFA FFT algorithms;
> I don't see how this constitutes a defense
> for rejecting conventional nomenclature.
>
> > > What do you have against poor Cooley and Tukey? =)
> >
> > They neglected to give their algorithm a short descriptive name
> > so that people would be able to identify it without confusion?
>
> Instead, you use a short descriptive name which
> (a) is already used and well-accepted for something else and
> (b) is not very accurate (mixed-radix Cooley-Tukey
> does not require that its factorization be into prime numbers).
>
> Do you rename all equations that are named after people?
> Newton's Equations? Kepler's Laws? Schrodinger's Equation?
> The Pythagorean Theorem?
> It's one thing to have personal names for your own use;
> it's a disservice to teach your names to others without explanation --
> especially if your personal names are normally used for something else!
>
> > The naming conventions are important
> > when searching through the literature
>
> Or when talking with anyone else who is familiar with these algorithms.
> Or when reading any textbook which describes these algorithms.
> Or when communicating with anyone who has done those things.
>
> Common language is important in mathematics just as in other fields.
>
> > It just strikes me as a bit pedantic
> > when you correct someone's usage
> > of the names of FFT algorithms
>
> When you respond to a person's questions
> by promoting your own terminology,
> at odds with what everyone else uses
> and what is used in all the textbooks and literature,
> I think you are inviting confusion.
> When this occurs in a thread that I am participating in,
> I try to correct the confusion.
>
> > especially when you don't bother with any explanation
> > of how they got the names that they have.
>
> Etymology is interesting, but it is not always relevant.
> If I am proofreading someone's paper and I see a spelling error,
> I don't respond with a discourse on the history of the word in English --
> I just correct the mistake. In this case, the only relevant point is that
> "Prime Factor Algorithm" is accepted to refer to
> a completely different FFT algorithm than mixed-radix Cooley-Tukey.
I'm not sure that I have ever met anyone
as passionate about nomenclature as you appear to be.
Yes. It is common practice to associate short descriptive names
with Newton's Equations, Maxwell's Equations, Schrodinger's Equation,
Kepler's Laws, etc.
Different Schools and University Departments may have different cultures.
Engineering Departments tend to be more pedantic than Physics Departments
about properly identifying these things.
It isn't my fault that
the conventional names for FFT algorithms were poorly chosen.
They are simply a sequence of historical accidents
that we are all obliged to live with the best we can.
And I don't think that I am contributing to the confusion
by calling the algorithm that I describe in .../svmt/doc/fft.pdf
a [least] prime factor algorithm.
There isn't the slightest bit of evidence
that anyone is confusing my [least] prime factor algorithm
with the conventional [mutually] prime factor algorithm.
The problem is that some people identify the FFT
exclusively with a radix 2 Cooley-Tukey algorithm.
All you can do is point out that Cooley-Tukey applies
to higher radix and mixed radix algorithms
as well as to radix 2 algorithms.
I describe a mixed radix algorithm where the radices
are the least prime factors of the transform length n.
It seems perfectly reasonable to me to call this
a [least] prime factor algorithm and leave it to people
like yourself to point out the conventional meaning of PFA.
And that seems to be working out rather well,
don't you agree?
Ray Andraka
Cooley-Tukey does not use the Prime factor algorithm, it uses the mixed
radix algorithm. The prime factor algorithm is a special case of the
mixed radix algorithm which requires that the DFT building blocks be
relatively prime. Cooley-Tukey is made up of radix two stages using the
mixed radix construction. As a result Cooley-Tukey FFTs are always a
power of two.
The mixed radix algorithm allows you to cascade smaller arbitrarily
sized FFT kernels to obtain a larger FFT. It generally requires a
complex multiply to effect a rotation of each element between each of
the small kernels, which are the 'twiddle factors' in the Cooley-Tukey
algorithm. The prime factor algorithm is a special case of the mixed
radix algorithm which requires each of the DFT kernels to be relatively
prime with respect to the all other kernels in the FFT. By meeting that
criterion, the complex multiplies between the kernels are eliminated
(the rotation is always 1,-1,j or -j). Clearly the Cooley-Tukey
algorithm is not like this.
To do a non-power of two FFT, you need non-power-of-two DFT kernels.
These generally are not as efficient computationally as the power of two
ones, so you do lose some efficiency. Some of the more efficient
implementations are those attributed to Singleton, Winograd, and Rader.
These are derived from decomposing the DFT in various ways, for example
by viewing it as a convolution. There are sometimes advantages to using
something other than the Cooley-Tukey. For example, the Winograd FFT is
an optimization to minimize the number of multiplications: the only
multiplies are in the middle butterfly, all the other butterflies have
no multiplications associated with them. For hardware implementations,
that can significantly reduce the gate count.
"E. Robert Tisdale" wrote:
> Rick Lyons wrote:
>
> > Meant that I'm *NOT* one of the people
> > who understand those algorithms!
>
> Take a look at
> The C++ Scalar, Vector, Matrix and Tensor class library
>
> http://www.netwood.net/~edwin/svmt/
>
> The portable document file in .../svmt/doc/fft.pdf
> derives a simple recursive formula
> for fast discrete Fourier transforms of any length n
> using a [least] prime factor algorithm --
> what some people call a Cooley-Tukey FFT algorithm.
> It shows how to estimate the computational complexity
> of the algorithm by counting multiplications.
> It also describes a digit reverse algorithm
> which can be used to transpose the elements in place
> either before or after the transform.
> Altogether, it's just six easy pages of reading.
>
> Hope this helps,
>
> E. Robert Tisdale <ed...@netwood.net>
--
-Ray Andraka, P.E.
President, the Andraka Consulting Group, Inc.
401/884-7930 Fax 401/884-7950
email rand...@ids.net
http://users.ids.net/~randraka
Gerald T Beauregard
> For example, on a 400MHz Pentium-II (Linux/gcc) with FFTW and
> single-precision complex transforms, a size-1024 FFT takes 172 us, a
> size-882 FFT takes 215 us, and a size-864 (2^5*3^3) FFT takes 167 us.
> (For purely real inputs, the numbers are 87/108/86 us, respectively.)
> YMMV.
Thanks for getting these numbers! Extrapolating from your data, here are a
couple of rules of thumb for FFTW (and other implementations on general
purpose processors?):
- Given two similar sizes, you'll get better performance using the size that
has smaller factors.
- It's hardly worth using a NonPof2 size for performance reasons alone,
since at best the NonPof2 FFT will be only slighly faster than the next
bigger Pof2 size.
Reasonable summary?
-Gerry
E. Robert Tisdale
> BZZZZZZT Wrong!
>
> Cooley-Tukey does not use the Prime factor algorithm,
> it uses the mixed radix algorithm. The prime factor algorithm
> is a special case of the mixed radix algorithm
> which requires that the DFT building blocks be relatively prime.
> Cooley-Tukey is made up of radix two stages
> using the mixed radix construction.
> As a result Cooley-Tukey FFTs are always a power of two.
Damn it Ray,
You just made a liar out of me.
I just told Steven G. Johnson <ste...@alum.mit.edu>
that no one was confusing my [least] prime factor algorithm
with the conventional [mutually] prime factor algorithm;-)
E. Robert Tisdale
> The moral of this story is that computers these days are a mess,
> and the only reliable way to optimize performance is by experimentation.
I assume that this applies to Analog Devices SHARC
http://www.analog.com/industry/dsp/
and Texas Instrument's TMS320C67X Family
http://www.ti.com/sc/docs/products/dsp/
> My generic advice, when speed is important,
> is to restrict yourself to products of 2, 3, and 5, as this gives good results
> under most circumstances (at least with FFTW).
I think it would be really nice to have a tool like FFTW
which would emit highly optimized C (or even Assembler)
for high end DSP chips for any transform length n.
Ray Andraka
sentence in Smith&Smith's chapter 9.6.1 Prime Factor Algorithm Introduction: "
The Prime Factor algorithm is a special form of the mixed-radix algorithm
presented in Section 9.7"
If you use PFA, you will get different decompositions depending on which
algorithm you use, as well as depending on the ordering of your factors. I'm
too rusty with the derivation to reproduce it (I tried, but without digging in
the books I can't remember). As I recall, modulo arithmetic on the twiddles
makes them fall out.
"Steven G. Johnson" wrote:
> Ray Andraka <rand...@ids.net> wrote:
> > algorithm. The prime factor algorithm is a special case of the mixed
> > radix algorithm which requires each of the DFT kernels to be relatively
> > prime with respect to the all other kernels in the FFT. By meeting that
>
> Hmm, we're talking about the same thing, but I'm not sure I would describe
> it in those terms. The Prime Factor Algorithm isn't really a "special
> case" of the mixed-radix algorithm (which usually refers to Cooley-Tukey),
> it is a different algorithm entirely. e.g. you can use both PFA and
> mixed-radix (C-T) to decompose 10=5*2, but the decompositions are
> different, and mixed-radix (i.e. Cooley-Tukey) has twiddle factors while
> PFA doesn't.
>
> Steven
--
Jerry Avins
>
> Rick Lyons wrote:
>
> > Meant that I'm *NOT* one of the people
> > who understand those algorithms!
>
> Take a look at
> The C++ Scalar, Vector, Matrix and Tensor class library
>
> http://www.netwood.net/~edwin/svmt/
>
> The portable document file in .../svmt/doc/fft.pdf
> derives a simple recursive formula
> for fast discrete Fourier transforms of any length n
> using a [least] prime factor algorithm --
> what some people call a Cooley-Tukey FFT algorithm.
> It shows how to estimate the computational complexity
> of the algorithm by counting multiplications.
> It also describes a digit reverse algorithm
> which can be used to transpose the elements in place
> either before or after the transform.
> Altogether, it's just six easy pages of reading.
>
> Hope this helps,
>
> E. Robert Tisdale <ed...@netwood.net>
Lord, how does he manage to confuse things so thoroughly with so few
words? He confounds prime and relatively prime; he endorses a count of
multiplications as a measure of complexity when on DSP chips,
multiplying, adding, and MAC have the same cost. And he utters his crap
with such confidence!
Jerry
--
Engineering is the art of making what you want from things you can get.
-----------------------------------------------------------------------
E. Robert Tisdale
> Just recently I accidentally tried to use
> my own implementation of a radix 2 FFT
> on something other than a power of two.
> At first I was concerned that the results were misscaled,
> but when I realized what I had done,
> I was stunned to see that it even had the right shape.
> I've been meaning to get around to figuring out why it works.
It doesn't always work.
The trick is to figure out when it does and doesn't work.
If you samples simply trail off to zero, you can almost always
pad with zeros out to the next power of 2 and use a radix 2 FFT.
This may require almost twice as much memory for the FFT
so it might not be so attractive.
If the samples don't trail off to zero,
you might not be able to get anything useful
by simply padding with zeros.
It depends upon what you need from your FFT.
Sometimes, nothing but an FFT of length n will do
and you will be obliged to use something other than
the radix 2 FFT if n is not a power of 2.
> I learned the FFT from a National Semiconductor databook
> (Application Note 487, Ashok Krishnamurthy).
> Their explanation made it pretty clear that non-radix-2 was a possibility.
> They showed that an FFT is the sum of two halves of an FFT,
> so you see right there how it could also be the sum of three thirds, etc.
Yes. But can you actually implement a non-radix-2 FFT?
How would you reorder the array for example?
Steven G. Johnson
>Sheesh, You made me go back and dig through my books again. From the first
>sentence in Smith&Smith's chapter 9.6.1 Prime Factor Algorithm Introduction: "
>The Prime Factor algorithm is a special form of the mixed-radix algorithm
>presented in Section 9.7"
>
>If you use PFA, you will get different decompositions depending on which
>algorithm you use, as well as depending on the ordering of your factors. I'm
>too rusty with the derivation to reproduce it (I tried, but without digging in
>the books I can't remember). As I recall, modulo arithmetic on the twiddles
>makes them fall out.
Hi Ray,
"Mixed-radix" normally refers to Cooley-Tukey, but you're right in that it
sometimes is used to refer to any FFT algorithm that divides the size by
various factors (e.g. Cooley-Tukey, PFA, or the generalized Bruun
algorithm). I just think that one should be careful in doing so, so as not
to confuse these other algorithms with Cooley-Tukey; I don't have S&S, but
presumably they are clear in their exposition.
Just for information:
For a DFT of size N = N1*N2, with the input array indexed by n=0..N-1 and
the output array indexed by k=0..N-1, the mixed-radix Cooley-Tukey
algorithm does the reindexing:
n = N2*n1 + n2
k = k1 + N1*k2
and PFA (for N1,N2 coprime) uses the entirely different reindexing:
n = N2*n1 + N1*n2
k = N2*(1/N2)*k1 + N1*(1/N2)*k2
where n1 and k1 go from 0..N1-1, n2 and k2 go from 0..N2-1, and (1/N1) and
(1/N2) are the inverses modulo N2 and N1, respectively (which exist
because the factors are coprime).
Recall that the original DFT had n*k, modulo N, in the exponent of the Nth
root of unity. In both cases, you have a 1d to 2d mapping of indices, and
1d sub-DFTs of size N1 and N2 (from the n1*k1 and n2*k2 terms). In C-T,
you get the twiddle factors from the k1*n2 cross-term in the exponent.
In PFA, the cross-terms vanish because they both contain N1*N2=N=0 mod N,
and thus there are no twiddle factors--it just looks like a 2D DFT in the
new indices.
(Of course, there are several (3?) variations of PFA, as you can freely
swap the coefficients between n and k.)
Regards,
Steven
Gordon Sande
<ed...@netwood.net> wrote:
>
>Yes. But can you actually implement a non-radix-2 FFT?
>How would you reorder the array for example?
>
By knowing the tricks of the trade and having a knowledge
of the literature. That is how the big boys do it E.Bob!
One way is to avoid the issue by using enough scratch memory
so that the results come out in the correct order. With
modern availability of memory this is often a very viable option.
If you are lucky the size with be either a square or a square times
a single factor. 36=6*6=2*3*3*2 for a "symmetric" factorization.
Four nested loops will do the job in the fairly common fashion.
210=36*5=2*3*5*3*2 is again a "symmetric" factorization. If it
is not "symmetric" then it takes more fuss. If there is enough memory
then just copy it out, real to permuted scratch and back and then
imaginary to permuted scratch and back except one gets curious if
there is half the available scratch but not all of it. This is
much more plausible for one dimension of a multiple dimensional
tranform. The next level of avoidance is to use scratch only
within the unsymmetric part and when all else fails one is down
to chasing permutations. All of this was worked out before 1970
except it does not neatly fit into the one page appendix of a
signal processing text so it has to be repeatedly rediscovered.
The permutation can also be factored and the factoring mixed
with the FT factors. Neither set of factors is done in the naive
order. These are the inplace-inorder variants. The point of the
symmetric factorization is that the permutation factoring is
then quite easy and obvious.
The PFA has its own collection of problems/tricks.
Been there, done that, bought the T-shirt, washed it so often
that it faded and wore out and so was thrown out.
E. Robert Tisdale
> On Wed, 29 Mar 2000 19:29:36 +0000,
> "E. Robert Tisdale" <ed...@netwood.net> wrote:
>
> > Yes. But can you actually implement a non-radix-2 FFT?
> > How would you reorder the array for example?
>
> By knowing the tricks of the trade
> and having a knowledge of the literature.
> That is how the big boys do it E.Bob!
>
> One way is to avoid the issue by using enough scratch memory
> so that the results come out in the correct order.
> With modern availability of memory this is often a very viable option.
>
> If you are lucky
> the size will be either a square or a square times a single factor.
> 36=6*6=2*3*3*2 for a "symmetric" factorization.
> Four nested loops will do the job in the fairly common fashion.
> 210=36*5=2*3*5*3*2 is again a "symmetric" factorization.
> If it is not "symmetric" then it takes more fuss.
> If there is enough memory then just copy it out,
> real to permuted scratch and back
> and then imaginary to permuted scratch and back
> except one gets curious
> if there is half the available scratch but not all of it.
> This is much more plausible
> for one dimension of a multiple dimensional transform.
> The next level of avoidance is to use scratch
> only within the unsymmetric part
> and when all else fails one is down to chasing permutations.
> All of this was worked out before 1970
> except it does not neatly fit into the one page appendix
> of a signal processing text so it has to be repeatedly rediscovered.
>
> The permutation can also be factored
> and the factoring mixed with the FT factors.
> Neither set of factors is done in the naive order.
> These are the inplace-inorder variants.
> The point of the symmetric factorization is that
> the permutation factoring is then quite easy and obvious.
>
> The PFA has its own collection of problems/tricks.
>
> Been there, done that, bought the T-shirt,
> washed it so often that it faded
> and wore out and so was thrown out.
Sounds like a "run-of-the-mill" hack
for a particular transform length n.
I suppose it is acceptable if you enjoy writing code
and can get your employer to pay you
for writing essentially the same code over and over again.
Have you ever tried to write anything
that would do a transform in place for any n?
You know, reusable write once and forget it code?
Gordon Sande
<ed...@netwood.net> wrote:
Since you asked, I did it in 1968 and the resulting code
has been heavily used in crystallography since then.
It is also the basis of the code in NAG, although they
did not use the full package.
>
Gumboot
> Many (most?) FFT implementations do assume power of 2 lengths.
Funny this should come up. Just recently I accidentally tried to use my own
implementation of a radix 2 FFT on something other than a power of two. At
first I was concerned that the results were misscaled, but when I realised
what I had done I was stunned to see that it even had the right shape. I've
been meaning to get around to figuring out why it works.
> PS: Are non power-of-2 length FFTs common knowledge, or did I just fall
> asleep during too many lectures in my DSP courses??
I learnt the FFT from a National Semiconductor databook (Application Note
487, Ashok Krishnamurthy). Their explanation made it pretty clear that
non-radix-2 was a possibility. They showed that an FFT is the sum of two
halves of an FFT, so you see right there how it could also be the sum of
three thirds, etc.
--
# Please try to quote no more than you need to show the context of your post.
# email me at "clear.net.nz", username "gumboot"
E. Robert Tisdale
> Since you asked, I did it in 1968 and the resulting code
> has been heavily used in crystallography since then.
> It is also the basis of the code in NAG,
> although they did not use the full package.
It's too bad that you can't share it.
Gordon Sande
If you had asked, and many did, in 1968 you would have been
provided with a copy. Times change. Some folks and polite and
professional.
E. Robert Tisdale
> If you had asked, and many did, in 1968,
> you would have been provided with a copy.
> Times change. Some folks are polite and professional.
If you speak C++, you are welcome to take a look at my FFTs.
There are both real to complex and complex to real FFTs
as well as complex to complex FFTs for any transform length n.
They aren't high performance implementations.
They are part of a reference library
so I tried to make them as simple as possible
so that they would be easy to read, understand and maintain.
The reference library also serves as sort of a tutorial
where I demonstrate how to use subvector and submatix objects
by implementing other functions them such as FFTs.
You can find the code in
.../svmt/src/vector/doubleComplexVector/doubleComplexVector.cc
if you download the library archive from
http://www.netwood.net/~edwin/svmt/
Gumboot
> If you samples simply trail off to zero, you can almost always pad with
> zeros out to the next power of 2 and use a radix 2 FFT.
Thinking about it, that's probably the net result of what happened, except
without actually padding the data. I tried out some prime numbers just to
be sure and got the same results, but since it wasn't what I was actually
trying to do at the time I didn't try data that didn't trail off to zero.
Now that you've made me think of it I probably won't bother looking into it
at all and assume that's what happened.
Gumboot:
> I learned the FFT from a National Semiconductor databook (Application Note
> 487, Ashok Krishnamurthy). [...]
E. Robert Tisdale:
> Yes. But can you actually implement a non-radix-2 FFT?
Not a clue. I might have thought about it had I cared to implement it, I
didn't at the time and haven't found cause to consider it yet. I'm not
going to do it on a dare, either.
> How would you reorder the array for example?
Talk about array reordering bothers me slightly. I think the concept is too
pervasive for something that doesn't have as much to do with the principle
so much as the implementation. Array reordering would be the result of the
way I had figured out how to implement the algorithm, not simply the result
of trying to generalise somebody else's implementation of a radix-2.
Gumboot
> By knowing the tricks of the trade and having a knowledge of the
> literature. That is how the big boys do it E.Bob!
Maybe that's how you'd do it, but I'm not a big boy. I'm just some kid with
too much graph paper and a mechanical pencil... oh, and an eraser, we
mustn't forget the eraser.
James P. Salsman
I rename any algorithm that is named after me, because there are
just too many to keep straight:
The Refined Puckette Frequency Synthesis
The Graphic Webliography
The Surprisingly Unreflexive Modified Hartley Transform
The Perceptual Lopar Cepstrum
The Segmenting Hidden Markov Model Rematch
Lightweight MIME Content Negotiation
... and those are just some of the serial algorithms.
And of course to be fair to my detractors, some of my unperfected
techniques include:
The Regrettable Listserv Pandorization
The Reliably Pointless Spelling Rebellion
The Needlessly Simple Explosive Synthesis
um, maybe I should stop there.
Cheers,
James
--
IMS Q&TI Editor project description: http://www.bovik.org/imsqtied.html
Open-source development: http://sourceforge.net/project/?group_id=3308