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}