diff options
Diffstat (limited to 'sndobj/rfftw/rfftwnd.c')
| -rw-r--r-- | sndobj/rfftw/rfftwnd.c | 570 |
1 files changed, 570 insertions, 0 deletions
diff --git a/sndobj/rfftw/rfftwnd.c b/sndobj/rfftw/rfftwnd.c new file mode 100644 index 0000000..91814f9 --- /dev/null +++ b/sndobj/rfftw/rfftwnd.c @@ -0,0 +1,570 @@ +
+/*
+ * Copyright (c) 1997-1999 Massachusetts Institute of Technology
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ */
+
+/* $Id: rfftwnd.c,v 1.1.1.1 2006/05/12 15:14:48 veplaini Exp $ */
+
+#include <fftw-int.h>
+#include <rfftw.h>
+
+/********************** prototypes for rexec2 routines **********************/
+
+extern void rfftw_real2c_aux(fftw_plan plan, int howmany,
+ fftw_real *in, int istride, int idist,
+ fftw_complex *out, int ostride, int odist,
+ fftw_real *work);
+extern void rfftw_c2real_aux(fftw_plan plan, int howmany,
+ fftw_complex *in, int istride, int idist,
+ fftw_real *out, int ostride, int odist,
+ fftw_real *work);
+extern void rfftw_real2c_overlap_aux(fftw_plan plan, int howmany,
+ fftw_real *in, int istride, int idist,
+ fftw_complex *out, int ostride, int odist,
+ fftw_real *work);
+extern void rfftw_c2real_overlap_aux(fftw_plan plan, int howmany,
+ fftw_complex *in, int istride, int idist,
+ fftw_real *out, int ostride, int odist,
+ fftw_real *work);
+
+/********************** Initializing the RFFTWND Plan ***********************/
+
+/*
+ * Create an fftwnd_plan specialized for specific arrays. (These
+ * arrays are ignored, however, if they are NULL or if the flags
+ * do not include FFTW_MEASURE.) The main advantage of being
+ * provided arrays like this is that we can do runtime timing
+ * measurements of our options, without worrying about allocating
+ * excessive scratch space.
+ */
+fftwnd_plan rfftwnd_create_plan_specific(int rank, const int *n,
+ fftw_direction dir, int flags,
+ fftw_real *in, int istride,
+ fftw_real *out, int ostride)
+{
+ fftwnd_plan p;
+ int i;
+ int rflags = flags & ~FFTW_IN_PLACE;
+ /* note that we always do rfftw transforms out-of-place in rexec2.c */
+
+ if (flags & FFTW_IN_PLACE) {
+ out = NULL;
+ ostride = istride;
+ }
+ istride = ostride = 1; /*
+ * strides don't work yet, since it is not
+ * clear whether they apply to real
+ * or complex data
+ */
+
+ if (!(p = fftwnd_create_plan_aux(rank, n, dir, flags)))
+ return 0;
+
+ for (i = 0; i < rank - 1; ++i)
+ p->n_after[i] = (n[rank - 1]/2 + 1) * (p->n_after[i] / n[rank - 1]);
+ if (rank > 0)
+ p->n[rank - 1] = n[rank - 1] / 2 + 1;
+
+ p->plans = fftwnd_new_plan_array(rank);
+ if (rank > 0 && !p->plans) {
+ rfftwnd_destroy_plan(p);
+ return 0;
+ }
+ if (rank > 0) {
+ p->plans[rank - 1] = rfftw_create_plan(n[rank - 1], dir, rflags);
+ if (!p->plans[rank - 1]) {
+ rfftwnd_destroy_plan(p);
+ return 0;
+ }
+ }
+ if (rank > 1) {
+ if (!(flags & FFTW_MEASURE) || in == 0
+ || (!p->is_in_place && out == 0)) {
+ if (!fftwnd_create_plans_generic(p->plans, rank - 1, n,
+ dir, flags | FFTW_IN_PLACE)) {
+ rfftwnd_destroy_plan(p);
+ return 0;
+ }
+ } else if (dir == FFTW_COMPLEX_TO_REAL || (flags & FFTW_IN_PLACE)) {
+ if (!fftwnd_create_plans_specific(p->plans, rank - 1, n,
+ p->n_after,
+ dir, flags | FFTW_IN_PLACE,
+ (fftw_complex *) in,
+ istride,
+ 0, 0)) {
+ rfftwnd_destroy_plan(p);
+ return 0;
+ }
+ } else {
+ if (!fftwnd_create_plans_specific(p->plans, rank - 1, n,
+ p->n_after,
+ dir, flags | FFTW_IN_PLACE,
+ (fftw_complex *) out,
+ ostride,
+ 0, 0)) {
+ rfftwnd_destroy_plan(p);
+ return 0;
+ }
+ }
+ }
+ p->nbuffers = 0;
+ p->nwork = fftwnd_work_size(rank, p->n, flags | FFTW_IN_PLACE,
+ p->nbuffers + 1);
+ if (p->nwork && !(flags & FFTW_THREADSAFE)) {
+ p->work = (fftw_complex *) fftw_malloc(p->nwork
+ * sizeof(fftw_complex));
+ if (!p->work) {
+ rfftwnd_destroy_plan(p);
+ return 0;
+ }
+ }
+ return p;
+}
+
+fftwnd_plan rfftw2d_create_plan_specific(int nx, int ny,
+ fftw_direction dir, int flags,
+ fftw_real *in, int istride,
+ fftw_real *out, int ostride)
+{
+ int n[2];
+
+ n[0] = nx;
+ n[1] = ny;
+
+ return rfftwnd_create_plan_specific(2, n, dir, flags,
+ in, istride, out, ostride);
+}
+
+fftwnd_plan rfftw3d_create_plan_specific(int nx, int ny, int nz,
+ fftw_direction dir, int flags,
+ fftw_real *in, int istride,
+ fftw_real *out, int ostride)
+{
+ int n[3];
+
+ n[0] = nx;
+ n[1] = ny;
+ n[2] = nz;
+
+ return rfftwnd_create_plan_specific(3, n, dir, flags,
+ in, istride, out, ostride);
+}
+
+/* Create a generic fftwnd plan: */
+
+fftwnd_plan rfftwnd_create_plan(int rank, const int *n,
+ fftw_direction dir, int flags)
+{
+ return rfftwnd_create_plan_specific(rank, n, dir, flags, 0, 1, 0, 1);
+}
+
+fftwnd_plan rfftw2d_create_plan(int nx, int ny,
+ fftw_direction dir, int flags)
+{
+ return rfftw2d_create_plan_specific(nx, ny, dir, flags, 0, 1, 0, 1);
+}
+
+fftwnd_plan rfftw3d_create_plan(int nx, int ny, int nz,
+ fftw_direction dir, int flags)
+{
+ return rfftw3d_create_plan_specific(nx, ny, nz, dir, flags, 0, 1, 0, 1);
+}
+
+/************************ Freeing the RFFTWND Plan ************************/
+
+void rfftwnd_destroy_plan(fftwnd_plan plan)
+{
+ fftwnd_destroy_plan(plan);
+}
+
+/************************ Printing the RFFTWND Plan ************************/
+
+void rfftwnd_fprint_plan(FILE *f, fftwnd_plan plan)
+{
+ fftwnd_fprint_plan(f, plan);
+}
+
+void rfftwnd_print_plan(fftwnd_plan plan)
+{
+ rfftwnd_fprint_plan(stdout, plan);
+}
+
+/*********** Computing the N-Dimensional FFT: Auxiliary Routines ************/
+
+void rfftwnd_real2c_aux(fftwnd_plan p, int cur_dim,
+ fftw_real *in, int istride,
+ fftw_complex *out, int ostride,
+ fftw_real *work)
+{
+ int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
+
+ if (cur_dim == p->rank - 2) {
+ /* just do the last dimension directly: */
+ if (p->is_in_place)
+ rfftw_real2c_aux(p->plans[p->rank - 1], n,
+ in, istride, (n_after * istride) * 2,
+ out, istride, n_after * istride,
+ work);
+ else
+ rfftw_real2c_aux(p->plans[p->rank - 1], n,
+ in, istride, p->plans[p->rank - 1]->n * istride,
+ out, ostride, n_after * ostride,
+ work);
+ } else { /* we have at least two dimensions to go */
+ int nr = p->plans[p->rank - 1]->n;
+ int n_after_r = p->is_in_place ? n_after * 2
+ : nr * (n_after / (nr/2 + 1));
+ int i;
+
+ /*
+ * process the subsequent dimensions recursively, in hyperslabs,
+ * to get maximum locality:
+ */
+ for (i = 0; i < n; ++i)
+ rfftwnd_real2c_aux(p, cur_dim + 1,
+ in + i * n_after_r * istride, istride,
+ out + i * n_after * ostride, ostride, work);
+ }
+
+ /* do the current dimension (in-place): */
+ fftw(p->plans[cur_dim], n_after,
+ out, n_after * ostride, ostride,
+ (fftw_complex *) work, 1, 0);
+ /* I hate this cast */
+}
+
+void rfftwnd_c2real_aux(fftwnd_plan p, int cur_dim,
+ fftw_complex *in, int istride,
+ fftw_real *out, int ostride,
+ fftw_real *work)
+{
+ int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
+
+ /* do the current dimension (in-place): */
+ fftw(p->plans[cur_dim], n_after,
+ in, n_after * istride, istride,
+ (fftw_complex *) work, 1, 0);
+
+ if (cur_dim == p->rank - 2) {
+ /* just do the last dimension directly: */
+ if (p->is_in_place)
+ rfftw_c2real_aux(p->plans[p->rank - 1], n,
+ in, istride, n_after * istride,
+ out, istride, (n_after * istride) * 2,
+ work);
+ else
+ rfftw_c2real_aux(p->plans[p->rank - 1], n,
+ in, istride, n_after * istride,
+ out, ostride, p->plans[p->rank - 1]->n * ostride,
+ work);
+ } else { /* we have at least two dimensions to go */
+ int nr = p->plans[p->rank - 1]->n;
+ int n_after_r = p->is_in_place ? n_after * 2 :
+ nr * (n_after / (nr/2 + 1));
+ int i;
+
+ /*
+ * process the subsequent dimensions recursively, in hyperslabs,
+ * to get maximum locality:
+ */
+ for (i = 0; i < n; ++i)
+ rfftwnd_c2real_aux(p, cur_dim + 1,
+ in + i * n_after * istride, istride,
+ out + i * n_after_r * ostride, ostride, work);
+ }
+}
+
+/*
+ * alternate version of rfftwnd_aux -- this version pushes the howmany
+ * loop down to the leaves of the computation, for greater locality
+ * in cases where dist < stride. It is also required for correctness
+ * if in==out, and we must call a special version of the executor.
+ * Note that work must point to 'howmany' copies of its data
+ * if in == out.
+ */
+
+void rfftwnd_real2c_aux_howmany(fftwnd_plan p, int cur_dim,
+ int howmany,
+ fftw_real *in, int istride, int idist,
+ fftw_complex *out, int ostride, int odist,
+ fftw_complex *work)
+{
+ int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
+ int k;
+
+ if (cur_dim == p->rank - 2) {
+ /* just do the last dimension directly: */
+ if (p->is_in_place)
+ for (k = 0; k < n; ++k)
+ rfftw_real2c_overlap_aux(p->plans[p->rank - 1], howmany,
+ in + (k * n_after * istride) * 2,
+ istride, idist,
+ out + (k * n_after * ostride),
+ ostride, odist,
+ (fftw_real *) work);
+ else {
+ int nlast = p->plans[p->rank - 1]->n;
+ for (k = 0; k < n; ++k)
+ rfftw_real2c_aux(p->plans[p->rank - 1], howmany,
+ in + k * nlast * istride,
+ istride, idist,
+ out + k * n_after * ostride,
+ ostride, odist,
+ (fftw_real *) work);
+ }
+ } else { /* we have at least two dimensions to go */
+ int nr = p->plans[p->rank - 1]->n;
+ int n_after_r = p->is_in_place ? n_after * 2 :
+ nr * (n_after / (nr/2 + 1));
+ int i;
+
+ /*
+ * process the subsequent dimensions recursively, in hyperslabs,
+ * to get maximum locality:
+ */
+ for (i = 0; i < n; ++i)
+ rfftwnd_real2c_aux_howmany(p, cur_dim + 1, howmany,
+ in + i * n_after_r * istride, istride, idist,
+ out + i * n_after * ostride, ostride, odist,
+ work);
+ }
+
+ /* do the current dimension (in-place): */
+ for (k = 0; k < n_after; ++k)
+ fftw(p->plans[cur_dim], howmany,
+ out + k * ostride, n_after * ostride, odist,
+ work, 1, 0);
+}
+
+void rfftwnd_c2real_aux_howmany(fftwnd_plan p, int cur_dim,
+ int howmany,
+ fftw_complex *in, int istride, int idist,
+ fftw_real *out, int ostride, int odist,
+ fftw_complex *work)
+{
+ int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
+ int k;
+
+ /* do the current dimension (in-place): */
+ for (k = 0; k < n_after; ++k)
+ fftw(p->plans[cur_dim], howmany,
+ in + k * istride, n_after * istride, idist,
+ work, 1, 0);
+
+ if (cur_dim == p->rank - 2) {
+ /* just do the last dimension directly: */
+ if (p->is_in_place)
+ for (k = 0; k < n; ++k)
+ rfftw_c2real_overlap_aux(p->plans[p->rank - 1], howmany,
+ in + (k * n_after * istride),
+ istride, idist,
+ out + (k * n_after * ostride) * 2,
+ ostride, odist,
+ (fftw_real *) work);
+ else {
+ int nlast = p->plans[p->rank - 1]->n;
+ for (k = 0; k < n; ++k)
+ rfftw_c2real_aux(p->plans[p->rank - 1], howmany,
+ in + k * n_after * istride,
+ istride, idist,
+ out + k * nlast * ostride,
+ ostride, odist,
+ (fftw_real *) work);
+ }
+ } else { /* we have at least two dimensions to go */
+ int nr = p->plans[p->rank - 1]->n;
+ int n_after_r = p->is_in_place ? n_after * 2
+ : nr * (n_after / (nr/2 + 1));
+ int i;
+
+ /*
+ * process the subsequent dimensions recursively, in hyperslabs,
+ * to get maximum locality:
+ */
+ for (i = 0; i < n; ++i)
+ rfftwnd_c2real_aux_howmany(p, cur_dim + 1, howmany,
+ in + i * n_after * istride, istride, idist,
+ out + i * n_after_r * ostride, ostride, odist,
+ work);
+ }
+}
+
+/********** Computing the N-Dimensional FFT: User-Visible Routines **********/
+
+void rfftwnd_real_to_complex(fftwnd_plan p, int howmany,
+ fftw_real *in, int istride, int idist,
+ fftw_complex *out, int ostride, int odist)
+{
+ fftw_complex *work = p->work;
+ int rank = p->rank;
+ int free_work = 0;
+
+ if (p->dir != FFTW_REAL_TO_COMPLEX)
+ fftw_die("rfftwnd_real_to_complex with complex-to-real plan");
+
+#ifdef FFTW_DEBUG
+ if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE)
+ && p->nwork && p->work)
+ fftw_die("bug with FFTW_THREADSAFE flag");
+#endif
+
+ if (p->is_in_place) {
+ ostride = istride;
+ odist = (idist == 1 && idist < istride) ? 1 : (idist / 2); /* ugh */
+ out = (fftw_complex *) in;
+ if (howmany > 1 && istride > idist && rank > 0) {
+ int new_nwork;
+
+ new_nwork = p->n[rank - 1] * howmany;
+ if (new_nwork > p->nwork) {
+ work = (fftw_complex *)
+ fftw_malloc(sizeof(fftw_complex) * new_nwork);
+ if (!work)
+ fftw_die("error allocating work array");
+ free_work = 1;
+ }
+ }
+ }
+ if (p->nwork && !work) {
+ work = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * p->nwork);
+ free_work = 1;
+ }
+ switch (rank) {
+ case 0:
+ break;
+ case 1:
+ if (p->is_in_place && howmany > 1 && istride > idist)
+ rfftw_real2c_overlap_aux(p->plans[0], howmany,
+ in, istride, idist,
+ out, ostride, odist,
+ (fftw_real *) work);
+ else
+ rfftw_real2c_aux(p->plans[0], howmany,
+ in, istride, idist,
+ out, ostride, odist,
+ (fftw_real *) work);
+ break;
+ default: /* rank >= 2 */
+ {
+ if (howmany > 1 && ostride > odist)
+ rfftwnd_real2c_aux_howmany(p, 0, howmany,
+ in, istride, idist,
+ out, ostride, odist,
+ work);
+ else {
+ int i;
+
+ for (i = 0; i < howmany; ++i)
+ rfftwnd_real2c_aux(p, 0,
+ in + i * idist, istride,
+ out + i * odist, ostride,
+ (fftw_real *) work);
+ }
+ }
+ }
+
+ if (free_work)
+ fftw_free(work);
+}
+
+void rfftwnd_complex_to_real(fftwnd_plan p, int howmany,
+ fftw_complex *in, int istride, int idist,
+ fftw_real *out, int ostride, int odist)
+{
+ fftw_complex *work = p->work;
+ int rank = p->rank;
+ int free_work = 0;
+
+ if (p->dir != FFTW_COMPLEX_TO_REAL)
+ fftw_die("rfftwnd_complex_to_real with real-to-complex plan");
+
+#ifdef FFTW_DEBUG
+ if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE)
+ && p->nwork && p->work)
+ fftw_die("bug with FFTW_THREADSAFE flag");
+#endif
+
+ if (p->is_in_place) {
+ ostride = istride;
+ odist = idist;
+ odist = (idist == 1 && idist < istride) ? 1 : (idist * 2); /* ugh */
+ out = (fftw_real *) in;
+ if (howmany > 1 && istride > idist && rank > 0) {
+ int new_nwork = p->n[rank - 1] * howmany;
+ if (new_nwork > p->nwork) {
+ work = (fftw_complex *)
+ fftw_malloc(sizeof(fftw_complex) * new_nwork);
+ if (!work)
+ fftw_die("error allocating work array");
+ free_work = 1;
+ }
+ }
+ }
+ if (p->nwork && !work) {
+ work = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * p->nwork);
+ free_work = 1;
+ }
+ switch (rank) {
+ case 0:
+ break;
+ case 1:
+ if (p->is_in_place && howmany > 1 && istride > idist)
+ rfftw_c2real_overlap_aux(p->plans[0], howmany,
+ in, istride, idist,
+ out, ostride, odist,
+ (fftw_real *) work);
+ else
+ rfftw_c2real_aux(p->plans[0], howmany,
+ in, istride, idist,
+ out, ostride, odist,
+ (fftw_real *) work);
+ break;
+ default: /* rank >= 2 */
+ {
+ if (howmany > 1 && ostride > odist)
+ rfftwnd_c2real_aux_howmany(p, 0, howmany,
+ in, istride, idist,
+ out, ostride, odist,
+ work);
+ else {
+ int i;
+
+ for (i = 0; i < howmany; ++i)
+ rfftwnd_c2real_aux(p, 0,
+ in + i * idist, istride,
+ out + i * odist, ostride,
+ (fftw_real *) work);
+ }
+ }
+ }
+
+ if (free_work)
+ fftw_free(work);
+}
+
+void rfftwnd_one_real_to_complex(fftwnd_plan p,
+ fftw_real *in, fftw_complex *out)
+{
+ rfftwnd_real_to_complex(p, 1, in, 1, 1, out, 1, 1);
+}
+
+void rfftwnd_one_complex_to_real(fftwnd_plan p,
+ fftw_complex *in, fftw_real *out)
+{
+ rfftwnd_complex_to_real(p, 1, in, 1, 1, out, 1, 1);
+}
|