package transforms;
public class FFT1d {
private final int FORWARD_FFT = -1;
private final int REVERSE_FFT = 1;
private int direction = FORWARD_FFT;
private static final float twoPI = (float) (2 * Math.PI);
private double realPart[] = null;
private double imaginaryPart[] = null;
public static final double gauss(
double x,
double xc, double sigma) {
double oneOnSigmaSquaredOn2 = 1 / (sigma * sigma) / 2;
return
Math.exp(-((x - xc) * (x - xc)) *
oneOnSigmaSquaredOn2) / Math.PI * oneOnSigmaSquaredOn2;
}
public static final double[]getGauss(int n){
double g[] = new double[n];
for (int x=0; x < n; x++)
g[x]=gauss(x,n/2,2);
return g;
}
public static final void testGauss() {
print(getGauss(10));
}
public static final void print(double d[]){
for (int i=0; i < d.length; i++)
System.out.println(d[i]);
}
public void visuallyTest() {
System.out.println("Testing FFT engine in 1D class...");
int N = 8;
double[] in_r = new double[N];
double[] in_i = new double[N];
double[] holder_r = new double[N];
double[] holder_i = new double[N];
for (int i = 0; i < N; i++) {
in_r[i] = (double) i / N;
in_i[i] = (double) 0;
}
forwardFFT(in_r, in_i);
for (int i = 0; i < N; i++) {
holder_r[i] = in_r[i];
holder_i[i] = in_i[i];
}
for (int i = 0; i < N; i++) {
in_r[i] = in_r[i];
in_i[i] = in_i[i];
}
reverseFFT(in_r, in_i);
double[] mag = magnitudeSpectrum(in_r, in_i);
System.out.println("i x1[i] re[i] im[i] tv[i]");
for (int i = 0; i < N; i++) {
System.out.println(i +
"\t" + i +
"\t" + holder_r[i] +
"\t\t" +
holder_i[i] +
"\t\t" + mag[i]);
}
}
public void timeFFT(int n) {
double realPart[] = new double[n];
double imaginaryPart[] = new double[n];
synthesizeRamp(realPart, imaginaryPart);
ip.gui.Timer t1 = new ip.gui.Timer();
runExperiment(t1, realPart, imaginaryPart);
}
public static final void synthesizeRamp(double[] realPart, double[] imaginaryPart) {
int n = realPart.length;
for (int i = 0; i < n; i++) {
realPart[i] = (double) i / n;
imaginaryPart[i] = (double) 0;
}
}
private void runExperiment(ip.gui.Timer t1,
double[] realPart,
double[] imaginaryPart) {
t1.start();
double w[] = getGauss(realPart.length);
for (int i = 0; i < 1000; i++) {
window(realPart,w);
forwardFFT(realPart, imaginaryPart);
getMaxPSDLocation(realPart, imaginaryPart);
}
t1.print(realPart.length +
" done 1000 times point 1D FFT (using double):");
}
private static final void window(double r[], double w[]){
for (int i=0; i < r.length; i++)
r[i]=r[i]*w[i];
}
public void swap(int i, int numBits) {
double tempr;
int j = bitr(i, numBits);
tempr = realPart[j];
realPart[j] = realPart[i];
realPart[i] = tempr;
tempr = imaginaryPart[j];
imaginaryPart[j] = imaginaryPart[i];
imaginaryPart[i] = tempr;
}
private void swapInt(int i, int j) {
double tempr;
int ti;
int tj;
ti = i - 1;
tj = j - 1;
tempr = realPart[tj];
realPart[tj] = realPart[ti];
realPart[ti] = tempr;
tempr = imaginaryPart[tj];
imaginaryPart[tj] = imaginaryPart[ti];
imaginaryPart[ti] = tempr;
}
public static final int getMaxPSDLocation(
double realPart[],
double imaginaryPart[]) {
double max;
max = -0.99e30;
int location = 0;
for (int i = 0; i < realPart.length; i++)
if (magnitude(realPart[i], imaginaryPart[i]) > max) {
max = realPart[i];
location = i;
}
return location;
}
public static final double magnitude(double r, double imag) {
return r * r + imag * imag;
}
double getMaxValue(double in[]) {
double max;
max = -0.99e30;
for (int i = 0; i < in.length; i++)
if (in[i] > max)
max = in[i];
return max;
}
void bitReverse2() {
int n = realPart.length;
int j = 1;
int k;
for (int i = 1; i < n; i++) {
if (i < j) {
swapInt(i, j);
}
k = n / 2;
while (k >= 1 && k < j) {
j = j - k;
k = k / 2;
}
j = j + k;
} }
int bitr(int j, int numBits) {
int ans = 0;
for (int i = 0; i < numBits; i++) {
ans = (ans << 1) + (j & 1);
j = j >> 1;
}
return ans;
}
public void forwardFFT(double[] in_r, double[] in_i) {
direction = FORWARD_FFT;
fft(in_r, in_i);
}
public void reverseFFT(double[] in_r, double[] in_i) {
direction = REVERSE_FFT;
fft(in_r, in_i);
}
public int log2(double d) {
return (int) (Math.log(d) / Math.log(2));
}
public void fft(double _realPart[], double _imaginaryPart[]) {
int id;
int localN;
double wtemp, Wjk_r, Wjk_i, Wj_r, Wj_i;
double theta, tempr, tempi;
int numBits = log2(_realPart.length);
int length = 1 << numBits; int n = length;
realPart = _realPart;
imaginaryPart = _imaginaryPart;
bitReverse2();
for (int m = 1; m <= numBits; m++) {
localN = 1 << m;
Wjk_r = 1;
Wjk_i = 0;
theta = twoPI / localN;
Wj_r = Math.cos(theta);
Wj_i = direction * Math.sin(theta);
int nby2 = localN / 2;
for (int j = 0; j < nby2; j++) {
for (int k = j; k < n; k += localN) {
id = k + nby2;
tempr = Wjk_r * realPart[id] - Wjk_i * imaginaryPart[id];
tempi = Wjk_r * imaginaryPart[id] + Wjk_i * realPart[id];
realPart[id] = realPart[k] - tempr;
imaginaryPart[id] = imaginaryPart[k] - tempi;
realPart[k] += tempr;
imaginaryPart[k] += tempi;
}
wtemp = Wjk_r;
Wjk_r = Wj_r * Wjk_r - Wj_i * Wjk_i;
Wjk_i = Wj_r * Wjk_i + Wj_i * wtemp;
}
}
}
public double[] getRealData() {
return realPart;
}
public double[] getImaginaryData() {
return imaginaryPart;
}
public double[] magnitudeSpectrum(double[] in_r, double[] in_i) {
int N = in_r.length;
double[] mag = new double[N];
for (int i = 0; i < in_r.length; i++)
mag[i] = Math.sqrt(in_r[i] * in_r[i] + in_i[i] * in_i[i]);
return (mag);
}
public static void testFFT() {
FFT1d f = new FFT1d();
f.timeFFT(1024);
}
public static void main(String args[]){
testFFT();
}
}