> 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.\n`):
> 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\n`):
> 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.\n`):
> fi:
> od:
> fi:
> if FLAG = 2 then
> print(`Has a text file been created with the `):
> print(`entries y(0),...,y(2m-1)\n`):
> print(`separated by a blank?\n`):
> 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:\\name.ext\n`):
> print(`for example:   A:\\DATA.DTA\n`):
> NAME := scanf(`%s`)[1]:
> INP := fopen(NAME,READ,TEXT):
> OK := FALSE:
> while OK = FALSE do
> print(`Input number m.\n`):
> 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.\n`):
> fi:
> od:
> else
> print(`The program will end so the input file can `):
> print(`be created.\n`):
> OK := FALSE:
> fi:
> fi:
> if FLAG = 3 then
> print(`Input the function F(x) in terms of x\n`):
> print(`for example:   cos(x)\n`):
> 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.\n`):
> 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.\n`):
> fi:
> od: 
> fi:
> if OK = TRUE then
> print(`Choice of output method:\n`):
> print(`1. Output to screen\n`):
> print(`2. Output to text file\n`):
> print(`Please enter 1 or 2.\n`):
> FLAG := scanf(`%d`)[1]:print(`Input is `):print(FLAG):
> if FLAG = 2 then
> print(`Input the file name in the form - drive:\\name.ext\n`):
> print(`for example:   A:\\OUTPUT.DTA\n`):
> NAME := scanf(`%s`)[1]:print(`Input file is `):print(NAME):
> OUP := fopen(NAME, WRITE,TEXT):
> else
> OUP := default:
> fi:
> fprintf(OUP, `FAST FOURIER TRANSFORM\n\n`):
> 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)\n\n`):
> 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\n`, K, C[JJ-1]):
> od:
> fprintf(OUP, `\nCoefficients a(0), ..., a(m)\n\n`):
> for JJ from 1 to M+1 do
> fprintf(OUP, `%.8f\n`, Re(C[JJ-1])):
> od:
> fprintf(OUP, `\nCoefficients b(1), ..., b(m-1)\n\n`):
> for JJ from 2 to M do fprintf(OUP, `%.8f\n`, Im(C[JJ-1])): od:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi: