<<>> <> <> <> 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] ~ { <> <<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]]; FetchWN: PROC [lgN, N: NAT, inv: BOOL] RETURNS [wN: 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]; }.