-
Notifications
You must be signed in to change notification settings - Fork 36
Expand file tree
/
Copy pathFourierTransform.cc
More file actions
154 lines (143 loc) · 4.62 KB
/
Copy pathFourierTransform.cc
File metadata and controls
154 lines (143 loc) · 4.62 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
// -*- c++ -*-
/*
* Copyright (c) 2010-2012, Jim Bosch
* All rights reserved.
*
* ndarray is distributed under a simple BSD-like license;
* see the LICENSE file that should be present in the root
* of the source distribution, or alternately available at:
* https://github.com/ndarray/ndarray
*/
#include "ndarray/fft/FFTWTraits.h"
#include "ndarray/fft/FourierTransform.h"
namespace ndarray {
template <typename T, int N>
template <int M>
Array<typename FourierTransform<T,N>::ElementX,M,M>
FourierTransform<T,N>::initializeX(Vector<Size,M> const & shape) {
OwnerX xOwner = detail::FFTWTraits<T>::allocateX(shape.product());
return Array<ElementX,M,M>(external(xOwner.get(), shape, ROW_MAJOR, xOwner));
}
template <typename T, int N>
template <int M>
Array<typename FourierTransform<T,N>::ElementK,M,M>
FourierTransform<T,N>::initializeK(Vector<Size,M> const & shape) {
Vector<Size,M> kShape(shape);
kShape[M-1] = detail::FourierTraits<T>::computeLastDimensionSize(shape[M-1]);
OwnerK kOwner = detail::FFTWTraits<T>::allocateK(kShape.product());
return Array<ElementK,M,M>(external(kOwner.get(), kShape, ROW_MAJOR, kOwner));
}
template <typename T, int N>
template <int M>
void
FourierTransform<T,N>::initialize(
Vector<Size,M> const & shape,
Array<ElementX,M,M> & x,
Array<ElementK,M,M> & k
) {
if (x.empty()) x = initializeX(shape);
if (k.empty()) k = initializeK(shape);
NDARRAY_ASSERT(x.getShape() == shape);
NDARRAY_ASSERT(std::equal(shape.begin(), shape.end()-1, k.getShape().begin()));
}
template <typename T, int N>
typename FourierTransform<T,N>::Ptr
FourierTransform<T,N>::planForward(
Index const & shape,
typename FourierTransform<T,N>::ArrayX & x,
typename FourierTransform<T,N>::ArrayK & k
) {
initialize(shape,x,k);
Vector<int,N> s = shape.template cast<int>();
return Ptr(
new FourierTransform(
detail::FFTWTraits<T>::forward(
N, s.begin(), 1,
x.getData(), NULL, 1, 0,
k.getData(), NULL, 1, 0,
FFTW_MEASURE | FFTW_DESTROY_INPUT
),
x.getManager(),
k.getManager()
)
);
}
template <typename T, int N>
typename FourierTransform<T,N>::Ptr
FourierTransform<T,N>::planInverse(
Index const & shape,
typename FourierTransform<T,N>::ArrayK & k,
typename FourierTransform<T,N>::ArrayX & x
) {
initialize(shape,x,k);
Vector<int,N> s = shape.template cast<int>();
return Ptr(
new FourierTransform(
detail::FFTWTraits<T>::inverse(
N, s.begin(), 1,
k.getData(), NULL, 1, 0,
x.getData(), NULL, 1, 0,
FFTW_MEASURE | FFTW_DESTROY_INPUT
),
x.getManager(),
k.getManager()
)
);
}
template <typename T, int N>
typename FourierTransform<T,N>::Ptr
FourierTransform<T,N>::planMultiplexForward(
MultiplexIndex const & shape,
typename FourierTransform<T,N>::MultiplexArrayX & x,
typename FourierTransform<T,N>::MultiplexArrayK & k
) {
initialize(shape,x,k);
Vector<int,N+1> s = shape.template cast<int>();
return Ptr(
new FourierTransform(
detail::FFTWTraits<T>::forward(
N, s.begin()+1, s[0],
x.getData(), NULL, 1, x.template getStride<0>(),
k.getData(), NULL, 1, k.template getStride<0>(),
FFTW_MEASURE | FFTW_DESTROY_INPUT
),
x.getManager(),
k.getManager()
)
);
}
template <typename T, int N>
typename FourierTransform<T,N>::Ptr
FourierTransform<T,N>::planMultiplexInverse(
MultiplexIndex const & shape,
typename FourierTransform<T,N>::MultiplexArrayK & k,
typename FourierTransform<T,N>::MultiplexArrayX & x
) {
initialize(shape,x,k);
Vector<int,N+1> s = shape.template cast<int>();
return Ptr(
new FourierTransform(
detail::FFTWTraits<T>::inverse(
N, s.begin()+1, s[0],
k.getData(), NULL, 1, k.template getStride<0>(),
x.getData(), NULL, 1, x.template getStride<0>(),
FFTW_MEASURE | FFTW_DESTROY_INPUT
),
x.getManager(),
k.getManager()
)
);
}
template <typename T, int N>
void FourierTransform<T,N>::execute() {
detail::FFTWTraits<T>::execute(
reinterpret_cast<typename detail::FFTWTraits<T>::Plan>(_plan)
);
}
template <typename T, int N>
FourierTransform<T,N>::~FourierTransform() {
detail::FFTWTraits<T>::destroy(
reinterpret_cast<typename detail::FFTWTraits<T>::Plan>(_plan)
);
}
} // namespace ndarray