Hello,
I found the VisAD project via a code excerpt for your FFT algorithm. I've
been looking at VisAD itself, and it looks very interesting, but in the
short term I just wanted to use your FFT algorithm. I've removed the VisAD
dependencies from the FFT code, just leaving the FT2D and FT1D methods, etc.
I just wanted to check whether it is acceptable (under the LGPL) to package
this as a (very small) library to use from my commercial code. I've attached
the modified code - really all I needed to do was remove the first few
methods and change the VisAD exceptions to IllegalArgumentExceptions.
Best Regards,
Ben.
--
Ben Webster
Director
IT-IS International Ltd
1 Wainstones Court
StokesleyBusiness Park
Stokesley
North Yorkshire
TS9 5JY
United Kingdom
T: +44 (0) 1642 714222
www.itisint.com
Company no: 5098475
Registered in the UK
The information contained in this e-mail and any files transmitted with
it is strictly Confidential and intended for the addressee only. If you have
received this e-mail in error please notify the originator. Reasonable
endeavours have been used to check for virus infection in this email,
including attachments; however, this does not warrant freedom from such
infection and you are responsible for conducting your own virus checking.
import java.text.DecimalFormat;
/*
VisAD system for interactive analysis and visualization of numerical
data. Copyright (C) 1996 - 2007 Bill Hibbard, Curtis Rueden, Tom
Rink, Dave Glowacki, Steve Emmerson, Tom Whittaker, Don Murray, and
Tommy Jasmin.
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA
*/
/**
* FFT is the VisAD class for Fourier Transforms, using the Fast Fourier
* Transform when the domain length is a power of two.
* <p>
*/
public class FFT {
/**
* compute 2-D Fourier transform, calling 1-D FT twice use FFT if rows
and
* cols are powers of 2
*
* @param rows
* first dimension for 2-D
* @param cols
* second dimension for 2-D
* @param x
* array for take Fourier transform of, dimensioned
[2][length]
* where length = rows * cols, and the first index (2) is
over
* real & imaginary parts
* @param forward
* true for forward and false for backward
* @return Fourier transform of x
* @throws VisADException
* a VisAD error occurred
*/
public static float[][] FT2D(int rows, int cols, float[][] x,
boolean forward) throws IllegalArgumentException {
if (x == null)
return null;
if (x.length != 2 || x[0].length != x[1].length) {
throw new IllegalArgumentException("bad x lengths");
}
int n = x[0].length;
if (rows * cols != n) {
throw new IllegalArgumentException(rows + " * " + cols
+ " must equal " + n);
}
float[][] y = new float[2][n];
float[][] z = new float[2][rows];
for (int c = 0; c < cols; c++) {
int i = c * rows;
for (int r = 0; r < rows; r++) {
z[0][r] = x[0][i + r];
z[1][r] = x[1][i + r];
}
z = FT1D(z, forward);
for (int r = 0; r < rows; r++) {
y[0][i + r] = z[0][r];
y[1][i + r] = z[1][r];
}
}
float[][] u = new float[2][n];
float[][] v = new float[2][cols];
for (int r = 0; r < rows; r++) {
int i = r;
for (int c = 0; c < cols; c++) {
v[0][c] = y[0][i + c * rows];
v[1][c] = y[1][i + c * rows];
}
v = FT1D(v, forward);
for (int c = 0; c < cols; c++) {
u[0][i + c * rows] = v[0][c];
u[1][i + c * rows] = v[1][c];
}
}
return u;
}
/**
* compute 2-D Fourier transform, calling 1-D FT twice use FFT if rows
and
* cols are powers of 2
*
* @param rows
* first dimension for 2-D
* @param cols
* second dimension for 2-D
* @param x
* array for take Fourier transform of, dimensioned
[2][length]
* where length = rows * cols, and the first index (2) is
over
* real & imaginary parts
* @param forward
* true for forward and false for backward
* @return Fourier transform of x
* @throws VisADException
* a VisAD error occurred
*/
public static double[][] FT2D(int rows, int cols, double[][] x,
boolean forward) throws IllegalArgumentException {
if (x == null)
return null;
if (x.length != 2 || x[0].length != x[1].length) {
throw new IllegalArgumentException("bad x lengths");
}
int n = x[0].length;
if (rows * cols != n) {
throw new IllegalArgumentException(rows + " * " + cols
+ " must equal " + n);
}
double[][] y = new double[2][n];
double[][] z = new double[2][rows];
for (int c = 0; c < cols; c++) {
int i = c * rows;
for (int r = 0; r < rows; r++) {
z[0][r] = x[0][i + r];
z[1][r] = x[1][i + r];
}
z = FT1D(z, forward);
for (int r = 0; r < rows; r++) {
y[0][i + r] = z[0][r];
y[1][i + r] = z[1][r];
}
}
double[][] u = new double[2][n];
double[][] v = new double[2][cols];
for (int r = 0; r < rows; r++) {
int i = r;
for (int c = 0; c < cols; c++) {
v[0][c] = y[0][i + c * rows];
v[1][c] = y[1][i + c * rows];
}
v = FT1D(v, forward);
for (int c = 0; c < cols; c++) {
u[0][i + c * rows] = v[0][c];
u[1][i + c * rows] = v[1][c];
}
}
return u;
}
/**
* compute 1-D Fourier transform use FFT if length (2nd dimension of x)
is a
* power of 2
*
* @param x
* array for take Fourier transform of, dimensioned
[2][length],
* the first index (2) is over real & imaginary parts
* @param forward
* true for forward and false for backward
* @return Fourier transform of x
* @throws VisADException
* a VisAD error occurred
*/
public static float[][] FT1D(float[][] x, boolean forward) {
if (x == null)
return null;
if (x.length != 2 || x[0].length != x[1].length) {
throw new IllegalArgumentException("bad x lengths");
}
int n = x[0].length;
int n2 = 1;
boolean fft = true;
while (n2 < n) {
n2 *= 2;
if (n2 > n) {
fft = false;
}
}
if (fft)
return FFT1D(x, forward);
float[][] temp = new float[2][n];
float angle = (float) (-2.0 * Math.PI / n);
if (!forward)
angle = -angle;
for (int i = 0; i < n; i++) {
temp[0][i] = (float) Math.cos(i * angle);
temp[1][i] = (float) Math.sin(i * angle);
}
float[][] y = new float[2][n];
for (int i = 0; i < n; i++) {
float re = 0.0f;
float im = 0.0f;
for (int j = 0; j < n; j++) {
int m = (i * j) % n;
re += x[0][j] * temp[0][m] - x[1][j] *
temp[1][m];
im += x[0][j] * temp[1][m] + x[1][j] *
temp[0][m];
}
y[0][i] = re;
y[1][i] = im;
}
if (!forward) {
for (int i = 0; i < n; i++) {
y[0][i] /= n;
y[1][i] /= n;
}
}
return y;
}
/**
* compute 1-D FFT transform length (2nd dimension of x) must be a
power of
* 2
*
* @param x
* array for take Fourier transform of, dimensioned
[2][length],
* the first index (2) is over real & imaginary parts
* @param forward
* true for forward and false for backward
* @return Fourier transform of x
* @throws VisADException
* a VisAD error occurred
*/
public static float[][] FFT1D(float[][] x, boolean forward) {
if (x == null)
return null;
if (x.length != 2 || x[0].length != x[1].length) {
throw new IllegalArgumentException("bad x lengths");
}
int n = x[0].length;
int n2 = 1;
while (n2 < n) {
n2 *= 2;
if (n2 > n) {
throw new IllegalArgumentException(
"x length must be power of 2");
}
}
n2 = n / 2;
float[][] temp = new float[2][n2];
float angle = (float) (-2.0 * Math.PI / n);
if (!forward)
angle = -angle;
for (int i = 0; i < n2; i++) {
temp[0][i] = (float) Math.cos(i * angle);
temp[1][i] = (float) Math.sin(i * angle);
}
float[][] y = FFT1D(x, temp);
if (!forward) {
for (int i = 0; i < n; i++) {
y[0][i] /= n;
y[1][i] /= n;
}
}
return y;
}
/** inner function for 1-D Fast Fourier Transform */
private static float[][] FFT1D(float[][] x, float[][] temp) {
int n = x[0].length;
int n2 = n / 2;
int k = 0;
int butterfly;
int buttered = 0;
if (n == 1) {
float[][] z1 = { { x[0][0] }, { x[1][0] } };
return z1;
}
butterfly = (temp[0].length / n2);
float[][] z = new float[2][n2];
float[][] w = new float[2][n2];
for (k = 0; k < n / 2; k++) {
int k2 = 2 * k;
z[0][k] = x[0][k2];
z[1][k] = x[1][k2];
w[0][k] = x[0][k2 + 1];
w[1][k] = x[1][k2 + 1];
}
z = FFT1D(z, temp);
w = FFT1D(w, temp);
float[][] y = new float[2][n];
for (k = 0; k < n2; k++) {
y[0][k] = z[0][k];
y[1][k] = z[1][k];
float re = w[0][k] * temp[0][buttered] - w[1][k]
* temp[1][buttered];
float im = w[0][k] * temp[1][buttered] + w[1][k]
* temp[0][buttered];
w[0][k] = re;
w[1][k] = im;
y[0][k] += w[0][k];
y[1][k] += w[1][k];
y[0][k + n2] = z[0][k] - w[0][k];
y[1][k + n2] = z[1][k] - w[1][k];
buttered += butterfly;
}
return y;
}
/**
* compute 1-D Fourier transform use FFT if length (2nd dimension of x)
is a
* power of 2
*
* @param x
* array for take Fourier transform of, dimensioned
[2][length],
* the first index (2) is over real & imaginary parts
* @param forward
* true for forward and false for backward
* @return Fourier transform of x
* @throws VisADException
* a VisAD error occurred
*/
public static double[][] FT1D(double[][] x, boolean forward)
throws IllegalArgumentException {
if (x == null)
return null;
if (x.length != 2 || x[0].length != x[1].length) {
throw new IllegalArgumentException("bad x lengths");
}
int n = x[0].length;
int n2 = 1;
boolean fft = true;
while (n2 < n) {
n2 *= 2;
if (n2 > n) {
fft = false;
}
}
if (fft)
return FFT1D(x, forward);
double[][] temp = new double[2][n];
double angle = -2.0 * Math.PI / n;
if (!forward)
angle = -angle;
for (int i = 0; i < n; i++) {
temp[0][i] = Math.cos(i * angle);
temp[1][i] = Math.sin(i * angle);
}
double[][] y = new double[2][n];
for (int i = 0; i < n; i++) {
double re = 0.0;
double im = 0.0;
for (int j = 0; j < n; j++) {
int m = (i * j) % n;
re += x[0][j] * temp[0][m] - x[1][j] *
temp[1][m];
im += x[0][j] * temp[1][m] + x[1][j] *
temp[0][m];
}
y[0][i] = re;
y[1][i] = im;
}
if (!forward) {
for (int i = 0; i < n; i++) {
y[0][i] /= n;
y[1][i] /= n;
}
}
return y;
}
/**
* compute 1-D FFT transform length (2nd dimension of x) must be a
power of
* 2
*
* @param x
* array for take Fourier transform of, dimensioned
[2][length],
* the first index (2) is over real & imaginary parts
* @param forward
* true for forward and false for backward
* @return Fourier transform of x
* @throws VisADException
* a VisAD error occurred
*/
public static double[][] FFT1D(double[][] x, boolean forward)
throws IllegalArgumentException {
if (x == null)
return null;
if (x.length != 2 || x[0].length != x[1].length) {
throw new IllegalArgumentException("bad x lengths");
}
int n = x[0].length;
int n2 = 1;
while (n2 < n) {
n2 *= 2;
if (n2 > n) {
throw new IllegalArgumentException(
"x length must be power of 2");
}
}
n2 = n / 2;
double[][] temp = new double[2][n2];
double angle = (double) (-2.0 * Math.PI / n);
if (!forward)
angle = -angle;
for (int i = 0; i < n2; i++) {
temp[0][i] = (double) Math.cos(i * angle);
temp[1][i] = (double) Math.sin(i * angle);
}
double[][] y = FFT1D(x, temp);
if (!forward) {
for (int i = 0; i < n; i++) {
y[0][i] /= n;
y[1][i] /= n;
}
}
return y;
}
/** inner function for 1-D Fast Fourier Transform */
private static double[][] FFT1D(double[][] x, double[][] temp) {
int n = x[0].length;
int n2 = n / 2;
int k = 0;
int butterfly;
int buttered = 0;
if (n == 1) {
double[][] z1 = { { x[0][0] }, { x[1][0] } };
return z1;
}
butterfly = (temp[0].length / n2);
double[][] z = new double[2][n2];
double[][] w = new double[2][n2];
for (k = 0; k < n2; k++) {
int k2 = 2 * k;
z[0][k] = x[0][k2];
z[1][k] = x[1][k2];
w[0][k] = x[0][k2 + 1];
w[1][k] = x[1][k2 + 1];
}
z = FFT1D(z, temp);
w = FFT1D(w, temp);
double[][] y = new double[2][n];
for (k = 0; k < n2; k++) {
y[0][k] = z[0][k];
y[1][k] = z[1][k];
double re = w[0][k] * temp[0][buttered] - w[1][k]
* temp[1][buttered];
double im = w[0][k] * temp[1][buttered] + w[1][k]
* temp[0][buttered];
w[0][k] = re;
w[1][k] = im;
y[0][k] += w[0][k];
y[1][k] += w[1][k];
y[0][k + n2] = z[0][k] - w[0][k];
y[1][k + n2] = z[1][k] - w[1][k];
buttered += butterfly;
}
return y;
}
/** test Fourier Transform methods */
public static void main(String args[]) throws IllegalArgumentException {
int n = 16;
int rows = 1, cols = 1;
boolean twod = false;
DecimalFormat format = new DecimalFormat("0.0000");
/*
* if (args.length > 0 && args[0].startsWith("AREA")) {
DisplayImpl
* display1 = new DisplayImplJ3D("display"); DisplayImpl
display1 = new
* DisplayImplJ3D("display"); AreaAdapter areaAdapter = new
* AreaAdapter(args[0]); Data image = areaAdapter.getData();
* FunctionType imageFunctionType = (FunctionType)
image.getType();
* RealType radianceType = (RealType) ((RealTupleType)
* imageFunctionType.getRange()).getComponent(0);
display.addMap(new
* ScalarMap(RealType.Latitude, Display.Latitude));
display.addMap(new
* ScalarMap(RealType.Longitude, Display.Longitude)); ScalarMap
rgbMap =
* new ScalarMap(radianceType, Display.RGB);
display.addMap(rgbMap); }
*/
if (args.length > 0)
n = Integer.valueOf(args[0]).intValue();
if (args.length > 1) {
rows = Integer.valueOf(args[0]).intValue();
cols = Integer.valueOf(args[1]).intValue();
n = rows * cols;
twod = true;
}
float[][] x = new float[2][n];
// double[][] x = new double[2][n];
System.out.println(" initial values");
if (twod) {
int i = 0;
for (int c = 0; c < cols; c++) {
for (int r = 0; r < rows; r++) {
x[0][i] = (float) (Math.sin(2 * Math.PI
* r / rows) * Math
.sin(2 * Math.PI * c /
cols));
x[1][i] = 0.0f;
// x[0][i] = (Math.sin(2 * Math.PI * r
/ rows) *
// Math.sin(2 * Math.PI * c / cols));
// x[1][i] = 0.0;
System.out.println("x[" + r + "][" + c
+ "] = "
+
format.format(x[0][i]) + " "
+
format.format(x[1][i]));
i++;
}
}
} else {
for (int i = 0; i < n; i++) {
x[0][i] = (float) Math.sin(2 * Math.PI * i / n);
x[1][i] = 0.0f;
// x[0][i] = Math.sin(2 * Math.PI * i / n);
// x[1][i] = 0.0;
System.out.println("x[" + i + "] = " +
format.format(x[0][i])
+ " " + format.format(x[1][i]));
}
}
x = twod ? FT2D(rows, cols, x, true) : FT1D(x, true);
System.out.println("\n fft");
if (twod) {
int i = 0;
for (int c = 0; c < cols; c++) {
for (int r = 0; r < rows; r++) {
System.out.println("x[" + r + "][" + c
+ "] = "
+
format.format(x[0][i]) + " "
+
format.format(x[1][i]));
i++;
}
}
} else {
for (int i = 0; i < n; i++) {
System.out.println("x[" + i + "] = " +
format.format(x[0][i])
+ " " + format.format(x[1][i]));
}
}
x = twod ? FT2D(rows, cols, x, false) : FT1D(x, false);
System.out.println("\n back fft");
if (twod) {
int i = 0;
for (int c = 0; c < cols; c++) {
for (int r = 0; r < rows; r++) {
System.out.println("x[" + r + "][" + c
+ "] = "
+
format.format(x[0][i]) + " "
+
format.format(x[1][i]));
i++;
}
}
} else {
for (int i = 0; i < n; i++) {
System.out.println("x[" + i + "] = " +
format.format(x[0][i])
+ " " + format.format(x[1][i]));
}
}
}
}
