Edited FFT for android.

FFT.java

/*
 * Copyright (c) 2007 - 2008 by Damien Di Fede <[email protected]>
 * 
 * This program is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
 * 
 * This program 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 Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public License along with this program; if not, write to the Free
 * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/**
 * FFT stands for Fast Fourier Transform. It is an efficient way to calculate the Complex Discrete Fourier Transform. There is not
 * much to say about this class other than the fact that when you want to analyze the spectrum of an audio buffer you will almost
 * always use this class. One restriction of this class is that the audio buffers you want to analyze must have a length that is a
 * power of two. If you try to construct an FFT with a <code>timeSize</code> that is not a power of two, an
 * IllegalArgumentException will be thrown.
 *
 * @author Damien Di Fede
 * @see FourierTransform
 * @see <a href="http://www.dspguide.com/ch12.htm">The Fast Fourier Transform</a>
 */
public class FFT extends FourierTransform {
    /**
     * Constructs an FFT that will accept sample buffers that are <code>timeSize</code> long and have been recorded with a sample
     * rate of <code>sampleRate</code>. <code>timeSize</code> <em>must</em> be a power of two. This will throw an exception if it
     * is not.
     *
     * @param timeSize   the length of the sample buffers you will be analyzing
     * @param sampleRate the sample rate of the audio you will be analyzing
     */
    public FFT(int timeSize, float sampleRate) {
        super(timeSize, sampleRate);
        if ((timeSize & (timeSize - 1)) != 0)
            throw new IllegalArgumentException("FFT: timeSize must be a power of two.");
        buildReverseTable();
        buildTrigTables();
    }

    protected void allocateArrays() {
        spectrum = new float[timeSize / 2 + 1];
        real = new float[timeSize];
        imag = new float[timeSize];
    }

    // performs an in-place fft on the data in the real and imag arrays
    // bit reversing is not necessary as the data will already be bit reversed
    private void fft() {
        for (int halfSize = 1; halfSize < real.length; halfSize *= 2) {
            // float k = -(float)Math.PI/halfSize;
            // phase shift step
            // float phaseShiftStepR = (float)Math.cos(k);
            // float phaseShiftStepI = (float)Math.sin(k);
            // using lookup table
            float phaseShiftStepR = cos(halfSize);
            float phaseShiftStepI = sin(halfSize);
            // current phase shift
            float currentPhaseShiftR = 1.0f;
            float currentPhaseShiftI = 0.0f;
            for (int fftStep = 0; fftStep < halfSize; fftStep++) {
                for (int i = fftStep; i < real.length; i += 2 * halfSize) {
                    int off = i + halfSize;
                    float tr = (currentPhaseShiftR * real[off]) - (currentPhaseShiftI * imag[off]);
                    float ti = (currentPhaseShiftR * imag[off]) + (currentPhaseShiftI * real[off]);
                    real[off] = real[i] - tr;
                    imag[off] = imag[i] - ti;
                    real[i] += tr;
                    imag[i] += ti;
                }
                float tmpR = currentPhaseShiftR;
                currentPhaseShiftR = (tmpR * phaseShiftStepR) - (currentPhaseShiftI * phaseShiftStepI);
                currentPhaseShiftI = (tmpR * phaseShiftStepI) + (currentPhaseShiftI * phaseShiftStepR);
            }
        }
    }

    public void forward(float[] buffer) {
        if (buffer.length != timeSize) {
            throw new IllegalArgumentException("FFT.forward: The length of the passed sample buffer must be equal to timeSize().");
        }
        doWindow(buffer);
        // copy samples to real/imag in bit-reversed order
        bitReverseSamples(buffer);
        // perform the fft
        fft();
        // fill the spectrum buffer with amplitudes
        fillSpectrum();
    }

    private int[] reverse;

    private void buildReverseTable() {
        int N = timeSize;
        reverse = new int[N];

        // set up the bit reversing table
        reverse[0] = 0;
        for (int limit = 1, bit = N / 2; limit < N; limit <<= 1, bit >>= 1)
            for (int i = 0; i < limit; i++)
                reverse[i + limit] = reverse[i] + bit;
    }

    // copies the values in the samples array into the real array
    // in bit reversed order. the imag array is filled with zeros.
    private void bitReverseSamples(float[] samples) {
        for (int i = 0; i < samples.length; i++) {
            real[i] = samples[reverse[i]];
            imag[i] = 0.0f;
        }
    }

    // lookup tables

    private float[] sinlookup;
    private float[] coslookup;

    private float sin(int i) {
        return sinlookup[i];
    }

    private float cos(int i) {
        return coslookup[i];
    }

    private void buildTrigTables() {
        int N = timeSize;
        sinlookup = new float[N];
        coslookup = new float[N];
        for (int i = 0; i < N; i++) {
            sinlookup[i] = (float) Math.sin(-(float) Math.PI / i);
            coslookup[i] = (float) Math.cos(-(float) Math.PI / i);
        }
    }

}

FourierTransform.java

/*
 * Copyright (c) 2007 - 2008 by Damien Di Fede <[email protected]>
 * 
 * This program is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
 * 
 * This program 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 Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public License along with this program; if not, write to the Free
 * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
/**
 * A Fourier Transform is an algorithm that transforms a signal in the time domain, such as a sample buffer, into a signal in the
 * frequency domain, often called the spectrum. The spectrum does not represent individual frequencies, but actually represents
 * frequency bands centered on particular frequencies. The center frequency of each band is usually expressed as a fraction of the
 * sampling rate of the time domain signal and is equal to the index of the frequency band divided by the total number of bands.
 * The total number of frequency bands is usually equal to the length of the time domain signal, but access is only provided to
 * frequency bands with indices less than half the length, because they correspond to frequencies below the <a
 * href="http://en.wikipedia.org/wiki/Nyquist_frequency">Nyquist frequency</a>. In other words, given a signal of length
 * <code>N</code>, there will be <code>N/2</code> frequency bands in the spectrum.
 * <p/>
 * As an example, if you construct a FourierTransform with a <code>timeSize</code> of 1024 and and a <code>sampleRate</code> of
 * 44100 Hz, then the spectrum will contain values for frequencies below 22010 Hz, which is the Nyquist frequency (half the sample
 * rate). If you ask for the value of band number 5, this will correspond to a frequency band centered on
 * <code>5/1024 * 44100 = 0.0048828125 * 44100 = 215 Hz</code>. The width of that frequency band is equal to <code>2/1024</code>,
 * expressed as a fraction of the total bandwidth of the spectrum. The total bandwith of the spectrum is equal to the Nyquist
 * frequency, which in this case is 22100, so the bandwidth is equal to about 50 Hz. It is not necessary for you to remember all
 * of these relationships, though it is good to be aware of them. The function <code>getFreq()</code> allows you to query the
 * spectrum with a frequency in Hz and the function <code>getBandWidth()</code> will return the bandwidth in Hz of each frequency
 * band in the spectrum.
 * <p/>
 * <b>Usage</b>
 * <p/>
 * A typical usage of a FourierTransform is to analyze a signal so that the frequency spectrum may be represented in some way,
 * typically with vertical lines. You could do this in Processing with the following code, where <code>audio</code> is an
 * AudioSource and <code>fft</code> is an FFT (one of the derived classes of FourierTransform).
 * <p/>
 * <pre>
 * fft.forward(audio.left);
 * for (int i = 0; i < fft.specSize(); i++) {
 *  // draw the line for frequency band i, scaling it by 4 so we can see it a bit better
 *  line(i, height, i, height - fft.getBand(i) * 4);
 * }
 * </pre>
 * <p/>
 * <b>Windowing</b>
 * <p/>
 * Windowing is the process of shaping the audio samples before transforming them to the frequency domain. If you call the
 * <code>window()</code> function with an appropriate constant, such as FourierTransform.HAMMING, the sample buffers passed to the
 * object for analysis will be shaped by the current window before being transformed. The result of using a window is to reduce
 * the noise in the spectrum somewhat.
 * <p/>
 * <b>Averages</b>
 * <p/>
 * FourierTransform also has functions that allow you to request the creation of an average spectrum. An average spectrum is
 * simply a spectrum with fewer bands than the full spectrum where each average band is the average of the amplitudes of some
 * number of contiguous frequency bands in the full spectrum.
 * <p/>
 * <code>linAverages()</code> allows you to specify the number of averages that you want and will group frequency bands into
 * groups of equal number. So if you have a spectrum with 512 frequency bands and you ask for 64 averages, each average will span
 * 8 bands of the full spectrum.
 * <p/>
 * <code>logAverages()</code> will group frequency bands by octave and allows you to specify the size of the smallest octave to
 * use (in Hz) and also how many bands to split each octave into. So you might ask for the smallest octave to be 60 Hz and to
 * split each octave into two bands. The result is that the bandwidth of each average is different. One frequency is an octave
 * above another when it's frequency is twice that of the lower frequency. So, 120 Hz is an octave above 60 Hz, 240 Hz is an
 * octave above 120 Hz, and so on. When octaves are split, they are split based on Hz, so if you split the octave 60-120 Hz in
 * half, you will get 60-90Hz and 90-120Hz. You can see how these bandwidths increase as your octave sizes grow. For instance, the
 * last octave will always span <code>sampleRate/4 - sampleRate/2</code>, which in the case of audio sampled at 44100 Hz is
 * 11025-22010 Hz. These logarithmically spaced averages are usually much more useful than the full spectrum or the linearly
 * spaced averages because they map more directly to how humans perceive sound.
 * <p/>
 * <code>calcAvg()</code> allows you to specify the frequency band you want an average calculated for. You might ask for 60-500Hz
 * and this function will group together the bands from the full spectrum that fall into that range and average their amplitudes
 * for you.
 * <p/>
 * If you don't want any averages calculated, then you can call <code>noAverages()</code>. This will not impact your ability to
 * use <code>calcAvg()</code>, it will merely prevent the object from calculating an average array every time you use
 * <code>forward()</code>.
 * <p/>
 * <b>Inverse Transform</b>
 * <p/>
 * FourierTransform also supports taking the inverse transform of a spectrum. This means that a frequency spectrum will be
 * transformed into a time domain signal and placed in a provided sample buffer. The length of the time domain signal will be
 * <code>timeSize()</code> long. The <code>set</code> and <code>scale</code> functions allow you the ability to shape the spectrum
 * already stored in the object before taking the inverse transform. You might use these to filter frequencies in a spectrum or
 * modify it in some other way.
 *
 * @author Damien Di Fede
 * @see <a href="http://www.dspguide.com/ch8.htm">The Discrete Fourier Transform</a>
 */
public abstract class FourierTransform {
    /**
     * A constant indicating no window should be used on sample buffers.
     */
    private static final int NONE = 0;
    /**
     * A constant indicating a Hamming window should be used on sample buffers.
     */
    private static final int HAMMING = 1;
    private static final int LINAVG = 2;
    private static final int LOGAVG = 3;
    private static final int NOAVG = 4;
    private static final float TWO_PI = (float) (2 * Math.PI);
    int timeSize;
    private int sampleRate;
    private float bandWidth;
    private int whichWindow;
    float[] real;
    float[] imag;
    float[] spectrum;
    private float[] averages;
    private int whichAverage;
    int octaves;
    int avgPerOctave;

    /**
     * Construct a FourierTransform that will analyze sample buffers that are <code>ts</code> samples long and contain samples with
     * a <code>sr</code> sample rate.
     *
     * @param ts the length of the buffers that will be analyzed
     * @param sr the sample rate of the samples that will be analyzed
     */
    FourierTransform(int ts, float sr) {
        timeSize = ts;
        sampleRate = (int) sr;
        bandWidth = (2f / timeSize) * ((float) sampleRate / 2f);
        noAverages();
        allocateArrays();
        whichWindow = NONE;
    }

    // allocating real, imag, and spectrum are the responsibility of derived
    // classes
    // because the size of the arrays will depend on the implementation being used
    // this enforces that responsibility
    protected abstract void allocateArrays();

    // fill the spectrum array with the amps of the data in real and imag
    // used so that this class can handle creating the average array
    // and also do spectrum shaping if necessary
    void fillSpectrum() {
        for (int i = 0; i < spectrum.length; i++) {
            spectrum[i] = (float) Math.sqrt(real[i] * real[i] + imag[i] * imag[i]);
        }

        if (whichAverage == LINAVG) {
            int avgWidth = spectrum.length / averages.length;
            for (int i = 0; i < averages.length; i++) {
                float avg = 0;
                int j;
                for (j = 0; j < avgWidth; j++) {
                    int offset = j + i * avgWidth;
                    if (offset < spectrum.length) {
                        avg += spectrum[offset];
                    } else {
                        break;
                    }
                }
                avg /= j + 1;
                averages[i] = avg;
            }
        } else if (whichAverage == LOGAVG) {
            for (int i = 0; i < octaves; i++) {
                float lowFreq, hiFreq, freqStep;
                if (i == 0) {
                    lowFreq = 0;
                } else {
                    lowFreq = (sampleRate / 2) / (float) Math.pow(2, octaves - i);
                }
                hiFreq = (sampleRate / 2) / (float) Math.pow(2, octaves - i - 1);
                freqStep = (hiFreq - lowFreq) / avgPerOctave;
                float f = lowFreq;
                for (int j = 0; j < avgPerOctave; j++) {
                    int offset = j + i * avgPerOctave;
                    averages[offset] = calcAvg(f, f + freqStep);
                    f += freqStep;
                }
            }
        }
    }

    /**
     * Sets the object to not compute averages.
     */
    private void noAverages() {
        averages = new float[0];
        whichAverage = NOAVG;
    }

    void doWindow(float[] samples) {
        switch (whichWindow) {
            case HAMMING:
                hamming(samples);
                break;
        }
    }

    // windows the data in samples with a Hamming window
    private void hamming(float[] samples) {
        for (int i = 0; i < samples.length; i++) {
            samples[i] *= (0.54f - 0.46f * Math.cos(TWO_PI * i / (samples.length - 1)));
        }
    }

    /**
     * Returns the amplitude of the requested frequency band.
     *
     * @param i the index of a frequency band
     * @return the amplitude of the requested frequency band
     */
    private float getBand(int i) {
        if (i < 0) i = 0;
        if (i > spectrum.length - 1) i = spectrum.length - 1;
        return spectrum[i];
    }

    /**
     * Returns the width of each frequency band in the spectrum (in Hz). It should be noted that the bandwidth of the first and
     * last frequency bands is half as large as the value returned by this function.
     *
     * @return the width of each frequency band in Hz.
     */
    private float getBandWidth() {
        return bandWidth;
    }

    /**
     * Returns the index of the frequency band that contains the requested frequency.
     *
     * @param freq the frequency you want the index for (in Hz)
     * @return the index of the frequency band that contains freq
     */
    private int freqToIndex(float freq) {
        // special case: freq is lower than the bandwidth of spectrum[0]
        if (freq < getBandWidth() / 2) return 0;
        // special case: freq is within the bandwidth of spectrum[spectrum.length - 1]
        if (freq > sampleRate / 2 - getBandWidth() / 2) return spectrum.length - 1;
        // all other cases
        float fraction = freq / (float) sampleRate;
        return Math.round(timeSize * fraction);
    }

    /**
     * Gets the amplitude of the requested frequency in the spectrum.
     *
     * @param freq the frequency in Hz
     * @return the amplitude of the frequency in the spectrum
     */
    public float getFreq(float freq) {
        return getBand(freqToIndex(freq));
    }

    /**
     * Calculate the average amplitude of the frequency band bounded by <code>lowFreq</code> and <code>hiFreq</code>, inclusive.
     *
     * @param lowFreq the lower bound of the band
     * @param hiFreq  the upper bound of the band
     * @return the average of all spectrum values within the bounds
     */
    private float calcAvg(float lowFreq, float hiFreq) {
        int lowBound = freqToIndex(lowFreq);
        int hiBound = freqToIndex(hiFreq);
        float avg = 0;
        for (int i = lowBound; i <= hiBound; i++) {
            avg += spectrum[i];
        }
        avg /= (hiBound - lowBound + 1);
        return avg;
    }

    /**
     * Performs a forward transform on <code>buffer</code>.
     *
     * @param buffer the buffer to analyze
     */
    public abstract void forward(float[] buffer);

}