From 416bd737074a287ea47106c73ea6bcfde40a75a8 Mon Sep 17 00:00:00 2001 From: John Glover Date: Fri, 24 Jun 2011 18:17:23 +0100 Subject: Change to using distutils. Currently only builds the simplsndobj module --- src/sndobj/rfftw/rgeneric.c | 188 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 188 insertions(+) create mode 100644 src/sndobj/rfftw/rgeneric.c (limited to 'src/sndobj/rfftw/rgeneric.c') diff --git a/src/sndobj/rfftw/rgeneric.c b/src/sndobj/rfftw/rgeneric.c new file mode 100644 index 0000000..7feb92d --- /dev/null +++ b/src/sndobj/rfftw/rgeneric.c @@ -0,0 +1,188 @@ + +/* + * 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 + * + */ + +/* + * rgeneric.c -- "generic" rfftw codelets. They work for all n (and + * they are slow) + */ +#include +#include + +/* this code assumes that r and m are both odd */ +void fftw_hc2hc_forward_generic(fftw_real *A, const fftw_complex *W, + int m, int r, int n, int dist) +{ + int i, j, k; + fftw_complex *tmp = (fftw_complex *) + fftw_malloc(r * sizeof(fftw_complex)); + fftw_real rsum, isum; + fftw_real *X, *YO, *YI; + int wp, wincr; + int iostride = m * dist; + X = A; + YO = A + r * iostride; + YI = A + iostride; + + /* compute the transform of the r 0th elements (which are real) */ + for (i = 0; i + i < r; ++i) { + rsum = 0.0; + isum = 0.0; + wincr = m * i; + for (j = 0, wp = 0; j < r; ++j) { + fftw_real tw_r = c_re(W[wp]); + fftw_real tw_i = c_im(W[wp]); + fftw_real re = X[j * iostride]; + rsum += re * tw_r; + isum += re * tw_i; + wp += wincr; + if (wp >= n) + wp -= n; + } + c_re(tmp[i]) = rsum; + c_im(tmp[i]) = isum; + } + + /* store the transform back onto the A array */ + X[0] = c_re(tmp[0]); + for (i = 1; i + i < r; ++i) { + X[i * iostride] = c_re(tmp[i]); + YO[-i * iostride] = c_im(tmp[i]); + } + + X += dist; + YI -= dist; + YO -= dist; + + /* compute the transform of the middle elements (which are complex) */ + for (k = 1; k + k < m; ++k, X += dist, YI -= dist, YO -= dist) { + for (i = 0; i < r; ++i) { + rsum = 0.0; + isum = 0.0; + wincr = k + m * i; + for (j = 0, wp = 0; j < r; ++j) { + fftw_real tw_r = c_re(W[wp]); + fftw_real tw_i = c_im(W[wp]); + fftw_real re = X[j * iostride]; + fftw_real im = YI[j * iostride]; + rsum += re * tw_r - im * tw_i; + isum += re * tw_i + im * tw_r; + wp += wincr; + if (wp >= n) + wp -= n; + } + c_re(tmp[i]) = rsum; + c_im(tmp[i]) = isum; + } + + /* store the transform back onto the A array */ + for (i = 0; i + i < r; ++i) { + X[i * iostride] = c_re(tmp[i]); + YO[-i * iostride] = c_im(tmp[i]); + } + for (; i < r; ++i) { + X[i * iostride] = -c_im(tmp[i]); + YO[-i * iostride] = c_re(tmp[i]); + } + } + + /* no final element, since m is odd */ + fftw_free(tmp); +} + +void fftw_hc2hc_backward_generic(fftw_real *A, const fftw_complex *W, + int m, int r, int n, int dist) +{ + int i, j, k; + int wp, wincr; + fftw_complex *tmp = (fftw_complex *) + fftw_malloc(r * sizeof(fftw_complex)); + fftw_real rsum, isum; + fftw_real *X, *YO, *YI; + int iostride = m * dist; + X = A; + YO = A + iostride; + YI = A + r * iostride; + + /* + * compute the transform of the r 0th elements (which are halfcomplex) + * yielding real numbers + */ + /* copy the input into the temporary array */ + c_re(tmp[0]) = X[0]; + for (i = 1; i + i < r; ++i) { + c_re(tmp[i]) = X[i * iostride]; + c_im(tmp[i]) = YI[-i * iostride]; + } + + for (i = 0; i < r; ++i) { + rsum = 0.0; + wincr = m * i; + for (j = 1, wp = wincr; j + j < r; ++j) { + fftw_real tw_r = c_re(W[wp]); + fftw_real tw_i = c_im(W[wp]); + fftw_real re = c_re(tmp[j]); + fftw_real im = c_im(tmp[j]); + rsum += re * tw_r + im * tw_i; + wp += wincr; + if (wp >= n) + wp -= n; + } + X[i * iostride] = 2.0 * rsum + c_re(tmp[0]); + } + + X += dist; + YI -= dist; + YO -= dist; + + /* compute the transform of the middle elements (which are complex) */ + for (k = 1; k + k < m; ++k, X += dist, YI -= dist, YO -= dist) { + /* copy the input into the temporary array */ + for (i = 0; i + i < r; ++i) { + c_re(tmp[i]) = X[i * iostride]; + c_im(tmp[i]) = YI[-i * iostride]; + } + for (; i < r; ++i) { + c_im(tmp[i]) = -X[i * iostride]; + c_re(tmp[i]) = YI[-i * iostride]; + } + + for (i = 0; i < r; ++i) { + rsum = 0.0; + isum = 0.0; + wincr = m * i; + for (j = 0, wp = k * i; j < r; ++j) { + fftw_real tw_r = c_re(W[wp]); + fftw_real tw_i = c_im(W[wp]); + fftw_real re = c_re(tmp[j]); + fftw_real im = c_im(tmp[j]); + rsum += re * tw_r + im * tw_i; + isum += im * tw_r - re * tw_i; + wp += wincr; + if (wp >= n) + wp -= n; + } + X[i * iostride] = rsum; + YO[i * iostride] = isum; + } + } + + /* no final element, since m is odd */ + fftw_free(tmp); +} -- cgit v1.2.3