Andreas Specovius / gammalib

Commit 79d739e69d3d817317458905a124a3aa416628c4

Authored by Jürgen Knödlseder
1 parent a103d9db

Use pointers instead of std::vector in GFft class

Inline Side-by-side
Showing 3 changed files with 558 additions and 1092 deletions
include/GFft.hpp
  View file @79d739e
... ... @@ -80,7 +80,7 @@ public:
80 80 void clear(void);
81 81 GFft* clone(void) const;
82 82 std::string classname(void) const;
83   - int dim() const;
  83 + int dim(void) const;
84 84 int size(void) const;
85 85 const std::vector<int>& shape(void) const;
86 86 const std::vector<int>& strides(void) const;
... ... @@ -96,81 +96,89 @@ protected:
96 96 void set_data(const GNdarray& array);
97 97  
98 98 // Low-level FFT methods
99   - void transform(std::vector<std::complex<double> >& data,
100   - const int& stride,
101   - const int& n,
102   - const GFftWavetable& wavetable,
103   - const bool& forward = true);
104   - void factor2(const std::vector<std::complex<double> >& in,
105   - const int& istride,
106   - std::vector<std::complex<double> >& out,
107   - const int& ostride,
108   - const GFftWavetable& wavetable,
109   - const int& sign,
110   - const int& product,
111   - const int& n,
112   - const int& index);
113   - void factor3(const std::vector<std::complex<double> >& in,
114   - const int& istride,
115   - std::vector<std::complex<double> >& out,
116   - const int& ostride,
117   - const GFftWavetable& wavetable,
118   - const int& sign,
119   - const int& product,
120   - const int& n,
121   - const int& index);
122   - void factor4(const std::vector<std::complex<double> >& in,
123   - const int& istride,
124   - std::vector<std::complex<double> >& out,
125   - const int& ostride,
126   - const GFftWavetable& wavetable,
127   - const int& sign,
128   - const int& product,
129   - const int& n,
130   - const int& index);
131   - void factor5(const std::vector<std::complex<double> >& in,
132   - const int& istride,
133   - std::vector<std::complex<double> >& out,
134   - const int& ostride,
135   - const GFftWavetable& wavetable,
136   - const int& sign,
137   - const int& product,
138   - const int& n,
139   - const int& index);
140   - void factor6(const std::vector<std::complex<double> >& in,
141   - const int& istride,
142   - std::vector<std::complex<double> >& out,
143   - const int& ostride,
144   - const GFftWavetable& wavetable,
145   - const int& sign,
146   - const int& product,
147   - const int& n,
148   - const int& index);
149   - void factor7(const std::vector<std::complex<double> >& in,
150   - const int& istride,
151   - std::vector<std::complex<double> >& out,
152   - const int& ostride,
153   - const GFftWavetable& wavetable,
154   - const int& sign,
155   - const int& product,
156   - const int& n,
157   - const int& index);
158   - void factorn(std::vector<std::complex<double> >& in,
159   - const int& istride,
160   - std::vector<std::complex<double> >& out,
161   - const int& ostride,
162   - const GFftWavetable& wavetable,
163   - const int& sign,
164   - const int& factor,
165   - const int& product,
166   - const int& n,
167   - const int& index);
  99 + void transform(std::complex<double>* data,
  100 + const int& stride,
  101 + const int& n,
  102 + const GFftWavetable& wavetable,
  103 + const bool& forward = true);
  104 + void factor2(const std::complex<double>* in,
  105 + const int& istride,
  106 + std::complex<double>* out,
  107 + const int& ostride,
  108 + const GFftWavetable& wavetable,
  109 + const int& sign,
  110 + const int& product,
  111 + const int& n,
  112 + const int& index);
  113 + void factor3(const std::complex<double>* in,
  114 + const int& istride,
  115 + std::complex<double>* out,
  116 + const int& ostride,
  117 + const GFftWavetable& wavetable,
  118 + const int& sign,
  119 + const int& product,
  120 + const int& n,
  121 + const int& index);
  122 + void factor4(const std::complex<double>* in,
  123 + const int& istride,
  124 + std::complex<double>* out,
  125 + const int& ostride,
  126 + const GFftWavetable& wavetable,
  127 + const int& sign,
  128 + const int& product,
  129 + const int& n,
  130 + const int& index);
  131 + void factor5(const std::complex<double>* in,
  132 + const int& istride,
  133 + std::complex<double>* out,
  134 + const int& ostride,
  135 + const GFftWavetable& wavetable,
  136 + const int& sign,
  137 + const int& product,
  138 + const int& n,
  139 + const int& index);
  140 + void factor6(const std::complex<double>* in,
  141 + const int& istride,
  142 + std::complex<double>* out,
  143 + const int& ostride,
  144 + const GFftWavetable& wavetable,
  145 + const int& sign,
  146 + const int& product,
  147 + const int& n,
  148 + const int& index);
  149 + void factor7(const std::complex<double>* in,
  150 + const int& istride,
  151 + std::complex<double>* out,
  152 + const int& ostride,
  153 + const GFftWavetable& wavetable,
  154 + const int& sign,
  155 + const int& product,
  156 + const int& n,
  157 + const int& index);
  158 + void factorn(std::complex<double>* in,
  159 + const int& istride,
  160 + std::complex<double>* out,
  161 + const int& ostride,
  162 + const GFftWavetable& wavetable,
  163 + const int& sign,
  164 + const int& factor,
  165 + const int& product,
  166 + const int& n,
  167 + const int& index);
  168 + std::vector<std::complex<double> > get_w(const GFftWavetable& wavetable,
  169 + const int& index,
  170 + const int& k,
  171 + const int& q,
  172 + const int& n,
  173 + const int& sign) const;
  174 + std::complex<double> timesi(const std::complex<double>& value) const;
168 175  
169 176 // Protected members
170   - std::vector<int> m_shape; //!< Array dimensions
171   - std::vector<int> m_strides; //!< Steps in each dimension
172   - std::vector<std::complex<double> > m_data; //!< Array data
173   - std::vector<GFftWavetable> m_wavetable; //!< Trigonometric coefficients
  177 + int m_size; //!< Size of data array
  178 + std::complex<double>* m_data; //!< Pointer on array data
  179 + std::vector<int> m_shape; //!< Array dimensions
  180 + std::vector<int> m_strides; //!< Steps in each dimension
  181 + std::vector<GFftWavetable> m_wavetable; //!< Trigonometric coefficients
174 182 };
175 183  
176 184  
... ... @@ -189,28 +197,28 @@ std::string GFft::classname(void) const
189 197 /***********************************************************************//**
190 198 * @brief 1-dimensional FFT element access operator
191 199 *
192   - * @param[in] ix Element index [0,...,shape(0)-1].
  200 + * @param[in] ix Element index [0,...,size()-1].
193 201 * @return Reference to FFT element.
194 202 ***************************************************************************/
195 203 inline
196 204 std::complex<double>& GFft::operator()(const int& ix)
197 205 {
198 206 // Return array element
199   - return m_data[ix];
  207 + return *(m_data + ix);
200 208 }
201 209  
202 210  
203 211 /***********************************************************************//**
204 212 * @brief 1-dimensional FFT element access operator (const variant)
205 213 *
206   - * @param[in] ix Element index [0,...,shape(0)-1].
  214 + * @param[in] ix Element index [0,...,size()-1].
207 215 * @return Const reference to FFT element.
208 216 ***************************************************************************/
209 217 inline
210 218 const std::complex<double>& GFft::operator()(const int& ix) const
211 219 {
212 220 // Return array element
213   - return m_data[ix];
  221 + return *(m_data + ix);
214 222 }
215 223  
216 224  
... ... @@ -222,7 +230,7 @@ const std::complex&lt;double&gt;&amp; GFft::operator()(const int&amp; ix) const
222 230 * Returns the dimension of the Fast Fourier Transformation.
223 231 ***************************************************************************/
224 232 inline
225   -int GFft::dim() const
  233 +int GFft::dim(void) const
226 234 {
227 235 return m_shape.size();
228 236 }
... ... @@ -236,9 +244,9 @@ int GFft::dim() const
236 244 * Returns the number of elements in the array.
237 245 ***************************************************************************/
238 246 inline
239   -int GFft::size() const
  247 +int GFft::size(void) const
240 248 {
241   - return m_data.size();
  249 + return m_size;
242 250 }
243 251  
244 252  
... ... @@ -341,4 +349,20 @@ GFft operator/(const GFft&amp; a, const GFft&amp; b)
341 349 return result;
342 350 }
343 351  
  352 +
  353 +/***********************************************************************//**
  354 + * @brief Return complex value times i
  355 + *
  356 + * @param[in] value Complex value.
  357 + * @return Complex value times i.
  358 + *
  359 + * Returns complex value times i.
  360 + ***************************************************************************/
  361 +inline
  362 +std::complex<double> GFft::timesi(const std::complex<double>& value) const
  363 +{
  364 + std::complex<double> result(-value.imag(), value.real());
  365 + return result;
  366 +}
  367 +
344 368 #endif /* GFFT_HPP */
... ...
include/GNdarray.hpp
  View file @79d739e
... ... @@ -102,7 +102,7 @@ public:
102 102 void clear(void);
103 103 GNdarray* clone(void) const;
104 104 std::string classname(void) const;
105   - int dim() const;
  105 + int dim(void) const;
106 106 int size(void) const;
107 107 const std::vector<int>& shape(void) const;
108 108 const std::vector<int>& strides(void) const;
... ... @@ -270,7 +270,7 @@ const double&amp; GNdarray::operator()(const std::vector&lt;int&gt;&amp; i) const
270 270 * Returns the dimension of the array.
271 271 ***************************************************************************/
272 272 inline
273   -int GNdarray::dim() const
  273 +int GNdarray::dim(void) const
274 274 {
275 275 return m_shape.size();
276 276 }
... ... @@ -284,7 +284,7 @@ int GNdarray::dim() const
284 284 * Returns the number of elements in the array.
285 285 ***************************************************************************/
286 286 inline
287   -int GNdarray::size() const
  287 +int GNdarray::size(void) const
288 288 {
289 289 return m_data.size();
290 290 }
... ...
src/numerics/GFft.cpp
  View file @79d739e
... ... @@ -229,8 +229,8 @@ GFft GFft::operator-(void) const
229 229 GFft fft = *this;
230 230  
231 231 // Negate all elements
232   - for (int i = 0; i < m_data.size(); ++i) {
233   - fft.m_data[i] = -fft.m_data[i];
  232 + for (int i = 0; i < fft.size(); ++i) {
  233 + fft(i) = -fft(i);
234 234 }
235 235  
236 236 // Return FFT
... ... @@ -285,7 +285,24 @@ void GFft::forward(const GNdarray&amp; array)
285 285 set_data(array);
286 286  
287 287 // Continue only if there are elements in the data member
288   - if (m_data.size() > 0) {
  288 + if (m_size > 0) {
  289 +/*
  290 + // Loop over all dimensions of array
  291 + for (int d = 0; d < dim(); ++d) {
  292 +
  293 + // Loop over all dimensions of array
  294 + for (int i = 0; i < dim(); ++i) {
  295 +
  296 + // Skip
  297 + if (d == i) {
  298 + continue;
  299 + }
  300 +
  301 + // Loop over all
  302 + int n = shape()[d];
  303 + int stride = strides()[d];
  304 + int istart =
  305 +*/
289 306  
290 307 // Perform FFT (is only correct for 1D for now)
291 308 transform(m_data, 1, m_shape[0], m_wavetable[0], true);
... ... @@ -313,7 +330,10 @@ GNdarray GFft::backward(void) const
313 330 if (array.size() > 0) {
314 331  
315 332 // Create copy of complex array
316   - std::vector<std::complex<double> > data = m_data;
  333 + std::complex<double>* data = new std::complex<double>[array.size()];
  334 + for (int i = 0; i < array.size(); ++i) {
  335 + *(data+i) = *(m_data+i);
  336 + }
317 337  
318 338 // Perform FFT (circumvent const correctness)
319 339 const_cast<GFft*>(this)->transform(data, 1, m_shape[0], m_wavetable[0],
... ... @@ -324,9 +344,12 @@ GNdarray GFft::backward(void) const
324 344  
325 345 // Extract N-dimensional array and normalise with 1/n
326 346 for (int i = 0; i < array.size(); ++i) {
327   - array(i) = data[i].real() * norm;
  347 + array(i) = (data+i)->real() * norm;
328 348 }
329 349  
  350 + // Delete copy of complex array
  351 + delete [] data;
  352 +
330 353 } // endif: there were elements in the array
331 354  
332 355 // Return N-dimensional array
... ... @@ -364,16 +387,16 @@ std::string GFft::print(const GChatter&amp; chatter) const
364 387 }
365 388 result.append(")");
366 389 result.append("\n"+gammalib::parformat("Size"));
367   - result.append(gammalib::str(size()));
  390 + result.append(gammalib::str(m_size));
368 391  
369 392 // VERBOSE: Put all FFT elements in string
370 393 if (chatter >= VERBOSE) {
371 394 result.append("\n"+gammalib::parformat("Elements"));
372   - for (int i = 0; i < m_data.size(); ++i) {
  395 + for (int i = 0; i < size(); ++i) {
373 396 if (i > 0) {
374 397 result.append(",");
375 398 }
376   - result.append(gammalib::str(m_data[i]));
  399 + result.append(gammalib::str((*this)(i)));
377 400 }
378 401 result.append(")");
379 402 }
... ... @@ -397,9 +420,10 @@ std::string GFft::print(const GChatter&amp; chatter) const
397 420 void GFft::init_members(void)
398 421 {
399 422 // Initialise members
  423 + m_size = 0;
  424 + m_data = NULL;
400 425 m_shape.clear();
401 426 m_strides.clear();
402   - m_data.clear();
403 427 m_wavetable.clear();
404 428  
405 429 // Return
... ... @@ -415,11 +439,19 @@ void GFft::init_members(void)
415 439 void GFft::copy_members(const GFft& fft)
416 440 {
417 441 // Copy members
  442 + m_size = fft.m_size;
418 443 m_shape = fft.m_shape;
419 444 m_strides = fft.m_strides;
420   - m_data = fft.m_data;
421 445 m_wavetable = fft.m_wavetable;
422 446  
  447 + // Copy data
  448 + if (m_size > 0) {
  449 + m_data = new std::complex<double>[m_size];
  450 + for (int i = 0; i < m_size; ++i) {
  451 + *(m_data+i) = *(fft.m_data+i);
  452 + }
  453 + }
  454 +
423 455 // Return
424 456 return;
425 457 }
... ... @@ -430,6 +462,15 @@ void GFft::copy_members(const GFft&amp; fft)
430 462 ***************************************************************************/
431 463 void GFft::free_members(void)
432 464 {
  465 + // Delete data
  466 + if (m_data != NULL) {
  467 + delete [] m_data;
  468 + }
  469 +
  470 + // Reset size and pointer
  471 + m_size = 0;
  472 + m_data = NULL;
  473 +
433 474 // Return
434 475 return;
435 476 }
... ... @@ -443,21 +484,21 @@ void GFft::free_members(void)
443 484 void GFft::set_data(const GNdarray& array)
444 485 {
445 486 // Clear information
446   - m_shape.clear();
447   - m_strides.clear();
448   - m_data.clear();
449   - m_wavetable.clear();
  487 + free_members();
  488 + init_members();
450 489  
451 490 // Continue only if there are elements in the array
452 491 if (array.size() > 0) {
453 492  
  493 + // Set size of array
  494 + m_size = array.size();
  495 +
454 496 // Allocate data and initialise data to zero
455   - std::complex<double> zero(0.0, 0.0);
456   - m_data.assign(array.size(), zero);
  497 + m_data = new std::complex<double>[m_size];
457 498  
458 499 // Set real part of data from n-dimensional array
459   - for (int i = 0; i < array.size(); ++i) {
460   - m_data[i] = std::complex<double>(array(i), 0.0);
  500 + for (int i = 0; i < m_size; ++i) {
  501 + *(m_data+i) = std::complex<double>(array(i), 0.0);
461 502 }
462 503  
463 504 // Save shape
... ... @@ -486,21 +527,24 @@ void GFft::set_data(const GNdarray&amp; array)
486 527 /***********************************************************************//**
487 528 * @brief Perform Fast Fourier Transform
488 529 *
489   - * @param[in] data FFT array.
490   - * @param[in] stride Step size when traversing FFT array in one dimension.
491   - * @param[in] n Number of elements in the dimension.
  530 + * @param[in] data Pointer to complex array to be transformed.
  531 + * @param[in] stride Step size when traversing complex array.
  532 + * @param[in] n Number of elements in complex array.
492 533 * @param[in] wavetable Trigonometric coefficients.
493 534 * @param[in] forward Forward transform if true, otherwise backward transform.
494 535 ***************************************************************************/
495   -void GFft::transform(std::vector<std::complex<double> >& data,
496   - const int& stride,
497   - const int& n,
498   - const GFftWavetable& wavetable,
499   - const bool& forward)
  536 +void GFft::transform(std::complex<double>* data,
  537 + const int& stride,
  538 + const int& n,
  539 + const GFftWavetable& wavetable,
  540 + const bool& forward)
500 541 {
501   - // Allocate scratch space
502   - std::complex<double> zero(0.0, 0.0);
503   - std::vector<std::complex<double> > scratch(n, zero);
  542 + // Allocate scratch space and initialize to zero
  543 + std::complex<double>* scratch = new std::complex<double>[n];
  544 + std::complex<double> zero(0.0, 0.0);
  545 + for (int i = 0; i < n; ++i) {
  546 + *(scratch+i) = zero;
  547 + }
504 548  
505 549 // Initialise transformation state
506 550 int state = 0;
... ... @@ -525,46 +569,46 @@ void GFft::transform(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; data,
525 569 int q = n / product;
526 570  
527 571 // Set state dependent stuff
528   - std::vector<std::complex<double> >* in = &data;
529   - std::vector<std::complex<double> >* out = &scratch;
530   - int istride;
531   - int ostride;
  572 + std::complex<double>* in;
  573 + std::complex<double>* out;
  574 + int istride;
  575 + int ostride;
532 576 if (state == 0) {
533   - in = &data;
  577 + in = data;
534 578 istride = stride;
535   - out = &scratch;
  579 + out = scratch;
536 580 ostride = 1;
537 581 state = 1;
538 582 }
539 583 else {
540   - in = &scratch;
  584 + in = scratch;
541 585 istride = 1;
542   - out = &data;
  586 + out = data;
543 587 ostride = stride;
544 588 state = 0;
545 589 }
546 590  
547 591 // Call factor dependent method
548 592 if (factor == 2) {
549   - factor2(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  593 + factor2(in, istride, out, ostride, wavetable, sign, product, n, index);
550 594 }
551 595 else if (factor == 3) {
552   - factor3(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  596 + factor3(in, istride, out, ostride, wavetable, sign, product, n, index);
553 597 }
554 598 else if (factor == 4) {
555   - factor4(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  599 + factor4(in, istride, out, ostride, wavetable, sign, product, n, index);
556 600 }
557 601 else if (factor == 5) {
558   - factor5(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  602 + factor5(in, istride, out, ostride, wavetable, sign, product, n, index);
559 603 }
560 604 else if (factor == 6) {
561   - factor6(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  605 + factor6(in, istride, out, ostride, wavetable, sign, product, n, index);
562 606 }
563 607 else if (factor == 7) {
564   - factor7(*in, istride, *out, ostride, wavetable, sign, product, n, index);
  608 + factor7(in, istride, out, ostride, wavetable, sign, product, n, index);
565 609 }
566 610 else {
567   - factorn(*in, istride, *out, ostride, wavetable, sign, factor, product, n, index);
  611 + factorn(in, istride, out, ostride, wavetable, sign, factor, product, n, index);
568 612 }
569 613  
570 614 } // endfor: looped over all factors
... ... @@ -574,10 +618,13 @@ void GFft::transform(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; data,
574 618 // the data array
575 619 if (state == 1) {
576 620 for (int i = 0; i < n; ++i) {
577   - data[stride*i] = scratch[i];
  621 + *(data+stride*i) = *(scratch+i);
578 622 }
579 623 }
580 624  
  625 + // Free scratch space
  626 + delete [] scratch;
  627 +
581 628 // Return
582 629 return;
583 630 }
... ... @@ -586,9 +633,9 @@ void GFft::transform(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; data,
586 633 /***********************************************************************//**
587 634 * @brief Compute FFT for factor 2
588 635 *
589   - * @param[in] in Input array.
  636 + * @param[in] in Pointer to input array.
590 637 * @param[in] istride Step size when traversing input array.
591   - * @param[in] out Output array.
  638 + * @param[in] out Pointer to output array.
592 639 * @param[in] ostride Step size when traversing output array.
593 640 * @param[in] wavetable Trigonometric coefficients.
594 641 * @param[in] sign Forward (-1) or backward (+1).
... ... @@ -596,15 +643,15 @@ void GFft::transform(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; data,
596 643 * @param[in] n Logical array length.
597 644 * @param[in] index Start index of trigonometric coefficients.
598 645 ***************************************************************************/
599   -void GFft::factor2(const std::vector<std::complex<double> >& in,
600   - const int& istride,
601   - std::vector<std::complex<double> >& out,
602   - const int& ostride,
603   - const GFftWavetable& wavetable,
604   - const int& sign,
605   - const int& product,
606   - const int& n,
607   - const int& index)
  646 +void GFft::factor2(const std::complex<double>* in,
  647 + const int& istride,
  648 + std::complex<double>* out,
  649 + const int& ostride,
  650 + const GFftWavetable& wavetable,
  651 + const int& sign,
  652 + const int& product,
  653 + const int& n,
  654 + const int& index)
608 655 {
609 656 // Compute ...
610 657 const int factor = 2;
... ... @@ -616,63 +663,23 @@ void GFft::factor2(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
616 663 // Loop over ...
617 664 for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
618 665  
619   - // Allocate variables
620   - double w1_real, w1_imag;
621   -
622   - // Set trigonometric coefficients for k=0 since they are not stored
623   - // in the wavetable object
624   - if (k == 0) {
625   - w1_real = 1.0;
626   - w1_imag = 0.0;
627   - }
628   - else {
629   -
630   - // Compute indices
631   - int twiddle = index + k - 1;
632   -
633   - // Set trigonometric coefficients for forward transform
634   - if (sign == -1) {
635   - w1_real = wavetable[twiddle].real();
636   - w1_imag = wavetable[twiddle].imag();
637   - }
638   -
639   - // ... otherwise set trigonometric coefficients for backward
640   - // tranform: w -> conjugate(w)
641   - else {
642   - w1_real = wavetable[twiddle].real();
643   - w1_imag = -wavetable[twiddle].imag();
644   - }
645   -
646   - } // endelse: k != 0
  666 + // Extract coefficients from wavetable
  667 + std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 1, sign);
647 668  
648   - // Loop over ...
  669 + // Compute x = W(2) z
649 670 for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
650 671  
651   - const double z0_real = in[istride*i].real();
652   - const double z0_imag = in[istride*i].imag();
653   - const double z1_real = in[istride*(i+m)].real();
654   - const double z1_imag = in[istride*(i+m)].imag();
655   -
656   - // compute x = W(2) z
657   -
658   - // x0 = z0 + z1
659   - const double x0_real = z0_real + z1_real;
660   - const double x0_imag = z0_imag + z1_imag;
  672 + // Get z
  673 + const std::complex<double> z0 = *(in+istride*i);
  674 + const std::complex<double> z1 = *(in+istride*(i+m));
661 675  
662   - // x1 = z0 - z1
663   - const double x1_real = z0_real - z1_real;
664   - const double x1_imag = z0_imag - z1_imag;
  676 + // Compute x
  677 + const std::complex<double> x0 = z0 + z1;
  678 + const std::complex<double> x1 = z0 - z1;
665 679  
666   - // out0 = 1 * x0
667   - //out[ostride*j].real(x0_real);
668   - //out[ostride*j].imag(x0_imag);
669   - out[ostride*j] = std::complex<double>(x0_real, x0_imag);
670   -
671   - // out1 = w1 * x1
672   - //out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
673   - //out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
674   - out[ostride*(j+p_1)] = std::complex<double>(w1_real * x1_real - w1_imag * x1_imag,
675   - w1_real * x1_imag + w1_imag * x1_real);
  680 + // out = w * x
  681 + *(out+ostride*j) = x0;
  682 + *(out+ostride*(j+p_1)) = w[0] * x1;
676 683  
677 684 } // endfor: k1
678 685  
... ... @@ -686,9 +693,9 @@ void GFft::factor2(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
686 693 /***********************************************************************//**
687 694 * @brief Compute FFT for factor 3
688 695 *
689   - * @param[in] in Input array.
  696 + * @param[in] in Pointer to input array.
690 697 * @param[in] istride Step size when traversing input array.
691   - * @param[in] out Output array.
  698 + * @param[in] out Pointer to output array.
692 699 * @param[in] ostride Step size when traversing output array.
693 700 * @param[in] wavetable Trigonometric coefficients.
694 701 * @param[in] sign Forward (-1) or backward (+1).
... ... @@ -696,15 +703,15 @@ void GFft::factor2(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
696 703 * @param[in] n Logical array length.
697 704 * @param[in] index Start index of trigonometric coefficients.
698 705 ***************************************************************************/
699   -void GFft::factor3(const std::vector<std::complex<double> >& in,
700   - const int& istride,
701   - std::vector<std::complex<double> >& out,
702   - const int& ostride,
703   - const GFftWavetable& wavetable,
704   - const int& sign,
705   - const int& product,
706   - const int& n,
707   - const int& index)
  706 +void GFft::factor3(const std::complex<double>* in,
  707 + const int& istride,
  708 + std::complex<double>* out,
  709 + const int& ostride,
  710 + const GFftWavetable& wavetable,
  711 + const int& sign,
  712 + const int& product,
  713 + const int& n,
  714 + const int& index)
708 715 {
709 716 // Compute ...
710 717 const int factor = 3;
... ... @@ -719,95 +726,32 @@ void GFft::factor3(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
719 726  
720 727 // Loop over ...
721 728 for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
722   -
723   - // Allocate variables
724   - double w1_real, w1_imag, w2_real, w2_imag;
725   -
726   - // Set trigonometric coefficients for k=0 since they are not stored
727   - // in the wavetable object
728   - if (k == 0) {
729   - w1_real = 1.0;
730   - w1_imag = 0.0;
731   - w2_real = 1.0;
732   - w2_imag = 0.0;
733   - }
734   - else {
735   -
736   - // Compute indices
737   - int twiddle1 = index + k - 1;
738   - int twiddle2 = twiddle1 + q;
739   -
740   - // Set trigonometric coefficients for forward transform
741   - if (sign == -1) {
742   - w1_real = wavetable[twiddle1].real();
743   - w1_imag = wavetable[twiddle1].imag();
744   - w2_real = wavetable[twiddle2].real();
745   - w2_imag = wavetable[twiddle2].imag();
746   - }
747 729  
748   - // ... otherwise set trigonometric coefficients for backward
749   - // tranform: w -> conjugate(w)
750   - else {
751   - w1_real = wavetable[twiddle1].real();
752   - w1_imag = -wavetable[twiddle1].imag();
753   - w2_real = wavetable[twiddle2].real();
754   - w2_imag = -wavetable[twiddle2].imag();
755   - }
756   - } // endelse: k != 0
  730 + // Extract coefficients from wavetable
  731 + std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 2, sign);
757 732  
758   - // ...
  733 + // Compute x = W(3) z
759 734 for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
760 735  
761 736 // Get z
762   - const double z0_real = in[istride*i].real();
763   - const double z0_imag = in[istride*i].imag();
764   - const double z1_real = in[istride*(i+m)].real();
765   - const double z1_imag = in[istride*(i+m)].imag();
766   - const double z2_real = in[istride*(i+2*m)].real();
767   - const double z2_imag = in[istride*(i+2*m)].imag();
768   -
769   - // compute x = W(3) z
770   -
771   - // t1 = z1 + z2
772   - const double t1_real = z1_real + z2_real;
773   - const double t1_imag = z1_imag + z2_imag;
774   -
775   - // t2 = z0 - t1/2
776   - const double t2_real = z0_real - t1_real / 2.0;
777   - const double t2_imag = z0_imag - t1_imag / 2.0;
778   -
779   - // t3 = (+/-) sin(pi/3)*(z1 - z2)
780   - const double t3_real = ((int) sign) * tau * (z1_real - z2_real);
781   - const double t3_imag = ((int) sign) * tau * (z1_imag - z2_imag);
782   -
783   - // x0 = z0 + t1
784   - const double x0_real = z0_real + t1_real;
785   - const double x0_imag = z0_imag + t1_imag;
786   -
787   - // x1 = t2 + i t3
788   - const double x1_real = t2_real - t3_imag;
789   - const double x1_imag = t2_imag + t3_real;
790   -
791   - // x2 = t2 - i t3
792   - const double x2_real = t2_real + t3_imag;
793   - const double x2_imag = t2_imag - t3_real;
794   -
795   - // out0 = 1 * x0
796   - //out[ostride*j].real(x0_real);
797   - //out[ostride*j].imag(x0_imag);
798   - out[ostride*j] = std::complex<double>(x0_real, x0_imag);
799   -
800   - // out1 = w1 * x1
801   - //out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
802   - //out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
803   - out[ostride*(j+p_1)] = std::complex<double>(w1_real * x1_real - w1_imag * x1_imag,
804   - w1_real * x1_imag + w1_imag * x1_real);
805   -
806   - // out2 = w2 * x2
807   - //out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
808   - //out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
809   - out[ostride*(j+2*p_1)] = std::complex<double>(w2_real * x2_real - w2_imag * x2_imag,
810   - w2_real * x2_imag + w2_imag * x2_real);
  737 + const std::complex<double> z0 = *(in+istride*i);
  738 + const std::complex<double> z1 = *(in+istride*(i+m));
  739 + const std::complex<double> z2 = *(in+istride*(i+2*m));
  740 +
  741 + // Compute t
  742 + const std::complex<double> t1 = z1 + z2;
  743 + const std::complex<double> t2 = z0 - t1/2.0;
  744 + const std::complex<double> t3 = double((int)sign) * tau * (z1 - z2);
  745 +
  746 + // Compute x
  747 + const std::complex<double> x0 = z0 + t1;
  748 + const std::complex<double> x1 = t2 + timesi(t3);
  749 + const std::complex<double> x2 = t2 - timesi(t3);
  750 +
  751 + // out = w * x
  752 + *(out+ostride*j) = x0;
  753 + *(out+ostride*(j+p_1)) = w[0] * x1;
  754 + *(out+ostride*(j+2*p_1)) = w[1] * x2;
811 755  
812 756 } // endfor: k1
813 757  
... ... @@ -821,9 +765,9 @@ void GFft::factor3(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
821 765 /***********************************************************************//**
822 766 * @brief Compute FFT for factor 4
823 767 *
824   - * @param[in] in Input array.
  768 + * @param[in] in Pointer to input array.
825 769 * @param[in] istride Step size when traversing input array.
826   - * @param[in] out Output array.
  770 + * @param[in] out Pointer to output array.
827 771 * @param[in] ostride Step size when traversing output array.
828 772 * @param[in] wavetable Trigonometric coefficients.
829 773 * @param[in] sign Forward (-1) or backward (+1).
... ... @@ -831,15 +775,15 @@ void GFft::factor3(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
831 775 * @param[in] n Logical array length.
832 776 * @param[in] index Start index of trigonometric coefficients.
833 777 ***************************************************************************/
834   -void GFft::factor4(const std::vector<std::complex<double> >& in,
835   - const int& istride,
836   - std::vector<std::complex<double> >& out,
837   - const int& ostride,
838   - const GFftWavetable& wavetable,
839   - const int& sign,
840   - const int& product,
841   - const int& n,
842   - const int& index)
  778 +void GFft::factor4(const std::complex<double>* in,
  779 + const int& istride,
  780 + std::complex<double>* out,
  781 + const int& ostride,
  782 + const GFftWavetable& wavetable,
  783 + const int& sign,
  784 + const int& product,
  785 + const int& n,
  786 + const int& index)
843 787 {
844 788 // Compute ...
845 789 const int factor = 4;
... ... @@ -850,121 +794,36 @@ void GFft::factor4(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
850 794  
851 795 // Loop over ...
852 796 for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
853   -
854   - // Allocate w
855   - double w1_real, w1_imag, w2_real, w2_imag, w3_real, w3_imag;
856   -
857   - // Set trigonometric coefficients for k=0 since they are not stored
858   - // in the wavetable object
859   - if (k == 0) {
860   - w1_real = 1.0;
861   - w1_imag = 0.0;
862   - w2_real = 1.0;
863   - w2_imag = 0.0;
864   - w3_real = 1.0;
865   - w3_imag = 0.0;
866   - }
867   - else
868   - {
869   -
870   - // Compute indices
871   - int twiddle1 = index + k - 1;
872   - int twiddle2 = twiddle1 + q;
873   - int twiddle3 = twiddle2 + q;
874   -
875   - // Set trigonometric coefficients for forward transform
876   - if (sign == -1) {
877   - w1_real = wavetable[twiddle1].real();
878   - w1_imag = wavetable[twiddle1].imag();
879   - w2_real = wavetable[twiddle2].real();
880   - w2_imag = wavetable[twiddle2].imag();
881   - w3_real = wavetable[twiddle3].real();
882   - w3_imag = wavetable[twiddle3].imag();
883   - }
884   -
885   - // ... otherwise set trigonometric coefficients for backward
886   - // tranform: w -> conjugate(w)
887   - else {
888   - w1_real = wavetable[twiddle1].real();
889   - w1_imag = -wavetable[twiddle1].imag();
890   - w2_real = wavetable[twiddle2].real();
891   - w2_imag = -wavetable[twiddle2].imag();
892   - w3_real = wavetable[twiddle3].real();
893   - w3_imag = -wavetable[twiddle3].imag();
894   - }
895 797  
896   - } // endelse: k != 0
  798 + // Extract coefficients from wavetable
  799 + std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 3, sign);
897 800  
898   - // Loop over ...
  801 + // Compute x = W(4) z
899 802 for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
900 803  
901   -
902 804 // Get z
903   - const double z0_real = in[istride*i].real();
904   - const double z0_imag = in[istride*i].imag();
905   - const double z1_real = in[istride*(i+m)].real();
906   - const double z1_imag = in[istride*(i+m)].imag();
907   - const double z2_real = in[istride*(i+2*m)].real();
908   - const double z2_imag = in[istride*(i+2*m)].imag();
909   - const double z3_real = in[istride*(i+3*m)].real();
910   - const double z3_imag = in[istride*(i+3*m)].imag();
911   -
912   - // compute x = W(4) z
  805 + const std::complex<double> z0 = *(in+istride*i);
  806 + const std::complex<double> z1 = *(in+istride*(i+m));
  807 + const std::complex<double> z2 = *(in+istride*(i+2*m));
  808 + const std::complex<double> z3 = *(in+istride*(i+3*m));
913 809  
914   - // t1 = z0 + z2
915   - const double t1_real = z0_real + z2_real;
916   - const double t1_imag = z0_imag + z2_imag;
917   -
918   - // t2 = z1 + z3
919   - const double t2_real = z1_real + z3_real;
920   - const double t2_imag = z1_imag + z3_imag;
921   -
922   - // t3 = z0 - z2
923   - const double t3_real = z0_real - z2_real;
924   - const double t3_imag = z0_imag - z2_imag;
925   -
926   - // t4 = (+/-) (z1 - z3)
927   - const double t4_real = ((int) sign) * (z1_real - z3_real);
928   - const double t4_imag = ((int) sign) * (z1_imag - z3_imag);
929   -
930   - // x0 = t1 + t2
931   - const double x0_real = t1_real + t2_real;
932   - const double x0_imag = t1_imag + t2_imag;
933   -
934   - // x1 = t3 + i t4
935   - const double x1_real = t3_real - t4_imag;
936   - const double x1_imag = t3_imag + t4_real;
937   -
938   - // x2 = t1 - t2
939   - const double x2_real = t1_real - t2_real;
940   - const double x2_imag = t1_imag - t2_imag;
941   -
942   - // x3 = t3 - i t4
943   - const double x3_real = t3_real + t4_imag;
944   - const double x3_imag = t3_imag - t4_real;
945   -
946   - // out0 = 1 * x0
947   - //out[ostride*j].real(x0_real);
948   - //out[ostride*j].imag(x0_imag);
949   - out[ostride*j] = std::complex<double>(x0_real, x0_imag);
950   -
951   - // out1 = w1 * x1
952   - //out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
953   - //out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
954   - out[ostride*(j+p_1)] = std::complex<double>(w1_real * x1_real - w1_imag * x1_imag,
955   - w1_real * x1_imag + w1_imag * x1_real);
956   -
957   - // out2 = w2 * x2
958   - //out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
959   - //out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
960   - out[ostride*(j+2*p_1)] = std::complex<double>(w2_real * x2_real - w2_imag * x2_imag,
961   - w2_real * x2_imag + w2_imag * x2_real);
962   -
963   - // out3 = w3 * x3
964   - //out[ostride*(j+3*p_1)].real(w3_real * x3_real - w3_imag * x3_imag);
965   - //out[ostride*(j+3*p_1)].imag(w3_real * x3_imag + w3_imag * x3_real);
966   - out[ostride*(j+3*p_1)] = std::complex<double>(w3_real * x3_real - w3_imag * x3_imag,
967   - w3_real * x3_imag + w3_imag * x3_real);
  810 + // Compute t
  811 + const std::complex<double> t1 = z0 + z2;
  812 + const std::complex<double> t2 = z1 + z3;
  813 + const std::complex<double> t3 = z0 - z2;
  814 + const std::complex<double> t4 = double((int)sign) * (z1 - z3);
  815 +
  816 + // Compute x
  817 + const std::complex<double> x0 = t1 + t2;
  818 + const std::complex<double> x1 = t3 + timesi(t4);
  819 + const std::complex<double> x2 = t1 - t2;
  820 + const std::complex<double> x3 = t3 - timesi(t4);
  821 +
  822 + // out = w * x
  823 + *(out+ostride*j) = x0;
  824 + *(out+ostride*(j+p_1)) = w[0] * x1;
  825 + *(out+ostride*(j+2*p_1)) = w[1] * x2;
  826 + *(out+ostride*(j+3*p_1)) = w[2] * x3;
968 827  
969 828 } // endfor: k1
970 829  
... ... @@ -978,9 +837,9 @@ void GFft::factor4(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
978 837 /***********************************************************************//**
979 838 * @brief Compute FFT for factor 5
980 839 *
981   - * @param[in] in Input array.
  840 + * @param[in] in Pointer to input array.
982 841 * @param[in] istride Step size when traversing input array.
983   - * @param[in] out Output array.
  842 + * @param[in] out Pointer to output array.
984 843 * @param[in] ostride Step size when traversing output array.
985 844 * @param[in] wavetable Trigonometric coefficients.
986 845 * @param[in] sign Forward (-1) or backward (+1).
... ... @@ -988,15 +847,15 @@ void GFft::factor4(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
988 847 * @param[in] n Logical array length.
989 848 * @param[in] index Start index of trigonometric coefficients.
990 849 ***************************************************************************/
991   -void GFft::factor5(const std::vector<std::complex<double> >& in,
992   - const int& istride,
993   - std::vector<std::complex<double> >& out,
994   - const int& ostride,
995   - const GFftWavetable& wavetable,
996   - const int& sign,
997   - const int& product,
998   - const int& n,
999   - const int& index)
  850 +void GFft::factor5(const std::complex<double>* in,
  851 + const int& istride,
  852 + std::complex<double>* out,
  853 + const int& ostride,
  854 + const GFftWavetable& wavetable,
  855 + const int& sign,
  856 + const int& product,
  857 + const int& n,
  858 + const int& index)
1000 859 {
1001 860 // Compute ...
1002 861 const int factor = 5;
... ... @@ -1013,170 +872,47 @@ void GFft::factor5(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1013 872 // Loop over ...
1014 873 for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
1015 874  
1016   - // Allocate w
1017   - double w1_real, w1_imag, w2_real, w2_imag, w3_real, w3_imag, w4_real,
1018   - w4_imag;
1019   -
1020   - // Set trigonometric coefficients for k=0 since they are not stored
1021   - // in the wavetable object
1022   - if (k == 0) {
1023   - w1_real = 1.0;
1024   - w1_imag = 0.0;
1025   - w2_real = 1.0;
1026   - w2_imag = 0.0;
1027   - w3_real = 1.0;
1028   - w3_imag = 0.0;
1029   - w4_real = 1.0;
1030   - w4_imag = 0.0;
1031   - }
1032   - else {
1033   -
1034   - // Compute indices
1035   - int twiddle1 = index + k - 1;
1036   - int twiddle2 = twiddle1 + q;
1037   - int twiddle3 = twiddle2 + q;
1038   - int twiddle4 = twiddle3 + q;
1039   -
1040   - // Set trigonometric coefficients for forward transform
1041   - if (sign == -1) {
1042   - w1_real = wavetable[twiddle1].real();
1043   - w1_imag = wavetable[twiddle1].imag();
1044   - w2_real = wavetable[twiddle2].real();
1045   - w2_imag = wavetable[twiddle2].imag();
1046   - w3_real = wavetable[twiddle3].real();
1047   - w3_imag = wavetable[twiddle3].imag();
1048   - w4_real = wavetable[twiddle4].real();
1049   - w4_imag = wavetable[twiddle4].imag();
1050   - }
1051   -
1052   - // ... otherwise set trigonometric coefficients for backward
1053   - // tranform: w -> conjugate(w)
1054   - else {
1055   - w1_real = wavetable[twiddle1].real();
1056   - w1_imag = -wavetable[twiddle1].imag();
1057   - w2_real = wavetable[twiddle2].real();
1058   - w2_imag = -wavetable[twiddle2].imag();
1059   - w3_real = wavetable[twiddle3].real();
1060   - w3_imag = -wavetable[twiddle3].imag();
1061   - w4_real = wavetable[twiddle4].real();
1062   - w4_imag = -wavetable[twiddle4].imag();
1063   - }
  875 + // Extract coefficients from wavetable
  876 + std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 4, sign);
1064 877  
1065   - } // endelse: k != 0
1066   -
1067   - // Loop over ...
  878 + // Compute x = W(5) z
1068 879 for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
1069 880  
1070 881 // Get z
1071   - const double z0_real = in[istride*i].real();
1072   - const double z0_imag = in[istride*i].imag();
1073   - const double z1_real = in[istride*(i+m)].real();
1074   - const double z1_imag = in[istride*(i+m)].imag();
1075   - const double z2_real = in[istride*(i+2*m)].real();
1076   - const double z2_imag = in[istride*(i+2*m)].imag();
1077   - const double z3_real = in[istride*(i+3*m)].real();
1078   - const double z3_imag = in[istride*(i+3*m)].imag();
1079   - const double z4_real = in[istride*(i+4*m)].real();
1080   - const double z4_imag = in[istride*(i+4*m)].imag();
1081   -
1082   - // compute x = W(5) z
1083   -
1084   - // t1 = z1 + z4
1085   - const double t1_real = z1_real + z4_real;
1086   - const double t1_imag = z1_imag + z4_imag;
1087   -
1088   - // t2 = z2 + z3
1089   - const double t2_real = z2_real + z3_real;
1090   - const double t2_imag = z2_imag + z3_imag;
1091   -
1092   - // t3 = z1 - z4
1093   - const double t3_real = z1_real - z4_real;
1094   - const double t3_imag = z1_imag - z4_imag;
1095   -
1096   - // t4 = z2 - z3
1097   - const double t4_real = z2_real - z3_real;
1098   - const double t4_imag = z2_imag - z3_imag;
1099   -
1100   - // t5 = t1 + t2
1101   - const double t5_real = t1_real + t2_real;
1102   - const double t5_imag = t1_imag + t2_imag;
1103   -
1104   - // t6 = (sqrt(5)/4)(t1 - t2)
1105   - const double t6_real = (tau) * (t1_real - t2_real);
1106   - const double t6_imag = (tau) * (t1_imag - t2_imag);
1107   -
1108   - // t7 = z0 - ((t5)/4)
1109   - const double t7_real = z0_real - t5_real / 4.0;
1110   - const double t7_imag = z0_imag - t5_imag / 4.0;
1111   -
1112   - // t8 = t7 + t6
1113   - const double t8_real = t7_real + t6_real;
1114   - const double t8_imag = t7_imag + t6_imag;
1115   -
1116   - // t9 = t7 - t6
1117   - const double t9_real = t7_real - t6_real;
1118   - const double t9_imag = t7_imag - t6_imag;
1119   -
1120   - // t10 = sin(2 pi/5) t3 + sin(2 pi/10) t4 */
1121   - const double t10_real = ((int) sign) * (sin_2pi_by_5 * t3_real +
1122   - sin_2pi_by_10 * t4_real);
1123   - const double t10_imag = ((int) sign) * (sin_2pi_by_5 * t3_imag +
1124   - sin_2pi_by_10 * t4_imag);
1125   -
1126   - // t11 = sin(2 pi/10) t3 - sin(2 pi/5) t4
1127   - const double t11_real = ((int) sign) * (sin_2pi_by_10 * t3_real -
1128   - sin_2pi_by_5 * t4_real);
1129   - const double t11_imag = ((int) sign) * (sin_2pi_by_10 * t3_imag -
1130   - sin_2pi_by_5 * t4_imag);
1131   -
1132   - // x0 = z0 + t5
1133   - const double x0_real = z0_real + t5_real;
1134   - const double x0_imag = z0_imag + t5_imag;
1135   -
1136   - // x1 = t8 + i t10
1137   - const double x1_real = t8_real - t10_imag;
1138   - const double x1_imag = t8_imag + t10_real;
1139   -
1140   - // x2 = t9 + i t11
1141   - const double x2_real = t9_real - t11_imag;
1142   - const double x2_imag = t9_imag + t11_real;
1143   -
1144   - // x3 = t9 - i t11
1145   - const double x3_real = t9_real + t11_imag;
1146   - const double x3_imag = t9_imag - t11_real;
1147   -
1148   - // x4 = t8 - i t10
1149   - const double x4_real = t8_real + t10_imag;
1150   - const double x4_imag = t8_imag - t10_real;
1151   -
1152   - // out0 = 1 * x0
1153   - //out[ostride*j].real(x0_real);
1154   - //out[ostride*j].imag(x0_imag);
1155   - out[ostride*j] = std::complex<double>(x0_real, x0_imag);
1156   -
1157   - // out1 = w1 * x1
1158   - //out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
1159   - //out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
1160   - out[ostride*(j+p_1)] = std::complex<double>(w1_real * x1_real - w1_imag * x1_imag,
1161   - w1_real * x1_imag + w1_imag * x1_real);
1162   -
1163   - // out2 = w2 * x2
1164   - //out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
1165   - //out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
1166   - out[ostride*(j+2*p_1)] = std::complex<double>(w2_real * x2_real - w2_imag * x2_imag,
1167   - w2_real * x2_imag + w2_imag * x2_real);
1168   -
1169   - // out3 = w3 * x3
1170   - //out[ostride*(j+3*p_1)].real(w3_real * x3_real - w3_imag * x3_imag);
1171   - //out[ostride*(j+3*p_1)].imag(w3_real * x3_imag + w3_imag * x3_real);
1172   - out[ostride*(j+3*p_1)] = std::complex<double>(w3_real * x3_real - w3_imag * x3_imag,
1173   - w3_real * x3_imag + w3_imag * x3_real);
1174   -
1175   - // out4 = w4 * x4
1176   - //out[ostride*(j+4*p_1)].real(w4_real * x4_real - w4_imag * x4_imag);
1177   - //out[ostride*(j+4*p_1)].imag(w4_real * x4_imag + w4_imag * x4_real);
1178   - out[ostride*(j+4*p_1)] = std::complex<double>(w4_real * x4_real - w4_imag * x4_imag,
1179   - w4_real * x4_imag + w4_imag * x4_real);
  882 + const std::complex<double> z0 = *(in+istride*i);
  883 + const std::complex<double> z1 = *(in+istride*(i+m));
  884 + const std::complex<double> z2 = *(in+istride*(i+2*m));
  885 + const std::complex<double> z3 = *(in+istride*(i+3*m));
  886 + const std::complex<double> z4 = *(in+istride*(i+4*m));
  887 +
  888 + // Compute t
  889 + const std::complex<double> t1 = z1 + z4;
  890 + const std::complex<double> t2 = z2 + z3;
  891 + const std::complex<double> t3 = z1 - z4;
  892 + const std::complex<double> t4 = z2 - z3;
  893 + const std::complex<double> t5 = t1 + t2;
  894 + const std::complex<double> t6 = tau * (t1 - t2);
  895 + const std::complex<double> t7 = z0 - t5/4.0;
  896 + const std::complex<double> t8 = t7 + t6;
  897 + const std::complex<double> t9 = t7 - t6;
  898 + const std::complex<double> t10 = double((int)sign) *
  899 + (sin_2pi_by_5*t3 + sin_2pi_by_10*t4);
  900 + const std::complex<double> t11 = double((int)sign) *
  901 + (sin_2pi_by_10*t3 - sin_2pi_by_5*t4);
  902 +
  903 + // Compute x
  904 + const std::complex<double> x0 = z0 + t5;
  905 + const std::complex<double> x1 = t8 + timesi(t10);
  906 + const std::complex<double> x2 = t9 + timesi(t11);
  907 + const std::complex<double> x3 = t9 - timesi(t11);
  908 + const std::complex<double> x4 = t8 - timesi(t10);
  909 +
  910 + // Compute out = w * x
  911 + *(out+ostride*j) = x0;
  912 + *(out+ostride*(j+p_1)) = w[0] * x1;
  913 + *(out+ostride*(j+2*p_1)) = w[1] * x2;
  914 + *(out+ostride*(j+3*p_1)) = w[2] * x3;
  915 + *(out+ostride*(j+4*p_1)) = w[3] * x4;
1180 916  
1181 917 } // endfor: k1
1182 918  
... ... @@ -1190,9 +926,9 @@ void GFft::factor5(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1190 926 /***********************************************************************//**
1191 927 * @brief Compute FFT for factor 6
1192 928 *
1193   - * @param[in] in Input array.
  929 + * @param[in] in Pointer to input array.
1194 930 * @param[in] istride Step size when traversing input array.
1195   - * @param[in] out Output array.
  931 + * @param[in] out Pointer to output array.
1196 932 * @param[in] ostride Step size when traversing output array.
1197 933 * @param[in] wavetable Trigonometric coefficients.
1198 934 * @param[in] sign Forward (-1) or backward (+1).
... ... @@ -1200,15 +936,15 @@ void GFft::factor5(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1200 936 * @param[in] n Logical array length.
1201 937 * @param[in] index Start index of trigonometric coefficients.
1202 938 ***************************************************************************/
1203   -void GFft::factor6(const std::vector<std::complex<double> >& in,
1204   - const int& istride,
1205   - std::vector<std::complex<double> >& out,
1206   - const int& ostride,
1207   - const GFftWavetable& wavetable,
1208   - const int& sign,
1209   - const int& product,
1210   - const int& n,
1211   - const int& index)
  939 +void GFft::factor6(const std::complex<double>* in,
  940 + const int& istride,
  941 + std::complex<double>* out,
  942 + const int& ostride,
  943 + const GFftWavetable& wavetable,
  944 + const int& sign,
  945 + const int& product,
  946 + const int& n,
  947 + const int& index)
1212 948 {
1213 949 // Compute ...
1214 950 const int factor = 6;
... ... @@ -1222,193 +958,57 @@ void GFft::factor6(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1222 958  
1223 959 // Loop over ...
1224 960 for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
1225   -
1226   - // Allocate w
1227   - double w1_real, w1_imag, w2_real, w2_imag, w3_real, w3_imag, w4_real,
1228   - w4_imag, w5_real, w5_imag;
1229   -
1230   - // Set trigonometric coefficients for k=0 since they are not stored
1231   - // in the wavetable object
1232   - if (k == 0) {
1233   - w1_real = 1.0;
1234   - w1_imag = 0.0;
1235   - w2_real = 1.0;
1236   - w2_imag = 0.0;
1237   - w3_real = 1.0;
1238   - w3_imag = 0.0;
1239   - w4_real = 1.0;
1240   - w4_imag = 0.0;
1241   - w5_real = 1.0;
1242   - w5_imag = 0.0;
1243   - }
1244   - else {
1245 961  
1246   - // Compute indices
1247   - int twiddle1 = index + k - 1;
1248   - int twiddle2 = twiddle1 + q;
1249   - int twiddle3 = twiddle2 + q;
1250   - int twiddle4 = twiddle3 + q;
1251   - int twiddle5 = twiddle4 + q;
1252   -
1253   - // Set trigonometric coefficients for forward transform
1254   - if (sign == -1) {
1255   - w1_real = wavetable[twiddle1].real();
1256   - w1_imag = wavetable[twiddle1].imag();
1257   - w2_real = wavetable[twiddle2].real();
1258   - w2_imag = wavetable[twiddle2].imag();
1259   - w3_real = wavetable[twiddle3].real();
1260   - w3_imag = wavetable[twiddle3].imag();
1261   - w4_real = wavetable[twiddle4].real();
1262   - w4_imag = wavetable[twiddle4].imag();
1263   - w5_real = wavetable[twiddle5].real();
1264   - w5_imag = wavetable[twiddle5].imag();
1265   - }
  962 + // Extract coefficients from wavetable
  963 + std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 5, sign);
1266 964  
1267   - // ... otherwise set trigonometric coefficients for backward
1268   - // tranform: w -> conjugate(w)
1269   - else {
1270   - w1_real = wavetable[twiddle1].real();
1271   - w1_imag = -wavetable[twiddle1].imag();
1272   - w2_real = wavetable[twiddle2].real();
1273   - w2_imag = -wavetable[twiddle2].imag();
1274   - w3_real = wavetable[twiddle3].real();
1275   - w3_imag = -wavetable[twiddle3].imag();
1276   - w4_real = wavetable[twiddle4].real();
1277   - w4_imag = -wavetable[twiddle4].imag();
1278   - w5_real = wavetable[twiddle5].real();
1279   - w5_imag = -wavetable[twiddle5].imag();
1280   - }
1281   -
1282   - } // endelse: k != 0
1283   -
1284   - // ...
  965 + // Compute x = W(6) z. W(6) is a combination of sums and differences of
  966 + // W(3) acting on the even and odd elements of z
1285 967 for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
1286 968  
1287 969 // Get z
1288   - const double z0_real = in[istride*i].real();
1289   - const double z0_imag = in[istride*i].imag();
1290   - const double z1_real = in[istride*(i+m)].real();
1291   - const double z1_imag = in[istride*(i+m)].imag();
1292   - const double z2_real = in[istride*(i+2*m)].real();
1293   - const double z2_imag = in[istride*(i+2*m)].imag();
1294   - const double z3_real = in[istride*(i+3*m)].real();
1295   - const double z3_imag = in[istride*(i+3*m)].imag();
1296   - const double z4_real = in[istride*(i+4*m)].real();
1297   - const double z4_imag = in[istride*(i+4*m)].imag();
1298   - const double z5_real = in[istride*(i+5*m)].real();
1299   - const double z5_imag = in[istride*(i+5*m)].imag();
1300   -
1301   - // compute x = W(6) z
1302   -
1303   - // W(6) is a combination of sums and differences of W(3) acting
1304   - // on the even and odd elements of z
1305   -
1306   - // ta1 = z2 + z4
1307   - const double ta1_real = z2_real + z4_real;
1308   - const double ta1_imag = z2_imag + z4_imag;
1309   -
1310   - // ta2 = z0 - ta1/2
1311   - const double ta2_real = z0_real - ta1_real / 2.0;
1312   - const double ta2_imag = z0_imag - ta1_imag / 2.0;
1313   -
1314   - // ta3 = (+/-) sin(pi/3)*(z2 - z4)
1315   - const double ta3_real = ((int) sign) * tau * (z2_real - z4_real);
1316   - const double ta3_imag = ((int) sign) * tau * (z2_imag - z4_imag);
1317   -
1318   - // a0 = z0 + ta1
1319   - const double a0_real = z0_real + ta1_real;
1320   - const double a0_imag = z0_imag + ta1_imag;
1321   -
1322   - // a1 = ta2 + i ta3
1323   - const double a1_real = ta2_real - ta3_imag;
1324   - const double a1_imag = ta2_imag + ta3_real;
1325   -
1326   - // a2 = ta2 - i ta3
1327   - const double a2_real = ta2_real + ta3_imag;
1328   - const double a2_imag = ta2_imag - ta3_real;
1329   -
1330   - // tb1 = z5 + z1
1331   - const double tb1_real = z5_real + z1_real;
1332   - const double tb1_imag = z5_imag + z1_imag;
1333   -
1334   - // tb2 = z3 - tb1/2
1335   - const double tb2_real = z3_real - tb1_real / 2.0;
1336   - const double tb2_imag = z3_imag - tb1_imag / 2.0;
1337   -
1338   - // tb3 = (+/-) sin(pi/3)*(z5 - z1)
1339   - const double tb3_real = ((int) sign) * tau * (z5_real - z1_real);
1340   - const double tb3_imag = ((int) sign) * tau * (z5_imag - z1_imag);
1341   -
1342   - // b0 = z3 + tb1
1343   - const double b0_real = z3_real + tb1_real;
1344   - const double b0_imag = z3_imag + tb1_imag;
1345   -
1346   - // b1 = tb2 + i tb3
1347   - const double b1_real = tb2_real - tb3_imag;
1348   - const double b1_imag = tb2_imag + tb3_real;
1349   -
1350   - // b2 = tb2 - i tb3
1351   - const double b2_real = tb2_real + tb3_imag;
1352   - const double b2_imag = tb2_imag - tb3_real;
1353   -
1354   - // x0 = a0 + b0
1355   - const double x0_real = a0_real + b0_real;
1356   - const double x0_imag = a0_imag + b0_imag;
1357   -
1358   - // x4 = a1 + b1
1359   - const double x4_real = a1_real + b1_real;
1360   - const double x4_imag = a1_imag + b1_imag;
1361   -
1362   - // x2 = a2 + b2
1363   - const double x2_real = a2_real + b2_real;
1364   - const double x2_imag = a2_imag + b2_imag;
1365   -
1366   - // x3 = a0 - b0
1367   - const double x3_real = a0_real - b0_real;
1368   - const double x3_imag = a0_imag - b0_imag;
1369   -
1370   - // x1 = a1 - b1
1371   - const double x1_real = a1_real - b1_real;
1372   - const double x1_imag = a1_imag - b1_imag;
1373   -
1374   - // x5 = a2 - b2
1375   - const double x5_real = a2_real - b2_real;
1376   - const double x5_imag = a2_imag - b2_imag;
1377   -
1378   - // out0 = 1 * x0
1379   - //out[ostride*j].real(x0_real);
1380   - //out[ostride*j].imag(x0_imag);
1381   - out[ostride*j] = std::complex<double>(x0_real, x0_imag);
1382   -
1383   - // out1 = w1 * x1
1384   - //out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
1385   - //out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
1386   - out[ostride*(j+p_1)] = std::complex<double>(w1_real * x1_real - w1_imag * x1_imag,
1387   - w1_real * x1_imag + w1_imag * x1_real);
1388   -
1389   - // out2 = w2 * x2
1390   - //out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
1391   - //out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
1392   - out[ostride*(j+2*p_1)] = std::complex<double>(w2_real * x2_real - w2_imag * x2_imag,
1393   - w2_real * x2_imag + w2_imag * x2_real);
1394   -
1395   - // out3 = w3 * x3
1396   - //out[ostride*(j+3*p_1)].real(w3_real * x3_real - w3_imag * x3_imag);
1397   - //out[ostride*(j+3*p_1)].imag(w3_real * x3_imag + w3_imag * x3_real);
1398   - out[ostride*(j+3*p_1)] = std::complex<double>(w3_real * x3_real - w3_imag * x3_imag,
1399   - w3_real * x3_imag + w3_imag * x3_real);
1400   -
1401   - // out4 = w4 * x4
1402   - //out[ostride*(j+4*p_1)].real(w4_real * x4_real - w4_imag * x4_imag);
1403   - //out[ostride*(j+4*p_1)].imag(w4_real * x4_imag + w4_imag * x4_real);
1404   - out[ostride*(j+4*p_1)] = std::complex<double>(w4_real * x4_real - w4_imag * x4_imag,
1405   - w4_real * x4_imag + w4_imag * x4_real);
  970 + const std::complex<double> z0 = *(in+istride*i);
  971 + const std::complex<double> z1 = *(in+istride*(i+m));
  972 + const std::complex<double> z2 = *(in+istride*(i+2*m));
  973 + const std::complex<double> z3 = *(in+istride*(i+3*m));
  974 + const std::complex<double> z4 = *(in+istride*(i+4*m));
  975 + const std::complex<double> z5 = *(in+istride*(i+5*m));
  976 +
  977 + // Compute ta
  978 + const std::complex<double> ta1 = z2 + z4;
  979 + const std::complex<double> ta2 = z0 - ta1/2.0;
  980 + const std::complex<double> ta3 = double((int)sign) * tau * (z2 - z4);
  981 +
  982 + // Compute a
  983 + const std::complex<double> a0 = z0 + ta1;
  984 + const std::complex<double> a1 = ta2 + timesi(ta3);
  985 + const std::complex<double> a2 = ta2 - timesi(ta3);
  986 +
  987 + // Compute tb
  988 + const std::complex<double> tb1 = z5 + z1;
  989 + const std::complex<double> tb2 = z3 - tb1/2.0;
  990 + const std::complex<double> tb3 = double((int)sign) * tau * (z5 - z1);
  991 +
  992 + // Compute b
  993 + const std::complex<double> b0 = z3 + tb1;
  994 + const std::complex<double> b1 = tb2 + timesi(tb3);
  995 + const std::complex<double> b2 = tb2 - timesi(tb3);
  996 +
  997 + // Compute x
  998 + const std::complex<double> x0 = a0 + b0;
  999 + const std::complex<double> x1 = a1 - b1;
  1000 + const std::complex<double> x2 = a2 + b2;
  1001 + const std::complex<double> x3 = a0 - b0;
  1002 + const std::complex<double> x4 = a1 + b1;
  1003 + const std::complex<double> x5 = a2 - b2;
1406 1004  
1407   - // out5 = w5 * x5
1408   - //out[ostride*(j+5*p_1)].real(w5_real * x5_real - w5_imag * x5_imag);
1409   - //out[ostride*(j+5*p_1)].imag(w5_real * x5_imag + w5_imag * x5_real);
1410   - out[ostride*(j+5*p_1)] = std::complex<double>(w5_real * x5_real - w5_imag * x5_imag,
1411   - w5_real * x5_imag + w5_imag * x5_real);
  1005 + // Compute out = w * x
  1006 + *(out+ostride*j) = x0;
  1007 + *(out+ostride*(j+p_1)) = w[0] * x1;
  1008 + *(out+ostride*(j+2*p_1)) = w[1] * x2;
  1009 + *(out+ostride*(j+3*p_1)) = w[2] * x3;
  1010 + *(out+ostride*(j+4*p_1)) = w[3] * x4;
  1011 + *(out+ostride*(j+5*p_1)) = w[4] * x5;
1412 1012  
1413 1013 } // endfor: k1
1414 1014  
... ... @@ -1422,9 +1022,9 @@ void GFft::factor6(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1422 1022 /***********************************************************************//**
1423 1023 * @brief Compute FFT for factor 7
1424 1024 *
1425   - * @param[in] in Input array.
  1025 + * @param[in] in Pointer to input array.
1426 1026 * @param[in] istride Step size when traversing input array.
1427   - * @param[in] out Output array.
  1027 + * @param[in] out Pointer to output array.
1428 1028 * @param[in] ostride Step size when traversing output array.
1429 1029 * @param[in] wavetable Trigonometric coefficients.
1430 1030 * @param[in] sign Forward (-1) or backward (+1).
... ... @@ -1432,15 +1032,15 @@ void GFft::factor6(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1432 1032 * @param[in] n Logical array length.
1433 1033 * @param[in] index Start index of trigonometric coefficients.
1434 1034 ***************************************************************************/
1435   -void GFft::factor7(const std::vector<std::complex<double> >& in,
1436   - const int& istride,
1437   - std::vector<std::complex<double> >& out,
1438   - const int& ostride,
1439   - const GFftWavetable& wavetable,
1440   - const int& sign,
1441   - const int& product,
1442   - const int& n,
1443   - const int& index)
  1035 +void GFft::factor7(const std::complex<double>* in,
  1036 + const int& istride,
  1037 + std::complex<double>* out,
  1038 + const int& ostride,
  1039 + const GFftWavetable& wavetable,
  1040 + const int& sign,
  1041 + const int& product,
  1042 + const int& n,
  1043 + const int& index)
1444 1044 {
1445 1045 // Compute ...
1446 1046 const int factor = 7;
... ... @@ -1468,282 +1068,76 @@ void GFft::factor7(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1468 1068  
1469 1069 // Loop over ...
1470 1070 for (int k = 0, i = 0, j = 0; k < q; ++k, j += jump) {
1471   -
1472   - //
1473   - double w1_real, w1_imag, w2_real, w2_imag, w3_real, w3_imag, w4_real,
1474   - w4_imag, w5_real, w5_imag, w6_real, w6_imag;
1475   -
1476   - // Set trigonometric coefficients for k=0 since they are not stored
1477   - // in the wavetable object
1478   - if (k == 0) {
1479   - w1_real = 1.0;
1480   - w1_imag = 0.0;
1481   - w2_real = 1.0;
1482   - w2_imag = 0.0;
1483   - w3_real = 1.0;
1484   - w3_imag = 0.0;
1485   - w4_real = 1.0;
1486   - w4_imag = 0.0;
1487   - w5_real = 1.0;
1488   - w5_imag = 0.0;
1489   - w6_real = 1.0;
1490   - w6_imag = 0.0;
1491   - }
1492   - else {
1493   -
1494   - // Compute indices
1495   - int twiddle1 = index + k - 1;
1496   - int twiddle2 = twiddle1 + q;
1497   - int twiddle3 = twiddle2 + q;
1498   - int twiddle4 = twiddle3 + q;
1499   - int twiddle5 = twiddle4 + q;
1500   - int twiddle6 = twiddle5 + q;
1501   -
1502   - // Set trigonometric coefficients for forward transform
1503   - if (sign == -1) {
1504   - w1_real = wavetable[twiddle1].real();
1505   - w1_imag = wavetable[twiddle1].imag();
1506   - w2_real = wavetable[twiddle2].real();
1507   - w2_imag = wavetable[twiddle2].imag();
1508   - w3_real = wavetable[twiddle3].real();
1509   - w3_imag = wavetable[twiddle3].imag();
1510   - w4_real = wavetable[twiddle4].real();
1511   - w4_imag = wavetable[twiddle4].imag();
1512   - w5_real = wavetable[twiddle5].real();
1513   - w5_imag = wavetable[twiddle5].imag();
1514   - w6_real = wavetable[twiddle6].real();
1515   - w6_imag = wavetable[twiddle6].imag();
1516   - }
1517   -
1518   - // ... otherwise set trigonometric coefficients for backward
1519   - // tranform: w -> conjugate(w)
1520   - else {
1521   - w1_real = wavetable[twiddle1].real();
1522   - w1_imag = -wavetable[twiddle1].imag();
1523   - w2_real = wavetable[twiddle2].real();
1524   - w2_imag = -wavetable[twiddle2].imag();
1525   - w3_real = wavetable[twiddle3].real();
1526   - w3_imag = -wavetable[twiddle3].imag();
1527   - w4_real = wavetable[twiddle4].real();
1528   - w4_imag = -wavetable[twiddle4].imag();
1529   - w5_real = wavetable[twiddle5].real();
1530   - w5_imag = -wavetable[twiddle5].imag();
1531   - w6_real = wavetable[twiddle6].real();
1532   - w6_imag = -wavetable[twiddle6].imag();
1533   - }
1534 1071  
1535   - } // endelse: k != 0
  1072 + // Extract coefficients from wavetable
  1073 + std::vector<std::complex<double> > w = get_w(wavetable, index, k, q, 6, sign);
1536 1074  
1537   - // ...
  1075 + // Compute x = W(7) z
1538 1076 for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
1539 1077  
1540 1078 // Get z
1541   - const double z0_real = in[istride*i].real();
1542   - const double z0_imag = in[istride*i].imag();
1543   - const double z1_real = in[istride*(i+m)].real();
1544   - const double z1_imag = in[istride*(i+m)].imag();
1545   - const double z2_real = in[istride*(i+2*m)].real();
1546   - const double z2_imag = in[istride*(i+2*m)].imag();
1547   - const double z3_real = in[istride*(i+3*m)].real();
1548   - const double z3_imag = in[istride*(i+3*m)].imag();
1549   - const double z4_real = in[istride*(i+4*m)].real();
1550   - const double z4_imag = in[istride*(i+4*m)].imag();
1551   - const double z5_real = in[istride*(i+5*m)].real();
1552   - const double z5_imag = in[istride*(i+5*m)].imag();
1553   - const double z6_real = in[istride*(i+6*m)].real();
1554   - const double z6_imag = in[istride*(i+6*m)].imag();
1555   -
1556   - // Compute x = W(7) z
1557   -
1558   - // t0 = z1 + z6
1559   - const double t0_real = z1_real + z6_real;
1560   - const double t0_imag = z1_imag + z6_imag;
1561   -
1562   - // t1 = z1 - z6
1563   - const double t1_real = z1_real - z6_real;
1564   - const double t1_imag = z1_imag - z6_imag;
1565   -
1566   - // t2 = z2 + z5
1567   - const double t2_real = z2_real + z5_real;
1568   - const double t2_imag = z2_imag + z5_imag;
1569   -
1570   - // t3 = z2 - z5
1571   - const double t3_real = z2_real - z5_real;
1572   - const double t3_imag = z2_imag - z5_imag;
1573   -
1574   - // t4 = z4 + z3
1575   - const double t4_real = z4_real + z3_real;
1576   - const double t4_imag = z4_imag + z3_imag;
1577   -
1578   - // t5 = z4 - z3
1579   - const double t5_real = z4_real - z3_real;
1580   - const double t5_imag = z4_imag - z3_imag;
1581   -
1582   - // t6 = t2 + t0
1583   - const double t6_real = t2_real + t0_real;
1584   - const double t6_imag = t2_imag + t0_imag;
1585   -
1586   - // t7 = t5 + t3
1587   - const double t7_real = t5_real + t3_real;
1588   - const double t7_imag = t5_imag + t3_imag;
1589   -
1590   - // b0 = z0 + t6 + t4
1591   - const double b0_real = z0_real + t6_real + t4_real;
1592   - const double b0_imag = z0_imag + t6_imag + t4_imag;
  1079 + const std::complex<double> z0 = *(in+istride*i);
  1080 + const std::complex<double> z1 = *(in+istride*(i+m));
  1081 + const std::complex<double> z2 = *(in+istride*(i+2*m));
  1082 + const std::complex<double> z3 = *(in+istride*(i+3*m));
  1083 + const std::complex<double> z4 = *(in+istride*(i+4*m));
  1084 + const std::complex<double> z5 = *(in+istride*(i+5*m));
  1085 + const std::complex<double> z6 = *(in+istride*(i+6*m));
  1086 +
  1087 + // Compute t
  1088 + const std::complex<double> t0 = z1 + z6;
  1089 + const std::complex<double> t1 = z1 - z6;
  1090 + const std::complex<double> t2 = z2 + z5;
  1091 + const std::complex<double> t3 = z2 - z5;
  1092 + const std::complex<double> t4 = z4 + z3;
  1093 + const std::complex<double> t5 = z4 - z3;
  1094 + const std::complex<double> t6 = t2 + t0;
  1095 + const std::complex<double> t7 = t5 + t3;
  1096 +
  1097 + // Compute b
  1098 + const std::complex<double> b0 = z0 + t6 + t4;
  1099 + const std::complex<double> b1 = tau1 * (t6 + t4);
  1100 + const std::complex<double> b2 = tau2 * (t0 - t4);
  1101 + const std::complex<double> b3 = tau3 * (t4 - t2);
  1102 + const std::complex<double> b4 = tau4 * (t2 - t0);
  1103 + const std::complex<double> b5 = double(-(int)sign) * tau5 * (t7 + t1);
  1104 + const std::complex<double> b6 = double(-(int)sign) * tau6 * (t1 - t5);
  1105 + const std::complex<double> b7 = double(-(int)sign) * tau7 * (t5 - t3);
  1106 + const std::complex<double> b8 = double(-(int)sign) * tau8 * (t3 - t1);
  1107 +
  1108 + // Compute T
  1109 + const std::complex<double> T0 = b0 + b1;
  1110 + const std::complex<double> T1 = b2 + b3;
  1111 + const std::complex<double> T2 = b4 - b3;
  1112 + const std::complex<double> T3 = -b2 - b4;
  1113 + const std::complex<double> T4 = b6 + b7;
  1114 + const std::complex<double> T5 = b8 - b7;
  1115 + const std::complex<double> T6 = -b8 - b6;
  1116 + const std::complex<double> T7 = T0 + T1;
  1117 + const std::complex<double> T8 = T0 + T2;
  1118 + const std::complex<double> T9 = T0 + T3;
  1119 + const std::complex<double> T10 = T4 + b5;
  1120 + const std::complex<double> T11 = T5 + b5;
  1121 + const std::complex<double> T12 = T6 + b5;
  1122 +
  1123 + // Compute x
  1124 + const std::complex<double> x0 = b0;
  1125 + const std::complex<double> x1 = T7 - timesi(T10);
  1126 + const std::complex<double> x2 = T9 - timesi(T12);
  1127 + const std::complex<double> x3 = T8 + timesi(T11);
  1128 + const std::complex<double> x4 = T8 - timesi(T11);
  1129 + const std::complex<double> x5 = T9 + timesi(T12);
  1130 + const std::complex<double> x6 = T7 + timesi(T10);
  1131 +
  1132 + // Compute out = w * x
  1133 + *(out+ostride*j) = x0;
  1134 + *(out+ostride*(j+p_1)) = w[0] * x1;
  1135 + *(out+ostride*(j+2*p_1)) = w[1] * x2;
  1136 + *(out+ostride*(j+3*p_1)) = w[2] * x3;
  1137 + *(out+ostride*(j+4*p_1)) = w[3] * x4;
  1138 + *(out+ostride*(j+5*p_1)) = w[4] * x5;
  1139 + *(out+ostride*(j+6*p_1)) = w[5] * x6;
1593 1140  
1594   - // b1 = ((cos(2pi/7) + cos(4pi/7) + cos(6pi/7))/3-1) (t6 + t4)
1595   - const double b1_real = tau1 * (t6_real + t4_real);
1596   - const double b1_imag = tau1 * (t6_imag + t4_imag);
1597   -
1598   - // b2 = ((2*cos(2pi/7) - cos(4pi/7) - cos(6pi/7))/3) (t0 - t4)
1599   - const double b2_real = tau2 * (t0_real - t4_real);
1600   - const double b2_imag = tau2 * (t0_imag - t4_imag);
1601   -
1602   - // b3 = ((cos(2pi/7) - 2*cos(4pi/7) + cos(6pi/7))/3) (t4 - t2)
1603   - const double b3_real = tau3 * (t4_real - t2_real);
1604   - const double b3_imag = tau3 * (t4_imag - t2_imag);
1605   -
1606   - // b4 = ((cos(2pi/7) + cos(4pi/7) - 2*cos(6pi/7))/3) (t2 - t0)
1607   - const double b4_real = tau4 * (t2_real - t0_real);
1608   - const double b4_imag = tau4 * (t2_imag - t0_imag);
1609   -
1610   - // b5 = sign * ((sin(2pi/7) + sin(4pi/7) - sin(6pi/7))/3) (t7 + t1)
1611   - const double b5_real = (-(int)sign) * tau5 * (t7_real + t1_real);
1612   - const double b5_imag = (-(int)sign) * tau5 * (t7_imag + t1_imag);
1613   -
1614   - // b6 = sign * ((2sin(2pi/7) - sin(4pi/7) + sin(6pi/7))/3) (t1 - t5)
1615   - const double b6_real = (-(int)sign) * tau6 * (t1_real - t5_real);
1616   - const double b6_imag = (-(int)sign) * tau6 * (t1_imag - t5_imag);
1617   -
1618   - // b7 = sign * ((sin(2pi/7) - 2sin(4pi/7) - sin(6pi/7))/3) (t5 - t3)
1619   - const double b7_real = (-(int)sign) * tau7 * (t5_real - t3_real);
1620   - const double b7_imag = (-(int)sign) * tau7 * (t5_imag - t3_imag);
1621   -
1622   - // b8 = sign * ((sin(2pi/7) + sin(4pi/7) + 2sin(6pi/7))/3) (t3 - t1) */
1623   - const double b8_real = (-(int)sign) * tau8 * (t3_real - t1_real);
1624   - const double b8_imag = (-(int)sign) * tau8 * (t3_imag - t1_imag);
1625   -
1626   - // T0 = b0 + b1
1627   - const double T0_real = b0_real + b1_real;
1628   - const double T0_imag = b0_imag + b1_imag;
1629   -
1630   - // T1 = b2 + b3
1631   - const double T1_real = b2_real + b3_real;
1632   - const double T1_imag = b2_imag + b3_imag;
1633   -
1634   - // T2 = b4 - b3
1635   - const double T2_real = b4_real - b3_real;
1636   - const double T2_imag = b4_imag - b3_imag;
1637   -
1638   - // T3 = -b2 - b4
1639   - const double T3_real = -b2_real - b4_real;
1640   - const double T3_imag = -b2_imag - b4_imag;
1641   -
1642   - // T4 = b6 + b7
1643   - const double T4_real = b6_real + b7_real;
1644   - const double T4_imag = b6_imag + b7_imag;
1645   -
1646   - // T5 = b8 - b7
1647   - const double T5_real = b8_real - b7_real;
1648   - const double T5_imag = b8_imag - b7_imag;
1649   -
1650   - // T6 = -b8 - b6
1651   - const double T6_real = -b8_real - b6_real;
1652   - const double T6_imag = -b8_imag - b6_imag;
1653   -
1654   - // T7 = T0 + T1
1655   - const double T7_real = T0_real + T1_real;
1656   - const double T7_imag = T0_imag + T1_imag;
1657   -
1658   - // T8 = T0 + T2
1659   - const double T8_real = T0_real + T2_real;
1660   - const double T8_imag = T0_imag + T2_imag;
1661   -
1662   - // T9 = T0 + T3
1663   - const double T9_real = T0_real + T3_real;
1664   - const double T9_imag = T0_imag + T3_imag;
1665   -
1666   - // T10 = T4 + b5
1667   - const double T10_real = T4_real + b5_real;
1668   - const double T10_imag = T4_imag + b5_imag;
1669   -
1670   - // T11 = T5 + b5
1671   - const double T11_real = T5_real + b5_real;
1672   - const double T11_imag = T5_imag + b5_imag;
1673   -
1674   - // T12 = T6 + b5
1675   - const double T12_real = T6_real + b5_real;
1676   - const double T12_imag = T6_imag + b5_imag;
1677   -
1678   - // x0 = b0
1679   - const double x0_real = b0_real;
1680   - const double x0_imag = b0_imag;
1681   -
1682   - // x1 = T7 - i T10
1683   - const double x1_real = T7_real + T10_imag;
1684   - const double x1_imag = T7_imag - T10_real;
1685   -
1686   - // x2 = T9 - i T12
1687   - const double x2_real = T9_real + T12_imag;
1688   - const double x2_imag = T9_imag - T12_real;
1689   -
1690   - // x3 = T8 + i T11
1691   - const double x3_real = T8_real - T11_imag;
1692   - const double x3_imag = T8_imag + T11_real;
1693   -
1694   - // x4 = T8 - i T11
1695   - const double x4_real = T8_real + T11_imag;
1696   - const double x4_imag = T8_imag - T11_real;
1697   -
1698   - // x5 = T9 + i T12
1699   - const double x5_real = T9_real - T12_imag;
1700   - const double x5_imag = T9_imag + T12_real;
1701   -
1702   - // x6 = T7 + i T10
1703   - const double x6_real = T7_real - T10_imag;
1704   - const double x6_imag = T7_imag + T10_real;
1705   -
1706   - // out0 = 1 * x0
1707   - //out[ostride*j].real(x0_real);
1708   - //out[ostride*j].imag(x0_imag);
1709   - out[ostride*j] = std::complex<double>(x0_real, x0_imag);
1710   -
1711   - // out1 = w1 * x1
1712   - //out[ostride*(j+p_1)].real(w1_real * x1_real - w1_imag * x1_imag);
1713   - //out[ostride*(j+p_1)].imag(w1_real * x1_imag + w1_imag * x1_real);
1714   - out[ostride*(j+p_1)] = std::complex<double>(w1_real * x1_real - w1_imag * x1_imag,
1715   - w1_real * x1_imag + w1_imag * x1_real);
1716   -
1717   - // out2 = w2 * x2
1718   - //out[ostride*(j+2*p_1)].real(w2_real * x2_real - w2_imag * x2_imag);
1719   - //out[ostride*(j+2*p_1)].imag(w2_real * x2_imag + w2_imag * x2_real);
1720   - out[ostride*(j+2*p_1)] = std::complex<double>(w2_real * x2_real - w2_imag * x2_imag,
1721   - w2_real * x2_imag + w2_imag * x2_real);
1722   -
1723   - // out3 = w3 * x3
1724   - //out[ostride*(j+3*p_1)].real(w3_real * x3_real - w3_imag * x3_imag);
1725   - //out[ostride*(j+3*p_1)].imag(w3_real * x3_imag + w3_imag * x3_real);
1726   - out[ostride*(j+3*p_1)] = std::complex<double>(w3_real * x3_real - w3_imag * x3_imag,
1727   - w3_real * x3_imag + w3_imag * x3_real);
1728   -
1729   - // out4 = w4 * x4
1730   - //out[ostride*(j+4*p_1)].real(w4_real * x4_real - w4_imag * x4_imag);
1731   - //out[ostride*(j+4*p_1)].imag(w4_real * x4_imag + w4_imag * x4_real);
1732   - out[ostride*(j+4*p_1)] = std::complex<double>(w4_real * x4_real - w4_imag * x4_imag,
1733   - w4_real * x4_imag + w4_imag * x4_real);
1734   -
1735   - // out5 = w5 * x5
1736   - //out[ostride*(j+5*p_1)].real(w5_real * x5_real - w5_imag * x5_imag);
1737   - //out[ostride*(j+5*p_1)].imag(w5_real * x5_imag + w5_imag * x5_real);
1738   - out[ostride*(j+5*p_1)] = std::complex<double>(w5_real * x5_real - w5_imag * x5_imag,
1739   - w5_real * x5_imag + w5_imag * x5_real);
1740   -
1741   - // out6 = w6 * x6
1742   - //out[ostride*(j+6*p_1)].real(w6_real * x6_real - w6_imag * x6_imag);
1743   - //out[ostride*(j+6*p_1)].imag(w6_real * x6_imag + w6_imag * x6_real);
1744   - out[ostride*(j+6*p_1)] = std::complex<double>(w6_real * x6_real - w6_imag * x6_imag,
1745   - w6_real * x6_imag + w6_imag * x6_real);
1746   -
1747 1141 } // endfor: k1
1748 1142  
1749 1143 } // endfor: k
... ... @@ -1756,9 +1150,9 @@ void GFft::factor7(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1756 1150 /***********************************************************************//**
1757 1151 * @brief Compute FFT for arbitrary factor
1758 1152 *
1759   - * @param[in] in Input array.
  1153 + * @param[in] in Pointer to input array.
1760 1154 * @param[in] istride Step size when traversing input array.
1761   - * @param[in] out Output array.
  1155 + * @param[in] out Pointer to output array.
1762 1156 * @param[in] ostride Step size when traversing output array.
1763 1157 * @param[in] wavetable Trigonometric coefficients.
1764 1158 * @param[in] sign Forward (-1) or backward (+1).
... ... @@ -1767,16 +1161,16 @@ void GFft::factor7(const std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1767 1161 * @param[in] n Logical array length.
1768 1162 * @param[in] index Start index of trigonometric coefficients.
1769 1163 ***************************************************************************/
1770   -void GFft::factorn(std::vector<std::complex<double> >& in,
1771   - const int& istride,
1772   - std::vector<std::complex<double> >& out,
1773   - const int& ostride,
1774   - const GFftWavetable& wavetable,
1775   - const int& sign,
1776   - const int& factor,
1777   - const int& product,
1778   - const int& n,
1779   - const int& index)
  1164 +void GFft::factorn(std::complex<double>* in,
  1165 + const int& istride,
  1166 + std::complex<double>* out,
  1167 + const int& ostride,
  1168 + const GFftWavetable& wavetable,
  1169 + const int& sign,
  1170 + const int& factor,
  1171 + const int& product,
  1172 + const int& n,
  1173 + const int& index)
1780 1174 {
1781 1175 // Compute ...
1782 1176 const int m = n / factor;
... ... @@ -1786,7 +1180,7 @@ void GFft::factorn(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1786 1180  
1787 1181 // ...
1788 1182 for (int i = 0; i < m; ++i) {
1789   - out[ostride*i] = in[istride*i];
  1183 + *(out+ostride*i) = *(in+istride*i);
1790 1184 }
1791 1185  
1792 1186 // ...
... ... @@ -1794,19 +1188,19 @@ void GFft::factorn(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1794 1188 for (int i = 0; i < m; ++i) {
1795 1189 const int idx = i + e * m;
1796 1190 const int idxc = i + (factor - e) * m;
1797   - out[ostride*idx] = in[istride*idx] + in[istride*idxc];
1798   - out[ostride*idxc] = in[istride*idx] - in[istride*idxc];
  1191 + *(out+ostride*idx) = *(in+istride*idx) + *(in+istride*idxc);
  1192 + *(out+ostride*idxc) = *(in+istride*idx) - *(in+istride*idxc);
1799 1193 }
1800 1194 }
1801 1195  
1802 1196 // e = 0
1803 1197 for (int i = 0; i < m; ++i) {
1804   - in[istride*i] = out[ostride*i];
  1198 + *(in+istride*i) = *(out+ostride*i);
1805 1199 }
1806 1200 for (int e = 1; e < (factor - 1) / 2 + 1; ++e) {
1807 1201 for (int i = 0; i < m; ++i) {
1808   - const int idx = i + e * m;
1809   - in[istride*i] += out[ostride*idx];
  1202 + const int idx = i + e * m;
  1203 + *(in+istride*i) += *(out+ostride*idx);
1810 1204 }
1811 1205 }
1812 1206  
... ... @@ -1826,8 +1220,8 @@ void GFft::factorn(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1826 1220  
1827 1221 // ...
1828 1222 for (int i = 0; i < m; ++i) {
1829   - in[istride*(i+em)] = out[ostride*i];
1830   - in[istride*(i+ecm)] = out[ostride*i];
  1223 + *(in+istride*(i+em)) = *(out+ostride*i);
  1224 + *(in+istride*(i+ecm)) = *(out+ostride*i);
1831 1225 }
1832 1226  
1833 1227 // ...
... ... @@ -1860,10 +1254,10 @@ void GFft::factorn(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1860 1254 for (int i = 0; i < m; ++i) {
1861 1255  
1862 1256 //
1863   - const double xp_real = out[ostride*(i + e1 * m)].real();
1864   - const double xp_imag = out[ostride*(i + e1 * m)].imag();
1865   - const double xm_real = out[ostride*(i + (factor-e1)*m)].real();
1866   - const double xm_imag = out[ostride*(i + (factor-e1)*m)].imag();
  1257 + const double xp_real = (out+ostride*(i + e1 * m))->real();
  1258 + const double xp_imag = (out+ostride*(i + e1 * m))->imag();
  1259 + const double xm_real = (out+ostride*(i + (factor-e1)*m))->real();
  1260 + const double xm_imag = (out+ostride*(i + (factor-e1)*m))->imag();
1867 1261  
1868 1262 //
1869 1263 const double ap = w_real * xp_real;
... ... @@ -1882,14 +1276,10 @@ void GFft::factorn(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1882 1276 double sumc_imag = bp - bm;
1883 1277  
1884 1278 //
1885   - //in[istride*(i+em)].real(in[istride*(i+em)].real() + sum_real);
1886   - //in[istride*(i+em)].imag(in[istride*(i+em)].imag() + sum_imag);
1887   - //in[istride*(i+ecm)].real(in[istride*(i+ecm)].real() + sumc_real);
1888   - //in[istride*(i+ecm)].imag(in[istride*(i+ecm)].imag() + sumc_imag);
1889   - in[istride*(i+em)] = std::complex<double>(in[istride*(i+em)].real() + sum_real,
1890   - in[istride*(i+em)].imag() + sum_imag);
1891   - in[istride*(i+ecm)] = std::complex<double>(in[istride*(i+ecm)].real() + sumc_real,
1892   - in[istride*(i+ecm)].imag() + sumc_imag);
  1279 + *(in+istride*(i+em)) = std::complex<double>((in+istride*(i+em))->real() + sum_real,
  1280 + (in+istride*(i+em))->imag() + sum_imag);
  1281 + *(in+istride*(i+ecm)) = std::complex<double>((in+istride*(i+ecm))->real() + sumc_real,
  1282 + (in+istride*(i+ecm))->imag() + sumc_imag);
1893 1283  
1894 1284  
1895 1285 } // endfor: i
... ... @@ -1904,20 +1294,20 @@ void GFft::factorn(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1904 1294  
1905 1295 // k = 0
1906 1296 for (int k1 = 0; k1 < p_1; ++k1) {
1907   - out[ostride*k1] = in[istride*k1];
  1297 + *(out+ostride*k1) = *(in+istride*k1);
1908 1298 }
1909 1299  
1910 1300 // k > 0
1911 1301 for (int e1 = 1; e1 < factor; ++e1) {
1912 1302 for (int k1 = 0; k1 < p_1; ++k1) {
1913   - out[ostride*(k1 + e1 * p_1)] = in[istride*(k1 + e1 * m)];
  1303 + *(out+ostride*(k1 + e1 * p_1)) = *(in+istride*(k1 + e1 * m));
1914 1304 }
1915 1305 }
1916 1306  
1917 1307 // e = 0
1918 1308 for (int k = 1, i = p_1, j = product; k < q; ++k, j += jump) {
1919 1309 for (int k1 = 0; k1 < p_1; ++k1, ++i, ++j) {
1920   - out[ostride*j] = in[istride*i];
  1310 + *(out+ostride*j) = *(in+istride*i);
1921 1311 }
1922 1312 }
1923 1313  
... ... @@ -1929,8 +1319,8 @@ void GFft::factorn(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1929 1319 for (int e1 = 1; e1 < factor; ++e1) {
1930 1320  
1931 1321 // Get x
1932   - double x_real = in[istride*(i + e1 * m)].real();
1933   - double x_imag = in[istride*(i + e1 * m)].imag();
  1322 + double x_real = (in+istride*(i + e1 * m))->real();
  1323 + double x_imag = (in+istride*(i + e1 * m))->imag();
1934 1324  
1935 1325 // Get w
1936 1326 double w_real, w_imag;
... ... @@ -1951,10 +1341,8 @@ void GFft::factorn(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1951 1341 }
1952 1342  
1953 1343 //
1954   - //out[ostride*(j + e1 * p_1)].real(w_real * x_real - w_imag * x_imag);
1955   - //out[ostride*(j + e1 * p_1)].imag(w_real * x_imag + w_imag * x_real);
1956   - out[ostride*(j + e1 * p_1)] = std::complex<double>(w_real * x_real - w_imag * x_imag,
1957   - w_real * x_imag + w_imag * x_real);
  1344 + *(out+ostride*(j + e1 * p_1)) = std::complex<double>(w_real * x_real - w_imag * x_imag,
  1345 + w_real * x_imag + w_imag * x_real);
1958 1346  
1959 1347 } // endfor: e1
1960 1348  
... ... @@ -1965,3 +1353,57 @@ void GFft::factorn(std::vector&lt;std::complex&lt;double&gt; &gt;&amp; in,
1965 1353 // Return
1966 1354 return;
1967 1355 }
  1356 +
  1357 +
  1358 +/***********************************************************************//**
  1359 + * @brief Extract coefficients from wavetable
  1360 + *
  1361 + * @param[in] index Start index of trigonometric coefficients.
  1362 + * @param[in] k ...
  1363 + * @param[in] q ...
  1364 + * @param[in] n Number of coefficients.
  1365 + * @param[in] sign Forward (-1) or backward (+1) transformation.
  1366 + * @param[in] wavetable Trigonometric coefficients.
  1367 + ***************************************************************************/
  1368 +std::vector<std::complex<double> > GFft::get_w(const GFftWavetable& wavetable,
  1369 + const int& index,
  1370 + const int& k,
  1371 + const int& q,
  1372 + const int& n,
  1373 + const int& sign) const
  1374 +{
  1375 + // Allocate w vectors
  1376 + std::vector<std::complex<double> > w(n);
  1377 +
  1378 + // Set trigonometric coefficients for k=0 since they are not stored in the
  1379 + // wavetable object
  1380 + if (k == 0) {
  1381 + w.assign(n, std::complex<double>(1.0, 0.0));
  1382 + }
  1383 +
  1384 + // ... otherwise use coefficients stored in wavetable
  1385 + else {
  1386 +
  1387 + // Compute indices
  1388 + int twiddle = index + k - 1;
  1389 +
  1390 + // If sign==-1 then set trigonometric coefficients for forward transform
  1391 + if (sign == -1) {
  1392 + for (int i = 0; i < n; ++i, twiddle += q) {
  1393 + w[i] = wavetable[twiddle];
  1394 + }
  1395 + }
  1396 +
  1397 + // ... otherwise set trigonometric coefficients for backward tranform
  1398 + // w -> conjugate(w)
  1399 + else {
  1400 + for (int i = 0; i < n; ++i, twiddle += q) {
  1401 + w[i] = std::conj(wavetable[twiddle]);
  1402 + }
  1403 + }
  1404 +
  1405 + } // endelse: k != 0
  1406 +
  1407 + // Return w vectors
  1408 + return w;
  1409 +}
... ...