alg083.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> 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).
> alg083 := proc() global IBR; local OK, FLAG, M, N, JJ, J, Y, A, NAME, INP, F, Z, OUP, TW, N2, C, NG, NU1, YY, WW, X, W, K, L, M1, NP, T1, T3;
> 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;
> printf(`This is the Fast Fourier Transform.\n\n`);
> printf(`The user must make provisions if the\n`);
> printf(`interval is not [-pi,pi].\n`);
> OK := FALSE;
> while OK = FALSE do
> printf(`Choice of input method:\n`);
> printf(`1. Input entry by entry from keyboard\n`);
> printf(`2. Input data from a text file\n`);
> printf(`3. Generate data using a function F\n`);
> printf(`Choose 1, 2, or 3 please\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 1 or FLAG = 2 or FLAG = 3 then
> OK := TRUE;
> fi;
> od;
> if FLAG = 1 then
> OK := FALSE;
> while OK = FALSE do
> printf(`Input m\n`);
> M := scanf(`%d`)[1];
> if M > 0 then
> OK := TRUE;
> N := 2*M;
> for JJ from 1 to N do
> J := JJ-1;
> printf(`Input y(%d).\n`, J);
> Y[JJ-1] := scanf(`%f`)[1];
> od;
> else
> printf(`Number must be a positive integer.\n`);
> fi;
> od;
> fi;
> if FLAG = 2 then
> printf(`Has a text file been created with the `);
> printf(`entries y(0),...,y(2m-1)\n`);
> printf(`separated by a blank?\n`);
> printf(`Enter Y or N\n`);
> A := scanf(`\n%c`)[1];
> if A = "Y" or A = "y" then
> printf(`Input the file name in the form - `);
> printf(`drive:\\name.ext\n`);
> printf(`for example:   A:\\DATA.DTA\n`);
> NAME := scanf(`%s`)[1];
> INP := fopen(NAME,READ,TEXT);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input number m.\n`);
> M := scanf(`%d`)[1];
> 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 printf(`Number must be a positive integer.\n`);
> fi;
> od;
> else
> printf(`The program will end so the input file can `);
> printf(`be created.\n`);
> OK := FALSE;
> fi;
> fi;
> if FLAG = 3 then
> printf(`Input the function F(x) in terms of x\n`);
> printf(`for example:   cos(x)\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,x);
> OK := FALSE;
> while OK = FALSE
> do printf(`Input the number m.\n`);
> M := scanf(`%d`)[1];
> 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 printf(`Number must be a postive integer.\n`);
> fi;
> od; 
> fi;
> if OK = TRUE then
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:   A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> 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):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg083();