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] ~ { 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]; }. ζ FFTSingleImpl.mesa Copyright Σ 1991, 1992 by Xerox Corporation. All rights reserved. Spreitze, January 25, 1991 10:09 am PST Do an FFT step on the N-element sub-array of cs at l, l+k, l+2k, ... 0 <= l < k; Nk = cs.length. ΚΖ–(cedarcode) style•NewlineDelimiter ™code™Kšœ Οeœ7™BK™'—K˜KšΟk œ žœ˜2K˜šΟn œžœž˜Kšžœž˜Kšžœ ˜—Kšœžœ ˜K˜š Ÿ œžœžœžœžœ˜CKšœžœ˜!K˜Kšžœ˜—K˜š Ÿœžœžœžœžœ˜KKšžœžœžœ˜MK˜Kšœžœ˜Kšžœ˜ —K˜š Ÿœžœžœžœžœ˜TKšœžœ˜"Kšžœžœ žœ˜MKš žœžœžœžœžœ˜>Kšžœ˜—K˜š Ÿœžœžœžœžœ˜PKšŸœžœ ˜Kšœžœ˜ Kšœžœ˜ Kšžœžœžœ˜K˜šžœ˜šžœ˜K˜+Kšœ"˜"Kšžœžœžœ˜—Kšžœ˜—Kšžœ˜ —K˜šŸœžœ%žœžœ˜CK™DK™Kšœžœ˜Kšœžœ˜Kšœžœ˜KšœΟdœ œ˜&KšœžœΟc˜Kšœžœ‘˜!Kš œ œ Οuœ œžœ˜šžœž˜˜Kšœžœ˜Kšœžœ ˜Kšœžœ˜Kšœžœ ˜Kšœžœ ˜Kšœžœ ˜Kšœžœ ˜K˜K˜K˜K˜—šœ˜Kšœžœ˜.Kšœžœ˜2šžœžœ žœž˜0K˜Kšžœ˜—Kšœ˜—K˜5Kšžœžœ˜—Kšœ œ ˜ Kšœ ’œ œ ˜Kšœ œ ’œ œ˜Kšœ œ ’œ œ˜&šžœžœžœž˜K˜K˜Kšœ œ ˜Kšœ ’œ œ œ˜*Kšœ œ ’œ œ˜!Kšœ œ ’œ œ˜(Kšžœ˜—Kšžœ˜—K˜Kšœžœ žœžœžœžœžœžœžœ˜=K˜šŸΠdnœžœŸœžœžœžœ œ˜Cš žœžœžœžœž˜)˜Kšœ œ˜Kšœ œ˜Kšœ œ ˜Kšœ œ˜—˜Kšœ œ˜Kšœ œ˜Kšœ œ˜Kšœ œ ˜Kšœ œ˜—˜Kšœ œ˜Kšœ œ˜Kš œ œ žœžœžœ˜)Kšœ œ˜Kš œ œ žœžœžœ˜)Kšœ œ ˜Kšœ œ˜—šžœ˜ Kšœ œ˜*Kšœ œ˜Kšœ œ˜Kšœ œ˜šžœžœžœ ž˜Kšœ œ˜Kšœ œ œ˜&Kšžœ˜—Kšœ œ ˜Kšœ œ˜——Kšžœ˜—K˜Kšœžœ˜"K˜šŸœžœžœ˜-K˜+KšŸœžœ˜Kš œžœžœžœžœ ˜0šžœžœžœ ž˜Kšœžœ ˜Kšœžœ ˜šžœžœžœ ž˜Kšœžœ˜Kšœ ’œžœ˜*Kšœ, ’œ˜2Kšžœ˜—K˜ Kšžœ˜—Kšžœžœžœ˜Kšžœ˜—K˜šŸ œžœ˜$KšŸœžœ ˜šžœžœžœž˜K˜#Kšžœ˜—Kšžœ˜—K˜šŸœžœžœžœžœžœžœ˜XKš žœžœžœžœžœ˜