1 // ******************************************************************************
2 //
3 // Title: Force Field X.
4 // Description: Force Field X - Software for Molecular Biophysics.
5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2024.
6 //
7 // This file is part of Force Field X.
8 //
9 // Force Field X is free software; you can redistribute it and/or modify it
10 // under the terms of the GNU General Public License version 3 as published by
11 // the Free Software Foundation.
12 //
13 // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16 // details.
17 //
18 // You should have received a copy of the GNU General Public License along with
19 // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20 // Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // Linking this library statically or dynamically with other modules is making a
23 // combined work based on this library. Thus, the terms and conditions of the
24 // GNU General Public License cover the whole combination.
25 //
26 // As a special exception, the copyright holders of this library give you
27 // permission to link this library with independent modules to produce an
28 // executable, regardless of the license terms of these independent modules, and
29 // to copy and distribute the resulting executable under terms of your choice,
30 // provided that you also meet, for each linked independent module, the terms
31 // and conditions of the license of that module. An independent module is a
32 // module which is not derived from or based on this library. If you modify this
33 // library, you may extend this exception to your version of the library, but
34 // you are not obligated to do so. If you do not wish to do so, delete this
35 // exception statement from your version.
36 //
37 // ******************************************************************************
38 package ffx.numerics.fft;
39
40 import static java.lang.System.arraycopy;
41 import static org.apache.commons.math3.util.FastMath.PI;
42 import static org.apache.commons.math3.util.FastMath.cos;
43 import static org.apache.commons.math3.util.FastMath.sin;
44
45 /**
46 * Compute the FFT of real, double precision data of arbitrary length n using a Complex transform.
47 *
48 * @author Michal J. Schnieders<br> Derived from: <br> Bruce R. Miller bruce.miller@nist.gov <br>
49 * Contribution of the National Institute of Standards and Technology, not subject to copyright.
50 * <br>
51 * Derived from:<br> GSL (Gnu Scientific Library) FFT Code by Brian Gough bjg@vvv.lanl.gov
52 * @see <ul>
53 * <li><a href="http://dx.doi.org/10.1109/TASSP.1987.1165220" target="_blank"> Henrik V.
54 * Sorenson, Douglas L. Jones, Michael T. Heideman, and C. Sidney Burrus. Real-valued fast
55 * fourier fft algorithms. IEEE Transactions on Acoustics, Speech, and Signal Processing,
56 * ASSP-35(6):849–863, 1987. </a>
57 * <li><a href="http://www.jstor.org/stable/2003354" target="_blank">J. W. Cooley and J. W.
58 * Tukey, Mathematics of Computation 19 (90), 297 (1965) </a>
59 * <li><a href="http://en.wikipedia.org/wiki/Fast_Fourier_transform" target="_blank">FFT at
60 * Wikipedia </a>
61 * </ul>
62 * @since 1.0
63 */
64 public class Real {
65
66 private final Complex complexFFT;
67 private final int n;
68 private final int halfN;
69 private final double cosTheta;
70 private final double sinTheta;
71 private final double[] work;
72
73 /**
74 * Constructs a Complex FFT of length (n / 2) for real data of length n.
75 *
76 * @param n the length of the real data.
77 */
78 public Real(int n) {
79 this.n = n;
80 halfN = n / 2;
81 double theta = PI / halfN;
82 cosTheta = cos(theta);
83 sinTheta = sin(theta);
84 complexFFT = new Complex(halfN);
85 work = new double[n + 2];
86 }
87
88 /**
89 * fft
90 *
91 * @param data Input data.
92 * @param offset Offset to the beginning of the data.
93 */
94 public void fft(double[] data, int offset) {
95 complexFFT.fft(data, offset, 2);
96 unpack(data, offset);
97 }
98
99 /**
100 * ifft
101 *
102 * @param data Input data.
103 * @param offset Offset to the beginning of the data.
104 */
105 public void ifft(double[] data, int offset) {
106 pack(data, offset);
107 complexFFT.ifft(data, offset, 2);
108 // Renormalize the 1/2 length fft.
109 int ii = offset;
110 for (int i = 0; i < n; i++) {
111 data[ii++] *= 2.0;
112 }
113 }
114
115 /**
116 * inverse
117 *
118 * @param data Input data.
119 * @param offset Offset to the beginning of the data.
120 */
121 public void inverse(double[] data, int offset) {
122 ifft(data, offset);
123
124 // normalize inverse FFT with 1/n.
125 double norm = normalization();
126 for (int i = 0; i < n; i++) {
127 final int index = offset + i;
128 data[index] *= norm;
129 }
130 }
131
132 /**
133 * Return the normalization factor. Multiply the elements of the back-transformed data to get the
134 * normalized inverse.
135 *
136 * @return a double.
137 */
138 private double normalization() {
139 return 1.0 / n;
140 }
141
142 /**
143 * Unpack following the forward Complex FFT.
144 *
145 * @param data Input data.
146 * @param offset Offset to the beginning of the data.
147 */
148 private void unpack(double[] data, int offset) {
149 arraycopy(data, offset, work, 0, n);
150 double wrs = cosTheta;
151 double wis = -sinTheta;
152 final double d0 = work[0];
153 final double d1 = work[1];
154 final int on = offset + n;
155 final int o1 = offset + 1;
156 data[offset] = d0 + d1;
157 data[o1] = 0.0;
158 data[on] = d0 - d1;
159 data[on + 1] = 0.0;
160 double wr = wrs;
161 double wi = wis;
162 for (int i = 1; i < halfN; i++) {
163 int i1 = 2 * i;
164 double rk = work[i1];
165 double ik = work[i1 + 1];
166 int i2 = n - 2 * i;
167 double rn = work[i2];
168 double in = work[i2 + 1];
169 double s1r = 0.5 * (rk + rn);
170 double s1i = 0.5 * (ik - in);
171 double s2r = 0.5 * (ik + in);
172 double s2i = 0.5 * (rn - rk);
173 double fr = s1r + wr * s2r - wi * s2i;
174 double fi = s1i + wi * s2r + wr * s2i;
175 data[i1 + offset] = fr;
176 data[i1 + 1 + offset] = fi;
177 double wrp = wr;
178 wr = wr * wrs - wi * wis;
179 wi = wrp * wis + wi * wrs;
180 }
181 }
182
183 /**
184 * Pack prior to inverse Complex FFT.
185 *
186 * @param data Input data.
187 * @param offset Offset to the beginning of the data.
188 */
189 private void pack(double[] data, int offset) {
190 arraycopy(data, offset, work, 0, n + 2);
191 double wrs = cosTheta;
192 double wis = sinTheta;
193 double wr = 1.0;
194 double wi = 0.0;
195 for (int i = 0; i < halfN; i++) {
196 int i1 = 2 * i;
197 double fkr = work[i1];
198 double fki = work[i1 + 1];
199 int i2 = n - 2 * i;
200 double fnr = work[i2];
201 double fni = -work[i2 + 1];
202 double s1r = fkr + fnr;
203 double s1i = fki + fni;
204 double s2r = fkr - fnr;
205 double s2i = fki - fni;
206 double fr = s2r * wr - s2i * wi;
207 double fi = s2r * wi + s2i * wr;
208 data[i1 + offset] = 0.5 * (s1r - fi);
209 data[i1 + 1 + offset] = 0.5 * (s1i + fr);
210 double wrp = wr;
211 wr = wr * wrs - wi * wis;
212 wi = wrp * wis + wi * wrs;
213 }
214 }
215 }