KJFFT.java

/**
 * --------------------------------------------------------------------------------
 * NoiseTube Mobile client (Java implementation; Android version)
 * <p>
 * Copyright (C) 2008-2010 SONY Computer Science Laboratory Paris
 * Portions contributed by Vrije Universiteit Brussel (BrusSense team), 2008-2011
 * Android port by Vrije Universiteit Brussel (BrusSense team), 2010-2011
 * --------------------------------------------------------------------------------
 * This library is free software; you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License, version 2.1, as published
 * by the Free Software Foundation.
 * <p>
 * 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 Lesser General Public License for more
 * details.
 * <p>
 * You should have received a copy of the GNU Lesser General Public License along
 * with this library; if not, write to:
 * Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor,
 * Boston, MA  02110-1301, USA.
 * <p>
 * Full GNU LGPL v2.1 text: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt
 * NoiseTube project source code repository: http://code.google.com/p/noisetube
 * --------------------------------------------------------------------------------
 * More information:
 * - NoiseTube project website: http://www.noisetube.net
 * - Sony Computer Science Laboratory Paris: http://csl.sony.fr
 * - VUB BrusSense team: http://www.brussense.be
 * --------------------------------------------------------------------------------
 */

package dosimetry.app.wikilek.filters;

import dosimetry.app.wikilek.tools.Utils;

/**
 * This class is based on the KJFFT.java class, which is part of the KJ DSS project.
 * <p>
 * Further information on the KJ DSS project :
 * Author : Kristofer Fudalewski
 * Email  : sirk_sytes@hotmail.com
 * Website: http://sirk.sytes.net
 * <p>
 * It is a Fast Fourier Transformation class.
 *
 * @author kristofer, sbarthol
 */
public class KJFFT {

  private float[] xre;
  private float[] xim;
  private float[] mag;

  private float[] fftSin;
  private float[] fftCos;
  private int[] fftBr;

  private float[] fc3;
  private float[] fco;
  private float[] fl3;
  private float[] fu3;
  private float[] flo;
  private float[] fuo;

  private int ss, ss2, nu;
  private float[] fTable;

  /**
   * @param pSampleSize The amount of the sample provided to the "calculate" method to use during
   *                    FFT calculations, this is used to prepare the calculation tables in advance.
   *                    This value is automatically rounded up to the nearest power of 2.
   */
  public KJFFT(int pSampleSize) {

    nu = (int) Math.ceil(Math.log(pSampleSize) / Math.log(2));

    // -- Calculate the nearest sample size to a power of 2.
    ss = (int) Math.pow(2, nu);

    ss2 = ss >> 1;

    // -- Allocate calculation buffers.
    xre = new float[ss];
    xim = new float[ss];
    mag = new float[ss2];

    // -- Allocate FFT SIN/COS tables.
    fftSin = new float[nu * ss2];
    fftCos = new float[nu * ss2];

    prepareTables();
    fTable = calculateFrequencyTable(Utils.LEQ_1SEC);
    computeFrequencyes();

  }

  /**
   * Compute all the base frequencies for
   * the Octaves bands and the 1/3 octaves bands
   */
  private void computeFrequencyes() {
    fc3 = new float[Utils.NB_3_BANDS];
    fco = new float[Utils.NB_OCTAVES];

    for (int i = -13; i < 11; i++) {
      fc3[i + 13] = (float) (1000.0 * Math.pow(2.0, i / 3.0));
    }
    for (int i = 0, j = 0; i < fc3.length; i++) {
      if (i % 3 == 0) {
        fco[j] = fc3[i + 1];
        j++;
      }
    }

    // upper and lower frequency for each central one
    fl3 = Utils.arr_op_const(fc3, (float) Math.pow(2.0, -1 / 6.0), Utils.Operator.MUL);
    fu3 = Utils.arr_op_const(fc3, (float) Math.pow(2.0, 1 / 6.0), Utils.Operator.MUL);
    flo = Utils.arr_op_const(fco, (float) Math.pow(2.0, -1 / 3.0), Utils.Operator.MUL);
    fuo = Utils.arr_op_const(fco, (float) Math.pow(2.0, 1 / 3.0), Utils.Operator.MUL);
  }

  /**
   * Bit swapping method.
   */
  private int bitrev(int pJ, int pNu) {

    int j1 = pJ;
    int j2;
    int k = 0;

    for (int i = 0; i < pNu; i++) {
      j2 = j1 >> 1;
      k = (k << 1) + j1 - (j2 << 1);
      j1 = j2;
    }

    return k;

  }


  /**
   * Converts sound data over time into pressure values. (FFT)
   *
   * @param samples The sample to compute FFT values on.
   * @return The results of the calculation, normalized between 0.0 and 1.0.
   */
  public float[] calculate(float[] samples) {

    int n2 = ss2;

    // -- Fill buffer.
    for (int a = 0; a < samples.length; a++) {
      xre[a] = samples[a];
      xim[a] = (float) 0.0;
    }

    // -- Clear the remainder of the buffer.
    for (int a = samples.length; a < ss; a++) {
      xre[a] = 0.0f;
      xim[a] = 0.0f;
    }

    float tr, ti, c, s;
    int k, kn2, x = 0;

    for (int l = 0; l < nu; l++) {
      k = 0;

      while (k < ss) {
        for (int i = 0; i < n2; i++) {
          // -- Tabled sin/cos
          c = fftCos[x];
          s = fftSin[x];
          kn2 = k + n2;

          tr = xre[kn2] * c + xim[kn2] * s;
          ti = xim[kn2] * c - xre[kn2] * s;

          xre[kn2] = xre[k] - tr;
          xim[kn2] = xim[k] - ti;
          xre[k] += tr;
          xim[k] += ti;

          k++;
          x++;
        }
        k += n2;
      }
      n2 >>= 1;
    }

    int r;

    // -- Reorder output.
    for (k = 0; k < ss; k++) {
      // -- Use tabled BR values.
      r = fftBr[k];

      if (r > k) {
        tr = xre[k];
        xre[k] = xre[r];
        xre[r] = tr;

        ti = xim[k];
        xim[k] = xim[r];
        xim[r] = ti;
      }
    }


    float[] dataFFT_RE = Utils.arr_op_const(xre, (float) (Math.sqrt(2) / (samples.length)), Utils.Operator.MUL);
    float[] dataFFT_IM = Utils.arr_op_const(xim, (float) (Math.sqrt(2) / (samples.length)), Utils.Operator.MUL);

    // -- Calculate magnitude.
    for (int i = 0; i < ss2; i++) {
      mag[i] = Math.abs((float) (dataFFT_RE[i] * dataFFT_RE[i]) + (dataFFT_IM[i] * dataFFT_IM[i]));
    }

    return mag;

  }

  /**
   * Calculates a table of frequencies represented by the amplitude data returned by the 'calculate' method.
   * Each element states the end of the frequency range of the corresponding FFT band (or bin). For example:
   * <p>
   * Range of band 0 =                 0.0 hz to frequencyTable[ 0 ] hz
   * Range of band 1 = frequencyTable[ 0 ] hz to frequencyTable[ 1 ] hz
   * Range of band 2 = frequencyTable[ 1 ] hz to frequencyTable[ 2 ] hz
   * ... and so on.
   * <p>
   * Calculation uses the sample size rounded to the nearest power of 2 of the FFT instance and the sample rate parameter
   * to build this table.
   *
   * @param pSampleRate The sample rate used to calculate the frequency table. Usually the sample rate of the input
   *                    to the FFT calculate method.
   * @return An array of frequency limits for each band.
   */
  public float[] calculateFrequencyTable(float pSampleRate) {

    float wFr = pSampleRate / 2.0f;

    // -- Calculate band width.
    float wBw = wFr / ss2;

    // -- Store for frequency table.
    float[] wFt = new float[ss2];

    // -- Build band range table.
    int b = 0;

    for (float wFp = (wBw / 2.0f); wFp < wFr - 1; wFp += wBw) {
      wFt[b] = wFp;
      b++;
    }

    return wFt;

  }

  private float[] ones(int sampleRate) {
    float[] wFt = new float[sampleRate];
    for (int i = 0; i < wFt.length; i++) {
      wFt[i] = 1;
    }
    return wFt;
  }

  /**
   * Returns the sample size this FFT instance uses for processing. It is automatically rounded to the nearest power of 2.
   *
   * @return The sample size used by the calculate method.
   */
  public int getInputSampleSize() {
    return ss;
  }

  /**
   * Returns the sample size this FFT instance returns after processing. It is automatically rounded to the nearest power of 2.
   *
   * @return The sample size returned by the calculate method.
   */
  public int getOutputSampleSize() {
    return ss2;
  }

  /**
   * Pre-calculates SIN/COS and bitrev tables in memory.
   */
  private void prepareTables() {

    int n2 = ss2;
    int nu1 = nu - 1;

    //       System.out.println( "bs: " + ( nu * n2 ) );

    float p, arg;
    int k = 0, x = 0;

    // -- Prepare SIN/COS tables.
    for (int l = 0; l < nu; l++) {
      k = 0;

      while (k < ss) {

        for (int i = 0; i < n2; i++) {
          p = bitrev(k >> nu1, nu);

          arg = 2 * (float) Math.PI * p / ss;
          fftSin[x] = (float) Math.sin(arg);
          fftCos[x] = (float) Math.cos(arg);

          k++;
          x++;
        }

        k += n2;
      }
      nu1--;
      n2 >>= 1;
    }

    // -- Prepare bitrev table.
    fftBr = new int[ss];

    for (k = 0; k < ss; k++) {
      fftBr[k] = bitrev(k, nu);
    }

  }

  /**
   * Computes the result in Octaves
   *
   * @param arr_samples
   * @return
   */
  public float[] frequencyAnalyzeOctaves(float[] arr_samples, boolean standard) {

    float[] msqt = computeMSQT(arr_samples, standard);

    // In Octaves bands now
//    	float[] msqo = new float[Utils.NB_OCTAVES];
//    	for(int i = 0; i < fco.length; i++) {
//    		msqo[i] = Utils.arr_sum(Utils.arr_trunk_index(msqt, i, i + 2));
//    	}

    // convert the pressure values to dB values
    float[] Lpo = new float[msqt.length];
    for (int i = 0; i < Lpo.length; i++) {
      Lpo[i] = (float) (10 * Math.log10(msqt[i]) + 94);
    }

    // will free some memory when the GC
    // decides to show up
    msqt = null;

    return Lpo;
  }

  /**
   * Computes the result in 1/3 Bands
   *
   * @param arr_samples
   * @return
   */
  public float[] frequencyAnalyze24Bands(float[] arr_samples, boolean standard) {

    float[] msqt = computeMSQT(arr_samples, standard);

    // convert the pressure values to dB values
    float[] Lpt = new float[msqt.length];
    for (int i = 0; i < Lpt.length; i++) {
      Lpt[i] = (float) (10 * Math.log10(msqt[i]) + 94);
    }

    // the mighty task of releasing some memory space
    msqt = null;

    return Lpt;
  }

  /**
   * Computes the data and put it together by frequencies
   *
   * @param arr_samples
   * @return
   */
  private float[] computeMSQT(float[] arr_samples, boolean standard) {

    float[][] filter = new float[Utils.NB_3_BANDS][];
    float[] p_band = new float[Utils.NB_3_BANDS];

    float[] fft = new float[Utils.BLOCK_SIZE];
    for (int i = 0; i < arr_samples.length / Utils.BLOCK_SIZE; i++) {
      float[] tmp_result = calculate(Utils.arr_trunk_index(Utils.outputToVoltage(arr_samples), i * Utils.BLOCK_SIZE, (i + 1) * Utils.BLOCK_SIZE));
      fft = Utils.arr_op(fft, tmp_result, Utils.Operator.ADD);
    }

    fft = Utils.arr_op_const(fft, (float) Math.floor(arr_samples.length / Utils.BLOCK_SIZE), Utils.Operator.DIV);

    float[] p_filtre = new float[fft.length];//ArrayPool.GetInstance().getArray(ArrayPool.ThatBig.NbSmaple);//new float[fft.length];

    // creates the filter matrix
    for (int i = 0; i < filter.length; i++) {
      // new float[fTable.length]
      filter[i] = new float[fTable.length];//ArrayPool.GetInstance().getArray(ArrayPool.ThatBig.Lower2);
    }

    // init p_band
    for (int i = 0; i < p_band.length; i++) {
      p_band[i] = 0.0F;

    }

    for (int i = 0; i < filter.length; i++) {
      for (int j = 0; j < filter[i].length; j++) {
        if (standard) {
          double numerator = fTable[j] * fTable[j] - fc3[i] * fc3[i];
          double denumerator = (fu3[i] - fl3[i]) * fTable[j];
          filter[i][j] = (float) (1 / (1 + Math.pow(numerator / denumerator, 6)));
        } else {
          int c1 = (fTable[j] >= fl3[i]) ? 1 : 0;
          int c2 = (fTable[j] < fu3[i]) ? 1 : 0;
          filter[i][j] = 1 * c1 * c2;
        }
      }
      for (int j = 0; j < filter[i].length; j++) {
        p_filtre[j] = fft[j] * filter[i][j] * filter[i][j];
      }
      for (int j = 0; j < filter[i].length; j++) {
        p_band[i] += p_filtre[j];
      }

    }

    // Big cleanup of all the arrays used in this function
    // the cleanup will take effect when the GC is activated
    //	for(int i = 0; i < filter.length; i++){
    //  		ArrayPool.GetInstance().releaseArray(filter[i]);
    // 	}
    //  filter = null;
    //  fft = null;
//    	ArrayPool.GetInstance().releaseArray(p_filtre);

    return p_band;
  }
}