Subversion Repositories shark

Rev

Rev 2 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 1
 
2
/*
3
 * Copyright (c) 1997-1999 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 <ports/fftw-int.h>
26
#include <ports/rfftw.h>
27
 
28
/* this code assumes that r and m are both odd */
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
 
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
}