NOTICE: This version of the NSF Unidata web site (archive.unidata.ucar.edu) is no longer being updated.
Current content can be found at unidata.ucar.edu.
To learn about what's going on, see About the Archive Site.
NOTICE: This version of the NSF Unidata web site (archive.unidata.ucar.edu) is no longer being updated.
Current content can be found at unidata.ucar.edu.
To learn about what's going on, see About the Archive Site.
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]));
}
}
}
}
visad archives: