unit Wiskunde; { This program has been designed and is CopyLefted by: * Han de Bruijn; Systems for Research "A little bit of Physics would be (===) * and Education (DTO/SOO), Mekelweg 6 NO Idleness in Mathematics"(HdB) @-O^O-@ * 2628 CD Delft, The Netherlands http://huizen.dto.tudelft.nl/deBruijn #/_\# * E-mail: Han.deBruijn@DTO.TUDelft.NL Tel: +31 15 27 82751. Fax: 81722 ### All I ask is to be credited when it is appropriate. To be used (with Delphi) as follows: program testje; uses Wiskunde; begin FFTest1D(false); FFTest2D; end. } INTERFACE procedure FFTest1D(blah : boolean); procedure FFTest2D; IMPLEMENTATION type complex = record Re : double; Im : double; end; type vektor = array of double; complex_vektor = array of complex; var beeld : array of array of complex; function plus(z1, z2 : complex) : complex; { Add } begin plus.Re := z1.Re + z2.Re; plus.Im := z1.Im + z2.Im; end; function min(z1, z2 : complex) : complex; { Substract } begin min.Re := z1.Re - z2.Re; min.Im := z1.Im - z2.Im; end; function maal(z1, z2 : complex) : complex; { Multiply } begin maal.Re := z1.Re * z2.Re - z1.Im * z2.Im; maal.Im := z1.Re * z2.Im + z1.Im * z2.Re; end; function deel(z1, z2 : complex) : complex; { Divide } var onder : double; begin onder := z2.Re * z2.Re + z2.Im * z2.Im; deel.Re := (z1.Re * z2.Re + z1.Im * z2.Im) / onder; deel.Im := (z1.Im * z2.Re - z1.Re * z2.Im) / onder; end; function absoluut(z : complex) : double; { Absolute value } begin absoluut := sqrt(sqr(z.Re) + sqr(z.Im)); end; procedure Drukaf(x : vektor); { Print } var orde,k : integer; begin orde := Length(x); for k := 0 to orde-1 do Write(' ',x[k]:9:5); Writeln; end; procedure Scherm(W,H : integer); { Print } var k,L : integer; begin for L := 0 to H-1 do begin for k := 0 to W-1 do Write(' ',Round(absoluut(beeld[k][L]))); Writeln; end; Writeln; end; function tweede(macht : integer) : integer; { Compute 2^macht } var k : integer; h : integer; begin h := 1; for k := 1 to macht do h := h shl 1; tweede := h; end; function log2(n : integer) : integer; { Logarithm base 2 } var k, L : integer; begin log2 := 0; if n <= 0 then Exit; for k := 0 to 32 do begin n := n div 2 ; L := k; if n = 0 then Break ; end; log2 := L; end; procedure Halen(rij : vektor; var row : complex_vektor); { Get Real and store into Complex } var i,orde,g : integer; w : complex; begin orde := Length(rij); for i := 1 to orde do begin { Real part symmetric around g Values at the Left stored } g := (orde div 2) + 1; if i <= g then w.Re := rij[i-1]; if i > g then w.Re := rij[orde-i+1]; { Imaginary part anti-symmetric Values at the Right stored } if i > g then w.Im := rij[i-1]; if i < g then w.Im := - rij[orde-i+1]; if (i = g) or (i = 0) then w.Im := 0; row[i-1] := w; end; end; procedure Brengen(row : complex_vektor; var rij : vektor); { Get Complex and store into Real } var i,orde,g : integer; w : complex; begin orde := Length(rij); for i := 1 to orde do begin w := row[i-1]; { Real part symmetric around g Store values at the Left } g := (orde div 2) + 1; if i <= g then rij[i-1] := w.Re; { Imaginary part anti-symmetric Store values at the Right } if i > g then rij[i-1] := w.Im; end; end; procedure Schudden(x : complex_vektor); { Bitreverse shuffling -------------------- } var t : complex; orde,i,j : integer; function kount(j : integer) : integer; { Reverse bit count } var k : integer; begin k := orde shr 1; while (k > 0) and (k < j) do begin j := j - k; k := k shr 1; end; j := j + k; kount := j; end; begin orde := Length(x); j := 1; for i := 2 to orde-1 do begin j := kount(j); if i > j then Continue; t := x[i-1]; x[i-1] := x[j-1]; x[j-1] := t; end; end; procedure FFT(var x : complex_vektor ; ifb : integer); { + 1 : forward Fast Fourier Transform - 1 : backward ---------------------- } var u,w,t : complex; le,lt, i,j,ip : integer; rang,orde : integer; arg : double; L : integer; begin Schudden(x); orde := Length(x); rang := log2(orde); lt := 1; for L := 1 to rang do begin le := lt; lt := lt * 2; u.Re := 1 ; u.Im := 0; arg := -ifb*Pi/le; w.Re := cos(arg); w.Im := sin(arg); for j := 1 to le do begin i := j; while i <= orde do begin ip := i+le; t := maal(x[ip-1],u); x[ip-1] := min(x[i-1],t); x[i-1] := plus(x[i-1],t); i := i + lt; end; u := maal(u,w); end; end; end; procedure FFTest1D; { Test Program 1-D } const test : integer = 27; var orde,k : integer; a,b : vektor; z : complex_vektor; faktor : double; begin orde := tweede(log2(test-1)+1); SetLength(a,orde); SetLength(b,orde); SetLength(z,orde); for k := 0 to test-1 do a[k] := Random; for k := test to orde-1 do a[k] := 0; Drukaf(a); { Fourier Transforms: } for k := 0 to orde-1 do begin z[k].Re := a[k]; z[k].Im := 0; end; FFT(z,1); for k := 0 to orde-1 do b[k] := z[k].Re; if blah then Drukaf(b); for k := 0 to orde-1 do b[k] := z[k].Im; if blah then Drukaf(b); Brengen(z,b); if blah then Drukaf(b); for k := 0 to orde-1 do z[k].Re := 0; for k := 0 to orde-1 do z[k].Im := 0; Halen(b,z); for k := 0 to orde-1 do b[k] := z[k].Re; if blah then Drukaf(b); for k := 0 to orde-1 do b[k] := z[k].Im; if blah then Drukaf(b); faktor := 1/orde; for k := 0 to orde-1 do begin z[k].Re := z[k].Re *faktor; z[k].Im := z[k].Im *faktor; end; FFT(z,-1); for k := 0 to orde-1 do b[k] := z[k].Re; Drukaf(b); end; procedure FFT2D(Ox,Oy,ibl : integer); { Two dimensional Fourier Transform } var k,L : integer; z : complex_vektor; begin SetLength(z,Oy); for k := 0 to Ox-1 do begin { In Y-direction } for L := 0 to Oy-1 do begin z[L].Re := beeld[k][L].Re; z[L].Im := beeld[k][L].Im; end; FFT(z,ibl); for L := 0 to Oy-1 do begin beeld[k][L].Re := z[L].Re; beeld[k][L].Im := z[L].Im; end; end; SetLength(z,Oy); for L := 0 to Oy-1 do begin { In X-direction } for k := 0 to Ox-1 do begin z[k].Re := beeld[k][L].Re; z[k].Im := beeld[k][L].Im; end; FFT(z,ibl); for k := 0 to Ox-1 do begin beeld[k][L].Re := z[k].Re; beeld[k][L].Im := z[k].Im; end; end; end; procedure FFTest2D; { Test Program 2-D } const tX : integer = 30; tY : integer = 20; var k,L,Ox,Oy : integer; faktor : double; begin Ox := tweede(log2(tX-1)+1); Oy := tweede(log2(tY-1)+1); SetLength(beeld,Ox,Oy); for k := 0 to Ox-1 do for L := 0 to Oy-1 do begin beeld[k][L].Re := 0; beeld[k][L].Im := 0; end; for k := 0 to tX-1 do for L := 0 to tY-1 do if Random > 0.5 then beeld[k][L].Re := 1; Scherm(tX,tY); FFT2D(Ox,Oy,1); { Scherm(Ox,Oy); } faktor := 1/(Ox*Oy); for k := 0 to Ox-1 do for L := 0 to Oy-1 do begin beeld[k][L].Re := beeld[k][L].Re *faktor; beeld[k][L].Im := beeld[k][L].Im *faktor; end; FFT2D(Ox,Oy,-1); Scherm(tX,tY); end; END.