← Back to branch summary

~ubuntu-branches/ubuntu/maverick/blender/maverick

« back to all changes in this revision

Viewing changes to extern/fftw/rdft/rdft2-radix2.c

  • Committer: Bazaar Package Importer
  • Author(s): Khashayar Naderehvandi, Khashayar Naderehvandi, Alessio Treglia
  • Date: 2009-01-22 16:53:59 UTC
  • mfrom: (14.1.1 experimental)
  • Revision ID: james.westby@ubuntu.com-20090122165359-v0996tn7fbit64ni
Tags: 2.48a+dfsg-1ubuntu1
https://launchpad.net/bugs/320045
[ Khashayar Naderehvandi ]
* Merge from debian experimental (LP: #320045), Ubuntu remaining changes:
  - Add patch correcting header file locations.
  - Add libvorbis-dev and libgsm1-dev to Build-Depends.
  - Use avcodec_decode_audio2() in source/blender/src/hddaudio.c

[ Alessio Treglia ]
* Add missing previous changelog entries.

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/fftw/rdft/rdft2-radix2.c
 
1
/*
 
2
 * Copyright (c) 2003, 2006 Matteo Frigo
 
3
 * Copyright (c) 2003, 2006 Massachusetts Institute of Technology
 
4
 *
 
5
 * This program is free software; you can redistribute it and/or modify
 
6
 * it under the terms of the GNU General Public License as published by
 
7
 * the Free Software Foundation; either version 2 of the License, or
 
8
 * (at your option) any later version.
 
9
 *
 
10
 * This program is distributed in the hope that it will be useful,
 
11
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
 * GNU General Public License for more details.
 
14
 *
 
15
 * You should have received a copy of the GNU General Public License
 
16
 * along with this program; if not, write to the Free Software
 
17
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
18
 *
 
19
 */
 
20
 
 
21
/* $Id: rdft2-radix2.c,v 1.32 2006-01-27 02:10:50 athena Exp $ */
 
22
 
 
23
/*
 
24
  Compute RDFT2 of even size via either a DFT or a vector RDFT of
 
25
  size n/2.
 
26
 
 
27
  This file is meant as a temporary hack until we do the right thing.
 
28
 
 
29
  The right thing is: 1) get rid of reduction to DFT, and 2) implement
 
30
  arbitrary even-radix reduction to RDFT.  We currently reduce to DFT
 
31
  so as to exploit the SIMD code.  We currently do only radix-2 in
 
32
  order to avoid generating yet another set of codelets.
 
33
*/
 
34
 
 
35
#include "rdft.h"
 
36
#include "dft.h"
 
37
 
 
38
typedef struct {
 
39
     int (*applicable) (const problem *p_, const planner *plnr);
 
40
     void (*apply) (const plan *ego_, R *r, R *rio, R *iio);
 
41
     problem *(*mkcld) (const problem_rdft2 *p);
 
42
     opcnt ops;
 
43
     const char *nam;
 
44
} madt;
 
45
 
 
46
typedef struct {
 
47
     solver super;
 
48
     const madt *adt;
 
49
} S;
 
50
 
 
51
typedef struct {
 
52
     plan_dft super;
 
53
     plan *cld;
 
54
     twid *td;
 
55
     INT is, os, ivs, ovs;
 
56
     INT n, vl;
 
57
     const S *slv;
 
58
} P;
 
59
 
 
60
/* common applicability function of forward problems */
 
61
static int applicable_f(const problem *p_, const planner *plnr)
 
62
{
 
63
     const problem_rdft2 *p = (const problem_rdft2 *) p_;
 
64
     UNUSED(plnr);
 
65
     return (1
 
66
             && p->kind == R2HC
 
67
             && p->vecsz->rnk <= 1
 
68
             && p->sz->rnk == 1
 
69
             && (p->sz->dims[0].n % 2) == 0
 
70
          );
 
71
}
 
72
 
 
73
static int applicable_f_dft(const problem *p_, const planner *plnr)
 
74
{
 
75
     UNUSED(plnr);
 
76
     if (applicable_f(p_, plnr)) {
 
77
          const problem_rdft2 *p = (const problem_rdft2 *) p_;
 
78
          return(p->r != p->rio
 
79
                 || (p->iio == p->rio + p->sz->dims[0].is
 
80
                     && p->sz->dims[0].os == 2 * p->sz->dims[0].is));
 
81
     }
 
82
     return 0;
 
83
}
 
84
 
 
85
/* common applicability function of backward problems */
 
86
static int applicable_b(const problem *p_, const planner *plnr)
 
87
{
 
88
     const problem_rdft2 *p = (const problem_rdft2 *) p_;
 
89
     return (1
 
90
             && p->kind == HC2R
 
91
             && (p->r == p->rio || !NO_DESTROY_INPUTP(plnr))
 
92
             && p->vecsz->rnk <= 1
 
93
             && p->sz->rnk == 1
 
94
             && (p->sz->dims[0].n % 2) == 0
 
95
          );
 
96
}
 
97
 
 
98
static int applicable_b_dft(const problem *p_, const planner *plnr)
 
99
{
 
100
     UNUSED(plnr);
 
101
     if (applicable_b(p_, plnr)) {
 
102
          const problem_rdft2 *p = (const problem_rdft2 *) p_;
 
103
          return(p->r != p->rio
 
104
                 || (p->iio == p->rio + p->sz->dims[0].os
 
105
                     && p->sz->dims[0].is == 2 * p->sz->dims[0].os));
 
106
     }
 
107
     return 0;
 
108
}
 
109
 
 
110
/*
 
111
 * forward rdft2 via dft
 
112
 */
 
113
static void k_f_dft(R *rio, R *iio, const R *W, INT n, INT dist)
 
114
{
 
115
     INT i;
 
116
     R *pp = rio, *pm = rio + n * dist;
 
117
     INT im = iio - rio;
 
118
 
 
119
     /* i = 0 and i = n */
 
120
     {
 
121
          E rop = pp[0], iop = pp[im];
 
122
          pp[0] = rop + iop;
 
123
          pm[0] = rop - iop;
 
124
          pp[im] = K(0.0);
 
125
          pm[im] = K(0.0);
 
126
          pp += dist; pm -= dist;
 
127
     }
 
128
 
 
129
     /* middle elements */
 
130
     for (W += 2, i = 2; i < n; i += 2, W += 2) {
 
131
          E rop = pp[0], iop = pp[im], rom = pm[0], iom = pm[im];
 
132
          E wr = W[0], wi = W[1];
 
133
          E re = rop + rom;
 
134
          E ie = iop - iom;
 
135
          E rd = rom - rop;
 
136
          E id = iop + iom;
 
137
          E tr = rd * wr - id * wi;
 
138
          E ti = id * wr + rd * wi;
 
139
          pp[0] = K(0.5) * (re + ti);
 
140
          pp[im] = K(0.5) * (ie + tr);
 
141
          pm[0] = K(0.5) * (re - ti);
 
142
          pm[im] = K(0.5) * (tr - ie);
 
143
          pp += dist; pm -= dist;
 
144
     }
 
145
 
 
146
     /* i = n/2 when n is even */
 
147
     if (!(n & 1)) pp[im] = -pp[im];
 
148
}
 
149
 
 
150
static void apply_f_dft(const plan *ego_, R *r, R *rio, R *iio)
 
151
{
 
152
     const P *ego = (const P *) ego_;
 
153
 
 
154
     {
 
155
          /* transform input as a vector of complex numbers */
 
156
          plan_dft *cld = (plan_dft *) ego->cld;
 
157
          cld->apply((plan *) cld, r, r + ego->is, rio, iio);
 
158
     }
 
159
 
 
160
     {
 
161
          INT i, vl = ego->vl, n2 = ego->n / 2;
 
162
          INT ovs = ego->ovs, os = ego->os;
 
163
          const R *W = ego->td->W;
 
164
          for (i = 0; i < vl; ++i, rio += ovs, iio += ovs)
 
165
               k_f_dft(rio, iio, W, n2, os);
 
166
     }
 
167
}
 
168
 
 
169
static problem *mkcld_f_dft(const problem_rdft2 *p)
 
170
{
 
171
     const iodim *d = p->sz->dims;
 
172
     return X(mkproblem_dft_d) (
 
173
          X(mktensor_1d)(d[0].n / 2, d[0].is * 2, d[0].os),
 
174
          X(tensor_copy)(p->vecsz),
 
175
          p->r, p->r + d[0].is, p->rio, p->iio);
 
176
}
 
177
 
 
178
static const madt adt_f_dft = {
 
179
     applicable_f_dft, apply_f_dft, mkcld_f_dft, {10, 8, 0, 0}, "r2hc2-dft"
 
180
};
 
181
 
 
182
/*
 
183
 * forward rdft2 via rdft
 
184
 */
 
185
static void k_f_rdft(R *rio, R *iio, const R *W, INT n, INT dist)
 
186
{
 
187
     INT i;
 
188
     R *pp = rio, *pm = rio + n * dist;
 
189
     INT im = iio - rio;
 
190
 
 
191
     /* i = 0 and i = n */
 
192
     {
 
193
          E rop = pp[0], iop = pp[im];
 
194
          pp[0] = rop + iop;
 
195
          pm[0] = rop - iop;
 
196
          pp[im] = K(0.0);
 
197
          pm[im] = K(0.0);
 
198
          pp += dist; pm -= dist;
 
199
     }
 
200
 
 
201
     /* middle elements */
 
202
     for (W += 2, i = 2; i < n; i += 2, W += 2) {
 
203
          E r0 = pp[0], r1 = pp[im], i0 = pm[0], i1 = pm[im];
 
204
          E wr = W[0], wi = W[1];
 
205
          E tr = r1 * wr + i1 * wi;
 
206
          E ti = i1 * wr - r1 * wi;
 
207
          pp[0] = r0 + tr;
 
208
          pp[im] = i0 + ti;
 
209
          pm[0] = r0 - tr;
 
210
          pm[im] = ti - i0;
 
211
          pp += dist; pm -= dist;
 
212
     }
 
213
 
 
214
     /* i = n/2 when n is even */
 
215
     if (!(n & 1)) pp[im] = -pp[im];
 
216
}
 
217
 
 
218
static void apply_f_rdft(const plan *ego_, R *r, R *rio, R *iio)
 
219
{
 
220
     const P *ego = (const P *) ego_;
 
221
 
 
222
     {
 
223
          plan_rdft *cld = (plan_rdft *) ego->cld;
 
224
          cld->apply((plan *) cld, r, rio);
 
225
     }
 
226
 
 
227
     {
 
228
          INT i, vl = ego->vl, n2 = ego->n / 2;
 
229
          INT ovs = ego->ovs, os = ego->os;
 
230
          const R *W = ego->td->W;
 
231
          for (i = 0; i < vl; ++i, rio += ovs, iio += ovs)
 
232
               k_f_rdft(rio, iio, W, n2, os);
 
233
     }
 
234
}
 
235
 
 
236
static problem *mkcld_f_rdft(const problem_rdft2 *p)
 
237
{
 
238
     const iodim *d = p->sz->dims;
 
239
 
 
240
     tensor *radix = X(mktensor_1d)(2, d[0].is, p->iio - p->rio);
 
241
     tensor *cld_vec = X(tensor_append)(radix, p->vecsz);
 
242
     X(tensor_destroy)(radix);
 
243
 
 
244
     return X(mkproblem_rdft_1_d) (
 
245
          X(mktensor_1d)(d[0].n / 2, 2 * d[0].is, d[0].os),
 
246
          cld_vec, p->r, p->rio, R2HC);
 
247
}
 
248
 
 
249
static const madt adt_f_rdft = {
 
250
     applicable_f, apply_f_rdft, mkcld_f_rdft, {6, 4, 0, 0}, "r2hc2-rdft"
 
251
};
 
252
 
 
253
 
 
254
/*
 
255
 * backward rdft2 via dft
 
256
 */
 
257
static void k_b_dft(R *rio, R *iio, const R *W, INT n, INT dist)
 
258
{
 
259
     INT i;
 
260
     R *pp = rio, *pm = rio + n * dist;
 
261
     INT im = iio - rio;
 
262
 
 
263
     /* i = 0 and i = n */
 
264
     {
 
265
          E rop = pp[0], iop = pm[0];
 
266
          pp[0] = rop + iop;
 
267
          pp[im] = rop - iop;
 
268
          pp += dist; pm -= dist;
 
269
     }
 
270
 
 
271
     /* middle elements */
 
272
     for (W += 2, i = 2; i < n; i += 2, W += 2) {
 
273
          E a = pp[0], b = pp[im], c = pm[0], d = pm[im];
 
274
          E wr = W[0], wi = W[1];
 
275
          E re = a + c, ti = a - c, ie = b - d, tr = b + d;
 
276
          E rd = tr * wr + ti * wi;
 
277
          E id = ti * wr - tr * wi;
 
278
          pp[0] = re - rd;
 
279
          pp[im] = ie + id;
 
280
          pm[0] = re + rd;
 
281
          pm[im] = id - ie;
 
282
          pp += dist; pm -= dist;
 
283
     }
 
284
 
 
285
     /* i = n/2 when n is even */
 
286
     if (!(n & 1)) { pp[0] *= K(2.0); pp[im] *= -K(2.0); }
 
287
}
 
288
 
 
289
static void apply_b_dft(const plan *ego_, R *r, R *rio, R *iio)
 
290
{
 
291
     const P *ego = (const P *) ego_;
 
292
     {
 
293
          INT i, vl = ego->vl, n2 = ego->n / 2;
 
294
          INT ivs = ego->ivs, is = ego->is;
 
295
          const R *W = ego->td->W;
 
296
          R *rio1 = rio, *iio1 = iio;
 
297
          for (i = 0; i < vl; ++i, rio1 += ivs, iio1 += ivs)
 
298
               k_b_dft(rio1, iio1, W, n2, is);
 
299
     }
 
300
 
 
301
     {
 
302
          plan_dft *cld = (plan_dft *) ego->cld;
 
303
          /* swap r/i because of backward transform */
 
304
          cld->apply((plan *) cld, iio, rio, r + ego->os, r);
 
305
     }
 
306
}
 
307
 
 
308
static problem *mkcld_b_dft(const problem_rdft2 *p)
 
309
{
 
310
     const iodim *d = p->sz->dims;
 
311
 
 
312
     return X(mkproblem_dft_d) (
 
313
          X(mktensor_1d)(d[0].n / 2, d[0].is, 2 * d[0].os),
 
314
          X(tensor_copy)(p->vecsz),
 
315
          p->iio, p->rio, p->r + d[0].os, p->r);
 
316
}
 
317
 
 
318
static const madt adt_b_dft = {
 
319
     applicable_b_dft, apply_b_dft, mkcld_b_dft, {10, 8, 0, 0}, "hc2r2-dft"
 
320
};
 
321
 
 
322
/*
 
323
 * backward rdft2 via backward rdft
 
324
 */
 
325
static void k_b_rdft(R *rio, R *iio, const R *W, INT n, INT dist)
 
326
{
 
327
     INT i;
 
328
     R *pp = rio, *pm = rio + n * dist;
 
329
     INT im = iio - rio;
 
330
 
 
331
     /* i = 0 and i = n */
 
332
     {
 
333
          E rop = pp[0], iop = pm[0];
 
334
          pp[0] = rop + iop;
 
335
          pp[im] = rop - iop;
 
336
          pp += dist; pm -= dist;
 
337
     }
 
338
 
 
339
     /* middle elements */
 
340
     for (W += 2, i = 2; i < n; i += 2, W += 2) {
 
341
          E a = pp[0], b = pp[im], c = pm[0], d = pm[im];
 
342
          E wr = W[0], wi = W[1];
 
343
          E r0 = a + c, r1 = a - c, i0 = b - d, i1 = b + d;
 
344
          pp[0] = r0;
 
345
          pm[0] = i0;
 
346
          pp[im] = r1 * wr - i1 * wi;
 
347
          pm[im] = i1 * wr + r1 * wi;
 
348
          pp += dist; pm -= dist;
 
349
     }
 
350
 
 
351
     /* i = n/2 when n is even */
 
352
     if (!(n & 1)) { pp[0] *= K(2.0); pp[im] *= -K(2.0); }
 
353
}
 
354
 
 
355
static void apply_b_rdft(const plan *ego_, R *r, R *rio, R *iio)
 
356
{
 
357
     const P *ego = (const P *) ego_;
 
358
 
 
359
     {
 
360
          INT i, vl = ego->vl, n2 = ego->n / 2;
 
361
          INT ivs = ego->ivs, is = ego->is;
 
362
          const R *W = ego->td->W;
 
363
          R *rio1 = rio, *iio1 = iio;
 
364
          for (i = 0; i < vl; ++i, rio1 += ivs, iio1 += ivs)
 
365
               k_b_rdft(rio1, iio1, W, n2, is);
 
366
     }
 
367
 
 
368
     {
 
369
          plan_rdft *cld = (plan_rdft *) ego->cld;
 
370
          cld->apply((plan *) cld, rio, r);
 
371
     }
 
372
}
 
373
 
 
374
static problem *mkcld_b_rdft(const problem_rdft2 *p)
 
375
{
 
376
     const iodim *d = p->sz->dims;
 
377
 
 
378
     tensor *radix = X(mktensor_1d)(2, p->iio - p->rio, d[0].os);
 
379
     tensor *cld_vec = X(tensor_append)(radix, p->vecsz);
 
380
     X(tensor_destroy)(radix);
 
381
 
 
382
     return X(mkproblem_rdft_1_d) (
 
383
          X(mktensor_1d)(d[0].n / 2, d[0].is, 2 * d[0].os),
 
384
          cld_vec, p->rio, p->r, HC2R);
 
385
}
 
386
 
 
387
static const madt adt_b_rdft = {
 
388
     applicable_b, apply_b_rdft, mkcld_b_rdft, {6, 4, 0, 0}, "hc2r2-rdft"
 
389
};
 
390
 
 
391
/*
 
392
 * common stuff
 
393
 */
 
394
static void awake(plan *ego_, enum wakefulness wakefulness)
 
395
{
 
396
     P *ego = (P *) ego_;
 
397
     static const tw_instr twinstr[] = { {TW_FULL, 0, 2}, {TW_NEXT, 1, 0} };
 
398
     X(plan_awake)(ego->cld, wakefulness);
 
399
     X(twiddle_awake)(wakefulness, &ego->td, twinstr, 
 
400
                      ego->n, 2, (ego->n / 2 + 1) / 2);
 
401
}
 
402
 
 
403
static void destroy(plan *ego_)
 
404
{
 
405
     P *ego = (P *) ego_;
 
406
     X(plan_destroy_internal) (ego->cld);
 
407
}
 
408
 
 
409
static void print(const plan *ego_, printer * p)
 
410
{
 
411
     const P *ego = (const P *) ego_;
 
412
     p->print(p, "(%s-%D%v%(%p%))", ego->slv->adt->nam,
 
413
              ego->n, ego->vl, ego->cld);
 
414
}
 
415
 
 
416
static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
 
417
{
 
418
     const S *ego = (const S *) ego_;
 
419
     P *pln;
 
420
     const problem_rdft2 *p;
 
421
     plan *cld;
 
422
     const iodim *d;
 
423
 
 
424
     static const plan_adt padt = {
 
425
          X(rdft2_solve), awake, print, destroy
 
426
     };
 
427
 
 
428
     if (!ego->adt->applicable(p_, plnr))
 
429
          return (plan *) 0;
 
430
 
 
431
     p = (const problem_rdft2 *) p_;
 
432
 
 
433
     cld = X(mkplan_d)(plnr, ego->adt->mkcld(p));
 
434
     if (!cld) return (plan *) 0;
 
435
 
 
436
     pln = MKPLAN_RDFT2(P, &padt, ego->adt->apply);
 
437
 
 
438
     d = p->sz->dims;
 
439
     pln->n = d[0].n;
 
440
     pln->os = d[0].os;
 
441
     pln->is = d[0].is;
 
442
     X(tensor_tornk1) (p->vecsz, &pln->vl, &pln->ivs, &pln->ovs);
 
443
     pln->cld = cld;
 
444
     pln->td = 0;
 
445
     pln->slv = ego;
 
446
 
 
447
     /* approximately */
 
448
     X(ops_madd)(pln->vl * ((pln->n/2 + 1) / 2), &ego->adt->ops,
 
449
                 &cld->ops, &pln->super.super.ops);
 
450
 
 
451
     return &(pln->super.super);
 
452
}
 
453
 
 
454
static solver *mksolver(const madt *adt)
 
455
{
 
456
     static const solver_adt sadt = { PROBLEM_RDFT2, mkplan };
 
457
     S *slv = MKSOLVER(S, &sadt);
 
458
     slv->adt = adt;
 
459
     return &(slv->super);
 
460
}
 
461
 
 
462
void X(rdft2_radix2_register)(planner *p)
 
463
{
 
464
     unsigned i;
 
465
     static const madt *const adts[] = {
 
466
          &adt_f_dft, &adt_f_rdft,
 
467
          &adt_b_dft, &adt_b_rdft
 
468
     };
 
469
 
 
470
     for (i = 0; i < sizeof(adts) / sizeof(adts[0]); ++i)
 
471
          REGISTER_SOLVER(p, mksolver(adts[i]));
 
472
}