/**
*
* Funf: Open Sensing Framework
* Copyright (C) 2010-2011 Nadav Aharony, Wei Pan, Alex Pentland.
* Acknowledgments: Alan Gardner
* Contact: nadav@media.mit.edu
*
* This file is part of Funf.
*
* Funf is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation, either version 3 of
* the License, or (at your option) any later version.
*
* Funf 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.
*
* You should have received a copy of the GNU Lesser General Public
* License along with Funf. If not, see <http://www.gnu.org/licenses/>.
*
*/
package edu.mit.media.funf.math;
public class FFT
{
int n, m;
// Lookup tables. Only need to recompute when size of FFT changes.
double[] cos;
double[] sin;
public FFT(int n)
{
this.n = n;
this.m = (int)(Math.log(n) / Math.log(2));
// Make sure n is a power of 2
if (n != (1<<m))
{
throw new RuntimeException("FFT length must be power of 2");
}
// precompute tables
cos = new double[n/2];
sin = new double[n/2];
for(int i=0; i<n/2; i++)
{
cos[i] = Math.cos(-2*Math.PI*i/n);
sin[i] = Math.sin(-2*Math.PI*i/n);
}
}
/***************************************************************
* fft.c
* Douglas L. Jones
* University of Illinois at Urbana-Champaign
* January 19, 1992
* http://cnx.rice.edu/content/m12016/latest/
*
* fft: in-place radix-2 DIT DFT of a complex input
*
* input:
* n: length of FFT: must be a power of two
* m: n = 2**m
* input/output
* x: double array of length n with real part of data
* y: double array of length n with imag part of data
*
* Permission to copy and use this program is granted
* as long as this header is included.
****************************************************************/
public void fft(double[] re, double[] im)
{
int i,j,k,n1,n2,a;
double c,s,t1,t2;
// Bit-reverse
j = 0;
n2 = n/2;
for (i=1; i < n - 1; i++)
{
n1 = n2;
while ( j >= n1 )
{
j = j - n1;
n1 = n1/2;
}
j = j + n1;
if (i < j)
{
t1 = re[i];
re[i] = re[j];
re[j] = t1;
t1 = im[i];
im[i] = im[j];
im[j] = t1;
}
}
// FFT
n1 = 0;
n2 = 1;
for (i=0; i < m; i++)
{
n1 = n2;
n2 = n2 + n2;
a = 0;
for (j=0; j < n1; j++)
{
c = cos[a];
s = sin[a];
a += 1 << (m-i-1);
for (k=j; k < n; k=k+n2)
{
t1 = c*re[k+n1] - s*im[k+n1];
t2 = s*re[k+n1] + c*im[k+n1];
re[k+n1] = re[k] - t1;
im[k+n1] = im[k] - t2;
re[k] = re[k] + t1;
im[k] = im[k] + t2;
}
}
}
}
}