FourierMath


1   package JSci.maths;
2   
3   /**
4   * The Fourier math library.
5   * This class cannot be subclassed or instantiated because all methods are static.
6   * Use <code>sort(transform(sort(...)))</code> for the discrete analogue of the continuous Fourier transform,
7   * and <code>sort(inverseTransform(sort(...)))</code> for the inverse transform.
8   * @jsci.planetmath FourierTransform
9   * @version 0.9
10  * @author Mark Hale
11  */
12  public final class FourierMath extends AbstractMath implements NumericalConstants {
13          private FourierMath() {}
14  
15          /**
16          * Fourier transform.
17          * @return an array containing positive frequencies in ascending order
18          * followed by negative frequencies in ascending order.
19          * @author Don Cross
20          */
21          public static Complex[] transform(final Complex data[]) {
22                  final int N=data.length;
23                  if(!isPowerOf2(N))
24                          throw new IllegalArgumentException  ("The number of samples must be a power of 2.");
25  
26                  final double arrayRe[]=new double[N];
27                  final double arrayIm[]=new double[N];
28  
29                  final int numBits=numberOfBitsNeeded(N);
30  // Simultaneous data copy and bit-reversal ordering into output
31                  for(int i=0;i<N;i++) {
32                          final int j=reverseBits(i,numBits);
33                          arrayRe[j]=data[i].real();
34                          arrayIm[j]=data[i].imag();
35                  }
36  // FFT
37                  fft(arrayRe, arrayIm, TWO_PI);
38  
39                  final Complex answer[]=new Complex[N];
40                  for(int i=0;i<N;i++)
41                          answer[i]=new Complex(arrayRe[i],arrayIm[i]);
42                  return answer;
43          }
44          /**
45          * Fourier transform.
46          * @return an array containing positive frequencies in ascending order
47          * followed by negative frequencies in ascending order.
48          * @author Don Cross
49          */
50          public static Complex[] transform(final double dataReal[], final double dataImag[]) {
51                  final int N=dataReal.length;
52                  if(!isPowerOf2(N))
53                          throw new IllegalArgumentException  ("The number of samples must be a power of 2.");
54  
55                  final double arrayRe[]=new double[N];
56                  final double arrayIm[]=new double[N];
57  
58                  final int numBits=numberOfBitsNeeded(N);
59  // Simultaneous data copy and bit-reversal ordering into output
60                  for(int i=0;i<N;i++) {
61                          final int j=reverseBits(i,numBits);
62                          arrayRe[j]=dataReal[i];
63                          arrayIm[j]=dataImag[i];
64                  }
65  // FFT
66                  fft(arrayRe, arrayIm, TWO_PI);
67  
68                  final Complex answer[]=new Complex[N];
69                  for(int i=0;i<N;i++)
70                          answer[i]=new Complex(arrayRe[i],arrayIm[i]);
71                  return answer;
72          }
73          /**
74          * Fourier transform.
75          * @return an array containing positive frequencies in ascending order
76          * followed by negative frequencies in ascending order.
77          */
78          public static Complex[] transform(final double data[]) {
79                  final int N=data.length;
80                  if(!isPowerOf2(N))
81                          throw new IllegalArgumentException  ("The number of samples must be a power of 2.");
82  
83                  final double arrayRe[]=new double[N];
84                  final double arrayIm[]=new double[N];
85  
86                  final int numBits=numberOfBitsNeeded(N);
87  // Simultaneous data copy and bit-reversal ordering into output
88                  for(int i=0;i<N;i++) {
89                          final int j=reverseBits(i,numBits);
90                          arrayRe[j]=data[i];
91                  }
92  // FFT
93                  fft(arrayRe, arrayIm, TWO_PI);
94  
95                  final Complex answer[]=new Complex[N];
96                  for(int i=0;i<N;i++)
97                          answer[i]=new Complex(arrayRe[i],arrayIm[i]);
98                  return answer;
99          }
100         /**
101         * Inverse Fourier transform.
102         * @return an array containing the positive time part of the signal
103         * followed by the negative time part.
104         * @author Don Cross
105         */
106         public static Complex[] inverseTransform(final Complex data[]) {
107                 final int N=data.length;
108                 if(!isPowerOf2(N))
109                         throw new IllegalArgumentException  ("Data length must be a power of 2.");
110 
111                 final double arrayRe[]=new double[N];
112                 final double arrayIm[]=new double[N];
113 
114                 final int numBits=numberOfBitsNeeded(N);
115 // Simultaneous data copy and bit-reversal ordering into output
116                 for(int i=0;i<N;i++) {
117                         final int j=reverseBits(i,numBits);
118                         arrayRe[j]=data[i].real();
119                         arrayIm[j]=data[i].imag();
120                 }
121 // inverse FFT
122                 fft(arrayRe, arrayIm, -TWO_PI);
123 // Normalize
124                 final Complex answer[]=new Complex[N];
125                 final double denom=N;
126                 for(int i=0;i<N;i++)
127                         answer[i]=new Complex(arrayRe[i]/denom,arrayIm[i]/denom);
128                 return answer;
129         }
130         /**
131         * Inverse Fourier transform.
132         * @return an array containing the positive time part of the signal
133         * followed by the negative time part.
134         * @author Don Cross
135         */
136         public static Complex[] inverseTransform(final double dataReal[], final double dataImag[]) {
137                 final int N=dataReal.length;
138                 if(!isPowerOf2(N))
139                         throw new IllegalArgumentException  ("Data length must be a power of 2.");
140 
141                 final double arrayRe[]=new double[N];
142                 final double arrayIm[]=new double[N];
143 
144                 final int numBits=numberOfBitsNeeded(N);
145 // Simultaneous data copy and bit-reversal ordering into output
146                 for(int i=0;i<N;i++) {
147                         final int j=reverseBits(i,numBits);
148                         arrayRe[j]=dataReal[i];
149                         arrayIm[j]=dataImag[i];
150                 }
151 // inverse FFT
152                 fft(arrayRe, arrayIm, -TWO_PI);
153 // Normalize
154                 final Complex answer[]=new Complex[N];
155                 final double denom=N;
156                 for(int i=0;i<N;i++)
157                         answer[i]=new Complex(arrayRe[i]/denom,arrayIm[i]/denom);
158                 return answer;
159         }
160         /**
161         * Inverse Fourier transform.
162         * @return an array containing the positive time part of the signal
163         * followed by the negative time part.
164         */
165         public static Complex[] inverseTransform(final double data[]) {
166                 final int N=data.length;
167                 if(!isPowerOf2(N))
168                         throw new IllegalArgumentException  ("Data length must be a power of 2.");
169 
170                 final double arrayRe[]=new double[N];
171                 final double arrayIm[]=new double[N];
172 
173                 final int numBits=numberOfBitsNeeded(N);
174 // Simultaneous data copy and bit-reversal ordering into output
175                 for(int i=0;i<N;i++) {
176                         final int j=reverseBits(i,numBits);
177                         arrayRe[j]=data[i];
178                 }
179 // inverse FFT
180                 fft(arrayRe, arrayIm, -TWO_PI);
181 // Normalize
182                 final Complex answer[]=new Complex[N];
183                 final double denom=N;
184                 for(int i=0;i<N;i++)
185                         answer[i]=new Complex(arrayRe[i]/denom,arrayIm[i]/denom);
186                 return answer;
187         }
188         /**
189         * Common FFT code.
190         * @param twoPi TWO_PI for transform, -TWO_PI for inverse transform.
191         */
192         private static void fft(double arrayRe[], double arrayIm[], final double twoPi) {
193                 final int N=arrayRe.length;
194                 int blockEnd=1;
195                 for(int blockSize=2; blockSize<=N; blockSize<<=1) {
196                         final double deltaAngle = twoPi/blockSize;
197                         double alpha = Math.sin(0.5*deltaAngle);
198                         alpha = 2.0*alpha*alpha;
199                         final double beta = Math.sin(deltaAngle);
200                         for(int i=0; i<N; i+=blockSize) {
201                                 double angleRe=1.0;
202                                 double angleIm=0.0;
203                                 for(int j=i,n=0; n<blockEnd; j++,n++) {
204                                         final int k = j+blockEnd;
205                                         // tmp = angle*array[k]
206                                         double tmpRe = angleRe*arrayRe[k]-angleIm*arrayIm[k];
207                                         double tmpIm = angleRe*arrayIm[k]+angleIm*arrayRe[k];
208                                         arrayRe[k] = arrayRe[j]-tmpRe;
209                                         arrayIm[k] = arrayIm[j]-tmpIm;
210                                         arrayRe[j] += tmpRe;
211                                         arrayIm[j] += tmpIm;
212                                         // angle = angle - (a-bi)*angle
213                                         tmpRe = alpha*angleRe + beta*angleIm;
214                                         tmpIm = alpha*angleIm - beta*angleRe;
215                     angleRe -= tmpRe;
216                     angleIm -= tmpIm;
217                                 }
218                         }
219                         blockEnd=blockSize;
220                 }
221         }
222         /**
223         * Returns true if x is a power of 2.
224         * @author Don Cross
225         */
226         private static boolean isPowerOf2(final int x) {
227                 final int BITS_PER_WORD=32;
228                 for(int i=1,y=2; i<BITS_PER_WORD; i++,y<<=1) {
229                         if(x==y)
230                                 return true;
231                 }
232                 return false;
233         }
234         /**
235         * Number of bits needed.
236         * @author Don Cross
237         */
238         private static int numberOfBitsNeeded(final int pwrOf2) {
239                 if(pwrOf2<2)
240                         throw new IllegalArgumentException  ();
241                 for(int i=0;;i++) {
242                         if((pwrOf2&(1<<i))>0)
243                                 return i;
244                 }
245         }
246         /**
247         * Reverse bits.
248         * @author Don Cross
249         */
250         private static int reverseBits(int index,final int numBits) {
251                 int i,rev;
252                 for(i=rev=0;i<numBits;i++) {
253                         rev=(rev<<1)|(index&1);
254                         index>>=1;
255                 }
256                 return rev;
257         }
258 
259         /**
260         * Sorts the output from the Fourier transfom methods into
261         * ascending frequency/time order.
262         */
263         public static Complex[] sort(final Complex output[]) {
264                 final Complex ret[]=new Complex[output.length];
265                 final int Nby2=output.length/2;
266                 for(int i=0;i<Nby2;i++) {
267                         ret[Nby2+i]=output[i];
268                         ret[i]=output[Nby2+i];
269                 }
270                 return ret;
271         }
272         public static double[] sort(final double input[]) {
273                 final double ret[]=new double[input.length];
274                 final int Nby2=input.length/2;
275                 for(int i=0;i<Nby2;i++) {
276                         ret[Nby2+i]=input[i];
277                         ret[i]=input[Nby2+i];
278                 }
279                 return ret;
280         }
281 }
282
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags