Mercurial > octave-nkf

comparison liboctave/oct-fftw.cc @ 4785:d3018a33c584

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-02-17 21:58:17 by jwe]
author jwe
date Tue, 17 Feb 2004 21:58:17 +0000
parents 743ef6154f8a
children fc316bde0053
comparison
equal deleted inserted replaced
4784:743ef6154f8a 4785:d3018a33c584
213 213
214 static inline void 214 static inline void
215 convert_packcomplex_1d (Complex *out, size_t nr, size_t nc, 215 convert_packcomplex_1d (Complex *out, size_t nr, size_t nc,
216 int stride, int dist) 216 int stride, int dist)
217 { 217 {
218 // Fill in the missing data 218 OCTAVE_QUIT;
219
220 // Fill in the missing data.
221
219 for (size_t i = 0; i < nr; i++) 222 for (size_t i = 0; i < nr; i++)
220 for (size_t j = nc/2+1; j < nc; j++) 223 for (size_t j = nc/2+1; j < nc; j++)
221 out[j*stride + i*dist] = conj(out[(nc - j)*stride + i*dist]); 224 out[j*stride + i*dist] = conj(out[(nc - j)*stride + i*dist]);
225
226 OCTAVE_QUIT;
222 } 227 }
223 228
224 static inline void 229 static inline void
225 convert_packcomplex_Nd (Complex *out, const dim_vector &dv) 230 convert_packcomplex_Nd (Complex *out, const dim_vector &dv)
226 { 231 {
228 size_t nr = dv(1); 233 size_t nr = dv(1);
229 size_t np = (dv.length() > 2 ? dv.numel () / nc / nr : 1); 234 size_t np = (dv.length() > 2 ? dv.numel () / nc / nr : 1);
230 size_t nrp = nr * np; 235 size_t nrp = nr * np;
231 Complex *ptr1, *ptr2; 236 Complex *ptr1, *ptr2;
232 237
233 // Create space for the missing elements 238 OCTAVE_QUIT;
239
240 // Create space for the missing elements.
241
234 for (size_t i = 0; i < nrp; i++) 242 for (size_t i = 0; i < nrp; i++)
235 { 243 {
236 ptr1 = out + i * (nc/2 + 1) + nrp*((nc-1)/2); 244 ptr1 = out + i * (nc/2 + 1) + nrp*((nc-1)/2);
237 ptr2 = out + i * nc; 245 ptr2 = out + i * nc;
238 for (size_t j = 0; j < nc/2+1; j++) 246 for (size_t j = 0; j < nc/2+1; j++)
239 *ptr2++ = *ptr1++; 247 *ptr2++ = *ptr1++;
240 } 248 }
241 249
242 // Fill in the missing data for the rank = 2 case directly for speed 250 OCTAVE_QUIT;
251
252 // Fill in the missing data for the rank = 2 case directly for speed.
253
243 for (size_t i = 0; i < np; i++) 254 for (size_t i = 0; i < np; i++)
244 { 255 {
245 for (size_t j = 1; j < nr; j++) 256 for (size_t j = 1; j < nr; j++)
246 for (size_t k = nc/2+1; k < nc; k++) 257 for (size_t k = nc/2+1; k < nc; k++)
247 out[k + (j + i*nr)*nc] = conj(out[nc - k + ((i+1)*nr - j)*nc]); 258 out[k + (j + i*nr)*nc] = conj(out[nc - k + ((i+1)*nr - j)*nc]);
248 259
249 for (size_t j = nc/2+1; j < nc; j++) 260 for (size_t j = nc/2+1; j < nc; j++)
250 out[j + i*nr*nc] = conj(out[(i*nr+1)*nc - j]); 261 out[j + i*nr*nc] = conj(out[(i*nr+1)*nc - j]);
251 } 262 }
252 263
253 // Now do the permutations needed for rank > 2 cases 264 OCTAVE_QUIT;
265
266 // Now do the permutations needed for rank > 2 cases.
267
254 size_t jstart = dv(0) * dv(1); 268 size_t jstart = dv(0) * dv(1);
255 size_t kstep = dv(0); 269 size_t kstep = dv(0);
256 size_t nel = dv.numel (); 270 size_t nel = dv.numel ();
271
257 for (int inner = 2; inner < dv.length(); inner++) 272 for (int inner = 2; inner < dv.length(); inner++)
258 { 273 {
259 size_t jmax = jstart * dv(inner); 274 size_t jmax = jstart * dv(inner);
260 for (size_t i = 0; i < nel; i+=jmax) 275 for (size_t i = 0; i < nel; i+=jmax)
261 for (size_t j = jstart, jj = jmax-jstart; j < jj; 276 for (size_t j = jstart, jj = jmax-jstart; j < jj;
267 out[i + j + k + l] = out[i + jj + k + l]; 282 out[i + j + k + l] = out[i + jj + k + l];
268 out[i + jj + k + l] = tmp; 283 out[i + jj + k + l] = tmp;
269 } 284 }
270 jstart = jmax; 285 jstart = jmax;
271 } 286 }
287
288 OCTAVE_QUIT;
272 } 289 }
273 290
274 int 291 int
275 octave_fftw::fft (const double *in, Complex *out, size_t npts, 292 octave_fftw::fft (const double *in, Complex *out, size_t npts,
276 size_t nsamples, int stride, int dist) 293 size_t nsamples, int stride, int dist)