Go to most recent revision |
Blame |
Last modification |
View Log
| RSS feed
/*
* 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 <ports/fftw-int.h>
#include <ports/rfftw.h>
/* 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);
}