/*
TiMidity++ -- MIDI to WAVE converter and player
Copyright (C) 1999-2002 Masanao Izumo <mo@goice.co.jp>
Copyright (C) 1995 Tuukka Toivonen <tt@cgs.fi>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif /* HAVE_CONFIG_H */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "timidity.h"
#include "common.h"
#include "fft.h"
#ifdef DEBUG
#ifdef sun
int fprintf(FILE *, const char *, ...);
#endif
#endif
#define swap(x, j, k) \
{ \
double tmp; \
tmp = x[j]; \
x[j] = x[k]; \
x[k] = tmp; \
}
#define swap2(x, j, k) \
{ \
double tmp; \
tmp = x[j]; \
x[j] = -x[k]; \
x[k] = tmp; \
}
static void make_table(double *trig, int *b, int n)
{
unsigned long i, j, k, bitmask;
/* check n */
for(i = n; !(i & 1); i >>= 1)
;
if(i != 1)
{
fprintf(stderr, "Invalid fft data size: %d\n", n);
exit(1);
}
/* make bitrev table */
for(i = 0; i < n; i++)
b[i] = 0;
bitmask = n / 2;
for(i = 1; i < n; i <<= 1, bitmask >>= 1)
for(j = 0; j < n; j += i * 2)
for(k = j + i; k < j + i * 2; k++)
b[k] = (int)((unsigned long)b[k] | bitmask);
/* make trig table */
for(i = 0; i < n; i++)
{
j = i * 2;
trig[j ] = cos((2 * M_PI) * i / (double)n);
trig[j+1] = sin((2 * M_PI) * i / (double)n);
}
for(i = 0; i < n; i++)
if(i < b[i])
{
j = b[i] * 2;
swap(trig, i*2, j);
swap(trig, i*2+1, j+1);
}
}
void realfft(double *x, int n_arg)
{
int n = n_arg;
static double *trig_table = NULL;
static int *bitrev_table;
int n1;
if(n == 0)
{
if(trig_table == NULL)
return;
free(trig_table);
free(bitrev_table);
trig_table = NULL;
return;
}
if(trig_table == NULL)
{
trig_table = (double *)safe_malloc((n * 2) * sizeof(double));
bitrev_table = (int *)safe_malloc(n * sizeof(int));
if(!(trig_table && bitrev_table))
{
fprintf(stderr, "fft: Can't allocate memroy.\n");
exit(1);
}
make_table(trig_table, bitrev_table, n);
if(x == NULL) /* initialize only */
return;
}
/* end of initialize tables */
/* first step: butterfly w^0 */
/* n1 = n/2, n/4, n/8, ... 1 */
for(n1 = n >> 1; n1 > 0; n1 >>= 1)
{
int i;
for(i = 0; i < n1; i++)
{
double x1, x2;
x[i] = (x1 = x[i]) + (x2 = x[i+n1]);
x[i+n1] = x1 - x2;
}
}
/* main loop */
/* n1 = n/8, n/16, n/32, ... 1 */
for(n1 = n >> 3; n1 > 0; n1 >>= 1)
{
int i;
int wi0;
wi0 = 8;
for(i = n1 << 2; i < n; i <<= 1, wi0 <<= 1)
{
double *imp; /* pointer to start of imaginary part */
double *rep; /* pointer to start of real part */
double *imp2; /* pointer to start of imaginary part */
double *rep2; /* pointer to start of real part */
int wi, j0;
int in;
in = i >> 1;
wi = wi0;
rep = x + i;
imp = rep + in;
rep2 = rep + n1;
imp2 = imp + n1;
for(j0 = 0; j0 < in; j0 += (n1 << 1), wi += 4)
{
int j, jn;
double c, s;
c = trig_table[wi];
s = trig_table[wi+1];
jn = j0 + n1;
/* butterfly */
for(j = j0; j < jn; j++)
{
double xi1, xr1, xi2, xr2, ti, tr;
rep[j] = (xr1 = rep[j]) + (tr = c * (xr2 = rep2[j])
- s * (xi2 = imp2[j]));
imp[j] = (xi1 = imp[j]) + (ti = s * xr2 + c * xi2);
rep2[j] = xr1 - tr;
imp2[j] = xi1 - ti;
}
}
}
}
/* move data */
{
int ii, i, j, i2, i4;
for(i = 4; i < n; i = ii)
{
i2 = i >> 1;
i4 = i2 >> 1;
ii = i << 1;
for(j = 0; j < i4; j++)
swap(x, j + i + i2, ii - j - 1);
for(j = 1; j < i2; j += 2)
swap2(x, j + i, ii - j - 1);
}
for(i = 0; i < n; i++)
if(i < bitrev_table[i])
swap(x, i, bitrev_table[i]);
for(i = n/2+1; i < n; i++)
x[i] = -x[i];
}
}