1 /*
2 Elmer, A Finite Element Software for Multiphysical Problems
3
4 Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
5
6 This library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
10
11 This library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with this library (in file ../../LGPL-2.1); if not, write
18 to the Free Software Foundation, Inc., 51 Franklin Street,
19 Fifth Floor, Boston, MA 02110-1301 USA
20 */
21
22 /*
23 * (Power of two) Fast Fourier Transform Subroutine Library
24 *
25 * Power of two algorithm to compute the transform:
26 *
27 * the basic idea is to compute partial f(k) sums for which exp(-i*2*pi*n*k/N)
28 * is equal and which are also equal to different output frequencies. for
29 * example for every (k) and (k+N/2) exp(-i*2*pi*n*k/N) is the same if (n) is
30 * even. the algebra is as follows:
31 *
32 * N even,
33 * F(n) = (1/N)*sum(f(k)*exp(-i*2*pi*n*k/N)), k=0..N-1,n=0..N-1
34 * take
35 * W = exp(-i*2*pi/N)
36 * that is
37 * F(n) = (1/N)*sum(f(k)*W^(n*k)), k=0..N-1, n=0..N-1
38 * = (1/N)*sum((f(k)*W^(k*n)+f(k+N/2)*W^((k+N/2)*n)), k=0..N/2-1,n=0..N-1
39 * = (1/N)*sum((f(k)+f(k+N/2)*W^(n*N/2))*W^(k*n)), k=0..N/2-1,n=0..N-1
40 * now
41 * W^(n*N/2) = exp(-i*2*pi*n*N/(2*N))
42 * = exp(-i*pi*n) = { 1: n even, -1: n odd }, n=0..N-1
43 * take
44 * h(k) = f(k)+f(l+N/2), k=0..N/2-1
45 * g(k) = (f(k)-f(l+N/2)))*W^(k), k=0..N/2-1
46 * then
47 * 2*H(l) = (1/(N/2))*sum(h(k)*exp(-i*2*pi*l*k/(N/2)), l=0..N/2-1, (n=2*l)
48 * 2*G(l) = (1/(N/2))*sum(g(k)*exp(-i*2*pi*l*k/(N/2)), l=0..N/2-1, (n=2*l+1)
49 * and
50 * F(2*l) = H(l), l=0..N/2-1
51 * F(2*l+1) = G(l), l=0..N/2-1
52 *
53 * that is you can compute one size (N) discrete fourier transform as two
54 * (M=N/2) size transforms. further, if (M) is even after each successive
55 * division (i.e. (N) is power of two) you can continue the division until
56 * you feel like writing the sum out explicitly.
57 *
58 * all the routines below at the end use cfftf (1D complex forward transform)
59 * to actually compute the transform and do not in other ways depend on the
60 * algorithm (which implies that input sequence lengths must be power of two,
61 * serious limitation with storage for multidimensional transforms).so you can
62 * substitute it for different algorithm if you like. for rfftf, rfftb, gfftf,
63 * and gfftb (the real 1D transforms) the sequence lengths must be even.
64 *
65 * BTW. zero extend your input to next largest power of two when using the
66 * library. it does not check for it....
67 *
68 * The following routines are available:
69 *
70 *------------------------------------------------------------------------------
71 *
72 * cfftf: forward complex FFT.
73 *
74 * Parameters:
75 *
76 * N : int / length of sequences.
77 * T : COMPLEX[N] / input sequence.
78 * F : COMPLEX[N] / transform sequence. may point to same memory as T.
79 *
80 * F(n) = sum(T(k)*exp(-i*2*pi*n*k/N)), k=0..N-1,n=0..N-1
81 *
82 *------------------------------------------------------------------------------
83 *
84 * cfftb: inverse complex FFT.
85 *
86 * Parameters:
87 *
88 * N : int / length of sequences.
89 * F : COMPLEX[N] / transform sequence.
90 * T : COMPLEX[N] / output sequence. may point to same memory as F.
91 *
92 * T(n) = sum(F(k)*exp(i*2*pi*n*k/N)), k=0..N-1,N=0..N-1
93 *
94 *------------------------------------------------------------------------------
95 *
96 * rfftf: forward real FFT. First (N/2+1) coefficients returned.
97 * F(n) = Real(F(N-n))-i*Imag(F(N-n)), n=N/2+1..N-1.
98 *
99 * Parameters:
100 *
101 * N : int / length of input sequence.
102 * T : double[N] / input sequence.
103 * F : COMPLEX[N/2+1] / transform sequence. may point to same memory
104 * / as T (which must then be of length (N+2)).
105 *
106 * F(n) = sum(T(k)*exp(-i*2*pi*n*k/N)), k=0..N-1,n=0..N/2
107 *
108 * real transform as half size complex transform:
109 *
110 * f(k) real
111 * F(n) = (1/N)*sum(f(k)*exp(-i*2*pi*n*k/N)), k=0..N-1,n=0..N-1
112 * Divide to two parts:
113 * g(k) = f(2*k), k=0..N/2-1
114 * h(k) = f(2*k+1), k=0..N/2-1
115 * then
116 * 2*G(n) = (1/(N/2))*sum(g(l)*exp(-i*2*pi*n*l/(N/2)),l=0..N/2-1,(k=2*l)
117 * 2*H(n) = (1/(N/2))*sum(h(l)*exp(-i*2*pi*n*l/(N/2)),l=0..N/2-1,(k=2*l+1)
118 * and
119 * F(n) = G(n) + exp(-i*2*pi*n/N)*H(n), n=0..N/2-1
120 * Now take
121 * y(k) = g(k) + i*h(k), k=0..N/2-1
122 * then
123 * Y(n) = (1/(N/2))*sum(y(k)*exp(-i*2*pi*n*k/(N/2))), k=0..N/2-1, n=0..N/2-1
124 * = G(n) + i*H(n)
125 * = G_even(n) + i*G_odd(n) + i*H_even(n) - H_odd(n)
126 * where
127 * G_even is even part of G, G_odd odd part of G,
128 * H_even is even part of H, and H_odd odd part of H.
129 * now
130 * 2*G_even(n) = Real( Y(n) + Y(N/2-n))
131 * 2*G_odd(n) = Imag( Y(n) - Y(N/2-n))
132 * 2*H_even(n) = Imag( Y(n) + Y(N/2-n))
133 * 2*H_odd(n) = Real(-Y(n) + Y(N/2-n))
134 * and (at last)
135 * F(n)=G_even(n)+i*G_odd(n)+exp(-i*2*pi*n/N)*(H_even(n)+i*H_odd(n)), n=0..N/2-1
136 * F(n)=Real(F(N-n))-i*Imag(F(N-n)), n=N/2+1..N-1
137 *
138 *------------------------------------------------------------------------------
139 *
140 * rfftb: inverse real FFT.
141 *
142 * Parameters:
143 *
144 * N : int / length of output sequence.
145 * F : COMPLEX[N/2+1] / transform sequence.
146 * T : double[N] / output sequence. may point to same memory as F.
147 *
148 * T(k) = sum(F(n)*exp(i*2*pi*k*n/N)), n=0..N-1,N=0..N-1
149 *
150 * Transform of complex input where the other half is complex conjugates of the
151 * other (that is the real part is even, complex part odd, and the result is
152 * real) as half size transform:
153 *
154 * F(N-n) = Real(F(n))-i*Imag(F(n))
155 * f(k) = sum(F(n)*exp(i*2*pi*k*n/N)), n=0..N-1, k=0..N-1
156 * Divide to two parts:
157 * G(n) = F(2*n), n=0..N/2-1
158 * H(n) = F(2*n+1)-F(2*n-1), n=0..N/2-1
159 * now
160 * h(k) = (exp(i*2*pi*k/N)-exp(-i*2*pi*k/N))*x(k)
161 * = i*2*sin(2*pi*k/N)*x(k)
162 * where
163 * X(n) = F(2*n+1), k=0..N/2-1
164 * x(k) = sum(X(l)*exp(i*2*pi*k*l/(N/2))), l=0..N/2-1,(n=2*l+1)
165 * g(k) = sum(G(l)*exp(i*2*pi*k*l/(N/2))), l=0..N/2-1,(n=2*l)
166 * take
167 * Y(n) = G(n) + H(n)
168 * so
169 * y(k) = g(k) + h(k)
170 * = g(k) + x(k) * (i*2*sin(2*pi*k/N)), k=0..N/2-1
171 * if
172 * 2*g_even(k) = Real( y(k) + y(N/2-k))
173 * 2*g_odd(k) = Real(-y(k) + y(N/2-k))
174 * 2*h_odd(k) = Imag( y(k) - y(N/2-k))
175 * 2*h_even(k) = Imag( y(k) + y(N/2-k))
176 * k=0..N/2-1
177 * then
178 * f(k) = g(k) + x(k)
179 * = g_even(k)+g_odd(k) -
180 * (h_odd(k)+h_evek(k))/(2*sin(2*pi*k/N)), k=0..N/2-1
181 * = g_even(N-k)-g_odd(N-k) -
182 * (-h_odd(N-k)+h_even(N-k))/(2*sin(2*pi*k/N)), k=N/2+1..N-1
183 *
184 *------------------------------------------------------------------------------
185 *
186 * cfftf2D: 2D forward complex FFT.
187 *
188 * Parameters:
189 *
190 * M,N : int / array dimensions.
191 * T : COMPLEX[M][N] / input array.
192 * F : COMPLEX[M][N] / transform array. may point to same memory as T.
193 *
194 * F(m,n) = sum(T(k,l)*exp(-i*2*pi*(m*k/M+n*l/N))), k,m=0..M-1; l,n=0..N-1
195 *
196 *------------------------------------------------------------------------------
197 *
198 * cfftb2D: 2D inverse complex FFT.
199 *
200 * Parameters:
201 *
202 * M,N : int / array dimensions.
203 * F : COMPLEX[M][N] / transform array.
204 * T : COMPLEX[M][N] / output array. may point to same memory as T.
205 *
206 * T(m,n) = sum(F(k,l)*exp(i*2*pi*(m*k/M+n*l/N))), k,m=0..M-1; l,n=0..N-1
207 *
208 *------------------------------------------------------------------------------
209 *
210 * cfftf3D: 3D forward complex FFT.
211 *
212 * Parameters:
213 *
214 * L,M,N : int / array dimensions.
215 * T : COMPLEX[L][M][N] / input array.
216 * F : COMPLEX[L][M][N] / transform array. may point to same memory as T.
217 *
218 * F(l,m,n) = sum(T(i,j,k)*exp(-i*2*pi*(l*i/L+m*j/M+n*k/N))),
219 * l,i=0..L-1; m,j=0..M-1; k,n=0..N-1
220 *
221 *------------------------------------------------------------------------------
222 *
223 * cfftb3D: 3D inverse complex FFT.
224 *
225 * Parameters:
226 *
227 * L,M,N : int / array dimensions.
228 * F : COMPLEX[L][M][N] / transform array.
229 * T : COMPLEX[L][M][N] / output array. may point to same memory as T.
230 *
231 * T(l,m,n) = sum(F(i,j,k)*exp(i*2*pi*(l*i/L+m*j/M+n*k/N))),
232 * l,i=0..L-1; m,j=0..M-1; k,n=0..N-1
233 *
234 *------------------------------------------------------------------------------
235 *
236 * cfftfND: multidimensional (max 32) forward complex FFT.
237 *
238 * Parameters:
239 *
240 * N : int / number of dimensions.
241 * D : int[N] / array dimensions.
242 * T : COMPLEX T[D[N-1]][D[N-2]]...[D[0]] / input array.
243 * F : COMPLEX F[D[N-1]][D[N-2]]...[D[0]] / transform array.
244 * / may point to same memory as T.
245 *
246 *------------------------------------------------------------------------------
247 *
248 * cfftbND: multidimensional inverse complex FFT.
249 *
250 * Parameters:
251 *
252 * N : int / number of dimensions.
253 * D : int[N] / array dimensions.
254 * F : COMPLEX T[D[N-1]][D[N-2]]...[D[0]] / transform array.
255 * T : COMPLEX F[D[N-1]][D[N-2]]...[D[0]] / output array.
256 * / may point to same memory as T.
257 *
258 *------------------------------------------------------------------------------
259 *
260 * in all the routines above F and T may point to same memory, if the memory
261 * is allocated according to larger of the arrays.
262 *
263 * The COMPLEX type is defined as:
264 *
265 * typedef
266 * struct {
267 * double Real;
268 * double Imag;
269 * } COMPLEX;
270 *
271 *------------------------------------------------------------------------------
272 *
273 * gfftf( nT, time, nF, freq ) - return largest magnitude coeff. (real input)
274 * int nT;
275 * double time[nT + 2];
276 * int nF;
277 * FREQ freq[nF];
278 *
279 *------------------------------------------------------------------------------
280 *
281 * gfftb( nF, freq, nT, time ) - given largest magnitude coeff. return
282 * int nF; approximation of original sequence.
283 * FREQ freq[nF];
284 * int nT;
285 * double time[nT + 2];
286 *
287 *------------------------------------------------------------------------------
288 *
289 * The FREQ type is defined as:
290 *
291 * typedef
292 * struct {
293 * double Real;
294 * double Imag;
295 * double Mag;
296 * int FBin;
297 * } FREQ;
298 *
299 * The magnitude entry is returned by gfftf(), not used by gfftb().
300 *
301 * Juha Ruokolainen / CSC
302 * LAST CHANGE: 19.3. -92
303 *
304 */
305
306 #include <stdlib.h>
307 #include <string.h>
308 #include <math.h>
309 #include <stdio.h>
310
311 #include "../config.h"
312
313 #define FALSE 0
314 #define TRUE 1
315
316 typedef
317 struct {
318 double Real;
319 double Imag;
320 } COMPLEX;
321
322 typedef
323 struct {
324 double Real;
325 double Imag;
326 double Mag;
327 int FBin;
328 } FREQ;
329
330 static double _FFT_PI = 3.14159265358979323844;
331
332 #define _FFT_MAX_LEVELS 30
333
334 static double _FFT_SIN[_FFT_MAX_LEVELS];
335 static double _FFT_COS[_FFT_MAX_LEVELS];
336
337 static int _FFT_I;
338 static int _FFT_Level;
339
340 /*
341 * log2: base two logarithm of an (power of two!!!) integer.
342 * the macro needs and uses variable k.
343 *
344 * Parameters:
345 *
346 * R : int / put result here.
347 * N : int / input number.
348 *
349 */
350 #define log2( R, N ) \
351 k = 1; \
352 for( R = 0; R < 32; R++ ) { \
353 if ( k & N ) break; \
354 k <<= 1; \
355 }
356
357 /*
358 * BitReverse: reverse the order of given number of lowest bits in an integer.
359 * the macro needs and uses variables i & j.
360 *
361 * Parameters:
362 *
363 * R: int / put result here.
364 * N: int / the input number.
365 * Nbits: int / number of bits to be reversed - 1.
366 *
367 */
368 #define BitReverse( R, N, Nbits ) \
369 R = 0; \
370 j = 1; \
371 for( i = 0; i <= Nbits; i++ ) { \
372 if ( j & N ) { \
373 R |= 1 << ( Nbits - i ); \
374 } \
375 j <<= 1; \
376 }
377
378 /*
379 * BitReverseArray: clear the messed up order in an array,
380 * the change is done in place.
381 *
382 * Parameters:
383 *
384 * N: int / number of elements in the array, and
385 * / k=log(N) (the logarithm is base two),
386 * / the number of bits.
387 * T: COMPLEX[N] / the array to be operated on.
388 */
BitReverseArray(N,T)389 void BitReverseArray( N, T )
390 int N;
391 COMPLEX *T;
392 {
393 COMPLEX swap;
394
395 int i;
396 int j;
397
398 int k;
399 int n;
400
401 int logN;
402
403 log2( logN, N );
404 logN--;
405
406 for( k = 0; k < N; k++ ) {
407 BitReverse( n, k, logN );
408 if ( n > k ) {
409 swap = T[k];
410 T[k] = T[n];
411 T[n] = swap;
412 }
413 }
414 }
415
416 /*
417 * sort: sort an (double) array to ascending order, and move the elements of
418 * another (integer) array accordingly. the latter can be used as track
419 * keeper of where an element in the sorted order at position (k) was in
420 * in the original order (Ord[k]), if it is initialized to contain
421 * numbers (0..N-1) before calling sort.
422 *
423 * Parameters:
424 *
425 * N: int / number of entries in the arrays.
426 * Key: double[N] / array to be sorted.
427 * Ord: int[N] / change this accordingly.
428 */
sort_swap(int i,int j,double * Key,int * Ord)429 void sort_swap(int i,int j,double *Key,int *Ord)
430 {
431 int ival;
432 double dval;
433
434 dval=Key[i];
435 Key[i]=Key[j];
436 Key[j]=dval;
437
438 ival=Ord[i];
439 Ord[i]=Ord[j];
440 Ord[j]=ival;
441 }
442
sort_shift(int lbeg,int lend,double * Key,int * Ord)443 void sort_shift(int lbeg,int lend,double *Key,int *Ord)
444 {
445 int i,j;
446
447 i = lbeg;
448 while( 2*i+1<=lend )
449 {
450 j=2*i+1;
451 if ( j<lend && Key[j]<Key[j+1]) j=j+1;
452 if ( Key[i]<Key[j]) {
453 sort_swap(i,j,Key,Ord);
454 i = j;
455 } else break;
456 }
457
458 }
sort(N,Key,Ord)459 void sort( N, Key, Ord )
460 int N;
461 int *Ord;
462 double *Key;
463 {
464 int lend,lbeg;
465
466 lend = N-1;
467 lbeg = (lend-1)/2;
468
469 while(lbeg>=0)
470 {
471 sort_shift( lbeg--,lend,Key,Ord );
472 }
473
474 while( lend>0) {
475 sort_swap( 0,lend,Key,Ord );
476 sort_shift( 0,--lend,Key,Ord );
477 }
478 }
479
480 /*
481 * FFTInit: initialize internal sin & cos tables for (pi / N)
482 * N being power of two up to _FFT_MAX_LEVELS.
483 *
484 * Parameters: NONE.
485 */
FFTInit()486 static int FFTInit( )
487 {
488 static int InitDone = FALSE;
489
490 int k;
491 int n;
492
493 if ( InitDone ) {
494 return 0;
495 }
496
497 n = ( 1 << _FFT_MAX_LEVELS );
498 for( k = 0; k < _FFT_MAX_LEVELS; k++ ) {
499 _FFT_COS[k] = cos( _FFT_PI / n );
500 _FFT_SIN[k] = -sin( _FFT_PI / n );
501 n /= 2;
502 }
503
504 InitDone = TRUE;
505 }
506
507
508 /*
509 * FFTKernel: discrete fourier transform. output is in bitreversed order.
510 *
511 * Parameters:
512 *
513 * N : int / length of sequences.
514 * F : COMPLEX[N] / input sequence.
515 * T : COMPLEX[N] / transform sequence.
516 * / may point to same memory as F.
517 *
518 * T(n) = sum(F(k)*exp(-i*2*pi*n*k/N)), k=0..N-1,n=0..N-1
519 *
520 */
FFTKernel(N,F,T)521 static int FFTKernel( N, F, T )
522 int N;
523 COMPLEX *F;
524 COMPLEX *T;
525 {
526 double ExpR;
527 double ExpI;
528
529 double CO;
530 double SI;
531
532 double TempR;
533 double TempI;
534
535 double t;
536
537 int k;
538
539 if ( N == 4 ) {
540 TempR = F[0].Real; TempI = F[0].Imag;
541
542 F[0].Real = TempR + F[2].Real;
543 F[0].Imag = TempI + F[2].Imag;
544
545 F[2].Real = TempR - F[2].Real;
546 F[2].Imag = TempI - F[2].Imag;
547
548 TempR = F[1].Real; TempI = F[1].Imag;
549
550 F[1].Real = TempR + F[3].Real;
551 F[1].Imag = TempI + F[3].Imag;
552
553 TempR -= F[3].Real;
554 TempI -= F[3].Imag;
555
556 F[3].Real = TempI;
557 F[3].Imag = -TempR;
558
559 TempR = F[0].Real; TempI = F[0].Imag;
560
561 T[_FFT_I].Real = TempR + F[1].Real;
562 T[_FFT_I].Imag = TempI + F[1].Imag;
563 _FFT_I++;
564
565 T[_FFT_I].Real = TempR - F[1].Real;
566 T[_FFT_I].Imag = TempI - F[1].Imag;
567 _FFT_I++;
568
569 TempR = F[2].Real; TempI = F[2].Imag;
570
571 T[_FFT_I].Real = TempR + F[3].Real;
572 T[_FFT_I].Imag = TempI + F[3].Imag;
573 _FFT_I++;
574
575 T[_FFT_I].Real = TempR - F[3].Real;
576 T[_FFT_I].Imag = TempI - F[3].Imag;
577 _FFT_I++;
578
579 return 0;
580 }
581
582 N /= 2;
583
584 /*
585 * ExpR + i*ExpI = exp(-i*2*pi*k/N ) = exp(-i*2*pi/N)^(k)
586 */
587 ExpR = CO = _FFT_COS[_FFT_Level];
588 ExpI = SI = _FFT_SIN[_FFT_Level];
589
590 TempR = F[0].Real;
591 TempI = F[0].Imag;
592
593 F[0].Real = TempR + F[N].Real;
594 F[0].Imag = TempI + F[N].Imag;
595
596 F[N].Real = TempR - F[N].Real;
597 F[N].Imag = TempI - F[N].Imag;
598
599 for( k = 1; k < N; k++ ) {
600 TempR = F[k].Real - F[k + N].Real;
601 TempI = F[k].Imag - F[k + N].Imag;
602
603 F[k].Real += F[k + N].Real;
604 F[k].Imag += F[k + N].Imag;
605
606 F[k + N].Real = TempR * ExpR - TempI * ExpI;
607 F[k + N].Imag = TempR * ExpI + TempI * ExpR;
608
609 t = ExpR;
610 ExpR = CO * t - SI * ExpI;
611 ExpI = CO * ExpI + SI * t;
612 }
613
614 _FFT_Level++;
615
616 FFTKernel( N, &F[0], T );
617 FFTKernel( N, &F[N], T );
618
619 _FFT_Level--;
620 }
621
622 /*
623 * cfftf: forward complex FFT.
624 *
625 * Parameters:
626 *
627 * N : int / length of sequences.
628 * T : COMPLEX[N] / input sequence.
629 * F : COMPLEX[N] / transform sequence. may point to same memory as T.
630 *
631 * F(n) = sum(T(k)*exp(-i*2*pi*n*k/N)), k=0..N-1,n=0..N-1
632 */
cfftf(N,T,F)633 void cfftf( N, T, F )
634 int N;
635 COMPLEX *T;
636 COMPLEX *F;
637 {
638 int logN;
639 int k;
640
641 FFTInit( );
642
643 log2( logN, N );
644 _FFT_Level = _FFT_MAX_LEVELS - logN + 1;
645
646 _FFT_I = 0;
647
648 if ( F != T ) {
649 for( k = 0; k < N; k++ ) F[k] = T[k];
650 }
651
652 FFTKernel( N, F, F );
653 BitReverseArray( N, F );
654 }
655
656 /*
657 * cfftb: inverse complex FFT.
658 *
659 * Parameters:
660 *
661 * N : int / length of sequences.
662 * F : COMPLEX[N] / transform sequence.
663 * T : COMPLEX[N] / output sequence. may point to same memory as F.
664 *
665 * T(n) = sum(F(k)*exp(i*2*pi*n*k/N)), k=0..N-1,N=0..N-1
666 *
667 */
cfftb(N,F,T)668 void cfftb( N, F, T )
669 int N;
670 COMPLEX *F;
671 COMPLEX *T;
672 {
673 int logN;
674 int k;
675
676 if ( F != T ) {
677 for( k = 0; k < N; k++ ) T[k].Real = F[k].Real;
678 }
679 for( k = 0; k < N; k++ ) T[k].Imag = -F[k].Imag;
680 cfftf( N, T, T );
681 for( k = 0; k < N; k++ ) T[k].Imag = -T[k].Imag;
682 }
683
684 #ifdef USE_ISO_C_BINDINGS
fcfftb(N,F,T)685 void fcfftb( N, F, T )
686 int *N;
687 COMPLEX *F;
688 COMPLEX *T;
689 {
690 cfftb( *N, F, T );
691 }
692 #else
693 void FC_FUNC(fcfftb,FCFFTB)( N, F, T )
694 int *N;
695 COMPLEX *F;
696 COMPLEX *T;
697 {
698 cfftb( *N, F, T );
699 }
700 #endif /* USE_ISO_C_BINDINGS */
701
702 /*
703 * rfftf: forward real FFT. First (N/2+1) coefficients returned.
704 * F(n) = Real(F(N-n))-i*Imag(F(N-n)), n=N/2+1..N-1.
705 *
706 * parameters:
707 *
708 * N : int / length of input sequence.
709 * T : double[N] / input sequence.
710 * F : COMPLEX[N/2+1] / transform sequence. may point to same memory
711 * / as T (which must then be of length (N+2)).
712 *
713 * F(n) = sum(T(k)*exp(-i*2*pi*n*k/N)), k=0..N-1,n=0..N/2
714 *
715 */
rfftf(N,T,F)716 void rfftf( N, T, F )
717 int N;
718 double *T;
719 COMPLEX *F;
720 {
721 COMPLEX *W;
722
723 double CO;
724 double SI;
725
726 double ExpR;
727 double ExpI;
728
729 double pi;
730 double t;
731
732 int k;
733
734 N /= 2;
735
736 W = (COMPLEX *)malloc( ( N + 1 ) * sizeof( COMPLEX ) );
737 /*
738 * if you like (or change the COMPLEX type) uncomment the following
739 *
740 * for( k = 0; k < N; k++ ) {
741 * W[k].Real = T[2*k];
742 * W[k].Imag = T[2*k+1];
743 * }
744 *
745 * and replace
746 * cfftf( N, (COMPLEX *)T, W );
747 * by
748 * cfftf( N, W, W );
749 */
750 cfftf( N, (COMPLEX *)T, W );
751 W[N] = W[0];
752
753 /*
754 * ExpR + i*ExpI = exp(-i*2*pi*k/N ) = exp(-i*2*pi/N)^(k)
755 */
756 ExpR = 1.0;
757 ExpI = 0.0;
758
759 pi = _FFT_PI / N;
760 CO = cos( pi );
761 SI = -sin( pi );
762
763 for( k = 0; k <= N; k++ ) {
764 F[k].Real = W[k].Imag + W[N - k].Imag;
765 F[k].Imag = -W[k].Real + W[N - k].Real;
766
767 t = F[k].Real;
768 F[k].Real = ExpR * t - ExpI * F[k].Imag;
769 F[k].Imag = ExpR * F[k].Imag + ExpI* t;
770
771 F[k].Real += W[k].Real + W[N - k].Real;
772 F[k].Imag += W[k].Imag - W[N - k].Imag;
773
774 F[k].Real /= 2; F[k].Imag /= 2;
775
776 t = ExpR;
777 ExpR = CO * t - SI * ExpI;
778 ExpI = CO * ExpI + SI * t;
779 }
780
781 free( W );
782 }
783
784 /*
785 * rfftb: inverse real FFT.
786 *
787 * Parameters:
788 *
789 * N : int / length of output sequence.
790 * F : COMPLEX[N/2+1] / transform sequence.
791 * T : double[N] / output sequence. may point to same memory as F.
792 *
793 * T(n) = sum(F(k)*exp(i*2*pi*n*k/N)), k=0..N-1,N=0..N-1
794 *
795 */
rfftb(N,F,T)796 void rfftb( N, F, T )
797 int N;
798 COMPLEX *F;
799 double *T;
800 {
801 COMPLEX *W;
802
803 double ExpR;
804 double ExpI;
805
806 double CO;
807 double SI;
808
809 double pi;
810
811 double t;
812
813 double Eves;
814 double Odds;
815
816 int k;
817 int n;
818
819 N /= 2;
820
821 W = (COMPLEX *)malloc( ( N + 1 ) * sizeof( COMPLEX ) );
822
823 W[0].Imag = F[0].Imag + 2 * F[1].Imag;
824 W[0].Real = F[0].Real;
825
826 W[N / 2].Real = F[N].Real;
827 W[N / 2].Imag = F[N].Imag - 2 * F[N - 1].Imag;
828
829 for( k = 1; k < N / 2; k++ ) {
830 n = 2 * k;
831 W[k].Real = F[n].Real + F[n + 1].Real - F[n - 1].Real;
832 W[k].Imag = F[n].Imag + F[n + 1].Imag - F[n - 1].Imag;
833 }
834
835 for( k = N / 2 + 1; k < N; k++ ) {
836 n = 2 * ( N - k );
837 W[k].Real = F[n].Real + F[n - 1].Real - F[n + 1].Real;
838 W[k].Imag = -( F[n].Imag + F[n - 1].Imag - F[n + 1].Imag );
839 }
840
841 Odds = F[1].Real;
842 Eves = 0.0;
843 for( k = 1; k < N / 2; k++ ) {
844 n = 2 * k;
845 Eves += F[n].Real;
846 Odds += F[n + 1].Real;
847 }
848 Odds *= 2;
849 Eves *= 2;
850 Eves += F[0].Real + F[N].Real;
851
852 cfftb( N, W, W );
853 W[N] = W[0];
854
855 /*
856 * ExpR + i*ExpI = exp(i*2*pi*k/N ) = exp(i*2*pi/N)^(k)
857 */
858 ExpR = 1.0;
859 ExpI = 0.0;
860
861 pi = _FFT_PI / N;
862 CO = cos( pi );
863 SI = sin( pi );
864
865 for( k = 1; k < N; k++ ) {
866 t = ExpR;
867 ExpR = CO * t - SI * ExpI;
868 ExpI = CO * ExpI + SI * t;
869
870 n = 2 * N - k;
871 T[n] = T[k] = 0.5 / ExpI;
872
873 T[k] *= -W[k].Imag;
874 T[k] += W[k].Real;
875
876 T[n] *= W[N - k].Imag;
877 T[n] += W[N - k].Real;
878 }
879 T[0] = Eves + Odds;
880 T[N] = Eves - Odds;
881
882 free( W );
883 }
884
885 #ifdef USE_ISO_C_BINDINGS
frfftb(N,F,T)886 void frfftb( N, F, T )
887 int *N;
888 COMPLEX *F;
889 double *T;
890 {
891 rfftb( *N, F, T );
892 }
893 #else
894 FC_FUNC(frfftb,FRFFTB)( N, F, T )
895 int *N;
896 COMPLEX *F;
897 double *T;
898 {
899 rfftb( *N, F, T );
900 }
901 #endif /* USE_ISO_C_BINDINGS */
902
903 /*
904 *
905 * gfftb: given group of frequency coefficients (usually those with largest
906 * magnitude), compute approximation of the original sequence.
907 *
908 * Parameters:
909 *
910 * nF : int / number of input coefficients.
911 * freq : FREQ[nF] / input coefficients.
912 * nT : int / length of output sequence.
913 * time : time[nT+2] / output sequence.
914 *
915 * FREQ :
916 *
917 * struct {
918 * double Real;
919 * double Imag;
920 * double Mag;
921 * int FBin;
922 * }
923 *
924 */
gfftb(nF,freq,nT,time)925 void gfftb( nF, freq, nT, time )
926 int nT;
927 int nF;
928 double *time;
929 FREQ *freq;
930 {
931 COMPLEX *F;
932
933 int k;
934
935 /*
936 * if you change the COMPLEX type replace
937 * F = (COMPLEX *)time;
938 * bzero( F, ( nT / 2 + 1 ) * sizeof( COMPLEX ) );
939 * by
940 * F = (COMPLEX *)calloc( ( nT / 2 + 1, sizeof( COMPLEX ) );
941 * and add
942 * free( F );
943 * to end of the routine.
944 */
945 F = (COMPLEX *)time;
946 /* bzero( F, ( nT / 2 + 1 ) * sizeof( COMPLEX ) ); */
947 memset( F, 0, ( nT / 2 + 1 ) * sizeof( COMPLEX ) );
948
949 for( k = 0; k < nF; k++ ) {
950 F[freq[k].FBin].Real = freq[k].Real;
951 F[freq[k].FBin].Imag = freq[k].Imag;
952 }
953
954 rfftb( nT, F, time );
955 }
956
957
958 /*
959 *
960 * gfftf: given real sequence compute nF frequency coefficients
961 * of largest magnitude.
962 *
963 * Parameters:
964 *
965 * nT : int / length of input sequence.
966 * time : double[nT] / input sequence.
967 * nF : int / number of coefficients wanted.
968 * freq : FREQ[nF] / output coefficients.
969 *
970 * FREQ :
971 *
972 * struct {
973 * double Real;
974 * double Imag;
975 * double Mag;
976 * int FBin;
977 * }
978 */
gfftf(nT,time,nF,freq)979 void gfftf( nT, time, nF, freq )
980 int nT;
981 int nF;
982 FREQ *freq;
983 double *time;
984 {
985 COMPLEX *F;
986
987 double *Magnitude;
988
989 int i;
990 int k;
991 int *Order;
992
993 nT /= 2;
994
995 F = (COMPLEX *)malloc( ( nT + 1 ) * sizeof( COMPLEX ) );
996
997 rfftf( 2 * nT, time, F );
998
999 /************************************************
1000 * search for largest magnitude coefficients.. *
1001 ************************************************/
1002 Magnitude = (double *)malloc( ( nT + 1 ) * sizeof( double ) );
1003 Order = (int *)malloc( ( nT + 1 ) * sizeof( int ) );
1004
1005 for( k = 0; k <= nT; k++ ) {
1006 Magnitude[k] = F[k].Real * F[k].Real + F[k].Imag * F[k].Imag;
1007 Order[k] = k;
1008 }
1009 sort( nT + 1, Magnitude, Order );
1010
1011 k = nT;
1012 for( i = 0; i < nF; i++, k-- ) {
1013 freq[i].Real = F[Order[k]].Real;
1014 freq[i].Imag = F[Order[k]].Imag;
1015 freq[i].Mag = Magnitude[k];
1016 freq[i].FBin = Order[k];
1017 }
1018
1019 free( F );
1020 free( Order );
1021 free( Magnitude );
1022 }
1023
1024 /*
1025 * cfftf2D: 2D forward complex FFT.
1026 *
1027 * Parameters:
1028 *
1029 * M,N : int / array dimensions.
1030 * T : COMPLEX[M][N] / input array.
1031 * F : COMPLEX[M][N] / transform array. may point to same memory as T.
1032 *
1033 * F(m,n) = sum(T(k,l)*exp(-i*2*pi*(m*k/M+n*l/N))), k,m=0..M-1; l,n=0..N-1
1034 */
cfftf2D(M,N,T,F)1035 void cfftf2D( M, N, T, F )
1036 int M,N;
1037 COMPLEX *T;
1038 COMPLEX *F;
1039 {
1040 COMPLEX *W;
1041
1042 int l;
1043 int k;
1044 int n;
1045
1046 W = (COMPLEX *)malloc( M * sizeof( COMPLEX ) );
1047
1048 for( k = 0, n = 0; k < M; k++, n += N )
1049 cfftf( N, &T[n], &F[n] );
1050
1051 for( l = 0; l < N; l++ ) {
1052 for( k = 0, n = l; k < M; k++, n += N ) W[k] = F[n];
1053
1054 cfftf( M, W, W );
1055
1056 for( k = 0, n = l; k < M; k++, n += N ) F[n] = W[k];
1057 }
1058
1059 free( W );
1060 }
1061
1062 /*
1063 * cfftb2D: 2D inverse complex FFT.
1064 *
1065 * Parameters:
1066 *
1067 * M,N : int / array dimensions.
1068 * F : COMPLEX[M][N] / transform array.
1069 * T : COMPLEX[M][N] / output array. may point to same memory as T.
1070 *
1071 * T(m,n) = sum(F(k,l)*exp(i*2*pi*(m*k/M+n*l/N))), k,m=0..M-1; l,n=0..N-1
1072 */
cfftb2D(M,N,F,T)1073 void cfftb2D( M, N, F, T )
1074 int M,N;
1075 COMPLEX *F;
1076 COMPLEX *T;
1077 {
1078 COMPLEX *W;
1079
1080 int k;
1081
1082 if ( T != F ) {
1083 for( k = 0; k < M * N; k++ ) T[k].Real = F[k].Real;
1084 }
1085
1086 for( k = 0; k < M * N; k++ ) T[k].Imag = -F[k].Imag;
1087
1088 cfftf2D( M, N, T, T );
1089
1090 for( k = 0; k < M * N; k++ ) T[k].Imag = -T[k].Imag;
1091 }
1092
1093 /*
1094 * cfftf3D: 3D forward complex FFT.
1095 *
1096 * Parameters:
1097 *
1098 * L,M,N : int / array dimensions.
1099 * T : COMPLEX[L][M][N] / input array.
1100 * F : COMPLEX[L][M][N] / transform array. may point to same memory as T.
1101 *
1102 * F(l,m,n) = sum(T(i,j,k)*exp(-i*2*pi*(l*i/L+m*j/M+n*k/N))),
1103 * l,i=0..L-1; m,j=0..M-1; k,n=0..N-1
1104 */
cfftf3D(L,M,N,T,F)1105 void cfftf3D( L, M, N, T, F )
1106 int L,M,N;
1107 COMPLEX *T;
1108 COMPLEX *F;
1109 {
1110 COMPLEX *W;
1111
1112 int j;
1113 int k;
1114 int n;
1115
1116 W = (COMPLEX *)malloc( L * sizeof( COMPLEX ) );
1117
1118 for( k = 0; k < L; k++ ) cfftf2D( M, N, &T[k*M*N], &F[k*M*N] );
1119
1120 for( j = 0; j < M * N; j++ ) {
1121 for( k = 0, n = j; k < L; k++, n += M * N ) W[k] = F[n];
1122
1123 cfftf( L, W, W );
1124
1125 for( k = 0, n = j; k < L; k++, n += M * N ) F[n] = W[k];
1126 }
1127
1128 free( W );
1129 }
1130
1131 /*
1132 * cfftb3D: 3D inverse complex FFT.
1133 *
1134 * Parameters:
1135 *
1136 * L,M,N : int / array dimensions.
1137 * F : COMPLEX[L][M][N] / transform array.
1138 * T : COMPLEX[L][M][N] / output array. may point to same memory as T.
1139 *
1140 * T(l,m,n) = sum(F(i,j,k)*exp(i*2*pi*(l*i/L+m*j/M+n*k/N))),
1141 * l,i=0..L-1; m,j=0..M-1; k,n=0..N-1
1142 */
cfftb3D(L,M,N,F,T)1143 void cfftb3D( L, M, N, F, T )
1144 int L,M,N;
1145 COMPLEX *F;
1146 COMPLEX *T;
1147 {
1148 int k;
1149
1150 if ( T != F ) {
1151 for( k = 0; k < L * M * N; k++ ) T[k].Real = F[k].Real;
1152 }
1153
1154 for( k = 0; k < L * M * N; k++ ) T[k].Imag = -F[k].Imag;
1155
1156 cfftf3D( L, M, N, T, T );
1157
1158 for( k = 0; k < L * M * N; k++ ) T[k].Imag = -T[k].Imag;
1159 }
1160
1161 /*
1162 * cfftfND: multidimensional (max 32) forward complex FFT.
1163 *
1164 * Parameters:
1165 *
1166 * N : int / number of dimensions.
1167 * D : int[N] / array dimensions.
1168 * T : COMPLEX T[D[N-1]][D[N-2]]...[D[0]] / input array.
1169 * F : COMPLEX F[D[N-1]][D[N-2]]...[D[0]] / transform array.
1170 * / may point to same memory as T.
1171 */
cfftfND(N,D,T,F)1172 void cfftfND( N, D, T, F )
1173 int N;
1174 int *D;
1175 COMPLEX *T;
1176 COMPLEX *F;
1177 {
1178 COMPLEX *W;
1179
1180 int TotN;
1181 int WN;
1182
1183 int i;
1184 int j;
1185 int k;
1186 int l;
1187 int m;
1188
1189 int S[32];
1190 int P[32];
1191
1192 WN = D[0];
1193 TotN = 1;
1194 for( i = 0; i < N; i++ ) {
1195 if ( D[i] > WN ) WN = D[i];
1196 S[i] = TotN;
1197 TotN *= D[i];
1198 }
1199
1200 W = (COMPLEX *)malloc( WN * sizeof( COMPLEX ) );
1201
1202 if ( F != T )
1203 for( i = 0; i < TotN; i++ )
1204 F[i] = T[i];
1205
1206 for( i = 0; i < N; i++ ) {
1207
1208 for( j = 0; j < N; j++ ) P[j] = 0;
1209
1210 k = 0;
1211 for( j = 0; j < TotN / D[i]; j++ ) {
1212
1213 if ( j != 0 )
1214 for( l = 0; l < N; l++ )
1215 if ( l != i ) {
1216 P[l]++;
1217 k += S[l];
1218 if ( P[l] == D[l] ) {
1219 k -= S[l + 1];
1220 P[l] = 0;
1221 } else {
1222 break;
1223 }
1224 }
1225
1226 for( l = 0, m = k; l < D[i]; l++, m += S[i] ) W[l] = F[m];
1227
1228 cfftf( D[i], W, W );
1229
1230 for( l = 0, m = k; l < D[i]; l++, m += S[i] ) F[m] = W[l];
1231 }
1232 }
1233
1234 free( W );
1235 }
1236
1237
1238 /*
1239 * cfftbND: multidimensional inverse complex FFT.
1240 *
1241 * Parameters:
1242 *
1243 * N : int / number of dimensions.
1244 * D : int[N] / array dimensions.
1245 * F : COMPLEX T[D[N-1]][D[N-2]]...[D[0]] / transform array.
1246 * T : COMPLEX F[D[N-1]][D[N-2]]...[D[0]] / output array.
1247 * / may point to same memory as T.
1248 */
cfftbND(N,D,F,T)1249 void cfftbND( N, D, F, T )
1250 int N;
1251 int *D;
1252 COMPLEX *F;
1253 COMPLEX *T;
1254 {
1255 int TotN;
1256 int k;
1257
1258 TotN = D[0];
1259 for( k = 1; k < N; k++ ) TotN *= D[k];
1260
1261 if ( T != F ) {
1262 for( k = 0; k < TotN; k++ ) T[k].Real = F[k].Real;
1263 }
1264
1265 for( k = 0; k < TotN; k++ ) T[k].Imag = -F[k].Imag;
1266
1267 cfftfND( N, D, T, T );
1268
1269 for( k = 0; k < TotN; k++ ) T[k].Imag = -T[k].Imag;
1270 }
1271