001/**
002 *
003 * Copyright (c) 2014, the Railo Company Ltd. All rights reserved.
004 *
005 * This library is free software; you can redistribute it and/or
006 * modify it under the terms of the GNU Lesser General Public
007 * License as published by the Free Software Foundation; either 
008 * version 2.1 of the License, or (at your option) any later version.
009 * 
010 * This library is distributed in the hope that it will be useful,
011 * but WITHOUT ANY WARRANTY; without even the implied warranty of
012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
013 * Lesser General Public License for more details.
014 * 
015 * You should have received a copy of the GNU Lesser General Public 
016 * License along with this library.  If not, see <http://www.gnu.org/licenses/>.
017 * 
018 **/
019/*
020*
021
022Licensed under the Apache License, Version 2.0 (the "License");
023you may not use this file except in compliance with the License.
024You may obtain a copy of the License at
025
026   http://www.apache.org/licenses/LICENSE-2.0
027
028Unless required by applicable law or agreed to in writing, software
029distributed under the License is distributed on an "AS IS" BASIS,
030WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
031See the License for the specific language governing permissions and
032limitations under the License.
033*/
034
035package lucee.runtime.img.math;
036
037public class FFT {
038
039    // Weighting factors
040    protected float[] w1;
041    protected float[] w2;
042    protected float[] w3;
043
044    public FFT( int logN ) {
045        // Prepare the weighting factors
046        w1 = new float[logN];
047        w2 = new float[logN];
048        w3 = new float[logN];
049        int N = 1;
050        for ( int k = 0; k < logN; k++ ) {
051            N <<= 1;
052            double angle = -2.0 * Math.PI / N;
053            w1[k] = (float)Math.sin(0.5 * angle);
054            w2[k] = -2.0f * w1[k] * w1[k];
055            w3[k] = (float)Math.sin(angle);
056        }
057    }
058
059    private void scramble( int n, float[] real, float[] imag ) {
060        int j = 0;
061
062        for ( int i = 0; i < n; i++ ) {
063            if ( i > j ) {
064                float t;
065                t = real[j];
066                real[j] = real[i];
067                real[i] = t;
068                t = imag[j];
069                imag[j] = imag[i];
070                imag[i] = t;
071            }
072            int m = n >> 1;
073            while (j >= m && m >= 2) {
074                j -= m;
075                m >>= 1;
076            }
077            j += m;
078        }
079    }
080
081    private void butterflies( int n, int logN, int direction, float[] real, float[] imag ) {
082        int N = 1;
083
084        for ( int k = 0; k < logN; k++ ) {
085            float w_re, w_im, wp_re, wp_im, temp_re, temp_im, wt;
086            int half_N = N;
087            N <<= 1;
088            wt = direction * w1[k];
089            wp_re = w2[k];
090            wp_im = direction * w3[k];
091            w_re = 1.0f;
092            w_im = 0.0f;
093            for ( int offset = 0; offset < half_N; offset++ ) {
094                for( int i = offset; i < n; i += N ) {
095                    int j = i + half_N;
096                    float re = real[j];
097                    float im = imag[j];
098                    temp_re = (w_re * re) - (w_im * im);
099                    temp_im = (w_im * re) + (w_re * im);
100                    real[j] = real[i] - temp_re;
101                    real[i] += temp_re;
102                    imag[j] = imag[i] - temp_im;
103                    imag[i] += temp_im;
104                }
105                wt = w_re;
106                w_re = wt * wp_re - w_im * wp_im + w_re;
107                w_im = w_im * wp_re + wt * wp_im + w_im;
108            }
109        }
110        if ( direction == -1 ) {
111            float nr = 1.0f / n;
112            for ( int i = 0; i < n; i++ ) {
113                real[i] *= nr;
114                imag[i] *= nr;
115            }
116        }
117    }
118
119    public void transform1D( float[] real, float[] imag, int logN, int n, boolean forward ) {
120        scramble( n, real, imag );
121        butterflies( n, logN, forward ? 1 : -1, real, imag );
122    }
123
124    public void transform2D( float[] real, float[] imag, int cols, int rows, boolean forward ) {
125        int log2cols = log2(cols);
126        int log2rows = log2(rows);
127        int n = Math.max(rows, cols);
128        float[] rtemp = new float[n];
129        float[] itemp = new float[n];
130
131        // FFT the rows
132        for ( int y = 0; y < rows; y++ ) {
133            int offset = y*cols;
134            System.arraycopy( real, offset, rtemp, 0, cols );
135            System.arraycopy( imag, offset, itemp, 0, cols );
136            transform1D(rtemp, itemp, log2cols, cols, forward);
137            System.arraycopy( rtemp, 0, real, offset, cols );
138            System.arraycopy( itemp, 0, imag, offset, cols );
139        }
140
141        // FFT the columns
142        for ( int x = 0; x < cols; x++ ) {
143            int index = x;
144            for ( int y = 0; y < rows; y++ ) {
145                rtemp[y] = real[index];
146                itemp[y] = imag[index];
147                index += cols;
148            }
149            transform1D(rtemp, itemp, log2rows, rows, forward);
150            index = x;
151            for ( int y = 0; y < rows; y++ ) {
152                real[index] = rtemp[y];
153                imag[index] = itemp[y];
154                index += cols;
155            }
156        }
157    }
158
159    private int log2( int n ) {
160        int m = 1;
161        int log2n = 0;
162
163        while (m < n) {
164            m *= 2;
165            log2n++;
166        }
167        return m == n ? log2n : -1;
168    }
169
170}