1 #ifndef vnl_fft_base_hxx_
2 #define vnl_fft_base_hxx_
3 /*
4 fsm
5 */
6 #include "vnl_fft_base.h"
7 #include <vnl/algo/vnl_fft.h>
8 #include <cassert>
9 #ifdef _MSC_VER
10 # include <vcl_msvc_warnings.h>
11 #endif
12
13 template <int D, class T>
transform(std::complex<T> * signal,int dir)14 void vnl_fft_base<D, T>::transform(std::complex<T> *signal, int dir)
15 {
16 assert((dir == +1) || (dir == -1));
17
18 // transform along each dimension, i, in turn.
19 for (int i=0; i<D; ++i) {
20 int N1 = 1; // n[0] n[1] ... n[i-1]
21 int N2 = 1; // n[i]
22 int N3 = 1; // n[i+1] n[i+2] ... n[D-1]
23 for (int j=0; j<D; ++j) {
24 int d = factors_[j].number();
25 if (j < i) N1 *= d;
26 if (j == i) N2 *= d;
27 if (j > i) N3 *= d;
28 }
29
30 // pretend the signal is N1xN2xN3. we want to transform
31 // along the second dimension.
32 for (int n1=0; n1<N1; ++n1) {
33 // FIXME: we could avoid one loop by using the LOT parameter
34 // but it's not entirely clear that would save us anything.
35
36 for (int n3=0; n3<N3; ++n3) {
37 // This relies on the assumption that std::complex<T> is layout
38 // compatible with "struct { T real; T imag; }". It is probably
39 // a valid assumption for all sane C++ libraries.
40 T *data = (T *) (signal + n1*N2*N3 + n3);
41
42 long info = 0;
43 vnl_fft_gpfa (/* A */ data,
44 /* B */ data + 1,
45 /* TRIGS */ factors_[i].trigs (),
46 /* INC */ 2*N3,
47 /* JUMP */ 0,
48 /* N */ N2,
49 /* LOT */ 1,
50 /* ISIGN */ dir,
51 /* NIPQ */ factors_[i].pqr (),
52 /* INFO */ &info);
53 assert(info != -1);
54 }
55 }
56 }
57 }
58
59 #undef VNL_FFT_BASE_INSTANTIATE
60 #define VNL_FFT_BASE_INSTANTIATE(D, T) \
61 template struct VNL_ALGO_EXPORT vnl_fft_base<D, T >
62
63 #endif
64