(***************************************************************************** DIGITAL SIGNAL PROCESSING TOOLS Version 1.03, 2001/06/15 (c) 1999 - Laurent de Soras FFTReal.h Fourier transformation of real number arrays. Portable ISO C++ ------------------------------------------------------------------------------ LEGAL Source code may be freely used for any purpose, including commercial applications. Programs must display in their "About" dialog-box (or documentation) a text telling they use these routines by Laurent de Soras. Modified source code can be distributed, but modifications must be clearly indicated. CONTACT Laurent de Soras 92 avenue Albert 1er 92500 Rueil-Malmaison France ldesoras@club-internet.fr ------------------------------------------------------------------------------ Translation to ObjectPascal by : Frederic Vanmol frederic@axiworld.be *****************************************************************************) unit FFTReal; interface uses Windows; (* Change this typedef to use a different floating point type in your FFTs (i.e. float, double or long double). *) type pflt_t = ^flt_t; flt_t = single; pflt_array = ^flt_array; flt_array = array[0..0] of flt_t; plongarray = ^longarray; longarray = array[0..0] of longint; const sizeof_flt : longint = SizeOf(flt_t); type // Bit reversed look-up table nested class TBitReversedLUT = class private _ptr : plongint; public constructor Create(const nbr_bits: integer); destructor Destroy; override; function get_ptr: plongint; end; // Trigonometric look-up table nested class TTrigoLUT = class private _ptr : pflt_t; public constructor Create(const nbr_bits: integer); destructor Destroy; override; function get_ptr(const level: integer): pflt_t; end; TFFTReal = class private _bit_rev_lut : TBitReversedLUT; _trigo_lut : TTrigoLUT; _sqrt2_2 : flt_t; _length : longint; _nbr_bits : integer; _buffer_ptr : pflt_t; public constructor Create(const length: longint); destructor Destroy; override; procedure do_fft(f: pflt_array; const x: pflt_array); procedure do_ifft(const f: pflt_array; x: pflt_array); procedure rescale(x: pflt_array); end; implementation uses Math; { TBitReversedLUT } constructor TBitReversedLUT.Create(const nbr_bits: integer); var length : longint; cnt : longint; br_index : longint; bit : longint; begin inherited Create; length := 1 shl nbr_bits; GetMem(_ptr, length*SizeOf(longint)); br_index := 0; plongarray(_ptr)^[0] := 0; for cnt := 1 to length-1 do begin // ++br_index (bit reversed) bit := length shr 1; br_index := br_index xor bit; while br_index and bit = 0 do begin bit := bit shr 1; br_index := br_index xor bit; end; plongarray(_ptr)^[cnt] := br_index; end; end; destructor TBitReversedLUT.Destroy; begin FreeMem(_ptr); _ptr := nil; inherited; end; function TBitReversedLUT.get_ptr: plongint; begin Result := _ptr; end; { TTrigLUT } constructor TTrigoLUT.Create(const nbr_bits: integer); var total_len : longint; PI : double; level : integer; level_len : longint; level_ptr : pflt_array; mul : double; i : longint; begin inherited Create; _ptr := nil; if (nbr_bits > 3) then begin total_len := (1 shl (nbr_bits - 1)) - 4; GetMem(_ptr, total_len * sizeof_flt); PI := ArcTan(1) * 4; for level := 3 to nbr_bits-1 do begin level_len := 1 shl (level - 1); level_ptr := pointer(get_ptr(level)); mul := PI / (level_len shl 1); for i := 0 to level_len-1 do level_ptr^[i] := cos(i * mul); end; end; end; destructor TTrigoLUT.Destroy; begin FreeMem(_ptr); _ptr := nil; inherited; end; function TTrigoLUT.get_ptr(const level: integer): pflt_t; var tempp : pflt_t; begin tempp := _ptr; inc(tempp, (1 shl (level-1)) - 4); Result := tempp; end; { TFFTReal } constructor TFFTReal.Create(const length: longint); begin inherited Create; _length := length; _nbr_bits := Floor(Ln(length) / Ln(2) + 0.5); _bit_rev_lut := TBitReversedLUT.Create(Floor(Ln(length) / Ln(2) + 0.5)); _trigo_lut := TTrigoLUT.Create(Floor(Ln(length) / Ln(2) + 0.05)); _sqrt2_2 := Sqrt(2) * 0.5; _buffer_ptr := nil; if _nbr_bits > 2 then GetMem(_buffer_ptr, _length * sizeof_flt); end; destructor TFFTReal.Destroy; begin if _buffer_ptr <> nil then begin FreeMem(_buffer_ptr); _buffer_ptr := nil; end; _bit_rev_lut.Free; _bit_rev_lut := nil; _trigo_lut.Free; _trigo_lut := nil; inherited; end; (*==========================================================================*/ /* Name: do_fft */ /* Description: Compute the FFT of the array. */ /* Input parameters: */ /* - x: pointer on the source array (time). */ /* Output parameters: */ /* - f: pointer on the destination array (frequencies). */ /* f [0...length(x)/2] = real values, */ /* f [length(x)/2+1...length(x)-1] = imaginary values of */ /* coefficents 1...length(x)/2-1. */ /*==========================================================================*) procedure TFFTReal.do_fft(f: pflt_array; const x: pflt_array); var sf, df : pflt_array; pass : integer; nbr_coef : longint; h_nbr_coef : longint; d_nbr_coef : longint; coef_index : longint; bit_rev_lut_ptr : plongarray; rev_index_0 : longint; rev_index_1 : longint; rev_index_2 : longint; rev_index_3 : longint; df2 : pflt_array; n1, n2, n3 : integer; sf_0, sf_2 : flt_t; sqrt2_2 : flt_t; v : flt_t; cos_ptr : pflt_array; i : longint; sf1r, sf2r : pflt_array; dfr, dfi : pflt_array; sf1i, sf2i : pflt_array; c, s : flt_t; temp_ptr : pflt_array; b_0, b_2 : flt_t; begin n1 := 1; n2 := 2; n3 := 3; (*______________________________________________ * * General case *______________________________________________ *) if _nbr_bits > 2 then begin if _nbr_bits and 1 <> 0 then begin df := pointer(_buffer_ptr); sf := f; end else begin df := f; sf := pointer(_buffer_ptr); end; // // Do the transformation in several passes // // First and second pass at once bit_rev_lut_ptr := pointer(_bit_rev_lut.get_ptr); coef_index := 0; repeat rev_index_0 := bit_rev_lut_ptr^[coef_index]; rev_index_1 := bit_rev_lut_ptr^[coef_index + 1]; rev_index_2 := bit_rev_lut_ptr^[coef_index + 2]; rev_index_3 := bit_rev_lut_ptr^[coef_index + 3]; df2 := pointer(longint(df) + (coef_index*sizeof_flt)); df2^[n1] := x^[rev_index_0] - x^[rev_index_1]; df2^[n3] := x^[rev_index_2] - x^[rev_index_3]; sf_0 := x^[rev_index_0] + x^[rev_index_1]; sf_2 := x^[rev_index_2] + x^[rev_index_3]; df2^[0] := sf_0 + sf_2; df2^[n2] := sf_0 - sf_2; inc(coef_index, 4); until (coef_index >= _length); // Third pass coef_index := 0; sqrt2_2 := _sqrt2_2; repeat sf^[coef_index] := df^[coef_index] + df^[coef_index + 4]; sf^[coef_index + 4] := df^[coef_index] - df^[coef_index + 4]; sf^[coef_index + 2] := df^[coef_index + 2]; sf^[coef_index + 6] := df^[coef_index + 6]; v := (df [coef_index + 5] - df^[coef_index + 7]) * sqrt2_2; sf^[coef_index + 1] := df^[coef_index + 1] + v; sf^[coef_index + 3] := df^[coef_index + 1] - v; v := (df^[coef_index + 5] + df^[coef_index + 7]) * sqrt2_2; sf [coef_index + 5] := v + df^[coef_index + 3]; sf [coef_index + 7] := v - df^[coef_index + 3]; inc(coef_index, 8); until (coef_index >= _length); // Next pass for pass := 3 to _nbr_bits-1 do begin coef_index := 0; nbr_coef := 1 shl pass; h_nbr_coef := nbr_coef shr 1; d_nbr_coef := nbr_coef shl 1; cos_ptr := pointer(_trigo_lut.get_ptr(pass)); repeat sf1r := pointer(longint(sf) + (coef_index * sizeof_flt)); sf2r := pointer(longint(sf1r) + (nbr_coef * sizeof_flt)); dfr := pointer(longint(df) + (coef_index * sizeof_flt)); dfi := pointer(longint(dfr) + (nbr_coef * sizeof_flt)); // Extreme coefficients are always real dfr^[0] := sf1r^[0] + sf2r^[0]; dfi^[0] := sf1r^[0] - sf2r^[0]; // dfr [nbr_coef] = dfr^[h_nbr_coef] := sf1r^[h_nbr_coef]; dfi^[h_nbr_coef] := sf2r^[h_nbr_coef]; // Others are conjugate complex numbers sf1i := pointer(longint(sf1r) + (h_nbr_coef * sizeof_flt)); sf2i := pointer(longint(sf1i) + (nbr_coef * sizeof_flt)); for i := 1 to h_nbr_coef-1 do begin c := cos_ptr^[i]; // cos (i*PI/nbr_coef); s := cos_ptr^[h_nbr_coef - i]; // sin (i*PI/nbr_coef); v := sf2r^[i] * c - sf2i^[i] * s; dfr^[i] := sf1r^[i] + v; dfi^[-i] := sf1r^[i] - v; // dfr [nbr_coef - i] = v := sf2r^[i] * s + sf2i^[i] * c; dfi^[i] := v + sf1i^[i]; dfi^[nbr_coef - i] := v - sf1i^[i]; end; inc(coef_index, d_nbr_coef); until (coef_index >= _length); // Prepare to the next pass temp_ptr := df; df := sf; sf := temp_ptr; end; end (*______________________________________________ * * Special cases *______________________________________________ *) // 4-point FFT else if _nbr_bits = 2 then begin f^[n1] := x^[0] - x^[n2]; f^[n3] := x^[n1] - x^[n3]; b_0 := x^[0] + x^[n2]; b_2 := x^[n1] + x^[n3]; f^[0] := b_0 + b_2; f^[n2] := b_0 - b_2; end // 2-point FFT else if _nbr_bits = 1 then begin f^[0] := x^[0] + x^[n1]; f^[n1] := x^[0] - x^[n1]; end // 1-point FFT else f^[0] := x^[0]; end; (*==========================================================================*/ /* Name: do_ifft */ /* Description: Compute the inverse FFT of the array. Notice that */ /* IFFT (FFT (x)) = x * length (x). Data must be */ /* post-scaled. */ /* Input parameters: */ /* - f: pointer on the source array (frequencies). */ /* f [0...length(x)/2] = real values, */ /* f [length(x)/2+1...length(x)-1] = imaginary values of */ /* coefficents 1...length(x)/2-1. */ /* Output parameters: */ /* - x: pointer on the destination array (time). */ /*==========================================================================*) procedure TFFTReal.do_ifft(const f: pflt_array; x: pflt_array); var n1, n2, n3 : integer; n4, n5, n6, n7 : integer; sf, df, df_temp : pflt_array; pass : integer; nbr_coef : longint; h_nbr_coef : longint; d_nbr_coef : longint; coef_index : longint; cos_ptr : pflt_array; i : longint; sfr, sfi : pflt_array; df1r, df2r : pflt_array; df1i, df2i : pflt_array; c, s, vr, vi : flt_t; temp_ptr : pflt_array; sqrt2_2 : flt_t; bit_rev_lut_ptr : plongarray; sf2 : pflt_array; b_0, b_1, b_2, b_3 : flt_t; begin n1 := 1; n2 := 2; n3 := 3; n4 := 4; n5 := 5; n6 := 6; n7 := 7; (*______________________________________________ * * General case *______________________________________________ *) if _nbr_bits > 2 then begin sf := f; if _nbr_bits and 1 <> 0 then begin df := pointer(_buffer_ptr); df_temp := x; end else begin df := x; df_temp := pointer(_buffer_ptr); end; // Do the transformation in several pass // First pass for pass := _nbr_bits-1 downto 3 do begin coef_index := 0; nbr_coef := 1 shl pass; h_nbr_coef := nbr_coef shr 1; d_nbr_coef := nbr_coef shl 1; cos_ptr := pointer(_trigo_lut.get_ptr(pass)); repeat sfr := pointer(longint(sf) + (coef_index*sizeof_flt)); sfi := pointer(longint(sfr) + (nbr_coef*sizeof_flt)); df1r := pointer(longint(df) + (coef_index*sizeof_flt)); df2r := pointer(longint(df1r) + (nbr_coef*sizeof_flt)); // Extreme coefficients are always real df1r^[0] := sfr^[0] + sfi^[0]; // + sfr [nbr_coef] df2r^[0] := sfr^[0] - sfi^[0]; // - sfr [nbr_coef] df1r^[h_nbr_coef] := sfr^[h_nbr_coef] * 2; df2r^[h_nbr_coef] := sfi^[h_nbr_coef] * 2; // Others are conjugate complex numbers df1i := pointer(longint(df1r) + (h_nbr_coef*sizeof_flt)); df2i := pointer(longint(df1i) + (nbr_coef*sizeof_flt)); for i := 1 to h_nbr_coef-1 do begin df1r^[i] := sfr^[i] + sfi^[-i]; // + sfr [nbr_coef - i] df1i^[i] := sfi^[i] - sfi^[nbr_coef - i]; c := cos_ptr^[i]; // cos (i*PI/nbr_coef); s := cos_ptr^[h_nbr_coef - i]; // sin (i*PI/nbr_coef); vr := sfr^[i] - sfi^[-i]; // - sfr [nbr_coef - i] vi := sfi^[i] + sfi^[nbr_coef - i]; df2r^[i] := vr * c + vi * s; df2i^[i] := vi * c - vr * s; end; inc(coef_index, d_nbr_coef); until (coef_index >= _length); // Prepare to the next pass if (pass < _nbr_bits - 1) then begin temp_ptr := df; df := sf; sf := temp_ptr; end else begin sf := df; df := df_temp; end end; // Antepenultimate pass sqrt2_2 := _sqrt2_2; coef_index := 0; repeat df^[coef_index] := sf^[coef_index] + sf^[coef_index + 4]; df^[coef_index + 4] := sf^[coef_index] - sf^[coef_index + 4]; df^[coef_index + 2] := sf^[coef_index + 2] * 2; df^[coef_index + 6] := sf^[coef_index + 6] * 2; df^[coef_index + 1] := sf^[coef_index + 1] + sf^[coef_index + 3]; df^[coef_index + 3] := sf^[coef_index + 5] - sf^[coef_index + 7]; vr := sf^[coef_index + 1] - sf^[coef_index + 3]; vi := sf^[coef_index + 5] + sf^[coef_index + 7]; df^[coef_index + 5] := (vr + vi) * sqrt2_2; df^[coef_index + 7] := (vi - vr) * sqrt2_2; inc(coef_index, 8); until (coef_index >= _length); // Penultimate and last pass at once coef_index := 0; bit_rev_lut_ptr := pointer(_bit_rev_lut.get_ptr); sf2 := df; repeat b_0 := sf2^[0] + sf2^[n2]; b_2 := sf2^[0] - sf2^[n2]; b_1 := sf2^[n1] * 2; b_3 := sf2^[n3] * 2; x^[bit_rev_lut_ptr^[0]] := b_0 + b_1; x^[bit_rev_lut_ptr^[n1]] := b_0 - b_1; x^[bit_rev_lut_ptr^[n2]] := b_2 + b_3; x^[bit_rev_lut_ptr^[n3]] := b_2 - b_3; b_0 := sf2^[n4] + sf2^[n6]; b_2 := sf2^[n4] - sf2^[n6]; b_1 := sf2^[n5] * 2; b_3 := sf2^[n7] * 2; x^[bit_rev_lut_ptr^[n4]] := b_0 + b_1; x^[bit_rev_lut_ptr^[n5]] := b_0 - b_1; x^[bit_rev_lut_ptr^[n6]] := b_2 + b_3; x^[bit_rev_lut_ptr^[n7]] := b_2 - b_3; inc(sf2, 8); inc(coef_index, 8); inc(bit_rev_lut_ptr, 8); until (coef_index >= _length); end (*______________________________________________ * * Special cases *______________________________________________ *) // 4-point IFFT else if _nbr_bits = 2 then begin b_0 := f^[0] + f [n2]; b_2 := f^[0] - f [n2]; x^[0] := b_0 + f [n1] * 2; x^[n2] := b_0 - f [n1] * 2; x^[n1] := b_2 + f [n3] * 2; x^[n3] := b_2 - f [n3] * 2; end // 2-point IFFT else if _nbr_bits = 1 then begin x^[0] := f^[0] + f^[n1]; x^[n1] := f^[0] - f^[n1]; end // 1-point IFFT else x^[0] := f^[0]; end; (*==========================================================================*/ /* Name: rescale */ /* Description: Scale an array by divide each element by its length. */ /* This function should be called after FFT + IFFT. */ /* Input/Output parameters: */ /* - x: pointer on array to rescale (time or frequency). */ /*==========================================================================*) procedure TFFTReal.rescale(x: pflt_array); var mul : flt_t; i : longint; begin mul := 1.0 / _length; i := _length - 1; repeat x^[i] := x^[i] * mul; dec(i); until (i < 0); end; end.