The Fast Fourier Transform (FFT)
This is an archived, read-only copy of the United-TI subforum , including posts and topic from May 2003 to April 2012. If you would like to discuss any of the topics in this forum, you can visit Cemetech's TI-BASIC subforum. Some of these topics may also be directly-linked to active Cemetech topics. If you are a Cemetech member with a linked United-TI account, you can link United-TI topics here with your current Cemetech topics.

This forum is locked: you cannot post, reply to, or edit topics. TI-Basic => TI-BASIC
The Fast Fourier Transform (FFT)
United-TI Archives -> TI-Basic
 
    »
» View previous topic :: View next topic  
Author Message
thornahawk
μολών λαβέ


Active Member


Joined: 27 Mar 2005
Posts: 569
Posted: 04 Oct 2007 10:38:00 pm    Post subject:
Here, now, is a fast Fourier transform program that works on lists of any length that will fit into calculator memory, not just powers of two. The algorithm, unfortunately, is not "in place", and thus has to use a scratch array to do its work. It also needs a list of prime numbers (with the list dimensions depending on your applications) stored in ∟PRIME, which can be generated by hand or through a sieving program.

This is based on a formulation of the FFT by C. de Boor.


Code:
PROGRAM:FFT
:L₁→L₂
:dim(L₁)→L
:1→F:L→P:1→D
:Repeat P=1
:∟PRIME(D)→Q
:1→T
:If gcd(P,Q)≠1
:Then:Q→Z:P/Z→P
:Else:D+1→D
:D>dim(∟PRIME)→T
:If T:Then
:P→Z:1→P
:End:End
:If T:Then
:℮^(2πiB/Z/F)→W
:1→S
:For(M,1,Z)
:For(K,1,F)
:For(J,1,P)
:0→Y
:For(I,1,Z)
:YS+L₁(K+F(J–1+P(Z–I)))→Y
:End
:Y→L₂(K+F(M–1+Z(J–1)))
:End
:SW→S
:End:End
:For(J,1,L)
:L₂(J)→Y
:L₁(J)→L₂(J)
:Y→L₁(J)
:End
:FZ→F
:End:End
:DelVar L₂
:L₁/√(L/L^A)→L₁


There is a slight wrinkle in the usage: I wanted the program to act like Mathematica's [font="Courier"]Fourier[] command with its [font="Courier"]FourierParameters option, so one must specify the values of A and B before running the program. The list to be FFT'd is stored in L1, and is replaced by the transformed list after the program has finished running.

The usual settings for A and B in most applications are A=1 and B=-1. If the inverse transform is desired, the values of A and B must be negated.

(edit 06/09/2010: I have attached the de Boor paper: [attachment=3189:fftmul.pdf])

thornahawk

Last edited by Guest on 09 Jun 2010 09:43:16 am; edited 1 time in total
Back to top
Pseudoprogrammer


Member


Joined: 12 Dec 2006
Posts: 121
Posted: 04 Oct 2007 10:52:34 pm    Post subject:
I see some optimizations, after your delvar, don't skip to a new line, and I'm assuming you used 0->Y instead of delvar for speed reasons.
Back to top
thornahawk
μολών λαβέ


Active Member


Joined: 27 Mar 2005
Posts: 569
Posted: 04 Oct 2007 11:09:48 pm    Post subject:
Ah, yes, the program can of course be optimized by lobbing off ending parentheses and the others you have mentioned. I only elected not to in the interest of neatness. ;)

In other words, optimizing is entirely up to you when you copy it to the calculator. :)

thornahawk

Last edited by Guest on 09 Jun 2010 09:37:25 am; edited 1 time in total
Back to top
Weregoose
Authentic INTJ


Super Elite (Last Title)


Joined: 25 Nov 2004
Posts: 3976
Posted: 04 Oct 2007 11:13:06 pm    Post subject:
Well, there are those who are crazy about optimizing, and then some people recognize code that has dropped parentheses, etc. as having sloppy form that makes it less official. The point here is more of what the code does, and anyone is welcome to adjust it as it's being written into the calculator.

[EDIT] – Crosspost with thornahawk.

Last edited by Guest on 04 Oct 2007 11:15:17 pm; edited 1 time in total
Back to top
thornahawk
μολών λαβέ


Active Member


Joined: 27 Mar 2005
Posts: 569
Posted: 08 Jan 2008 08:46:08 pm    Post subject:
For reference, here is a TI-Basic implementation of the usual power-of-two in-place Cooley-Tukey FFT. Usage is the same as the program previously posted (including the setting of A and B to appropriate values).


Code:
PROGRAM:FFTCT
dim(L₁)→L
1→J
For(I,1,L−1) \\ bit-reverse permutation
If J>I:Then
L₁(J)→T
L₁(I)→L₁(J)
T→L₁(I)
End
L/2→K
While K<J
J−K→J:K/2→K
End
J+K→J
End
1→M \\ the actual FFT
While L>M
2M→V:π^r(B/M)→θ \\ ^r - radian sign
isin(θ)−2sin(θ/2)²→P
1→W
For(K,1,M)
For(I,K,L,V)
I+M→J
WL₁(J)→T
L₁(I)−T→L₁(J)
L₁(I)+T→L₁(I)
End
W+WP→W
End
V→M
End
L₁/√(L/L^A)→L₁


thornahawk

Last edited by Guest on 09 Jun 2010 09:48:20 am; edited 1 time in total
Back to top
thornahawk
μολών λαβέ


Active Member


Joined: 27 Mar 2005
Posts: 569
Posted: 09 Jun 2010 09:35:15 am    Post subject:
One can do better than the Cooley-Tukey FFT if the input is real-valued, which is often the case in applications. Edson modified the original algorithm to get rid of the redundancies in the calculations. To wit:


Code:
PROGRAM:FFTRE
:dim(L₁)→L
:1→J
:For(I,1,L-1) \\ bit-reversal permutation
:If J>I:Then
:L₁(J)→T
:L₁(I)→L₁(J)
:T→L₁(I)
:End
:L/2→K
:While K<J
:J-K→J:K/2→K
:End
:J+K→J
:End
:DelVar H:1→M \\ main FFT routine
:While M<L
:H→Q:M→H:2M→M
:For(I,1,L,M)
:I+H→K
:L₁(I)→U
:L₁(K)→V
:U+V→L₁(I)
:U-V→L₁(K)
:If Q:Then
:K+Q→K
:-L₁(K)→L₁(K)
:End:End
:π^r/H→θ \\ ^r - radian sign
:For(J,1,Q-1)
:cos(Jθ)→C
:sin(Jθ)→S
:For(I,1,L,M)
:I+J→R:I-J+H→T
:L₁(R+H)→V
:L₁(T+H)→W
:CV+SW→U
:CW-SV→V
:L₁(T)→W
:V+W→L₁(T+H)
:V-W→L₁(R+H)
:L₁(R)→W
:W-U→L₁(T)
:W+U→L₁(R)
:End:End:End
:L₁


As with FFTCT, FFTRE can only handle list lengths that are powers of two. The output is in the so-called "conjugate even" format. Unfortunately, FFTRE has no support for the A and B parameters as in FFT and FFTCT. Here is an example of how it's used, and a comparison with FFTCT (output rounded to five significant digits after the decimal point):


Code:

e^(-seq(X,X,0,15)/8)→L₁:prgmFFTRE
{7.35865, 1.0778, 0.61251, 0.519005, 0.486094, 0.471298, 0.463927, 0.460381, 0.459318, -0.085648, -0.178261, -0.28725, -0.428977, -0.638934, -1.01659, -1.97093}
e^(-seq(X,X,0,15)/8)→L₁:1→A:-1→B:prgmFFTCT
{7.35865, 1.0778-1.97093i, 0.61251-1.01659i, 0.519005-0.638934i, 0.486094-0.428977i, 0.471298-0.28725i, 0.463927-0.178261i, 0.460381-0.085648i, 0.459318, 0.460381+0.085648i, 0.463927+0.178261i, 0.471298+0.28725i, 0.486094+0.428977i, 0.519005+0.638934i, 0.61251+1.01659i, 1.0778+1.97093i}


To explain the "conjugate even" format a bit, note that for the length-16 list L₁, the first and ninth (16/2+1) list elements are real and identical in both outputs. For the second element of FFTCT's output, its real part is the second element of FFTRE's output, and its imaginary part is the 16th element of FFTRE's output (correspondingly, the 16th element of FFTCT's output holds the conjugate of the second element). The redundancy inherent in FFTCT's output is eliminated in FFTRE; instead of two list elements holding a complex number and its conjugate, the two elements hold two real numbers corresponding to the real and imaginary parts. There is a similar correspondence between the other complex elements of FFTCT and the elements of FFTRE.

The algorithm is adapted from this paper.

thornahawk
Back to top
Display posts from previous:   
Register to Join the Conversation
Have your own thoughts to add to this or any other topic? Want to ask a question, offer a suggestion, share your own programs and projects, upload a file to the file archives, get help with calculator and computer programming, or simply chat with like-minded coders and tech and calculator enthusiasts via the site-wide AJAX SAX widget? Registration for a free Cemetech account only takes a minute.

» Go to Registration page
The Fast Fourier Transform (FFT)
    »
» View previous topic :: View next topic  
Page 1 of 1 » All times are UTC - 5 Hours

 

Advertisement