/*
* Open Source Physics software is free software as described near the bottom of this code file.
*
* For additional information and documentation on Open Source Physics please see:
* <http://www.opensourcephysics.org/>
*/
package org.opensourcephysics.numerics;
/**
* FFT2D computes the FFT of 2 dimensional complex, double precision data.
*
* This class has been copied from Bruce Miller's FFT package for use in the
* Open Source Physics Project. The original package contains code for other transformations
* and other data types.
*
* The data is stored in a 1-dimensional array in Row-Major order.
* The physical layout in the array data, of the mathematical data d[i,j] is as follows:
* <PRE>
* Re(d[i,j]) = data[i*rowspan + 2*j]
* Im(d[i,j]) = data[i*rowspan + 2*j + 1]
* </PRE>
* where <code>rowspan</code> must be at least 2*ncols (it defaults to 2*ncols).
* The transformed data is returned in the original data array in
* <a href="package-summary.html#wraparound">wrap-around</A> order along each dimension.
*
* @author Bruce R. Miller bruce.miller@nist.gov
* @author Contribution of the National Institute of Standards and Technology,
* @author not subject to copyright.
*/
public class FFT2D {
static final double PI2 = 2*Math.PI;
int nrows;
int ncols;
FFT rowFFT, colFFT;
double[] acol, ccol;
/**
* Create an FFT for transforming nrows*ncols points of Complex, double precision
* data.
* @param nrows
* @param ncols
*/
public FFT2D(int nrows, int ncols) {
if((nrows<=0)||(ncols<=0)) {
throw new IllegalArgumentException("The array dimensions >=0 : "+nrows+","+ncols); //$NON-NLS-1$ //$NON-NLS-2$
}
this.nrows = nrows;
this.ncols = ncols;
acol = new double[2*nrows];
if(nrows%2==1) {
ccol = new double[2*nrows]; // temp storage for center column if nrows is odd
}
rowFFT = new FFT(ncols);
colFFT = ((nrows==ncols) ? rowFFT : new FFT(nrows));
}
protected void checkData(double data[], int rowspan) {
if(rowspan<2*ncols) {
throw new IllegalArgumentException("The row span "+rowspan+"is shorter than the row length "+2*ncols); //$NON-NLS-1$ //$NON-NLS-2$
}
if(nrows*rowspan>data.length) {
throw new IllegalArgumentException("The data array is too small for "+nrows+"x"+rowspan+" data.length="+data.length); //$NON-NLS-1$ //$NON-NLS-2$ //$NON-NLS-3$
}
}
/**
* Compute the Fast Fourier Transform of data leaving the result in data.
* The array data must be dimensioned (at least) 2*nrows*ncols, consisting of
* alternating real and imaginary parts.
* @param data
*/
public void transform(double data[]) {
transform_internal(data, 2*ncols);
}
/**
* Compute the Fast Fourier Transform of data leaving the result in data.
* The array data must be dimensioned (at least) 2*nrows*ncols, consisting of
* alternating real and imaginary parts.
* @param data
* @param rowspan
*/
void transform_internal(double data[], int rowspan) {
checkData(data, rowspan);
for(int i = 0; i<nrows; i++) {
rowFFT.transform_internal(data, i*rowspan, 2, FFT.FORWARD);
}
for(int j = 0; j<ncols; j++) {
colFFT.transform_internal(data, 2*j, rowspan, FFT.FORWARD);
}
}
/**
* Compute the (unnomalized) inverse FFT of data, leaving it in place.
* @param data
*/
public void backtransform(double data[]) {
backtransform_internal(data, 2*ncols);
}
/**
* Compute the (unnomalized) inverse FFT of data, leaving it in place.
* @param data
* @param rowspan
*/
void backtransform_internal(double data[], int rowspan) {
checkData(data, rowspan);
for(int j = 0; j<ncols; j++) {
colFFT.transform_internal(data, 2*j, rowspan, FFT.BACKWARD);
}
for(int i = 0; i<nrows; i++) {
rowFFT.transform_internal(data, i*rowspan, 2, FFT.BACKWARD);
}
}
/**
* Compute the (nomalized) inverse FFT of data, leaving it in place.
* @param data
*/
public void inverse(double data[]) {
inverse_internal(data, 2*ncols);
}
/**
* Compute the (nomalized) inverse FFT of data, leaving it in place.
* @param data
* @param rowspan
*/
void inverse_internal(double data[], int rowspan) {
backtransform_internal(data, rowspan);
double norm = 1.0/((double) nrows*ncols);
for(int i = 0; i<nrows; i++) {
data[i] *= norm;
}
}
/**
* Gets an array containing the frequencies in natural order.
* Data are separated by delta and there are n data points.
*
* @param delta
* @param n size of frequency array
* @return array of frequencies
*/
double[] getNaturalFreq(double delta, int n) {
double[] freq = new double[n];
double f = -0.5/delta, df = -2*f/(n-n%2);
for(int i = 0; i<n; i++) {
freq[i] = f;
f += df;
}
return freq;
}
/**
* Gets the minimum frequency given the domain and the number of points.
*
* @param min double
* @param max double
* @param n int
* @return double
*/
public double getFreqMin(double min, double max, int n) {
return -(n/2)/(max-min);
}
/**
* Gets the maximum frequency given the domain and the number of points.
*
* @param min double
* @param max double
* @param n int
* @return double
*/
public double getFreqMax(double min, double max, int n) {
return((n+1)/2-1)/(max-min);
}
/**
* Gets an array containing the mode numbers in natural order.
*
* @return array of mode numbers
*/
public double[] getNaturalModes(int n) {
double[] bins = new double[n];
double w = -(n-n%2)/2;
for(int i = 0; i<n; i++) {
bins[i] = w;
w++;
}
return bins;
}
/**
* Gets an array containing the mode numbers in wrap-around order.
*
* @return array of mode numbers
*/
public double[] getWrappedModes(int n) {
double[] bins = new double[n];
for(int i = 0; i<n; i++) {
bins[i] = (i<(n+1)/2) ? i : (i-n);
}
return bins;
}
/**
* Gets an array containing the angular frequencies (wavenumbers) in natural order.
* The first data point is at xmin (tmin) and the last data point is at xmax (tmax).
*
* @param xmin
* @param xmax
* @return the array of frequencies
*/
public double[] getWrappedOmegaX(double xmin, double xmax) {
return getWrappedFreq((xmax-xmin)/(nrows-nrows%2)/PI2, nrows);
}
/**
* Gets an array containing the angular frequencies (wavenumbers) in natural order.
* The first data point is at ymin and the last data point is at ymax.
*
* @param ymin
* @param ymax
* @return the array of frequencies
*/
public double[] getWrappedOmegaY(double ymin, double ymax) {
return getWrappedFreq((ymax-ymin)/(ncols-ncols%2)/PI2, ncols);
}
/**
* Gets an array containing the frequencies in wrap-around order.
* Samples in the data are separated by delta.
*
* @param delta
* @return the array of frequencies
*/
public double[] getWrappedFreq(double delta, int n) {
double[] freq = new double[n];
double f = -0.5/delta, df = -2*f/(n-n%2);
for(int i = 0; i<n; i++) {
freq[i] = (i<(n+1)/2) ? i*df : (i-n)*df;
}
return freq;
}
/**
* Gets an array containing the frequencies in natural order.
* Data are separated by delta.
*
* @param delta
* @return the array of frequencies
*/
public double[] getNaturalFreqX(double delta) {
return getNaturalFreq(delta, nrows);
}
/**
* Gets an array containing the frequencies in natural order.
* The first data point is at xmin (tmin) and the last data point is at xmax (tmax).
*
* @param xmin
* @param xmax
* @return the array of frequencies
*/
public double[] getNaturalFreqX(double xmin, double xmax) {
return getNaturalFreq((xmax-xmin)/(nrows-nrows%2), nrows);
}
/**
* Gets an array containing the angular frequencies (wavenumbers) in natural order.
* Data are separated by delta.
*
* @param delta
* @return the array of frequencies
*/
public double[] getNaturalOmegaX(double delta) {
return getNaturalFreq(delta/PI2, nrows);
}
/**
* Gets an array containing the angular frequencies (wavenumbers) in natural order.
* The first data point is at xmin (tmin) and the last data point is at xmax (tmax).
*
* @param xmin
* @param xmax
* @return the array of frequencies
*/
public double[] getNaturalOmegaX(double xmin, double xmax) {
return getNaturalFreq((xmax-xmin)/(nrows-nrows%2)/PI2, nrows);
}
/**
* Gets an array containing the frequencies in natural order if samples in the orginal data are
* separated by delta in y.
*
* @param delta
* @return the array of frequencies
*/
public double[] getNaturalFreqY(double delta) {
return getNaturalFreq(delta, ncols);
}
/**
* Gets an array containing the frequencies in natural order.
* The first data point is at ymin and the last data point is at ymax.
*
* @param ymin
* @param ymax
* @return the array of frequencies
*/
public double[] getNaturalFreqY(double ymin, double ymax) {
return getNaturalFreq((ymax-ymin)/(ncols-ncols%2), ncols);
}
/**
* Gets an array containing the frequencies in natural order if samples in the orginal data are
* separated by delta in y.
*
* @param delta
* @return the array of frequencies
*/
public double[] getNaturalOmegaY(double delta) {
return getNaturalFreq(delta/PI2, ncols);
}
/**
* Gets an array containing the frequencies in natural order.
* The first data point is at ymin and the last data point is at ymax.
*
* @param ymin
* @param ymax
* @return the array of frequencies
*/
public double[] getNaturalOmegaY(double ymin, double ymax) {
return getNaturalFreq((ymax-ymin)/(ncols-ncols%2)/PI2, ncols);
}
/**
* Reorder and normalize the transformed data from most negative frequency
* to most positive frequency leaving the result in data.
* @param data
*/
public void toNaturalOrder(double data[]) {
if(ccol!=null) {
System.arraycopy(data, (ncols/2)*acol.length, ccol, 0, ccol.length); // save center column if ncols is odd
}
for(int i = 0; i<ncols/2; i++) { // swap columns
int offset = i*acol.length;
System.arraycopy(data, offset, acol, 0, acol.length); // save the i-th column data
System.arraycopy(data, (ncols/2+i+ncols%2)*acol.length, data, offset, acol.length); // replace the i-th column data
System.arraycopy(acol, 0, data, (ncols/2+i)*acol.length, acol.length);
}
if(ccol!=null) {
System.arraycopy(ccol, 0, data, data.length-ccol.length, ccol.length);
}
for(int i = 0; i<ncols; i++) { // swap rows
int n = acol.length/2;
int offset = i*acol.length;
System.arraycopy(data, offset, acol, 0, acol.length); // save a column of data
System.arraycopy(acol, n+n%2, data, offset, n-n%2); // copy the second half of the data into the first half
System.arraycopy(acol, 0, data, offset+n-n%2, n+n%2); // copy the first half of the data into the second half
}
double norm = 1.0/((double) nrows*ncols);
for(int i = 0, n = data.length; i<n; i++) { // normalize
data[i] *= norm;
}
}
}
/*
* Open Source Physics software is free software; you can redistribute
* it and/or modify it under the terms of the GNU General Public License (GPL) as
* published by the Free Software Foundation; either version 2 of the License,
* or(at your option) any later version.
* Code that uses any portion of the code in the org.opensourcephysics package
* or any subpackage (subdirectory) of this package must must also be be released
* under the GNU GPL license.
*
* This software 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
* or view the license online at http://www.gnu.org/copyleft/gpl.html
*
* Copyright (c) 2007 The Open Source Physics project
* http://www.opensourcephysics.org
*/