[PD-cvs] pd/src d_fft_fftsg.c, 1.1, 1.2 d_fft_fftw.c, 1.1, 1.2 d_fftsg_h.c, 1.1, 1.2

Miller Puckette millerpuckette at users.sourceforge.net
Sun Jul 29 22:17:09 CEST 2007


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 at 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");
}






More information about the Pd-cvs mailing list