fft.java

dump3 morpheus 0.2.9 src

/**
 * DuMP3 version morpheus_0.2.9 - a duplicate/similar file finder in Java<BR>
 * Copyright 2005 Alexander Gr&auml;sser<BR>
 * All Rights Reserved, http://dump3.sourceforge.net/<BR>
 * <BR>
 * This file is part of DuMP3.<BR>
 * <BR>
 * DuMP3 is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software
 * Foundation; either version 2 of the License, or (at your option) any later version.<BR>
 * <BR>
 * DuMP3 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.<BR>
 * <BR>
 * You should have received a copy of the GNU General Public License along with DuMP3; if not, write to the Free Software Foundation, Inc., 51 Franklin St,
 * Fifth Floor, Boston, MA 02110-1301 USA
 */
package net.za.grasser.duplicate.util.fft;

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

/**
 * Compilation: javac FFT.java<BR>
 * Execution: java FFT N<BR>
 * Dependencies: Complex.java<BR>
 * Compute the FFT and inverse FFT of a length N complex sequence. Bare bones implementation that runs in O(N log N) time.<BR>
 * Limitations: assumes N is a power of 2 - not the most memory efficient algorithm<BR>
 * Windowing functions from <a href="http://astronomy.swin.edu.au/~pbourke/analysis/windows/">Swinn University</a> and <a
 * href="http://ccrma.stanford.edu/CCRMA/Courses/422/projects/kbd/">Stanford University</a>.
 * 
 * @author <a href="http://www.cs.princeton.edu/introcs/97data/FFT.java.html">Robert Sedgewick and Kevin Wayne</a>
 * @version $Revision: 1.7 $
 */
public class FFT {
  /**
   * <code>table</code> FFT - table contains precalculated plan
   */
  private static final List<Complex[]> table = Collections.synchronizedList(new ArrayList<Complex[]>(16));
  static {
    for (int i = 0; i < 17; i++) {
      fillTable(i);
    }
  }

  /**
   * fill the plan table
   * 
   * @param i
   */
  private static final void fillTable(final int i) {
    final int N = (int)Math.pow(2.0, i);
    final Complex[] te = new Complex[N];
    for (int k = 0; k < N; k++) {
      final double kth = -2 * k * Math.PI / (N << 1);
      te[k] = new Complex(Math.cos(kth), Math.sin(kth));
    }
    table.add(i, te);
  }

  /**
   * @param N
   * @return log N / log 2
   */
  private static final int lookup(final int N) {
    switch (N) {
      case 1:
        return 0;
      case 2:
        return 1;
      case 4:
        return 2;
      case 8:
        return 3;
      case 16:
        return 4;
      case 32:
        return 5;
      case 64:
        return 6;
      case 128:
        return 7;
      case 256:
        return 8;
      case 512:
        return 9;
      case 1024:
        return 10;
      case 2048:
        return 11;
      case 4096:
        return 12;
      case 8192:
        return 13;
      case 16384:
        return 14;
      case 32768:
        return 15;
      case 65536:
        return 16;
      default:
        return (int)(Math.log(N) / Math.log(2.0));
    }
  }

  /**
   * private constructor for utility class
   */
  private FFT() {
    super();
  }

  /**
   * compute the FFT of x[], assuming its length is a power of 2
   * 
   * @param x
   * @return the FFT of the Complex array x
   */
  public static Complex[] fft(final Complex[] x) {
    final int N = x.length;
    final Complex[] y = new Complex[N];
    // base case
    if (N == 1) {
      y[0] = x[0];
      return y;
    }
    // radix 2 Cooley-Tukey FFT
    if (N % 2 != 0) {
      throw new RuntimeException("N is not a power of 2");
    }
    final Complex[] even = new Complex[(N >> 1)];
    final Complex[] odd = new Complex[(N >> 1)];
    for (int k = 0; k < N >> 1; k++) {
      even[k] = x[(k << 1)];
    }
    for (int k = 0; k < N >> 1; k++) {
      odd[k] = x[((k << 1) + 1)];
    }
    final Complex[] q = fft(even);
    final Complex[] r = fft(odd);
    for (int k = 0; k < N >> 1; k++) {
      final double kth = -2 * k * Math.PI / N;
      final Complex wk = new Complex(Math.cos(kth), Math.sin(kth));
      y[k] = q[k].plus(wk.times(r[k]));
      y[(k + (N >> 1))] = q[k].minus(wk.times(r[k]));
    }
    return y;
  }

  /**
   * compute the FFT of x[], assuming its length is a power of 2 (faster than fft)
   * 
   * @param x
   * @return the FFT of the Complex array x
   */
  public static Complex[] fft2(final Complex[] x) {
    int N = x.length;
    final Complex[] y = new Complex[N];
    // base case
    if (N == 1) {
      y[0] = x[0];
      return y;
    }
    // radix 2 Cooley-Tukey FFT
    if (N % 2 != 0) {
      throw new RuntimeException("N is not a power of 2");
    }
    N = N >> 1;
    final Complex[] even = new Complex[N];
    final Complex[] odd = new Complex[N];
    for (int k = 0; k < N; k++) {
      even[k] = x[(k << 1)];
    }
    for (int k = 0; k < N; k++) {
      odd[k] = x[((k << 1) + 1)];
    }
    final Complex[] q = fft2(even);
    final Complex[] r = fft2(odd);
    Complex[] te = null;
    try {
      te = table.get(lookup(N));
    } catch (final Exception e) {
      // ignore - te will be null
    }
    if (te == null) {
      final int l = lookup(N);
      fillTable(l);
      te = table.get(l);
    }
    for (int k = 0; k < N; k++) {
      final Complex wk = te[k].times(r[k]);
      y[k] = q[k].plus(wk);
      y[(k + N)] = q[k].minus(wk);
    }
    return y;
  }

  /**
   * compute the inverse FFT of x[], assuming its length is a power of 2
   * 
   * @param x
   * @return the inverse FFT of the Complex array x
   */
  public static Complex[] ifft(final Complex[] x) {
    final int N = x.length;
    final double one_over_N = 1.0 / N;
    // take conjugate
    for (int i = 0; i < N; i++) {
      x[i] = x[i].conjugate();
    }
    final Complex[] y = fft(x);
    // take conjugate again
    for (int i = 0; i < N; i++) {
      y[i] = y[i].conjugate();
    }
    // divide by N
    for (int i = 0; i < N; i++) {
      y[i] = y[i].times(one_over_N);
    }
    return y;
  }

  /**
   * compute the convolution of x and y
   * 
   * @param x
   * @param y
   * @return convolved Complex array
   */
  public static Complex[] convolve(final Complex[] x, final Complex[] y) {
    if (x.length != y.length) {
      throw new RuntimeException("Dimensions don't agree");
    }
    final int N = x.length;
    final Complex[] a = fft(x);
    final Complex[] b = fft(y);
    final Complex[] c = new Complex[N];
    for (int i = 0; i < N; i++) {
      c[i] = a[i].times(b[i]);
    }
    // compute inverse FFT
    return ifft(c);
  }

  /**
   * This routine allocates the memory for the FFT window function lookup table, and it initializes the table with calculated values of the standard Raised
   * Cosine window. This is the general case for Hamming and Hanning windows.
   * 
   * @param len
   * @param alpha must be between 0 and 1 exclusive
   * @return a Raised Cosine window table
   */
  public static double[] setupRaisedCosineWindow(final int len, final double alpha) {
    if (alpha <= 0.0 || alpha >= 1.0) {
      throw new IllegalArgumentException("alpha is invalid. Must be between 0 and 1 exclusive.");
    }
    int i = 0;
    int j = 0;
    final int l = len >> 1;
    final double[] win_tbl = new double[len];
    final double two_x_pi = 2.0 * Math.PI;
    final double one_min_alpha = 1.0 - alpha;
    for (i = -l; i < l; i++) {
      win_tbl[j++] = alpha + one_min_alpha * Math.cos(two_x_pi * i / len);
    }
    return win_tbl;
  }

  /**
   * This routine allocates the memory for the FFT window function lookup table, and it initializes the table with calculated values of the standard Hanning
   * window.
   * 
   * @param len
   * @return a Hanning window table
   */
  public static double[] setupHanningWindow(final int len) {
    return setupRaisedCosineWindow(len, 0.5);
  }

  /**
   * This routine allocates the memory for the FFT window function lookup table, and it initializes the table with calculated values of the standard Hamming
   * window.
   * 
   * @param len
   * @return a Hamming window table
   */
  public static double[] setupHammingWindow(final int len) {
    return setupRaisedCosineWindow(len, 0.54);
  }

  /**
   * This routine allocates the memory for the FFT window function lookup table, and it initializes the table with calculated values of the standard Welch
   * window.
   * 
   * @param len
   * @return a Welch window table
   */
  public static double[] setupWelchWindow(final int len) {
    final double[] win_tbl = new double[len];
    for (int i = 0; i < len; i++) {
      win_tbl[i] = 1.0 - Math.pow((((double)i - (len >> 1)) / (len >> 1)), 2.0);
    }
    return win_tbl;
  }

  /**
   * This routine allocates the memory for the FFT window function lookup table, and it initializes the table with calculated values of the standard Bartlett or
   * triangular window.
   * 
   * @param len
   * @return a Bartlett window table
   */
  public static double[] setupBartlettWindow(final int len) {
    int j = 0;
    final double[] win_tbl = new double[len];
    for (double i = -(len >> 1); i < len >> 1; i++) {
      if (i < 0.0) {
        win_tbl[j++] = 1.0 + 2.0 * i / len;
      } else {
        win_tbl[j++] = 1.0 - 2.0 * i / len;
      }
    }
    return win_tbl;
  }

  /**
   * This routine allocates the memory for the FFT window function lookup table, and it initializes the table with calculated values of the standard Kaiser
   * window.<BR>
   * Translated from <a href="http://ccrma.stanford.edu/CCRMA/Courses/422/projects/kbd/">Kaiser-Bessel Derived Window</a>.
   * 
   * @param len
   * @param alpha
   * @return a Kaiser window table
   */
  public static double[] setupKaiserWindow(final int len, final double alpha) {
    double sumvalue = 0.0;
    int i = 0;
    final double[] window = new double[len];
    for (i = 0; i < len >> 1; i++) {
      sumvalue += BesselI0((Math.PI * alpha * Math.sqrt((1.0 - Math.pow((4.0 * i / len - 1.0), 2.0)))));
      window[i] = sumvalue;
    }
    sumvalue += BesselI0((Math.PI * alpha * Math.sqrt((1.0 - Math.pow((4.0 * (len >> 1) / len - 1.0), 2.0)))));
    // normalize the window and fill in the righthand side of the window:
    for (i = 0; i < len >> 1; i++) {
      window[i] = Math.sqrt(window[i] / sumvalue);
      window[len - 1 - i] = window[i];
    }
    return window;
  }

  /**
   * BesselI0 -- Regular Modified Cylindrical Bessel Function (Bessel I).<BR>
   * Translated from <a href="http://ccrma.stanford.edu/CCRMA/Courses/422/projects/kbd/">Kaiser-Bessel Derived Window</a>
   * 
   * @param x alpha
   * @return Bessel value
   */
  private static double BesselI0(final double x) {
    final double denominator;
    final double numerator;
    final double z;
    if (x == 0.0) {
      return 1.0;
    }
    z = x * x;
    numerator = z
        * (z
            * (z
                * (z
                    * (z
                        * (z
                            * (z
                                * (z
                                    * (z
                                        * (z
                                            * (z
                                                * (z * (z * (z * 0.210580722890567e-22 + 0.380715242345326e-19) + 0.479440257548300e-16) + 0.435125971262668e-13) + 0.300931127112960e-10) + 0.160224679395361e-7) + 0.654858370096785e-5) + 0.202591084143397e-2) + 0.463076284721000e0) + 0.754337328948189e2) + 0.830792541809429e4) + 0.571661130563785e6) + 0.216415572361227e8) + 0.356644482244025e9)
        + 0.144048298227235e10;
    denominator = z * (z * (z - 0.307646912682801e4) + 0.347626332405882e7) - 0.144048298227235e10;
    return -numerator / denominator;
  }

  /**
   * The routine applies the window function to the data set in preparation for the FFT.<BR>
   * This routine operates destructively on the data structure passed.
   * 
   * @param buf
   * @param win_tbl
   * @return the windowed data
   */
  public static Complex[] window(final Complex[] buf, final double[] win_tbl) {
    if (win_tbl != null) {
      int i = 0;
      final int len = buf.length;
      for (i = 0; i < len; i++) {
        // use 'dirty' NaN check - otherwise use !Double.isNaN(win_tbl[i])
        if (win_tbl[i] == win_tbl[i]) {
          buf[i] = buf[i].times(win_tbl[i]);
        } else {
          buf[i] = new Complex(Double.MIN_VALUE, 0.0);
        }
      }
    }
    return buf;
  }
}