1 /*
2 TiMidity++ -- MIDI to WAVE converter and player
3 Copyright (C) 1999-2002 Masanao Izumo <mo@goice.co.jp>
4 Copyright (C) 1995 Tuukka Toivonen <tt@cgs.fi>
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program 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
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 */
20
21 #ifdef HAVE_CONFIG_H
22 #include "config.h"
23 #endif /* HAVE_CONFIG_H */
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <math.h>
27 #include "timidity.h"
28 #include "common.h"
29 #include "fft.h"
30
31 #ifdef DEBUG
32 #ifdef sun
33 int fprintf(FILE *, const char *, ...);
34 #endif
35 #endif
36
37
38 #define swap(x, j, k) \
39 { \
40 double tmp; \
41 tmp = x[j]; \
42 x[j] = x[k]; \
43 x[k] = tmp; \
44 }
45
46 #define swap2(x, j, k) \
47 { \
48 double tmp; \
49 tmp = x[j]; \
50 x[j] = -x[k]; \
51 x[k] = tmp; \
52 }
53
make_table(double * trig,int * b,int n)54 static void make_table(double *trig, int *b, int n)
55 {
56 unsigned long i, j, k, bitmask;
57
58 /* check n */
59 for(i = n; !(i & 1); i >>= 1)
60 ;
61 if(i != 1)
62 {
63 fprintf(stderr, "Invalid fft data size: %d\n", n);
64 exit(1);
65 }
66
67 /* make bitrev table */
68 for(i = 0; i < n; i++)
69 b[i] = 0;
70 bitmask = n / 2;
71 for(i = 1; i < n; i <<= 1, bitmask >>= 1)
72 for(j = 0; j < n; j += i * 2)
73 for(k = j + i; k < j + i * 2; k++)
74 b[k] = (int)((unsigned long)b[k] | bitmask);
75
76 /* make trig table */
77 for(i = 0; i < n; i++)
78 {
79 j = i * 2;
80 trig[j ] = cos((2 * M_PI) * i / (double)n);
81 trig[j+1] = sin((2 * M_PI) * i / (double)n);
82 }
83
84 for(i = 0; i < n; i++)
85 if(i < b[i])
86 {
87 j = b[i] * 2;
88 swap(trig, i*2, j);
89 swap(trig, i*2+1, j+1);
90 }
91
92 }
93
realfft(double * x,int n_arg)94 void realfft(double *x, int n_arg)
95 {
96 int n = n_arg;
97 static double *trig_table = NULL;
98 static int *bitrev_table;
99 int n1;
100
101 if(n == 0)
102 {
103 if(trig_table == NULL)
104 return;
105 free(trig_table);
106 free(bitrev_table);
107 trig_table = NULL;
108 return;
109 }
110
111 if(trig_table == NULL)
112 {
113 trig_table = (double *)safe_malloc((n * 2) * sizeof(double));
114 bitrev_table = (int *)safe_malloc(n * sizeof(int));
115
116 if(!(trig_table && bitrev_table))
117 {
118 fprintf(stderr, "fft: Can't allocate memroy.\n");
119 exit(1);
120 }
121
122 make_table(trig_table, bitrev_table, n);
123
124 if(x == NULL) /* initialize only */
125 return;
126 }
127 /* end of initialize tables */
128
129
130 /* first step: butterfly w^0 */
131 /* n1 = n/2, n/4, n/8, ... 1 */
132 for(n1 = n >> 1; n1 > 0; n1 >>= 1)
133 {
134 int i;
135
136 for(i = 0; i < n1; i++)
137 {
138 double x1, x2;
139
140 x[i] = (x1 = x[i]) + (x2 = x[i+n1]);
141 x[i+n1] = x1 - x2;
142 }
143 }
144
145 /* main loop */
146 /* n1 = n/8, n/16, n/32, ... 1 */
147 for(n1 = n >> 3; n1 > 0; n1 >>= 1)
148 {
149 int i;
150 int wi0;
151
152 wi0 = 8;
153 for(i = n1 << 2; i < n; i <<= 1, wi0 <<= 1)
154 {
155 double *imp; /* pointer to start of imaginary part */
156 double *rep; /* pointer to start of real part */
157 double *imp2; /* pointer to start of imaginary part */
158 double *rep2; /* pointer to start of real part */
159 int wi, j0;
160 int in;
161
162 in = i >> 1;
163 wi = wi0;
164 rep = x + i;
165 imp = rep + in;
166 rep2 = rep + n1;
167 imp2 = imp + n1;
168
169 for(j0 = 0; j0 < in; j0 += (n1 << 1), wi += 4)
170 {
171 int j, jn;
172 double c, s;
173
174 c = trig_table[wi];
175 s = trig_table[wi+1];
176 jn = j0 + n1;
177
178 /* butterfly */
179 for(j = j0; j < jn; j++)
180 {
181 double xi1, xr1, xi2, xr2, ti, tr;
182
183 rep[j] = (xr1 = rep[j]) + (tr = c * (xr2 = rep2[j])
184 - s * (xi2 = imp2[j]));
185 imp[j] = (xi1 = imp[j]) + (ti = s * xr2 + c * xi2);
186 rep2[j] = xr1 - tr;
187 imp2[j] = xi1 - ti;
188 }
189 }
190 }
191 }
192
193 /* move data */
194 {
195 int ii, i, j, i2, i4;
196
197 for(i = 4; i < n; i = ii)
198 {
199 i2 = i >> 1;
200 i4 = i2 >> 1;
201 ii = i << 1;
202 for(j = 0; j < i4; j++)
203 swap(x, j + i + i2, ii - j - 1);
204 for(j = 1; j < i2; j += 2)
205 swap2(x, j + i, ii - j - 1);
206 }
207
208 for(i = 0; i < n; i++)
209 if(i < bitrev_table[i])
210 swap(x, i, bitrev_table[i]);
211 for(i = n/2+1; i < n; i++)
212 x[i] = -x[i];
213 }
214 }
215