MathgeomGLS

UFourier.pas
223 строки · 8.6 Кб
Перенос по словам
1
unit UFourier;
2

3
(**************************************************************************
4
 ************************* UFourier ***************************************
5
 * Purpose: Implements Fast Fourier Transform and InverseFastFourierTrans *
6
 *                                                                        *
7
 *                                                                        *
8
 * procedures                                                             *
9
 *   FFT(_NumSamples           : Integer;                                 *
10
 *       _InValues, _OutValues : TComplex1DArray);                        *
11
 *                                                                        *
12
 *   IFFT(_NumSamples            : Integer;                               *
13
 *        _InValues, _OutValues : TComplex1DArray);                       *
14
 *                                                                        *
15
 * Note:                                                                  *
16
 *  Based on the DMath Package by Jean Debord                             *
17
 *  (http://www.unilim.fr/pages_perso/jean.debord/index.htm)              *
18
 *                                                                        *
19
 *  Adapted to use TComplex and TComplex1DArrays - 06.09.2004             *
20
 * ---------------------------------------------------------------------- *
21
 * (c) W.Blecher 2004 - ITMC RWTH Aachen                                  *
22
 **************************************************************************
23
 **************************************************************************)
24

25
interface
26

27
uses  System.SysUtils,
28
      UGlobalTools;
29

30
procedure FFT(_NumSamples    : Integer;
31
              _InValues      : TComplex1DArray;
32
              var _OutValues : TComplex1DArray);
33

34
procedure IFFT(_NumSamples            : Integer;
35
               _InValues      : TComplex1DArray;
36
               var _OutValues : TComplex1DArray);
37

38
implementation
39

40
const MaxPowerOfTwo : integer = 25;
41

42
(***************************************************************
43
 ************************** IsPowerOfTwo ***********************
44
 * Purpose: Checks if an Integer is a Power of Two             *
45
 *                                                             *
46
 * Input:   _X : Integer                                       *
47
 *                                                             *
48
 *        Number to check                                      *
49
 * Output: True/False                                          *
50
 * ----------------------------------------------------------- *
51
 * http://www.unilim.fr/pages_perso/jean.debord/index.htm      *
52
 ***************************************************************
53
 ***************************************************************)
54
function IsPowerOfTwo(_X : Integer) : Boolean;
55
var
56
  I, Y : Integer;
57
begin
58
  Y := 2;
59
  for I := 1 to Pred(MaxPowerOfTwo) do
60
    begin
61
      if _X = Y then
62
        begin
63
          IsPowerOfTwo := True;
64
          Exit;
65
        end;
66
      Y := Y shl 1;
67
    end;
68
  IsPowerOfTwo := False;
69
end;
70

71
(***************************************************************
72
 ************************** NumberOfBitsNeeded *****************
73
 * Purpose: How much Bits are needed to store the max value    *
74
 *                                                             *
75
 * Input:   _PowerOfTwo : Integer                              *
76
 *        Max Array Index                                      *
77

78
 * Output: Integer                                             *
79
 * ----------------------------------------------------------- *
80
 * http://www.unilim.fr/pages_perso/jean.debord/index.htm      *
81
 ***************************************************************
82
 ***************************************************************)
83
function NumberOfBitsNeeded(_PowerOfTwo : Integer) : Integer;
84
var
85
  I : Integer;
86
begin
87
  for I := 0 to MaxPowerOfTwo do
88
    begin
89
      if (_PowerOfTwo and (1 shl I)) <> 0 then
90
        begin
91
          NumberOfBitsNeeded := I;
92
          Exit;
93
        end;
94
    end;
95
  NumberOfBitsNeeded := 0;
96
end;
97

98
(***************************************************************
99
 ************************ ReverseBits **************************
100
 * ----------------------------------------------------------- *
101
 * http://www.unilim.fr/pages_perso/jean.debord/index.htm      *
102
 ***************************************************************
103
 ***************************************************************)
104
function ReverseBits(_Index, _NumBits : Integer) : Integer;
105
var
106
  I, Rev : Integer;
107
begin
108
  Rev := 0;
109
  for I := 0 to _NumBits - 1 do
110
    begin
111
      Rev := (Rev shl 1) or (_Index and 1);
112
      _Index := _Index shr 1;
113
    end;
114
  ReverseBits := Rev;
115
end;
116

117
(***************************************************************
118
 ************************* FourierTransform ********************
119
 * Purpose: Calculate the Fourier Transform of _ValsIn         *
120
 * ----------------------------------------------------------- *
121
 * http://www.unilim.fr/pages_perso/jean.debord/index.htm      *
122
 ***************************************************************
123
 ***************************************************************)
124
procedure FourierTransform(_AngleNumerator : Single;
125
                           _NumSamples     : Integer;
126
                           _ValsIn         : TComplex1DArray;
127
                           var _ValsOut    : TComplex1DArray);
128
var
129
  NumBits, I, J, K, N, BlockSize, BlockEnd : Integer;
130
  Delta_angle, Delta_ar                    : Single;
131
  Alpha, Beta                              : Single;
132
  Tr, Ti, Ar, Ai                           : Single;
133
begin
134
  if not IsPowerOfTwo(_NumSamples) or (_NumSamples < 2) then
135
    begin
136
      raise Exception.Create('Error in procedure Fourier: NumSamples=' + IntToStr(_NumSamples)+ ' is not a positive integer power of 2');
137
    end;
138

139
  NumBits := NumberOfBitsNeeded(_NumSamples);
140
  for I := 0 to _NumSamples - 1 do
141
    begin
142
      J := ReverseBits(I, NumBits);
143
      _ValsOut[J].Re := _ValsIn[I].Re;
144
      _ValsOut[J].img := _ValsIn[I].Img;
145
    end;
146

147
  BlockEnd := 1;
148
  BlockSize := 2;
149
  while BlockSize <= _NumSamples do
150
    begin
151
      Delta_angle := _AngleNumerator / BlockSize;
152
      Alpha := Sin(0.5 * Delta_angle);
153
      Alpha := 2.0 * Alpha * Alpha;
154
      Beta := Sin(Delta_angle);
155

156
      I := 0;
157
      while I < _NumSamples do
158
        begin
159
          Ar := 1.0;    (* cos(0) *)
160
          Ai := 0.0;    (* sin(0) *)
161

162
          J := I;
163
          for N := 0 to BlockEnd - 1 do
164
            begin
165
              K := J + BlockEnd;
166
              Tr := Ar * _ValsOut[K].Re - Ai * _ValsOut[K].Img;
167
              Ti := Ar * _ValsOut[K].Img + Ai * _ValsOut[K].Re;
168
              _ValsOut[K].Re := _ValsOut[J].Re - Tr;
169
              _ValsOut[K].Img := _ValsOut[J].Img - Ti;
170
              _ValsOut[J].Re := _ValsOut[J].Re + Tr;
171
              _ValsOut[J].Img := _ValsOut[J].Img + Ti;
172
              Delta_ar := Alpha * Ar + Beta * Ai;
173
              Ai := Ai - (Alpha * Ai - Beta * Ar);
174
              Ar := Ar - Delta_ar;
175
              Inc(J);
176
            end;
177

178
          I := I + BlockSize;
179
        end;
180

181
      BlockEnd := BlockSize;
182
      BlockSize := BlockSize shl 1;
183
    end;
184
end;
185

186
(***************************************************************
187
 ************************* FFT *********************************
188
 * Purpose: Caller for calculation of Fast Fourier Transform   *
189
 * ----------------------------------------------------------- *
190
 * http://www.unilim.fr/pages_perso/jean.debord/index.htm      *
191
 ***************************************************************
192
 ***************************************************************)
193
procedure FFT(_NumSamples    : Integer;
194
              _InValues      : TComplex1DArray;
195
              var _OutValues : TComplex1DArray);
196
begin
197
  FourierTransform(2 * PI, _NumSamples, _InValues, _OutValues);
198
end;
199

200
(***************************************************************
201
 ************************ IFFT *********************************
202
 * Purpose: Caller of Iverse Fast Fourier Transform            *
203
 * ----------------------------------------------------------- *
204
 * http://www.unilim.fr/pages_perso/jean.debord/index.htm      *
205
 ***************************************************************
206
 ***************************************************************)
207
procedure IFFT(_NumSamples    : Integer;
208
               _InValues      : TComplex1DArray;
209
               var _OutValues : TComplex1DArray);
210
var
211
  I : Integer;
212
begin
213
  FourierTransform(- 2 * PI, _NumSamples, _InValues, _OutValues);
214

215
  { Normalize the resulting time samples }
216
  for I := 0 to _NumSamples - 1 do
217
    begin
218
      _OutValues[I].Re := _OutValues[I].Re / _NumSamples;
219
      _OutValues[I].Img := _OutValues[I].Img / _NumSamples;
220
    end;
221
end;
222

223
end.
224
MathgeomGLS

Использование cookies
