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 bar()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 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 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