d_fft_fftw.c 4.02 KB
Newer Older
Miller Puckette's avatar
Miller Puckette committed
1
2
3
4
5
6
/* 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 ---------- */

Ivica Bukvic's avatar
Ivica Bukvic committed
7
/* changes and additions for FFTW3 by Thomas Grill                      */
Miller Puckette's avatar
Miller Puckette committed
8

Ivica Bukvic's avatar
Ivica Bukvic committed
9
10
#include "m_pd.h"
#include <fftw3.h>
Miller Puckette's avatar
Miller Puckette committed
11
12
13

int ilog2(int n);

Ivica Bukvic's avatar
Ivica Bukvic committed
14
#define MINFFT 0
Miller Puckette's avatar
Miller Puckette committed
15
16
17
#define MAXFFT 30

/* from the FFTW website:
Ivica Bukvic's avatar
Ivica Bukvic committed
18
19
20
21
 #include <fftw3.h>
     ...
     {
         fftw_complex *in, *out;
Miller Puckette's avatar
Miller Puckette committed
22
23
     fftw_plan p;
     ...
Ivica Bukvic's avatar
Ivica Bukvic committed
24
25
26
         in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
         out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
         p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
Miller Puckette's avatar
Miller Puckette committed
27
     ...
Ivica Bukvic's avatar
Ivica Bukvic committed
28
         fftw_execute(p); 
Miller Puckette's avatar
Miller Puckette committed
29
30
     ...
     fftw_destroy_plan(p);  
Ivica Bukvic's avatar
Ivica Bukvic committed
31
32
         fftw_free(in); fftw_free(out);
     }
Miller Puckette's avatar
Miller Puckette committed
33
34
35
36
37
38
39
40

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

*/

Ivica Bukvic's avatar
Ivica Bukvic committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
/* complex stuff */

typedef struct {
    fftwf_plan plan;
    fftwf_complex *in,*out;
} cfftw_info;

static cfftw_info cfftw_fwd[MAXFFT+1 - MINFFT],cfftw_bwd[MAXFFT+1 - MINFFT];

static cfftw_info *cfftw_getplan(int n,int fwd)
{
    cfftw_info *info;
    int logn = ilog2(n);
    if (logn < MINFFT || logn > MAXFFT)
        return (0);
    info = (fwd?cfftw_fwd:cfftw_bwd)+(logn-MINFFT);
    if (!info->plan) 
    {
        info->in = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * n);
        info->out = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * n);
        info->plan = fftwf_plan_dft_1d(n, info->in, info->out, fwd?FFTW_FORWARD:FFTW_BACKWARD, FFTW_MEASURE);        
    }
    return info;
}

Miller Puckette's avatar
Miller Puckette committed
66

Ivica Bukvic's avatar
Ivica Bukvic committed
67
68
69
70
/* real stuff */

typedef struct {
    fftwf_plan plan;
71
    t_float *in,*out;
Ivica Bukvic's avatar
Ivica Bukvic committed
72
73
74
75
76
} rfftw_info;

static rfftw_info rfftw_fwd[MAXFFT+1 - MINFFT],rfftw_bwd[MAXFFT+1 - MINFFT];

static rfftw_info *rfftw_getplan(int n,int fwd)
Miller Puckette's avatar
Miller Puckette committed
77
{
Ivica Bukvic's avatar
Ivica Bukvic committed
78
79
    rfftw_info *info;
    int logn = ilog2(n);
Miller Puckette's avatar
Miller Puckette committed
80
81
    if (logn < MINFFT || logn > MAXFFT)
        return (0);
Ivica Bukvic's avatar
Ivica Bukvic committed
82
83
84
    info = (fwd?rfftw_fwd:rfftw_bwd)+(logn-MINFFT);
    if (!info->plan) 
    {
85
86
        info->in = (t_float*) fftwf_malloc(sizeof(t_float) * n);
        info->out = (t_float*) fftwf_malloc(sizeof(t_float) * n);
Ivica Bukvic's avatar
Ivica Bukvic committed
87
88
89
        info->plan = fftwf_plan_r2r_1d(n, info->in, info->out, fwd?FFTW_R2HC:FFTW_HC2R, FFTW_MEASURE);
    }
    return info;
Miller Puckette's avatar
Miller Puckette committed
90
91
}

Ivica Bukvic's avatar
Ivica Bukvic committed
92
93


94
EXTERN void mayer_fht(t_float *fz, int n)
Miller Puckette's avatar
Miller Puckette committed
95
96
97
98
{
    post("FHT: not yet implemented");
}

99
static void mayer_do_cfft(int n, t_float *fz1, t_float *fz2, int fwd)
Miller Puckette's avatar
Miller Puckette committed
100
101
{
    int i;
102
    t_float *fz;
Ivica Bukvic's avatar
Ivica Bukvic committed
103
    cfftw_info *p = cfftw_getplan(n, fwd);
Miller Puckette's avatar
Miller Puckette committed
104
105
    if (!p)
        return;
Ivica Bukvic's avatar
Ivica Bukvic committed
106
        
107
    for (i = 0, fz = (t_float *)p->in; i < n; i++)
Ivica Bukvic's avatar
Ivica Bukvic committed
108
109
110
111
        fz[i*2] = fz1[i], fz[i*2+1] = fz2[i];
        
    fftwf_execute(p->plan);
    
112
    for (i = 0, fz = (t_float *)p->out; i < n; i++)
Ivica Bukvic's avatar
Ivica Bukvic committed
113
        fz1[i] = fz[i*2], fz2[i] = fz[i*2+1];
Miller Puckette's avatar
Miller Puckette committed
114
115
}

116
EXTERN void mayer_fft(int n, t_float *fz1, t_float *fz2)
Miller Puckette's avatar
Miller Puckette committed
117
{
Ivica Bukvic's avatar
Ivica Bukvic committed
118
    mayer_do_cfft(n, fz1, fz2, 1);
Miller Puckette's avatar
Miller Puckette committed
119
120
}

121
EXTERN void mayer_ifft(int n, t_float *fz1, t_float *fz2)
Miller Puckette's avatar
Miller Puckette committed
122
{
Ivica Bukvic's avatar
Ivica Bukvic committed
123
    mayer_do_cfft(n, fz1, fz2, 0);
Miller Puckette's avatar
Miller Puckette committed
124
125
}

Ivica Bukvic's avatar
Ivica Bukvic committed
126
127
128
129
130
131
/* 
    in the following the sign flips are done to 
    be compatible with the mayer_fft implementation,
    but it's probably the mayer_fft that should be corrected...
*/ 

132
EXTERN void mayer_realfft(int n, t_float *fz)
Miller Puckette's avatar
Miller Puckette committed
133
{
Ivica Bukvic's avatar
Ivica Bukvic committed
134
135
136
137
138
139
140
141
142
143
144
145
    int i;
    rfftw_info *p = rfftw_getplan(n, 1);
    if (!p)
        return;
        
    for (i = 0; i < n; i++)
        p->in[i] = fz[i];
    fftwf_execute(p->plan);
    for (i = 0; i < n/2+1; i++)
        fz[i] = p->out[i];
    for (; i < n; i++)
        fz[i] = -p->out[i];
Miller Puckette's avatar
Miller Puckette committed
146
147
}

148
EXTERN void mayer_realifft(int n, t_float *fz)
Miller Puckette's avatar
Miller Puckette committed
149
{
Ivica Bukvic's avatar
Ivica Bukvic committed
150
151
152
153
154
155
156
157
158
159
160
161
    int i;
    rfftw_info *p = rfftw_getplan(n, 0);
    if (!p)
        return;
        
    for (i = 0; i < n/2+1; i++)
        p->in[i] = fz[i];
    for (; i < n; i++)
        p->in[i] = -fz[i];    
    fftwf_execute(p->plan);
    for (i = 0; i < n; i++)
        fz[i] = p->out[i];
Miller Puckette's avatar
Miller Puckette committed
162
163
}