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}