|
1 (***************************************************************************** |
|
2 |
|
3 DIGITAL SIGNAL PROCESSING TOOLS |
|
4 Version 1.03, 2001/06/15 |
|
5 (c) 1999 - Laurent de Soras |
|
6 |
|
7 FFTReal.h |
|
8 Fourier transformation of real number arrays. |
|
9 Portable ISO C++ |
|
10 |
|
11 ------------------------------------------------------------------------------ |
|
12 |
|
13 LEGAL |
|
14 |
|
15 Source code may be freely used for any purpose, including commercial |
|
16 applications. Programs must display in their "About" dialog-box (or |
|
17 documentation) a text telling they use these routines by Laurent de Soras. |
|
18 Modified source code can be distributed, but modifications must be clearly |
|
19 indicated. |
|
20 |
|
21 CONTACT |
|
22 |
|
23 Laurent de Soras |
|
24 92 avenue Albert 1er |
|
25 92500 Rueil-Malmaison |
|
26 France |
|
27 |
|
28 ldesoras@club-internet.fr |
|
29 |
|
30 ------------------------------------------------------------------------------ |
|
31 |
|
32 Translation to ObjectPascal by : |
|
33 Frederic Vanmol |
|
34 frederic@axiworld.be |
|
35 |
|
36 *****************************************************************************) |
|
37 |
|
38 |
|
39 unit |
|
40 FFTReal; |
|
41 |
|
42 interface |
|
43 |
|
44 uses |
|
45 Windows; |
|
46 |
|
47 (* Change this typedef to use a different floating point type in your FFTs |
|
48 (i.e. float, double or long double). *) |
|
49 type |
|
50 pflt_t = ^flt_t; |
|
51 flt_t = single; |
|
52 |
|
53 pflt_array = ^flt_array; |
|
54 flt_array = array[0..0] of flt_t; |
|
55 |
|
56 plongarray = ^longarray; |
|
57 longarray = array[0..0] of longint; |
|
58 |
|
59 const |
|
60 sizeof_flt : longint = SizeOf(flt_t); |
|
61 |
|
62 |
|
63 |
|
64 type |
|
65 // Bit reversed look-up table nested class |
|
66 TBitReversedLUT = class |
|
67 private |
|
68 _ptr : plongint; |
|
69 public |
|
70 constructor Create(const nbr_bits: integer); |
|
71 destructor Destroy; override; |
|
72 function get_ptr: plongint; |
|
73 end; |
|
74 |
|
75 // Trigonometric look-up table nested class |
|
76 TTrigoLUT = class |
|
77 private |
|
78 _ptr : pflt_t; |
|
79 public |
|
80 constructor Create(const nbr_bits: integer); |
|
81 destructor Destroy; override; |
|
82 function get_ptr(const level: integer): pflt_t; |
|
83 end; |
|
84 |
|
85 TFFTReal = class |
|
86 private |
|
87 _bit_rev_lut : TBitReversedLUT; |
|
88 _trigo_lut : TTrigoLUT; |
|
89 _sqrt2_2 : flt_t; |
|
90 _length : longint; |
|
91 _nbr_bits : integer; |
|
92 _buffer_ptr : pflt_t; |
|
93 public |
|
94 constructor Create(const length: longint); |
|
95 destructor Destroy; override; |
|
96 |
|
97 procedure do_fft(f: pflt_array; const x: pflt_array); |
|
98 procedure do_ifft(const f: pflt_array; x: pflt_array); |
|
99 procedure rescale(x: pflt_array); |
|
100 end; |
|
101 |
|
102 |
|
103 |
|
104 |
|
105 |
|
106 |
|
107 |
|
108 implementation |
|
109 |
|
110 uses |
|
111 Math; |
|
112 |
|
113 { TBitReversedLUT } |
|
114 |
|
115 constructor TBitReversedLUT.Create(const nbr_bits: integer); |
|
116 var |
|
117 length : longint; |
|
118 cnt : longint; |
|
119 br_index : longint; |
|
120 bit : longint; |
|
121 begin |
|
122 inherited Create; |
|
123 |
|
124 length := 1 shl nbr_bits; |
|
125 GetMem(_ptr, length*SizeOf(longint)); |
|
126 |
|
127 br_index := 0; |
|
128 plongarray(_ptr)^[0] := 0; |
|
129 for cnt := 1 to length-1 do |
|
130 begin |
|
131 // ++br_index (bit reversed) |
|
132 bit := length shr 1; |
|
133 br_index := br_index xor bit; |
|
134 while br_index and bit = 0 do |
|
135 begin |
|
136 bit := bit shr 1; |
|
137 br_index := br_index xor bit; |
|
138 end; |
|
139 |
|
140 plongarray(_ptr)^[cnt] := br_index; |
|
141 end; |
|
142 end; |
|
143 |
|
144 destructor TBitReversedLUT.Destroy; |
|
145 begin |
|
146 FreeMem(_ptr); |
|
147 _ptr := nil; |
|
148 inherited; |
|
149 end; |
|
150 |
|
151 function TBitReversedLUT.get_ptr: plongint; |
|
152 begin |
|
153 Result := _ptr; |
|
154 end; |
|
155 |
|
156 { TTrigLUT } |
|
157 |
|
158 constructor TTrigoLUT.Create(const nbr_bits: integer); |
|
159 var |
|
160 total_len : longint; |
|
161 PI : double; |
|
162 level : integer; |
|
163 level_len : longint; |
|
164 level_ptr : pflt_array; |
|
165 mul : double; |
|
166 i : longint; |
|
167 begin |
|
168 inherited Create; |
|
169 |
|
170 _ptr := nil; |
|
171 |
|
172 if (nbr_bits > 3) then |
|
173 begin |
|
174 total_len := (1 shl (nbr_bits - 1)) - 4; |
|
175 GetMem(_ptr, total_len * sizeof_flt); |
|
176 |
|
177 PI := ArcTan(1) * 4; |
|
178 for level := 3 to nbr_bits-1 do |
|
179 begin |
|
180 level_len := 1 shl (level - 1); |
|
181 level_ptr := pointer(get_ptr(level)); |
|
182 mul := PI / (level_len shl 1); |
|
183 |
|
184 for i := 0 to level_len-1 do |
|
185 level_ptr^[i] := cos(i * mul); |
|
186 end; |
|
187 end; |
|
188 end; |
|
189 |
|
190 destructor TTrigoLUT.Destroy; |
|
191 begin |
|
192 FreeMem(_ptr); |
|
193 _ptr := nil; |
|
194 inherited; |
|
195 end; |
|
196 |
|
197 function TTrigoLUT.get_ptr(const level: integer): pflt_t; |
|
198 var |
|
199 tempp : pflt_t; |
|
200 begin |
|
201 tempp := _ptr; |
|
202 inc(tempp, (1 shl (level-1)) - 4); |
|
203 Result := tempp; |
|
204 end; |
|
205 |
|
206 { TFFTReal } |
|
207 |
|
208 constructor TFFTReal.Create(const length: longint); |
|
209 begin |
|
210 inherited Create; |
|
211 |
|
212 _length := length; |
|
213 _nbr_bits := Floor(Ln(length) / Ln(2) + 0.5); |
|
214 _bit_rev_lut := TBitReversedLUT.Create(Floor(Ln(length) / Ln(2) + 0.5)); |
|
215 _trigo_lut := TTrigoLUT.Create(Floor(Ln(length) / Ln(2) + 0.05)); |
|
216 _sqrt2_2 := Sqrt(2) * 0.5; |
|
217 |
|
218 _buffer_ptr := nil; |
|
219 if _nbr_bits > 2 then |
|
220 GetMem(_buffer_ptr, _length * sizeof_flt); |
|
221 end; |
|
222 |
|
223 destructor TFFTReal.Destroy; |
|
224 begin |
|
225 if _buffer_ptr <> nil then |
|
226 begin |
|
227 FreeMem(_buffer_ptr); |
|
228 _buffer_ptr := nil; |
|
229 end; |
|
230 |
|
231 _bit_rev_lut.Free; |
|
232 _bit_rev_lut := nil; |
|
233 _trigo_lut.Free; |
|
234 _trigo_lut := nil; |
|
235 |
|
236 inherited; |
|
237 end; |
|
238 |
|
239 (*==========================================================================*/ |
|
240 /* Name: do_fft */ |
|
241 /* Description: Compute the FFT of the array. */ |
|
242 /* Input parameters: */ |
|
243 /* - x: pointer on the source array (time). */ |
|
244 /* Output parameters: */ |
|
245 /* - f: pointer on the destination array (frequencies). */ |
|
246 /* f [0...length(x)/2] = real values, */ |
|
247 /* f [length(x)/2+1...length(x)-1] = imaginary values of */ |
|
248 /* coefficents 1...length(x)/2-1. */ |
|
249 /*==========================================================================*) |
|
250 procedure TFFTReal.do_fft(f: pflt_array; const x: pflt_array); |
|
251 var |
|
252 sf, df : pflt_array; |
|
253 pass : integer; |
|
254 nbr_coef : longint; |
|
255 h_nbr_coef : longint; |
|
256 d_nbr_coef : longint; |
|
257 coef_index : longint; |
|
258 bit_rev_lut_ptr : plongarray; |
|
259 rev_index_0 : longint; |
|
260 rev_index_1 : longint; |
|
261 rev_index_2 : longint; |
|
262 rev_index_3 : longint; |
|
263 df2 : pflt_array; |
|
264 n1, n2, n3 : integer; |
|
265 sf_0, sf_2 : flt_t; |
|
266 sqrt2_2 : flt_t; |
|
267 v : flt_t; |
|
268 cos_ptr : pflt_array; |
|
269 i : longint; |
|
270 sf1r, sf2r : pflt_array; |
|
271 dfr, dfi : pflt_array; |
|
272 sf1i, sf2i : pflt_array; |
|
273 c, s : flt_t; |
|
274 temp_ptr : pflt_array; |
|
275 b_0, b_2 : flt_t; |
|
276 begin |
|
277 n1 := 1; |
|
278 n2 := 2; |
|
279 n3 := 3; |
|
280 |
|
281 (*______________________________________________ |
|
282 * |
|
283 * General case |
|
284 *______________________________________________ |
|
285 *) |
|
286 |
|
287 if _nbr_bits > 2 then |
|
288 begin |
|
289 if _nbr_bits and 1 <> 0 then |
|
290 begin |
|
291 df := pointer(_buffer_ptr); |
|
292 sf := f; |
|
293 end |
|
294 else |
|
295 begin |
|
296 df := f; |
|
297 sf := pointer(_buffer_ptr); |
|
298 end; |
|
299 |
|
300 // |
|
301 // Do the transformation in several passes |
|
302 // |
|
303 |
|
304 // First and second pass at once |
|
305 bit_rev_lut_ptr := pointer(_bit_rev_lut.get_ptr); |
|
306 coef_index := 0; |
|
307 |
|
308 repeat |
|
309 rev_index_0 := bit_rev_lut_ptr^[coef_index]; |
|
310 rev_index_1 := bit_rev_lut_ptr^[coef_index + 1]; |
|
311 rev_index_2 := bit_rev_lut_ptr^[coef_index + 2]; |
|
312 rev_index_3 := bit_rev_lut_ptr^[coef_index + 3]; |
|
313 |
|
314 df2 := pointer(longint(df) + (coef_index*sizeof_flt)); |
|
315 df2^[n1] := x^[rev_index_0] - x^[rev_index_1]; |
|
316 df2^[n3] := x^[rev_index_2] - x^[rev_index_3]; |
|
317 |
|
318 sf_0 := x^[rev_index_0] + x^[rev_index_1]; |
|
319 sf_2 := x^[rev_index_2] + x^[rev_index_3]; |
|
320 |
|
321 df2^[0] := sf_0 + sf_2; |
|
322 df2^[n2] := sf_0 - sf_2; |
|
323 |
|
324 inc(coef_index, 4); |
|
325 until (coef_index >= _length); |
|
326 |
|
327 |
|
328 // Third pass |
|
329 coef_index := 0; |
|
330 sqrt2_2 := _sqrt2_2; |
|
331 |
|
332 repeat |
|
333 sf^[coef_index] := df^[coef_index] + df^[coef_index + 4]; |
|
334 sf^[coef_index + 4] := df^[coef_index] - df^[coef_index + 4]; |
|
335 sf^[coef_index + 2] := df^[coef_index + 2]; |
|
336 sf^[coef_index + 6] := df^[coef_index + 6]; |
|
337 |
|
338 v := (df [coef_index + 5] - df^[coef_index + 7]) * sqrt2_2; |
|
339 sf^[coef_index + 1] := df^[coef_index + 1] + v; |
|
340 sf^[coef_index + 3] := df^[coef_index + 1] - v; |
|
341 |
|
342 v := (df^[coef_index + 5] + df^[coef_index + 7]) * sqrt2_2; |
|
343 sf [coef_index + 5] := v + df^[coef_index + 3]; |
|
344 sf [coef_index + 7] := v - df^[coef_index + 3]; |
|
345 |
|
346 inc(coef_index, 8); |
|
347 until (coef_index >= _length); |
|
348 |
|
349 |
|
350 // Next pass |
|
351 for pass := 3 to _nbr_bits-1 do |
|
352 begin |
|
353 coef_index := 0; |
|
354 nbr_coef := 1 shl pass; |
|
355 h_nbr_coef := nbr_coef shr 1; |
|
356 d_nbr_coef := nbr_coef shl 1; |
|
357 |
|
358 cos_ptr := pointer(_trigo_lut.get_ptr(pass)); |
|
359 repeat |
|
360 sf1r := pointer(longint(sf) + (coef_index * sizeof_flt)); |
|
361 sf2r := pointer(longint(sf1r) + (nbr_coef * sizeof_flt)); |
|
362 dfr := pointer(longint(df) + (coef_index * sizeof_flt)); |
|
363 dfi := pointer(longint(dfr) + (nbr_coef * sizeof_flt)); |
|
364 |
|
365 // Extreme coefficients are always real |
|
366 dfr^[0] := sf1r^[0] + sf2r^[0]; |
|
367 dfi^[0] := sf1r^[0] - sf2r^[0]; // dfr [nbr_coef] = |
|
368 dfr^[h_nbr_coef] := sf1r^[h_nbr_coef]; |
|
369 dfi^[h_nbr_coef] := sf2r^[h_nbr_coef]; |
|
370 |
|
371 // Others are conjugate complex numbers |
|
372 sf1i := pointer(longint(sf1r) + (h_nbr_coef * sizeof_flt)); |
|
373 sf2i := pointer(longint(sf1i) + (nbr_coef * sizeof_flt)); |
|
374 |
|
375 for i := 1 to h_nbr_coef-1 do |
|
376 begin |
|
377 c := cos_ptr^[i]; // cos (i*PI/nbr_coef); |
|
378 s := cos_ptr^[h_nbr_coef - i]; // sin (i*PI/nbr_coef); |
|
379 |
|
380 v := sf2r^[i] * c - sf2i^[i] * s; |
|
381 dfr^[i] := sf1r^[i] + v; |
|
382 dfi^[-i] := sf1r^[i] - v; // dfr [nbr_coef - i] = |
|
383 |
|
384 v := sf2r^[i] * s + sf2i^[i] * c; |
|
385 dfi^[i] := v + sf1i^[i]; |
|
386 dfi^[nbr_coef - i] := v - sf1i^[i]; |
|
387 end; |
|
388 |
|
389 inc(coef_index, d_nbr_coef); |
|
390 until (coef_index >= _length); |
|
391 |
|
392 // Prepare to the next pass |
|
393 temp_ptr := df; |
|
394 df := sf; |
|
395 sf := temp_ptr; |
|
396 end; |
|
397 end |
|
398 |
|
399 (*______________________________________________ |
|
400 * |
|
401 * Special cases |
|
402 *______________________________________________ |
|
403 *) |
|
404 |
|
405 // 4-point FFT |
|
406 else if _nbr_bits = 2 then |
|
407 begin |
|
408 f^[n1] := x^[0] - x^[n2]; |
|
409 f^[n3] := x^[n1] - x^[n3]; |
|
410 |
|
411 b_0 := x^[0] + x^[n2]; |
|
412 b_2 := x^[n1] + x^[n3]; |
|
413 |
|
414 f^[0] := b_0 + b_2; |
|
415 f^[n2] := b_0 - b_2; |
|
416 end |
|
417 |
|
418 // 2-point FFT |
|
419 else if _nbr_bits = 1 then |
|
420 begin |
|
421 f^[0] := x^[0] + x^[n1]; |
|
422 f^[n1] := x^[0] - x^[n1]; |
|
423 end |
|
424 |
|
425 // 1-point FFT |
|
426 else |
|
427 f^[0] := x^[0]; |
|
428 end; |
|
429 |
|
430 |
|
431 (*==========================================================================*/ |
|
432 /* Name: do_ifft */ |
|
433 /* Description: Compute the inverse FFT of the array. Notice that */ |
|
434 /* IFFT (FFT (x)) = x * length (x). Data must be */ |
|
435 /* post-scaled. */ |
|
436 /* Input parameters: */ |
|
437 /* - f: pointer on the source array (frequencies). */ |
|
438 /* f [0...length(x)/2] = real values, */ |
|
439 /* f [length(x)/2+1...length(x)-1] = imaginary values of */ |
|
440 /* coefficents 1...length(x)/2-1. */ |
|
441 /* Output parameters: */ |
|
442 /* - x: pointer on the destination array (time). */ |
|
443 /*==========================================================================*) |
|
444 procedure TFFTReal.do_ifft(const f: pflt_array; x: pflt_array); |
|
445 var |
|
446 n1, n2, n3 : integer; |
|
447 n4, n5, n6, n7 : integer; |
|
448 sf, df, df_temp : pflt_array; |
|
449 pass : integer; |
|
450 nbr_coef : longint; |
|
451 h_nbr_coef : longint; |
|
452 d_nbr_coef : longint; |
|
453 coef_index : longint; |
|
454 cos_ptr : pflt_array; |
|
455 i : longint; |
|
456 sfr, sfi : pflt_array; |
|
457 df1r, df2r : pflt_array; |
|
458 df1i, df2i : pflt_array; |
|
459 c, s, vr, vi : flt_t; |
|
460 temp_ptr : pflt_array; |
|
461 sqrt2_2 : flt_t; |
|
462 bit_rev_lut_ptr : plongarray; |
|
463 sf2 : pflt_array; |
|
464 b_0, b_1, b_2, b_3 : flt_t; |
|
465 begin |
|
466 n1 := 1; |
|
467 n2 := 2; |
|
468 n3 := 3; |
|
469 n4 := 4; |
|
470 n5 := 5; |
|
471 n6 := 6; |
|
472 n7 := 7; |
|
473 |
|
474 (*______________________________________________ |
|
475 * |
|
476 * General case |
|
477 *______________________________________________ |
|
478 *) |
|
479 |
|
480 if _nbr_bits > 2 then |
|
481 begin |
|
482 sf := f; |
|
483 |
|
484 if _nbr_bits and 1 <> 0 then |
|
485 begin |
|
486 df := pointer(_buffer_ptr); |
|
487 df_temp := x; |
|
488 end |
|
489 else |
|
490 begin |
|
491 df := x; |
|
492 df_temp := pointer(_buffer_ptr); |
|
493 end; |
|
494 |
|
495 // Do the transformation in several pass |
|
496 |
|
497 // First pass |
|
498 for pass := _nbr_bits-1 downto 3 do |
|
499 begin |
|
500 coef_index := 0; |
|
501 nbr_coef := 1 shl pass; |
|
502 h_nbr_coef := nbr_coef shr 1; |
|
503 d_nbr_coef := nbr_coef shl 1; |
|
504 |
|
505 cos_ptr := pointer(_trigo_lut.get_ptr(pass)); |
|
506 |
|
507 repeat |
|
508 sfr := pointer(longint(sf) + (coef_index*sizeof_flt)); |
|
509 sfi := pointer(longint(sfr) + (nbr_coef*sizeof_flt)); |
|
510 df1r := pointer(longint(df) + (coef_index*sizeof_flt)); |
|
511 df2r := pointer(longint(df1r) + (nbr_coef*sizeof_flt)); |
|
512 |
|
513 // Extreme coefficients are always real |
|
514 df1r^[0] := sfr^[0] + sfi^[0]; // + sfr [nbr_coef] |
|
515 df2r^[0] := sfr^[0] - sfi^[0]; // - sfr [nbr_coef] |
|
516 df1r^[h_nbr_coef] := sfr^[h_nbr_coef] * 2; |
|
517 df2r^[h_nbr_coef] := sfi^[h_nbr_coef] * 2; |
|
518 |
|
519 // Others are conjugate complex numbers |
|
520 df1i := pointer(longint(df1r) + (h_nbr_coef*sizeof_flt)); |
|
521 df2i := pointer(longint(df1i) + (nbr_coef*sizeof_flt)); |
|
522 |
|
523 for i := 1 to h_nbr_coef-1 do |
|
524 begin |
|
525 df1r^[i] := sfr^[i] + sfi^[-i]; // + sfr [nbr_coef - i] |
|
526 df1i^[i] := sfi^[i] - sfi^[nbr_coef - i]; |
|
527 |
|
528 c := cos_ptr^[i]; // cos (i*PI/nbr_coef); |
|
529 s := cos_ptr^[h_nbr_coef - i]; // sin (i*PI/nbr_coef); |
|
530 vr := sfr^[i] - sfi^[-i]; // - sfr [nbr_coef - i] |
|
531 vi := sfi^[i] + sfi^[nbr_coef - i]; |
|
532 |
|
533 df2r^[i] := vr * c + vi * s; |
|
534 df2i^[i] := vi * c - vr * s; |
|
535 end; |
|
536 |
|
537 inc(coef_index, d_nbr_coef); |
|
538 until (coef_index >= _length); |
|
539 |
|
540 |
|
541 // Prepare to the next pass |
|
542 if (pass < _nbr_bits - 1) then |
|
543 begin |
|
544 temp_ptr := df; |
|
545 df := sf; |
|
546 sf := temp_ptr; |
|
547 end |
|
548 else |
|
549 begin |
|
550 sf := df; |
|
551 df := df_temp; |
|
552 end |
|
553 end; |
|
554 |
|
555 // Antepenultimate pass |
|
556 sqrt2_2 := _sqrt2_2; |
|
557 coef_index := 0; |
|
558 |
|
559 repeat |
|
560 df^[coef_index] := sf^[coef_index] + sf^[coef_index + 4]; |
|
561 df^[coef_index + 4] := sf^[coef_index] - sf^[coef_index + 4]; |
|
562 df^[coef_index + 2] := sf^[coef_index + 2] * 2; |
|
563 df^[coef_index + 6] := sf^[coef_index + 6] * 2; |
|
564 |
|
565 df^[coef_index + 1] := sf^[coef_index + 1] + sf^[coef_index + 3]; |
|
566 df^[coef_index + 3] := sf^[coef_index + 5] - sf^[coef_index + 7]; |
|
567 |
|
568 vr := sf^[coef_index + 1] - sf^[coef_index + 3]; |
|
569 vi := sf^[coef_index + 5] + sf^[coef_index + 7]; |
|
570 |
|
571 df^[coef_index + 5] := (vr + vi) * sqrt2_2; |
|
572 df^[coef_index + 7] := (vi - vr) * sqrt2_2; |
|
573 |
|
574 inc(coef_index, 8); |
|
575 until (coef_index >= _length); |
|
576 |
|
577 |
|
578 // Penultimate and last pass at once |
|
579 coef_index := 0; |
|
580 bit_rev_lut_ptr := pointer(_bit_rev_lut.get_ptr); |
|
581 sf2 := df; |
|
582 |
|
583 repeat |
|
584 b_0 := sf2^[0] + sf2^[n2]; |
|
585 b_2 := sf2^[0] - sf2^[n2]; |
|
586 b_1 := sf2^[n1] * 2; |
|
587 b_3 := sf2^[n3] * 2; |
|
588 |
|
589 x^[bit_rev_lut_ptr^[0]] := b_0 + b_1; |
|
590 x^[bit_rev_lut_ptr^[n1]] := b_0 - b_1; |
|
591 x^[bit_rev_lut_ptr^[n2]] := b_2 + b_3; |
|
592 x^[bit_rev_lut_ptr^[n3]] := b_2 - b_3; |
|
593 |
|
594 b_0 := sf2^[n4] + sf2^[n6]; |
|
595 b_2 := sf2^[n4] - sf2^[n6]; |
|
596 b_1 := sf2^[n5] * 2; |
|
597 b_3 := sf2^[n7] * 2; |
|
598 |
|
599 x^[bit_rev_lut_ptr^[n4]] := b_0 + b_1; |
|
600 x^[bit_rev_lut_ptr^[n5]] := b_0 - b_1; |
|
601 x^[bit_rev_lut_ptr^[n6]] := b_2 + b_3; |
|
602 x^[bit_rev_lut_ptr^[n7]] := b_2 - b_3; |
|
603 |
|
604 inc(sf2, 8); |
|
605 inc(coef_index, 8); |
|
606 inc(bit_rev_lut_ptr, 8); |
|
607 until (coef_index >= _length); |
|
608 end |
|
609 |
|
610 (*______________________________________________ |
|
611 * |
|
612 * Special cases |
|
613 *______________________________________________ |
|
614 *) |
|
615 |
|
616 // 4-point IFFT |
|
617 else if _nbr_bits = 2 then |
|
618 begin |
|
619 b_0 := f^[0] + f [n2]; |
|
620 b_2 := f^[0] - f [n2]; |
|
621 |
|
622 x^[0] := b_0 + f [n1] * 2; |
|
623 x^[n2] := b_0 - f [n1] * 2; |
|
624 x^[n1] := b_2 + f [n3] * 2; |
|
625 x^[n3] := b_2 - f [n3] * 2; |
|
626 end |
|
627 |
|
628 // 2-point IFFT |
|
629 else if _nbr_bits = 1 then |
|
630 begin |
|
631 x^[0] := f^[0] + f^[n1]; |
|
632 x^[n1] := f^[0] - f^[n1]; |
|
633 end |
|
634 |
|
635 // 1-point IFFT |
|
636 else |
|
637 x^[0] := f^[0]; |
|
638 end; |
|
639 |
|
640 (*==========================================================================*/ |
|
641 /* Name: rescale */ |
|
642 /* Description: Scale an array by divide each element by its length. */ |
|
643 /* This function should be called after FFT + IFFT. */ |
|
644 /* Input/Output parameters: */ |
|
645 /* - x: pointer on array to rescale (time or frequency). */ |
|
646 /*==========================================================================*) |
|
647 procedure TFFTReal.rescale(x: pflt_array); |
|
648 var |
|
649 mul : flt_t; |
|
650 i : longint; |
|
651 begin |
|
652 mul := 1.0 / _length; |
|
653 i := _length - 1; |
|
654 |
|
655 repeat |
|
656 x^[i] := x^[i] * mul; |
|
657 dec(i); |
|
658 until (i < 0); |
|
659 end; |
|
660 |
|
661 end. |