FFTSingleImpl.mesa
Copyright Ó 1991, 1992 by Xerox Corporation. All rights reserved.
Spreitze, January 25, 1991 10:09 am PST
DIRECTORY Commander, Complex, FFTSingle, IO, Rope;
FFTSingleImpl:
CEDAR
PROGRAM
IMPORTS Commander, Complex, IO
EXPORTS FFTSingle
= {OPEN FFTSingle;
NewComplexSeq:
PUBLIC
PROC [size:
NAT]
RETURNS [cs: ComplexSeq] ~ {
cs ¬ NEW [ComplexSequence[size]];
cs.length ¬ 0;
RETURN};
CsAppend:
PUBLIC
PROC [cs: ComplexSeq, c:
COMPLEX]
RETURNS [ComplexSeq] ~ {
IF cs.length >= cs.size THEN cs ¬ CsCopy[cs, 0, cs.size + MAX[cs.size/2, 1]];
cs[cs.length] ¬ c;
cs.length ¬ cs.length.SUCC;
RETURN [cs]};
CsCopy:
PUBLIC
PROC [cs: ComplexSeq, start, size:
INT]
RETURNS [new: ComplexSeq] ~ {
new ¬ NEW [ComplexSequence[size]];
IF start < cs.length THEN new.length ¬ cs.length - start ELSE new.length ¬ 0;
FOR i: NAT IN [0..new.length) DO new[i] ¬ cs[start+i] ENDLOOP;
RETURN};
DestructiveFft:
PUBLIC
PROC [cs: ComplexSeq, inv:
BOOL]
RETURNS [ComplexSeq] ~ {
N: NAT ~ cs.length;
lgN: NAT;
exact: BOOL;
IF N<2 THEN RETURN [cs];
[lgN, exact] ¬ Lg[N];
IF exact
THEN {
alt: ComplexSeq ¬ NewComplexSeq[cs.length];
DStep[cs, alt, lgN, N, 1, 0, inv];
IF NOT inv THEN DivideByN[cs]}
ELSE DumbDft[cs, inv];
RETURN [cs]};
DStep:
PROC [cs, alt: ComplexSeq, lgN, N, k, l:
NAT, inv:
BOOL] ~ {
Do an FFT step on the N-element sub-array of cs at l, l+k, l+2k, ...
0 <= l < k; Nk = cs.length.
halfN: NAT ~ N/2;
halfKN: NAT ~ cs.length/2;
twoK: NAT ~ k+k;
wN: ComplexSeq ~ FetchWN[lgN, N, inv];
jkl: NAT ¬ l; --will be j*k + l
j2kl: NAT ¬ l; --will be j*2k + l
Bj, wNjCj: COMPLEX;
SELECT halfN
FROM
=2 => {
i1: NAT ~ l+k;
i2: NAT ~ l+twoK;
i3: NAT ~ i2+k;
c0: COMPLEX ¬ cs[l];
c1: COMPLEX ¬ cs[i1];
c2: COMPLEX ¬ cs[i2];
c3: COMPLEX ¬ cs[i3];
alt[l] ¬ Complex.Add[c0, c2];
alt[i2] ¬ Complex.Sub[c0, c2];
alt[i1] ¬ Complex.Add[c1, c3];
alt[i3] ¬ Complex.Sub[c1, c3]};
>2 => {
DStep[cs, alt, lgN.PRED, halfN, twoK, l, inv];
DStep[cs, alt, lgN.PRED, halfN, twoK, l + k, inv];
FOR idx:
NAT ¬ l, idx+k
WHILE idx < cs.length
DO
alt[idx] ¬ cs[idx];
ENDLOOP;
};
=1 => {alt[l] ¬ cs[l]; alt[l+halfKN] ¬ cs[l+halfKN]};
ENDCASE => ERROR;
Bj ¬ alt[l];
wNjCj ¬ alt[l+k];
cs[l] ¬ Complex.Add[Bj, wNjCj];
cs[l+halfKN] ¬ Complex.Sub[Bj, wNjCj];
FOR j:
NAT
IN (0 .. halfN)
DO
jkl ¬ jkl + k;
j2kl ¬ j2kl + twoK;
Bj ¬ alt[j2kl];
wNjCj ¬ Complex.Mul[wN[j], alt[j2kl + k]];
cs[jkl] ¬ Complex.Add[Bj, wNjCj];
cs[jkl+halfKN] ¬ Complex.Sub[Bj, wNjCj];
ENDLOOP;
RETURN};
w: ARRAY [0..32) OF ARRAY BOOL OF ComplexSeq ¬ ALL[ALL[NIL]];
FetchW
N:
PROC [lgN,
N:
NAT, inv:
BOOL]
RETURNS [w
N: ComplexSeq] ~ {
IF w[lgN][inv] =
NIL
THEN
SELECT lgN
FROM
0 => {
wN ¬ NewComplexSeq[1];
wN[0] ¬ [1.0, 0.0];
wN.length ¬ 1;
w[lgN][inv] ¬ wN};
1 => {
wN ¬ NewComplexSeq[2];
wN[0] ¬ [1.0, 0.0];
wN[1] ¬ [-1.0, 0.0];
wN.length ¬ 2;
w[lgN][inv] ¬ wN};
2 => {
wN ¬ NewComplexSeq[4];
wN[0] ¬ [1.0, 0.0];
wN[1] ¬ [0.0, IF inv THEN 1.0 ELSE -1.0];
wN[2] ¬ [-1.0, 0.0];
wN[3] ¬ [0.0, IF inv THEN -1.0 ELSE 1.0];
wN.length ¬ 4;
w[lgN][inv] ¬ wN};
ENDCASE => {
wh: ComplexSeq ~ FetchWN[lgN-1, N/2, inv];
wN ¬ NewComplexSeq[N];
wN[0] ¬ [1.0, 0.0];
wN[1] ¬ Complex.SqRt[wh[1]];
FOR j:
NAT
IN (0..N/2)
DO
wN[2*j] ¬ wh[j];
wN[2*j+1] ¬ Complex.Mul[wh[j], wN[1]];
ENDLOOP;
wN.length ¬ N;
w[lgN][inv] ¬ wN};
RETURN [w[lgN][inv]]};
twoPi: REAL ¬ 2*3.141592653589793;
DumbDft:
PROC [cs: ComplexSeq, inv:
BOOL] ~ {
org: ComplexSeq ~ CsCopy[cs, 0, cs.length];
N: NAT ~ org.length;
dda: REAL ~ (IF inv THEN twoPi ELSE -twoPi) / N;
FOR n:
NAT
IN [0 .. N)
DO
sum: COMPLEX ¬ [0, 0];
da: REAL ~ n*dda;
FOR m:
NAT
IN [0 .. N)
DO
a: REAL ~ da*m;
wNnm: COMPLEX ~ Complex.FromPolar[1.0, a];
sum ¬ Complex.Add[sum, Complex.Mul[org[m], wNnm]];
ENDLOOP;
cs[n] ¬ sum;
ENDLOOP;
IF NOT inv THEN DivideByN[cs];
RETURN};
DivideByN:
PROC [cs: ComplexSeq] ~ {
N: NAT ~ cs.length;
FOR j:
NAT
IN [0..N)
DO
cs[j] ¬ Complex.Div[cs[j], [N, 0]];
ENDLOOP;
RETURN};
DestructiveMap:
PUBLIC
PROC [cs: ComplexSeq, proc:
PROC [
COMPLEX]
RETURNS [
COMPLEX]] ~ {
FOR i: NAT IN [0..cs.length) DO cs[i] ¬ proc[cs[i]] ENDLOOP;
RETURN};
Lg:
PROC [
N:
NAT]
RETURNS [lg:
NAT, exact:
BOOL] ~ {
IF N=0 THEN ERROR;
IF N=1 THEN RETURN [0, TRUE];
{halfN: NAT ~ N/2;
[lg, exact] ¬ Lg[halfN];
exact ¬ exact AND N=2*halfN;
lg ¬ lg.SUCC;
RETURN}};
Test:
PROC [cmd: Commander.Handle]
RETURNS [result:
REF
ANY ¬
NIL, msg: Rope.
ROPE ¬
NIL] ~ {
N: NAT ~ IO.RIS[cmd.commandLine].GetInt[];
cs: ComplexSeq ¬ NewComplexSeq[N];
FOR i:
NAT
IN [0..N)
DO
x: REAL ¬ cmd.in.GetReal[];
y: REAL ¬ cmd.in.GetReal[];
IF cs # CsAppend[cs, [x, y]] THEN ERROR;
ENDLOOP;
SELECT cmd.procData.clientData
FROM
$FFT => cs ¬ DestructiveFft[cs, FALSE];
$InvFFT => cs ¬ DestructiveFft[cs, TRUE];
$Dumb => DumbDft[cs, FALSE];
$InvDumb => DumbDft[cs, TRUE];
ENDCASE => ERROR;
cs ¬ cs;
FOR i:
NAT
IN [0 .. N)
DO
cmd.out.PutF["%g %g\n", [real[cs[i].x]], [real[cs[i].y]] ];
ENDLOOP;
RETURN};
Commander.Register["DFT", Test, "N - order N DFT (filter)", $FFT];
Commander.Register["InvDFT", Test, "N - order N DFT (filter)", $InvFFT];
Commander.Register["DumbDFT", Test, "N - order N dumb DFT (filter)", $Dumb];
Commander.Register["InvDumbDFT", Test, "N - order N dumb DFT (filter)", $InvDumb];
}.