Update of /cvsroot/pure-data/pd/src In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv29112/pd/src
Added Files: d_fft_fftsg.c d_fft_fftw.c d_fftsg_h.c Log Message: add files that didn't seem to be uploaded
--- NEW FILE: d_fftsg_h.c --- /* Fast Fourier/Cosine/Sine Transform dimension :one data length :power of 2 decimation :frequency radix :split-radix data :inplace table :not use functions cdft: Complex Discrete Fourier Transform rdft: Real Discrete Fourier Transform ddct: Discrete Cosine Transform ddst: Discrete Sine Transform dfct: Cosine Transform of RDFT (Real Symmetric DFT) dfst: Sine Transform of RDFT (Real Anti-symmetric DFT) function prototypes void cdft(int, int, double *); void rdft(int, int, double *); void ddct(int, int, double *); [...3408 lines suppressed...] }
void dstsub4(int n, double *a) { int m; double wki, wdr, wdi, xr;
wki = WR5000; m = n >> 1; if (m == 2) { wdr = wki * WI2500; wdi = wki * WR2500; xr = wdi * a[3] - wdr * a[1]; a[3] = wdr * a[3] + wdi * a[1]; a[1] = xr; } a[m] *= wki; }
--- NEW FILE: d_fft_fftsg.c --- /****************** begin Pd-specific prologue ***********************/
/* This is the file, "fftsg.c" from the OOURA FFT package. The copyright notice from Ooura's README is:
Copyright: Copyright(C) 1996-2001 Takuya OOURA email: ooura@mmm.t.u-tokyo.ac.jp download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html You may use, copy, modify this code for any purpose and without fee. You may distribute this ORIGINAL package.
After the following prologue, the code is exactly Ooura's original, which I believe the above notice permits me to redistribute. See Ooura's website for another, more permissive-sounding copyright notice. -MSP */
/* ---------- Pd interface to OOURA FFT; imitate Mayer API ---------- */ [...3423 lines suppressed...] void dstsub(int n, double *a, int nc, double *c) { int j, k, kk, ks, m; double wkr, wki, xr;
m = n >> 1; ks = nc / n; kk = 0; for (j = 1; j < m; j++) { k = n - j; kk += ks; wkr = c[kk] - c[nc - kk]; wki = c[kk] + c[nc - kk]; xr = wki * a[k] - wkr * a[j]; a[k] = wkr * a[k] + wki * a[j]; a[j] = xr; } a[m] *= c[0]; }
--- NEW FILE: d_fft_fftw.c --- /* Copyright (c) 1997- Miller Puckette and others. * For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution. */
/* --------- Pd interface to FFTW library; imitate Mayer API ---------- */
#include "m_pd.h" #ifdef MSW #include <malloc.h> #else #include <alloca.h> #endif
#error oops -- I'm talking to the old fftw. Apparently it's still changing.
#include <fftw.h>
int ilog2(int n);
#define MINFFT 5 #define MAXFFT 30
/* from the FFTW website: fftw_complex in[N], out[N]; fftw_plan p; ... p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE); ... fftw_one(p, in, out); ... fftw_destroy_plan(p);
FFTW_FORWARD or FFTW_BACKWARD, and indicates the direction of the transform you are interested in. Alternatively, you can use the sign of the exponent in the transform, -1 or +1, which corresponds to FFTW_FORWARD or FFTW_BACKWARD respectively. The flags argument is either FFTW_MEASURE
*/
static fftw_plan fftw_pvec[2 * (MAXFFT+1 - MINFFT)];
static fftw_plan fftw_getplan(int n, int dir) { logn = ilog2(n); if (logn < MINFFT || logn > MAXFFT) return (0); int indx = 2*(logn-MINFFT) + inverse); if (!fftw_pvec[indx] fftw_pvec[indx] = fftw_create_plan(N, dir, FFTW_MEASURE); return (fftw_pvec[indx]); }
EXTERN void mayer_fht(float *fz, int n) { post("FHT: not yet implemented"); }
static void mayer_dofft(int n, float *fz1, float *fz2, int dir) { float *inbuf, *outbuf, *fp1, *fp2, *fp3; int i; fftw_plan p = fftw_getplan(n, dir); inbuf = alloca(n * (4 * sizeof(float))); outbuf = inbuf + 2*n; if (!p) return; for (i = 0, fp1 = fz1, fp2 = fz2, fp3 = inbuf; i < n; i++) fp3[0] = *fp1++, fp3[1] = *fp2++, fp3 += 2; fftw_one(p, inbuf, outbuf); for (i = 0, fp1 = fz1, fp2 = fz2, fp3 = outbuf; i < n; i++) *fp1++ = fp3[0], *fp2++ = fp3[1], fp3 += 2; }
EXTERN void mayer_fft(int n, float *fz1, float *fz2) { mayer_dofft(n, fz1, fz2, FFTW_FORWARD); }
EXTERN void mayer_ifft(int n, float *fz1, float *fz2) { mayer_dofft(n, fz1, fz2, FFTW_BACKWARD); }
EXTERN void mayer_realfft(int n, float *fz) { post("rfft: not yet implemented"); }
EXTERN void mayer_realifft(int n, float *fz) { post("rifft: not yet implemented"); }