source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/core/fft.c

Last change on this file was 10371, checked in by aboeckma, 11 years ago

updated mafft version. Added extensions (no svn ignore, yet)

File size: 3.7 KB
Line 
1#include "mltaln.h"
2#include "mtxutl.h"
3
4/*
5  from "C gengo niyoru saishin algorithm jiten" ISBN4-87408-414-1 Haruhiko Okumura
6*/
7static void make_sintbl(int n, float sintbl[])
8{
9        int i, n2, n4, n8;
10        double c, s, dc, ds, t;
11
12        n2 = n / 2;  n4 = n / 4;  n8 = n / 8;
13        t = sin(PI / n);
14        dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
15        t = 2 * dc;  c = sintbl[n4] = 1;  s = sintbl[0] = 0;
16        for (i = 1; i < n8; i++) {
17                c -= dc;  dc += t * c;
18                s += ds;  ds -= t * s;
19                sintbl[i] = s;  sintbl[n4 - i] = c;
20        }
21        if (n8 != 0) sintbl[n8] = sqrt(0.5);
22        for (i = 0; i < n4; i++)
23                sintbl[n2 - i] = sintbl[i];
24        for (i = 0; i < n2 + n4; i++)
25                sintbl[i + n2] = - sintbl[i];
26}
27/*
28  {\tt fft()}.
29*/
30static void make_bitrev(int n, int bitrev[])
31{
32        int i, j, k, n2;
33
34        n2 = n / 2;  i = j = 0;
35        for ( ; ; ) {
36                bitrev[i] = j;
37                if (++i >= n) break;
38                k = n2;
39                while (k <= j) {  j -= k;  k /= 2;  }
40                j += k;
41        }
42}
43/*
44*/
45int fft(int n, Fukusosuu *x, int freeflag)
46{
47        static TLS int    last_n = 0;    /*  {\tt n} */
48        static TLS int   *bitrev = NULL; /*  */
49        static TLS float *sintbl = NULL; /*  */
50        int i, j, k, ik, h, d, k2, n4, inverse;
51        float t, s, c, dR, dI;
52
53                if (freeflag)
54                {
55                        if (bitrev) free(bitrev); bitrev = NULL;
56                        if (sintbl) free(sintbl); sintbl = NULL;
57                        return( 0 );
58                }
59
60        /*  */
61        if (n < 0) {
62                n = -n;  inverse = 1;  /*  */
63        } else inverse = 0;
64        n4 = n / 4;
65        if (n != last_n || n == 0) {
66                last_n = n;
67#if 0
68                if (sintbl != NULL) {
69                                        free(sintbl);
70                                        sintbl = NULL;
71                                }
72                if (bitrev != NULL) {
73                                        free(bitrev);
74                                        bitrev = NULL;
75                                }
76                if (n == 0) return 0;  /*  */
77                sintbl = (float *)malloc((n + n4) * sizeof(float));
78                bitrev = (int *)malloc(n * sizeof(int));
79#else /* by T. Nishiyama */
80                                sintbl = realloc(sintbl, (n + n4) * sizeof(float));
81                                bitrev = realloc(bitrev, n * sizeof(int));
82#endif
83                if (sintbl == NULL || bitrev == NULL) {
84                        fprintf(stderr, "\n");  return 1;
85                }
86                make_sintbl(n, sintbl);
87                make_bitrev(n, bitrev);
88        }
89        for (i = 0; i < n; i++) {    /*  */
90                j = bitrev[i];
91                if (i < j) {
92                        t = x[i].R;  x[i].R = x[j].R;  x[j].R = t;
93                        t = x[i].I;  x[i].I = x[j].I;  x[j].I = t;
94                }
95        }
96        for (k = 1; k < n; k = k2) {    /*  */
97#if 0
98                                fprintf( stderr, "%d / %d\n", k, n );
99#endif
100                h = 0;  k2 = k + k;  d = n / k2;
101                for (j = 0; j < k; j++) {
102#if 0
103                                        if( j % 1 == 0 )
104                                                fprintf( stderr, "%d / %d\r", j, k );
105#endif
106                        c = sintbl[h + n4];
107                        if (inverse) s = - sintbl[h];
108                        else         s =   sintbl[h];
109                        for (i = j; i < n; i += k2) {
110#if 0
111                                                                if( k>=4194000 ) fprintf( stderr, "in loop %d -  %d < %d, k2=%d\r", j, i, n, k2 );
112#endif
113                                ik = i + k;
114                                dR = s * x[ik].I + c * x[ik].R;
115                                dI = c * x[ik].I - s * x[ik].R;
116                                x[ik].R = x[i].R - dR;  x[i].R += dR;
117                                x[ik].I = x[i].I - dI;  x[i].I += dI;
118                        }
119                        h += d;
120                }
121        }
122        if (! inverse)    /* n */
123                for (i = 0; i < n; i++) {  x[i].R /= n;  x[i].I /= n;  }
124        return 0;  /*  */
125}
Note: See TracBrowser for help on using the repository browser.