--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/demos/spectrum/3rdparty/fftreal/fftreal.pas Fri Jun 11 14:24:45 2010 +0300
@@ -0,0 +1,661 @@
+(*****************************************************************************
+
+ 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.