001/* 002Copyright 2006 Jerry Huxtable 003 004Licensed under the Apache License, Version 2.0 (the "License"); 005you may not use this file except in compliance with the License. 006You may obtain a copy of the License at 007 008 http://www.apache.org/licenses/LICENSE-2.0 009 010Unless required by applicable law or agreed to in writing, software 011distributed under the License is distributed on an "AS IS" BASIS, 012WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 013See the License for the specific language governing permissions and 014limitations under the License. 015*/ 016 017package com.jhlabs.math; 018 019public class FFT { 020 021 // Weighting factors 022 protected float[] w1; 023 protected float[] w2; 024 protected float[] w3; 025 026 public FFT( int logN ) { 027 // Prepare the weighting factors 028 w1 = new float[logN]; 029 w2 = new float[logN]; 030 w3 = new float[logN]; 031 int N = 1; 032 for ( int k = 0; k < logN; k++ ) { 033 N <<= 1; 034 double angle = -2.0 * Math.PI / N; 035 w1[k] = (float)Math.sin(0.5 * angle); 036 w2[k] = -2.0f * w1[k] * w1[k]; 037 w3[k] = (float)Math.sin(angle); 038 } 039 } 040 041 private void scramble( int n, float[] real, float[] imag ) { 042 int j = 0; 043 044 for ( int i = 0; i < n; i++ ) { 045 if ( i > j ) { 046 float t; 047 t = real[j]; 048 real[j] = real[i]; 049 real[i] = t; 050 t = imag[j]; 051 imag[j] = imag[i]; 052 imag[i] = t; 053 } 054 int m = n >> 1; 055 while (j >= m && m >= 2) { 056 j -= m; 057 m >>= 1; 058 } 059 j += m; 060 } 061 } 062 063 private void butterflies( int n, int logN, int direction, float[] real, float[] imag ) { 064 int N = 1; 065 066 for ( int k = 0; k < logN; k++ ) { 067 float w_re, w_im, wp_re, wp_im, temp_re, temp_im, wt; 068 int half_N = N; 069 N <<= 1; 070 wt = direction * w1[k]; 071 wp_re = w2[k]; 072 wp_im = direction * w3[k]; 073 w_re = 1.0f; 074 w_im = 0.0f; 075 for ( int offset = 0; offset < half_N; offset++ ) { 076 for( int i = offset; i < n; i += N ) { 077 int j = i + half_N; 078 float re = real[j]; 079 float im = imag[j]; 080 temp_re = (w_re * re) - (w_im * im); 081 temp_im = (w_im * re) + (w_re * im); 082 real[j] = real[i] - temp_re; 083 real[i] += temp_re; 084 imag[j] = imag[i] - temp_im; 085 imag[i] += temp_im; 086 } 087 wt = w_re; 088 w_re = wt * wp_re - w_im * wp_im + w_re; 089 w_im = w_im * wp_re + wt * wp_im + w_im; 090 } 091 } 092 if ( direction == -1 ) { 093 float nr = 1.0f / n; 094 for ( int i = 0; i < n; i++ ) { 095 real[i] *= nr; 096 imag[i] *= nr; 097 } 098 } 099 } 100 101 public void transform1D( float[] real, float[] imag, int logN, int n, boolean forward ) { 102 scramble( n, real, imag ); 103 butterflies( n, logN, forward ? 1 : -1, real, imag ); 104 } 105 106 public void transform2D( float[] real, float[] imag, int cols, int rows, boolean forward ) { 107 int log2cols = log2(cols); 108 int log2rows = log2(rows); 109 int n = Math.max(rows, cols); 110 float[] rtemp = new float[n]; 111 float[] itemp = new float[n]; 112 113 // FFT the rows 114 for ( int y = 0; y < rows; y++ ) { 115 int offset = y*cols; 116 System.arraycopy( real, offset, rtemp, 0, cols ); 117 System.arraycopy( imag, offset, itemp, 0, cols ); 118 transform1D(rtemp, itemp, log2cols, cols, forward); 119 System.arraycopy( rtemp, 0, real, offset, cols ); 120 System.arraycopy( itemp, 0, imag, offset, cols ); 121 } 122 123 // FFT the columns 124 for ( int x = 0; x < cols; x++ ) { 125 int index = x; 126 for ( int y = 0; y < rows; y++ ) { 127 rtemp[y] = real[index]; 128 itemp[y] = imag[index]; 129 index += cols; 130 } 131 transform1D(rtemp, itemp, log2rows, rows, forward); 132 index = x; 133 for ( int y = 0; y < rows; y++ ) { 134 real[index] = rtemp[y]; 135 imag[index] = itemp[y]; 136 index += cols; 137 } 138 } 139 } 140 141 private int log2( int n ) { 142 int m = 1; 143 int log2n = 0; 144 145 while (m < n) { 146 m *= 2; 147 log2n++; 148 } 149 return m == n ? log2n : -1; 150 } 151 152}