/*
* Copyright (c) 2007 by Damien Di Fede <ddf@compartmental.net>
*
* 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.
*/
package ddf.minim.analysis;
import ddf.minim.Minim;
/**
* 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.
*
* @see FourierTransform
* @see <a href="http://www.dspguide.com/ch12.htm">The Fast Fourier Transform</a>
*
* @author Damien Di Fede
*
*/
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];
}
public void scaleBand(int i, float s)
{
if (s < 0)
{
Minim.error("Can't scale a frequency band by a negative value.");
return;
}
if (spectrum[i] != 0)
{
real[i] /= spectrum[i];
imag[i] /= spectrum[i];
spectrum[i] *= s;
real[i] *= spectrum[i];
imag[i] *= spectrum[i];
}
if (i != 0 && i != timeSize / 2)
{
real[timeSize - i] = real[i];
imag[timeSize - i] = -imag[i];
}
}
public void setBand(int i, float a)
{
if (a < 0)
{
Minim.error("Can't set a frequency band to a negative value.");
return;
}
if (real[i] == 0 && imag[i] == 0)
{
real[i] = a;
spectrum[i] = a;
}
else
{
real[i] /= spectrum[i];
imag[i] /= spectrum[i];
spectrum[i] = a;
real[i] *= spectrum[i];
imag[i] *= spectrum[i];
}
if (i != 0 && i != timeSize / 2)
{
real[timeSize - i] = real[i];
imag[timeSize - i] = -imag[i];
}
}
// 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)
{
Minim
.error("FFT.forward: The length of the passed sample buffer must be equal to timeSize().");
return;
}
doWindow(buffer);
// copy samples to real/imag in bit-reversed order
bitReverseSamples(buffer);
// perform the fft
fft();
// fill the spectrum buffer with amplitudes
fillSpectrum();
}
/**
* Performs a forward transform on the passed buffers.
*
* @param buffReal the real part of the time domain signal to transform
* @param buffImag the imaginary part of the time domain signal to transform
*/
public void forward(float[] buffReal, float[] buffImag)
{
if (buffReal.length != timeSize || buffImag.length != timeSize)
{
Minim
.error("FFT.forward: The length of the passed buffers must be equal to timeSize().");
return;
}
setComplex(buffReal, buffImag);
bitReverseComplex();
fft();
fillSpectrum();
}
public void inverse(float[] buffer)
{
if (buffer.length > real.length)
{
Minim
.error("FFT.inverse: the passed array's length must equal FFT.timeSize().");
return;
}
// conjugate
for (int i = 0; i < timeSize; i++)
{
imag[i] *= -1;
}
bitReverseComplex();
fft();
// copy the result in real into buffer, scaling as we do
for (int i = 0; i < buffer.length; i++)
{
buffer[i] = real[i] / real.length;
}
}
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;
}
}
// bit reverse real[] and imag[]
private void bitReverseComplex()
{
float[] revReal = new float[real.length];
float[] revImag = new float[imag.length];
for (int i = 0; i < real.length; i++)
{
revReal[i] = real[reverse[i]];
revImag[i] = imag[reverse[i]];
}
real = revReal;
imag = revImag;
}
// 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);
}
}
}