restart:

# FAST FOURIER TRANSFORM ALGORITHM 8.3

# To compute the coefficients in the discrete approximation

# for the data (x(J),y(J)), 0<=J<=2m-1 where m=2^p and

# x(J)=-pi+J*pi/m for 0<=J<=2m-1.

# INPUT: m: y(0),y(1),...y(2m-1).

# OUTPUT: complex numbers c(0),...,c(2m-1): real numbers

# a(0),...,a(m): b(1),...,b(m-1).

IBR := proc(J,NU) local J1, K, JJ, J2:

J1 := J:

K := 0:

for JJ from 1 to NU do

J2 := floor(J1/2):

K := 2*K+J1-2*J2:

J1 := J2:

od:

RETURN(K):

end:

print(`This is the Fast Fourier Transform.`):

print(`The user must make provisions if the`):

print(`interval is not [-pi,pi].`):

OK := FALSE:

while OK = FALSE do

print(`Choice of input method:`):

print(`1. Input entry by entry from keyboard`):

print(`2. Input data from a text file`):

print(`3. Generate data using a function F`):

print(`Choose 1, 2, or 3 please`):

FLAG := scanf(`%d`)[1]:print(`Input is `):print(FLAG):

if FLAG = 1 or FLAG = 2 or FLAG = 3 then

OK := TRUE:

fi:

od:

if FLAG = 1 then

OK := FALSE:

while OK = FALSE do

print(`Input m\134n`):

M := scanf(`%d`)[1]:print(`m = `):print(M):

if M > 0 then

OK := TRUE:

N := 2*M:

for JJ from 1 to N do

J := JJ-1: print(`Input y(I) for I = `):print(J):

Y[JJ-1] := scanf(`%f`)[1]: print(`Input is `):print(Y[JJ-1]):

od:

else

print(`Number must be a positive integer.\134n`):

fi:

od:

fi:

if FLAG = 2 then

print(`Has a text file been created with the `):

print(`entries y(0),...,y(2m-1)\134n`):

print(`separated by a blank?\134n`):

print(`Enter 1 for yes and 2 for no`):

A := scanf(`%d`)[1]:print(`Input is `):print(A):

if A = 1 then

print(`Input the file name in the form - `):

print(`drive:\134\134name.ext\134n`):

print(`for example: A:\134\134DATA.DTA\134n`):

NAME := scanf(`%s`)[1]:

INP := fopen(NAME,READ,TEXT):

OK := FALSE:

while OK = FALSE do

print(`Input number m.\134n`):

M := scanf(`%d`)[1]:print(`Input is `):print(M):

N := 2*M:

if N > 0 then

for JJ from 1 to N do

Y[JJ-1] := fscanf(INP, `%f`)[1]:

od:

fclose(INP):

OK := TRUE:

else print(`Number must be a positive integer.\134n`):

fi:

od:

else

print(`The program will end so the input file can `):

print(`be created.\134n`):

OK := FALSE:

fi:

if FLAG = 3 then

print(`Input the function F(x) in terms of x\134n`):

print(`for example: cos(x)\134n`):

F := scanf(`%a`)[1]:print(`Input is F(x) = `):print(F):

F := unapply(F,x):

OK := FALSE:

while OK = FALSE

do print(`Input the number m.\134n`):

M := scanf(`%d`)[1]:print(`Input is `):print(M):

N := 2*M:

if N > 0 then

for JJ from 1 to N do

Z := -Pi+(JJ-1)*Pi/M:

Y[JJ-1] := evalf(F(Z)):

od:

OK := TRUE:

else print(`Number must be a postive integer.\134n`):

fi:

od:

fi:

if OK = TRUE then

print(`Choice of output method:\134n`):

print(`1. Output to screen\134n`):

print(`2. Output to text file\134n`):

print(`Please enter 1 or 2.\134n`):

FLAG := scanf(`%d`)[1]:print(`Input is `):print(FLAG):

if FLAG = 2 then

print(`Input the file name in the form - drive:\134\134name.ext\134n`):

print(`for example: A:\134\134OUTPUT.DTA\134n`):

NAME := scanf(`%s`)[1]:print(`Input file is `):print(NAME):

OUP := fopen(NAME, WRITE,TEXT):

else

OUP := default:

fi:

fprintf(OUP, `FAST FOURIER TRANSFORM\134n\134n`):

TW := ln(2):

# Step 1

# Use N2 for m, NG for p, NU1 for q, WW for zeta

N2 := floor(N/2):

# Step 2

for JJ from 1 to N do

C[JJ-1] := Y[JJ-1]:

od:

Z := N:

NG := round(ln(Z)/TW):

NU1 := NG-1:

YY := 2*Pi*I/N:

WW := exp(YY):

# Step 3

for JJ from 1 to N2 do

X := 1:

YY := 1:

for J from 1 to JJ do

YY := X*WW:

X := YY:

od:

W[JJ-1] := X:

W[N2+JJ-1] := -X:

od:

# Step 4

K := 0:

W[-1] := 1:

# Step 5

for L from 1 to NG do

# Step 6

while K < N-1 do

# Step 7

for JJ from 1 to N2 do

# Step 8

Z := exp(NU1*TW):

M1 := round(Z):

# IBR does the bit reversal

M1 := floor(K/M1):

NP := IBR(M1,NG):

# T1 is used in place of eta

T1 := C[K+N2]:

# Step 9

if NP <> 0 then

X := T1:

T1 := X*W[NP-1]:

fi:

C[K+N2] := C[K]-T1:

C[K] := C[K]+T1:

# Step 10

K := K+1:

od:

# Step 11

K := K+N2:

od:

# Step 12

K := 0:

N2 := floor(N2/2):

NU1 := NU1-1:

od:

# Step 13

while K < N-1 do

# Step 14

JJ := IBR(K,NG):

# Step 15

if JJ > K then

T3 := C[K]:

C[K] := C[JJ]:

C[JJ] := T3:

fi:

# Step 16

K := K+1:

od:

# Steps 17 and 18

fprintf(OUP, `Coefficients c(0), ... , c(2m-1)\134n\134n`):

for JJ from 1 to N do

YY := -(JJ-1)*Pi*I:

X := exp(YY):

YY := X*C[JJ-1]:

C[JJ-1] := evalc(evalf(YY/(0.5*N))):

K := JJ-1:

fprintf(OUP, `%3d %a\134n`, K, C[JJ-1]):

od:

fprintf(OUP, `\134nCoefficients a(0), ..., a(m)\134n\134n`):

for JJ from 1 to M+1 do

fprintf(OUP, `%.8f\134n`, Re(C[JJ-1])):

od:

fprintf(OUP, `\134nCoefficients b(1), ..., b(m-1)\134n\134n`):

for JJ from 2 to M do fprintf(OUP, `%.8f\134n`, Im(C[JJ-1])): od:

if OUP <> default then

fclose(OUP):

print(`Output file `,NAME,` created successfully`):

fi:

SURUaGlzfmlzfnRoZX5GYXN0fkZvdXJpZXJ+VHJhbnNmb3JtLkc2Ig==

SUVUaGV+dXNlcn5tdXN0fm1ha2V+cHJvdmlzaW9uc35pZn50aGVHNiI=

STppbnRlcnZhbH5pc35ub3R+Wy1waSxwaV0uRzYi

SThDaG9pY2V+b2Z+aW5wdXR+bWV0aG9kOkc2Ig==

SUYxLn5JbnB1dH5lbnRyeX5ieX5lbnRyeX5mcm9tfmtleWJvYXJkRzYi

ST8yLn5JbnB1dH5kYXRhfmZyb21+YX50ZXh0fmZpbGVHNiI=

SUQzLn5HZW5lcmF0ZX5kYXRhfnVzaW5nfmF+ZnVuY3Rpb25+Rkc2Ig==

STlDaG9vc2V+MSx+Mix+b3J+M35wbGVhc2VHNiI=

SSpJbnB1dH5pc35HNiI=

IiIi

SSlJbnB1dH5tfCtHNiI=

SSVtfj1+RzYi

IiIj

STRJbnB1dH55KEkpfmZvcn5Jfj1+RzYi

IiIh

SSpJbnB1dH5pc35HNiI=

JCErISkzI1IoPiEiKQ==

STRJbnB1dH55KEkpfmZvcn5Jfj1+RzYi

IiIi

SSpJbnB1dH5pc35HNiI=

JCErLEw/LXUhIio=

STRJbnB1dH55KEkpfmZvcn5Jfj1+RzYi

IiIj

SSpJbnB1dH5pc35HNiI=

JCIiIUYj

STRJbnB1dH55KEkpfmZvcn5Jfj1+RzYi

IiIk

SSpJbnB1dH5pc35HNiI=

JCIrKzZTbkMhIio=

STpDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDp8K0c2Ig==

STUxLn5PdXRwdXR+dG9+c2NyZWVufCtHNiI=

STgyLn5PdXRwdXR+dG9+dGV4dH5maWxlfCtHNiI=

STZQbGVhc2V+ZW50ZXJ+MX5vcn4yLnwrRzYi

SSpJbnB1dH5pc35HNiI=

IiIi

FAST FOURIER TRANSFORM Coefficients c(0), ... , c(2m-1) 0 -12.33700550 1 9.869604400+4.934802200*I 2 -7.402203300 3 9.869604400-4.934802200*I Coefficients a(0), ..., a(m) -12.33700550 9.86960440 -7.40220330 Coefficients b(1), ..., b(m-1) 4.93480220

JSFH