/*
* 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;
/**
* FFTReal computes the discrete Fourier coefficients
* a[0], ...., a[N/2] and b[1], ...., b[N/2 - 1]
* of the discrete partial Fourier sum
* a[0] + a[1]*cos(N/2*omega*x)
* + Sum (k=1,2,...,N/2-1) (a[2*k] * cos(k * omega * x) + a[2*k+1] * sin(k * omega * x))
* given real functional values y[0], ...., y[N-1].
*
* Adapted by W. Christian for use in the OSP project.
*
* @author Bruce R. Miller bruce.miller@nist.gov
* @author Contribution of the National Institute of Standards and Technology,
* @author Derived from GSL (Gnu Scientific Library)
* @author GSL's FFT Code by Brian Gough bjg@vvv.lanl.gov
* @author Since GSL is released under
* @author <H HREF="http://www.gnu.org/copyleft/gpl.html">GPL</A>,
* @author this class must also be.
*
* @version 1.0
*/
public class FFTReal {
static final double PI2 = 2*Math.PI;
int n; // number of data points
FFT fft = new FFT(); // complex fft that does the computation
/**
* Constructs a real FFT transformation for n data points.
*/
public FFTReal() {
setN(2); // avoids null pointer exceptions
}
/**
* Constructs a real FFT transformation for n data points.
*
* @param n the number of data points
*/
public FFTReal(int n) {
setN(n);
}
/**
* Sets the number of data points.
*
* @param n int
*/
public void setN(int n) {
if(n%2!=0) {
throw new IllegalArgumentException(n+" is not even"); //$NON-NLS-1$
}
this.n = n;
fft.setN(n/2);
}
/**
* Gets the number of data points.
*
* @return int
*/
public int getN() {
return n;
}
/**
* Computes the Fast Fourier Transform of the data leaving the result in data.
*
* The given array is returned after it has been transformed.
*
* @param data double[] the data to be transformed
* @return double[] the data after the FFT
*/
public double[] transform(double data[]) {
if(data.length!=n) {
setN(data.length);
}
fft.transform(data);
shuffle(data, +1);
return data;
}
/**
* Computes the (unnomalized) inverse FFT of data, leaving it in place.
*
* The given array is returned after it has been transformed.
*
* @param data double[] the data to be transformed
* @return double[] the data after the FFT
*/
public double[] backtransform(double data[]) {
if(data.length!=n) {
setN(data.length);
}
shuffle(data, -1);
fft.backtransform(data);
return data;
}
/**
* Computes the (nomalized) inverse FFT of data, leaving it in place.
*
* The given array is returned after it has been transformed.
*
* @param data double[] the data to be transformed
* @return double[] the data after the FFT
*/
public double[] inverse(double data[]) {
backtransform(data);
/* normalize inverse fft with 2/n */
double norm = 2.0/(n);
for(int i = 0; i<n; i++) {
data[i] *= norm;
}
return data;
}
/**
* Gets an array containing the frequencies in natural order.
* Data are separated by delta.
*
* @param delta
* @return the array of frequencies
*/
public double[] getNaturalFreq(double delta) {
int n = this.n/2;
double[] freq = new double[n];
double f = 0, df = 0.5/n/delta;
for(int i = 0; i<n; i++) {
freq[i] = f;
f += df;
}
return freq;
}
/**
* 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[] getNaturalFreq(double xmin, double xmax) {
return getNaturalFreq((xmax-xmin)/(n-n%2));
}
/**
* Gets an array containing the frequencies in natural order.
* Data are separated by delta.
*
* @param delta
* @return the array of frequencies
*/
public double[] getNaturalOmega(double delta) {
return getNaturalFreq(delta/PI2);
}
/**
* 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[] getNaturalOmega(double xmin, double xmax) {
return getNaturalFreq((xmax-xmin)/(n-n%2)/PI2);
}
/**
* Rearrage the coefficients.
* @param data double[]
* @param sign int
*/
private void shuffle(double data[], int sign) {
int nh = n/2;
int nq = n/4;
if(n==6) {
nq = 2;
}
double c1 = 0.5, c2 = -0.5*sign;
double theta = sign*Math.PI/nh;
double wtemp = Math.sin(0.5*theta);
double wpr = -2.0*wtemp*wtemp;
double wpi = -Math.sin(theta);
double wr = 1.0+wpr;
double wi = wpi;
for(int i = 1; i<nq; i++) {
int i1 = 2*i;
int i3 = n-i1;
double h1r = c1*(data[i1]+data[i3]);
double h1i = c1*(data[i1+1]-data[i3+1]);
double h2r = -c2*(data[i1+1]+data[i3+1]);
double h2i = c2*(data[i1]-data[i3]);
data[i1] = h1r+wr*h2r-wi*h2i;
data[i1+1] = h1i+wr*h2i+wi*h2r;
data[i3] = h1r-wr*h2r+wi*h2i;
data[i3+1] = -h1i+wr*h2i+wi*h2r;
wtemp = wr;
wr += wtemp*wpr-wi*wpi;
wi += wtemp*wpi+wi*wpr;
}
double d0 = data[0];
if(sign==1) {
data[0] = d0+data[1];
data[1] = d0-data[1];
} else {
data[0] = c1*(d0+data[1]);
data[1] = c1*(d0-data[1]);
}
if(n%4==0) {
data[nh+1] *= -1;
}
}
}
/*
* 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
*/