FFTImpl.mesa
McCreight, November 14, 1984 10:59:51 am PST
Russ Atkinson (RRA) May 24, 1985 5:54:01 pm PDT
FFTImpl:
CEDAR
PROGRAM
IMPORTS Basics, Complex, IO, Random, RealFns
EXPORTS FFT = {
OPEN FFT;
pi: REAL = 3.14159276;
radiansPerDeg: REAL ← pi/180.;
TimeToFreq:
PUBLIC
PROC [ v: CVecRef ] = {
p: NAT ← 1;
exp: NAT ← 2; -- 2**p
WHILE exp<v.size
DO
exp ← 2*exp;
p ← p+1;
ENDLOOP;
IF v.size#exp THEN ERROR; -- not a power of 2
In-place bit-reversing shuffle the data
FOR i:
NAT
IN [0..v.size)
DO
j: NAT = BitRev[i, p];
IF i<j
THEN
{ t: Complex.VEC = v[i]; v[i] ← v[j]; v[j] ← t };
ENDLOOP;
Do butterfly stages
FOR dist:
NAT ← 1, 2*dist
WHILE dist<v.size
DO
dist is the number of points in this stage of DFT ... 1,2,3,..,v.size/2
w: Complex.VEC = Complex.FromPolar[r: 1.0, radians: pi/dist];
u: Complex.VEC ← [x: 1.0, y: 0.0];
FOR j:
NAT
IN [0..dist)
DO
t: Complex.VEC;
FOR i:
NAT ← j, i+2*dist
WHILE i<v.size
DO
k: NAT = i+dist;
t: Complex.VEC = Complex.Mul[v[k], u];
v[k] ← Complex.Sub[v[i], t];
v[i] ← Complex.Add[v[i], t];
ENDLOOP;
t ← Complex.Mul[u, w];
u ← t; -- eliminates premature changes in u due to compiler bug
ENDLOOP;
ENDLOOP;
};
FreqToTime:
PUBLIC
PROC [v: CVecRef ] = {
factor: REAL = 1.0/v.size;
FOR i:
NAT
IN [0..v.size)
DO
v[i] ← Complex.Conjugate[v[i]];
ENDLOOP;
TimeToFreq[v];
FOR i:
NAT
IN [0..v.size)
DO
v[i] ← Complex.Mul[[x: factor, y: 0], Complex.Conjugate[v[i]]];
ENDLOOP;
};
BitRev:
PROC [ i:
CARDINAL, width: [0..16] ← 16 ]
RETURNS [
CARDINAL ] = {
reverse the low-order width bits of i
OPEN Basics;
masks:
ARRAY [0..4)
OF
WORD =
[5555h, 3333h, 0f0fh, 00ffh];
shift: [1..16] ← 1;
c: WORD ← i;
FOR i: [0..4)
IN [0..4)
DO
c ←
BITOR[
BITSHIFT[
BITAND[c,masks[i]], shift],
BITAND[BITSHIFT[c,-shift],masks[i]]];
shift ← shift+shift;
ENDLOOP;
RETURN[BITSHIFT[c,width-16]];
};
NewCVec:
PUBLIC
PROC [ size:
NAT ← 32, totalTime: Seconds ← 1.0, like: CVecRef ←
NIL ]
RETURNS [ v: CVecRef ] = {
IF like#NIL THEN {size ← like.size; totalTime ← like.totalTime};
v ← NEW[CVec[size]];
totalTime ← totalTime;
FOR i: NAT IN [0..v.size) DO v[i] ← [x: 0, y: 0] ENDLOOP;
};
CopyCVec:
PUBLIC
PROC [ from: CVecRef ]
RETURNS [ v: CVecRef ] = {
v ← NEW[CVec[from.size]];
v.totalTime ← from.totalTime;
FOR i: NAT IN [0..v.size) DO v[i] ← from[i] ENDLOOP;
};
AddCosWave:
PUBLIC
PROC [v: CVecRef, freq: Hz, magnitude: Volts ← 1.0, phase: Degrees ← 0.0 ] = {
FOR i:
NAT
IN [0..v.size)
DO
v[i].x ← v[i].x+magnitude*RealFns.CosDeg[360.0*i*freq*v.totalTime/v.size+phase];
ENDLOOP;
};
AddNoise:
PUBLIC
PROC [v: CVecRef, powerPerHz:
DBM ← 1.0, ohms:
REAL ← 50.0 ] = {
noise: CVecRef = NewCVec[like: v];
bwPerStep: Hz = noise.size/noise.totalTime;
voltagePerStep: Volts = RealFns.SqRt[ohms*bwPerStep*DBMToWatts[powerPerHz]];
FOR i:
NAT
IN [0..noise.size)
DO
noise[i] ← Complex.FromPolar[
r: voltagePerStep,
radians: radiansPerDeg*Random.ChooseInt[rs, 0, 359]];
ENDLOOP;
FreqToTime[noise];
FOR i:
NAT
IN [0..v.size)
DO
v[i].x ← v[i].x+noise[i].x;
ENDLOOP;
};
DBMToWatts:
PUBLIC
PROC [ dbm:
DBM ]
RETURNS [ watts: Watts ] = {
RETURN[RealFns.Power[base: 10.0, exponent: (dbm-30)/10.0]];
};
WattsToDBM:
PUBLIC
PROC [ watts: Watts ]
RETURNS [ dbm:
DBM ] = {
RETURN[10.0*RealFns.Log[base: 10.0, arg: watts]+30.0];
};
Print:
PUBLIC
PROC [ v: CVecRef ]
RETURNS [ Rope.
ROPE ] = {
s: IO.STREAM = IO.ROS[];
FOR i:
NAT
IN [0..v.size)
DO
IF i>0 THEN s.PutRope[", "];
s.PutF[format: "%d: %g+j%g", v1: IO.int[i], v2: IO.real[v[i].x], v3: IO.real[v[i].y]];
ENDLOOP;
RETURN[IO.RopeFromROS[s]];
};
rs: Random.RandomStream ← Random.Create[];