diff options
Diffstat (limited to 'demos/spectrum/fftreal/fftreal.pas')
| -rw-r--r-- | demos/spectrum/fftreal/fftreal.pas | 661 |
1 files changed, 661 insertions, 0 deletions
diff --git a/demos/spectrum/fftreal/fftreal.pas b/demos/spectrum/fftreal/fftreal.pas new file mode 100644 index 0000000..ea63754 --- /dev/null +++ b/demos/spectrum/fftreal/fftreal.pas @@ -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. |