| (***************************************************************************** |
| |
| 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. |