001    /*
002    *
003    
004    Licensed under the Apache License, Version 2.0 (the "License");
005    you may not use this file except in compliance with the License.
006    You may obtain a copy of the License at
007    
008       http://www.apache.org/licenses/LICENSE-2.0
009    
010    Unless required by applicable law or agreed to in writing, software
011    distributed under the License is distributed on an "AS IS" BASIS,
012    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013    See the License for the specific language governing permissions and
014    limitations under the License.
015    */
016    
017    package railo.runtime.img.math;
018    
019    public 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    }