diff options
Diffstat (limited to 'tests/manual/spectrum/3rdparty/fftreal/fftreal.pas')
| -rw-r--r-- | tests/manual/spectrum/3rdparty/fftreal/fftreal.pas | 661 |
1 files changed, 0 insertions, 661 deletions
diff --git a/tests/manual/spectrum/3rdparty/fftreal/fftreal.pas b/tests/manual/spectrum/3rdparty/fftreal/fftreal.pas deleted file mode 100644 index ea637545..00000000 --- a/tests/manual/spectrum/3rdparty/fftreal/fftreal.pas +++ /dev/null @@ -1,661 +0,0 @@ -(***************************************************************************** - - 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. |