1
2 /*
3 * Copyright (c) 1997-1999, 2003 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 /*
22 * rgeneric.c -- "generic" rfftw codelets. They work for all n (and
23 * they are slow)
24 */
25 #include "fftw-int.h"
26 #include "rfftw.h"
27
28 /* this code assumes that r and m are both odd */
fftw_hc2hc_forward_generic(fftw_real * A,const fftw_complex * W,int m,int r,int n,int dist)29 void fftw_hc2hc_forward_generic(fftw_real *A, const fftw_complex *W,
30 int m, int r, int n, int dist)
31 {
32 int i, j, k;
33 fftw_complex *tmp = (fftw_complex *)
34 fftw_malloc(r * sizeof(fftw_complex));
35 fftw_real rsum, isum;
36 fftw_real *X, *YO, *YI;
37 int wp, wincr;
38 int iostride = m * dist;
39 X = A;
40 YO = A + r * iostride;
41 YI = A + iostride;
42
43 /* compute the transform of the r 0th elements (which are real) */
44 for (i = 0; i + i < r; ++i) {
45 rsum = 0.0;
46 isum = 0.0;
47 wincr = m * i;
48 for (j = 0, wp = 0; j < r; ++j) {
49 fftw_real tw_r = c_re(W[wp]);
50 fftw_real tw_i = c_im(W[wp]);
51 fftw_real re = X[j * iostride];
52 rsum += re * tw_r;
53 isum += re * tw_i;
54 wp += wincr;
55 if (wp >= n)
56 wp -= n;
57 }
58 c_re(tmp[i]) = rsum;
59 c_im(tmp[i]) = isum;
60 }
61
62 /* store the transform back onto the A array */
63 X[0] = c_re(tmp[0]);
64 for (i = 1; i + i < r; ++i) {
65 X[i * iostride] = c_re(tmp[i]);
66 YO[-i * iostride] = c_im(tmp[i]);
67 }
68
69 X += dist;
70 YI -= dist;
71 YO -= dist;
72
73 /* compute the transform of the middle elements (which are complex) */
74 for (k = 1; k + k < m; ++k, X += dist, YI -= dist, YO -= dist) {
75 for (i = 0; i < r; ++i) {
76 rsum = 0.0;
77 isum = 0.0;
78 wincr = k + m * i;
79 for (j = 0, wp = 0; j < r; ++j) {
80 fftw_real tw_r = c_re(W[wp]);
81 fftw_real tw_i = c_im(W[wp]);
82 fftw_real re = X[j * iostride];
83 fftw_real im = YI[j * iostride];
84 rsum += re * tw_r - im * tw_i;
85 isum += re * tw_i + im * tw_r;
86 wp += wincr;
87 if (wp >= n)
88 wp -= n;
89 }
90 c_re(tmp[i]) = rsum;
91 c_im(tmp[i]) = isum;
92 }
93
94 /* store the transform back onto the A array */
95 for (i = 0; i + i < r; ++i) {
96 X[i * iostride] = c_re(tmp[i]);
97 YO[-i * iostride] = c_im(tmp[i]);
98 }
99 for (; i < r; ++i) {
100 X[i * iostride] = -c_im(tmp[i]);
101 YO[-i * iostride] = c_re(tmp[i]);
102 }
103 }
104
105 /* no final element, since m is odd */
106 fftw_free(tmp);
107 }
108
fftw_hc2hc_backward_generic(fftw_real * A,const fftw_complex * W,int m,int r,int n,int dist)109 void fftw_hc2hc_backward_generic(fftw_real *A, const fftw_complex *W,
110 int m, int r, int n, int dist)
111 {
112 int i, j, k;
113 int wp, wincr;
114 fftw_complex *tmp = (fftw_complex *)
115 fftw_malloc(r * sizeof(fftw_complex));
116 fftw_real rsum, isum;
117 fftw_real *X, *YO, *YI;
118 int iostride = m * dist;
119 X = A;
120 YO = A + iostride;
121 YI = A + r * iostride;
122
123 /*
124 * compute the transform of the r 0th elements (which are halfcomplex)
125 * yielding real numbers
126 */
127 /* copy the input into the temporary array */
128 c_re(tmp[0]) = X[0];
129 for (i = 1; i + i < r; ++i) {
130 c_re(tmp[i]) = X[i * iostride];
131 c_im(tmp[i]) = YI[-i * iostride];
132 }
133
134 for (i = 0; i < r; ++i) {
135 rsum = 0.0;
136 wincr = m * i;
137 for (j = 1, wp = wincr; j + j < r; ++j) {
138 fftw_real tw_r = c_re(W[wp]);
139 fftw_real tw_i = c_im(W[wp]);
140 fftw_real re = c_re(tmp[j]);
141 fftw_real im = c_im(tmp[j]);
142 rsum += re * tw_r + im * tw_i;
143 wp += wincr;
144 if (wp >= n)
145 wp -= n;
146 }
147 X[i * iostride] = 2.0 * rsum + c_re(tmp[0]);
148 }
149
150 X += dist;
151 YI -= dist;
152 YO -= dist;
153
154 /* compute the transform of the middle elements (which are complex) */
155 for (k = 1; k + k < m; ++k, X += dist, YI -= dist, YO -= dist) {
156 /* copy the input into the temporary array */
157 for (i = 0; i + i < r; ++i) {
158 c_re(tmp[i]) = X[i * iostride];
159 c_im(tmp[i]) = YI[-i * iostride];
160 }
161 for (; i < r; ++i) {
162 c_im(tmp[i]) = -X[i * iostride];
163 c_re(tmp[i]) = YI[-i * iostride];
164 }
165
166 for (i = 0; i < r; ++i) {
167 rsum = 0.0;
168 isum = 0.0;
169 wincr = m * i;
170 for (j = 0, wp = k * i; j < r; ++j) {
171 fftw_real tw_r = c_re(W[wp]);
172 fftw_real tw_i = c_im(W[wp]);
173 fftw_real re = c_re(tmp[j]);
174 fftw_real im = c_im(tmp[j]);
175 rsum += re * tw_r + im * tw_i;
176 isum += im * tw_r - re * tw_i;
177 wp += wincr;
178 if (wp >= n)
179 wp -= n;
180 }
181 X[i * iostride] = rsum;
182 YO[i * iostride] = isum;
183 }
184 }
185
186 /* no final element, since m is odd */
187 fftw_free(tmp);
188 }
189