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 package ffx.numerics.fft;
39
40
41
42
43
44
45
46
47 public class Real3D {
48
49 private final int nextX, nextY, nextZ;
50 private final int nX, nY, nZ;
51 private final int nX1, nZ2;
52 private final double[] work;
53 private final double[] recip;
54 private final Real fftX;
55 private final Complex fftY, fftZ;
56
57
58
59
60
61
62
63
64 public Real3D(int nX, int nY, int nZ) {
65 this.nX = nX / 2;
66 this.nY = nY;
67 this.nZ = nZ;
68 nX1 = this.nX + 1;
69 nZ2 = 2 * nZ;
70 nextX = 2;
71 nextY = nX + 2;
72 nextZ = nextY * nY;
73 work = new double[nZ2];
74 recip = new double[nX1 * nY * nZ];
75 fftX = new Real(nX);
76 fftY = new Complex(nY);
77 fftZ = new Complex(nZ);
78 }
79
80
81
82
83
84
85
86
87
88
89
90 public static int iReal3D(int i, int j, int k, int nX, int nY) {
91 int xSide = nX + 2;
92 int xySlice = xSide * nY;
93 return i + j * xSide + k * xySlice;
94 }
95
96
97
98
99
100
101 public void convolution(final double[] input) {
102 for (int z = 0; z < nZ; z++) {
103 for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
104 fftX.fft(input, offset);
105 }
106 for (int offset = z * nextZ, x = 0; x < nX1; x++, offset += nextX) {
107 fftY.fft(input, offset, nextY);
108 }
109 }
110 for (int index = 0, offset = 0, y = 0; y < nY; y++) {
111 for (int x = 0; x < nX1; x++, offset += nextX) {
112 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
113 work[i] = input[z];
114 work[i + 1] = input[z + 1];
115 }
116 fftZ.fft(work, 0, 2);
117 for (int i = 0; i < nZ2; i += 2) {
118 double r = recip[index++];
119 work[i] *= r;
120 work[i + 1] *= r;
121 }
122 fftZ.ifft(work, 0, 2);
123 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
124 input[z] = work[i];
125 input[z + 1] = work[i + 1];
126 }
127 }
128 }
129 for (int z = 0; z < nZ; z++) {
130 for (int offset = z * nextZ, x = 0; x < nX1; x++, offset += nextX) {
131 fftY.ifft(input, offset, nextY);
132 }
133 for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
134 fftX.ifft(input, offset);
135 }
136 }
137 }
138
139
140
141
142
143
144 public void fft(final double[] input) {
145 for (int z = 0; z < nZ; z++) {
146 for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
147 fftX.fft(input, offset);
148 }
149 for (int offset = z * nextZ, x = 0; x < nX1; x++, offset += nextX) {
150 fftY.fft(input, offset, nextY);
151 }
152 }
153 for (int x = 0; x < nX1; x++) {
154 for (int offset = x * 2, y = 0; y < nY; y++, offset += nextY) {
155 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
156 work[i] = input[z];
157 work[i + 1] = input[z + 1];
158 }
159 fftZ.fft(work, 0, 2);
160 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
161 input[z] = work[i];
162 input[z + 1] = work[i + 1];
163 }
164 }
165 }
166 }
167
168
169
170
171
172
173 public void ifft(final double[] input) {
174 for (int x = 0; x < nX1; x++) {
175 for (int offset = x * 2, y = 0; y < nY; y++, offset += nextY) {
176 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
177 work[i] = input[z];
178 work[i + 1] = input[z + 1];
179 }
180 fftZ.ifft(work, 0, 2);
181 for (int i = 0, z = offset; i < nZ2; i += 2, z += nextZ) {
182 input[z] = work[i];
183 input[z + 1] = work[i + 1];
184 }
185 }
186 }
187 for (int z = 0; z < nZ; z++) {
188 for (int offset = z * nextZ, x = 0; x < nX1; x++, offset += nextX) {
189 fftY.ifft(input, offset, nextY);
190 }
191 for (int offset = z * nextZ, y = 0; y < nY; y++, offset += nextY) {
192 fftX.ifft(input, offset);
193 }
194 }
195 }
196
197
198
199
200
201
202 public void setRecip(double[] recip) {
203
204 for (int index = 0, offset = 0, y = 0; y < nY; y++) {
205 for (int x = 0; x < nX1; x++, offset += 1) {
206 for (int i = 0, z = offset; i < nZ; i++, z += nX1 * nY) {
207 this.recip[index++] = recip[z];
208 }
209 }
210 }
211 }
212 }