Subversion Repositories shark

Rev

Rev 2 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 1
/*
2
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
3
 *
4
 * This program is free software; you can redistribute it and/or modify
5
 * it under the terms of the GNU General Public License as published by
6
 * the Free Software Foundation; either version 2 of the License, or
7
 * (at your option) any later version.
8
 *
9
 * This program is distributed in the hope that it will be useful,
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12
 * GNU General Public License for more details.
13
 *
14
 * You should have received a copy of the GNU General Public License
15
 * along with this program; if not, write to the Free Software
16
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
17
 *
18
 */
19
 
20
/*
21
 * twiddle.c -- compute twiddle factors
22
 * These are the twiddle factors for *direct* fft.  Flip sign to get
23
 * the inverse
24
 */
25
 
26
/* $Id: twiddle.c,v 1.1.1.1 2002-03-29 14:12:57 pj Exp $ */
27
#ifdef FFTW_USING_CILK
28
#include <ports/cilk.h>
29
#include <ports/cilk-compat.h>
30
#endif
31
 
32
#include <ports/fftw-int.h>
33
#include <math.h>
34
#include <stdlib.h>
35
 
36
#ifndef TRUE
37
#define TRUE (1 == 1)
38
#endif
39
 
40
#ifndef FALSE
41
#define FALSE (1 == 0)
42
#endif
43
 
44
static fftw_complex *fftw_compute_rader_twiddle(int n, int r, int g)
45
{
46
     FFTW_TRIG_REAL twoPiOverN;
47
     int m = n / r;
48
     int i, j, gpower;
49
     fftw_complex *W;
50
 
51
     twoPiOverN = FFTW_K2PI / (FFTW_TRIG_REAL) n;
52
     W = (fftw_complex *) fftw_malloc((r - 1) * m * sizeof(fftw_complex));
53
     for (i = 0; i < m; ++i)
54
          for (gpower = 1, j = 0; j < r - 1; ++j, gpower = (gpower * g) % r) {
55
               int k = i * (r - 1) + j;
56
               FFTW_TRIG_REAL
57
                   ij = (FFTW_TRIG_REAL) (i * gpower);
58
               c_re(W[k]) = FFTW_TRIG_COS(twoPiOverN * ij);
59
               c_im(W[k]) = FFTW_FORWARD * FFTW_TRIG_SIN(twoPiOverN * ij);
60
          }
61
 
62
     return W;
63
}
64
 
65
/*
66
 * compute the W coefficients (that is, powers of the root of 1)
67
 * and store them into an array.
68
 */
69
static fftw_complex *fftw_compute_twiddle(int n, const fftw_codelet_desc *d)
70
{
71
     FFTW_TRIG_REAL twoPiOverN;
72
     int i, j;
73
     fftw_complex *W;
74
 
75
     twoPiOverN = FFTW_K2PI / (FFTW_TRIG_REAL) n;
76
 
77
     if (!d) {
78
          /* generic codelet, needs all twiddles in order */
79
          W = (fftw_complex *) fftw_malloc(n * sizeof(fftw_complex));
80
          for (i = 0; i < n; ++i) {
81
               c_re(W[i]) = FFTW_TRIG_COS(twoPiOverN * (FFTW_TRIG_REAL) i);
82
               c_im(W[i]) = FFTW_FORWARD * FFTW_TRIG_SIN(twoPiOverN * (FFTW_TRIG_REAL) i);
83
          }
84
     } else if (d->type == FFTW_RADER)
85
          W = fftw_compute_rader_twiddle(n, d->size, d->signature);
86
     else {
87
          int r = d->size;
88
          int m = n / r, m_alloc;
89
          int r1 = d->ntwiddle;
90
          int istart;
91
 
92
          if (d->type == FFTW_TWIDDLE) {
93
               istart = 0;
94
               m_alloc = m;
95
          } else if (d->type == FFTW_HC2HC) {
96
               /*
97
                * This is tricky, do not change lightly.
98
                */
99
               m = (m + 1) / 2;
100
               m_alloc = m - 1;
101
               istart = 1;
102
          } else {
103
               fftw_die("compute_twiddle: invalid argument\n");
104
               /* paranoia for gcc */
105
               m_alloc = 0;
106
               istart = 0;
107
          }
108
 
109
          W = (fftw_complex *) fftw_malloc(r1 * m_alloc * sizeof(fftw_complex));
110
          for (i = istart; i < m; ++i)
111
               for (j = 0; j < r1; ++j) {
112
                    int k = (i - istart) * r1 + j;
113
                    FFTW_TRIG_REAL
114
                        ij = (FFTW_TRIG_REAL) (i * d->twiddle_order[j]);
115
                    c_re(W[k]) = FFTW_TRIG_COS(twoPiOverN * ij);
116
                    c_im(W[k]) = FFTW_FORWARD * FFTW_TRIG_SIN(twoPiOverN * ij);
117
               }
118
     }
119
 
120
     return W;
121
}
122
 
123
/*
124
 * these routines implement a simple reference-count-based
125
 * management of twiddle structures
126
 */
127
static fftw_twiddle *twlist = (fftw_twiddle *) 0;
128
int fftw_twiddle_size = 0;      /* total allocated size, for debugging */
129
 
130
/* true if the two codelets can share the same twiddle factors */
131
static int compatible(const fftw_codelet_desc *d1, const fftw_codelet_desc *d2)
132
{
133
     int i;
134
 
135
     /* true if they are the same codelet */
136
     if (d1 == d2)
137
          return TRUE;
138
 
139
     /* false if one is null and the other is not */
140
     if (!d1 || !d2)
141
          return FALSE;
142
 
143
     /* false if size is different */
144
     if (d1->size != d2->size)
145
          return FALSE;
146
 
147
     /* false if different types (FFTW_TWIDDLE/FFTW_HC2HC/FFTW_RADER) */
148
     if (d1->type != d2->type)
149
          return FALSE;
150
 
151
     /* false if they need different # of twiddles */
152
     if (d1->ntwiddle != d2->ntwiddle)
153
          return FALSE;
154
 
155
     /* false if the twiddle orders are different */
156
     for (i = 0; i < d1->ntwiddle; ++i)
157
          if (d1->twiddle_order[i] != d2->twiddle_order[i])
158
               return FALSE;
159
 
160
     return TRUE;
161
}
162
 
163
fftw_twiddle *fftw_create_twiddle(int n, const fftw_codelet_desc *d)
164
{
165
     fftw_twiddle *tw;
166
 
167
     /* lookup this n in the twiddle list */
168
     for (tw = twlist; tw; tw = tw->next)
169
          if (n == tw->n && compatible(d, tw->cdesc)) {
170
               ++tw->refcnt;
171
               return tw;
172
          }
173
     /* not found --- allocate a new struct twiddle */
174
     tw = (fftw_twiddle *) fftw_malloc(sizeof(fftw_twiddle));
175
     fftw_twiddle_size += n;
176
 
177
     tw->n = n;
178
     tw->cdesc = d;
179
     tw->twarray = fftw_compute_twiddle(n, d);
180
     tw->refcnt = 1;
181
 
182
     /* enqueue the new struct */
183
     tw->next = twlist;
184
     twlist = tw;
185
 
186
     return tw;
187
}
188
 
189
void fftw_destroy_twiddle(fftw_twiddle * tw)
190
{
191
     fftw_twiddle **p;
192
     --tw->refcnt;
193
 
194
     if (tw->refcnt == 0) {
195
          /* remove from the list of known twiddle factors */
196
          for (p = &twlist; p; p = &((*p)->next))
197
               if (*p == tw) {
198
                    *p = tw->next;
199
                    fftw_twiddle_size -= tw->n;
200
                    fftw_free(tw->twarray);
201
                    fftw_free(tw);
202
                    return;
203
               }
204
          fftw_die("BUG in fftw_destroy_twiddle\n");
205
     }
206
}