6 replies to this topic

[HelloSkitty]

I've been trying for a few days to program the algorithm for a Fast Fourier Transform, but it's been quite unsuccessful. I use 2 classes in java to implement it:

public class FFT { //Uses Cooley-Tukey Algorithm
	
	public static ComplexNumber findFFT(double[] data,int k) { //recursive function
		if (data.length==1) {  //base case
			return new ComplexNumber(data[0]);
		}
		double[] dataE=new double[data.length/2];
		double[] dataO=new double[data.length/2];
		for (int i=0;i<data.length;i++) {  //radix-2
			if (i%2==0) {
				dataE[i/2]=data[i];
			} else dataO[i/2]=data[i];
		}
		if (k<data.length/2) { //recurse
			return findFFT(dataE,k).add(findFFT(dataO,k).mult(new ComplexNumber(Math.cos(-2*k*Math.PI/data.length),Math.sin(-2*k*Math.PI/data.length))));
		}
		k-=data.length/2;
		return findFFT(dataE,k).add(findFFT(dataO,k).mult(new ComplexNumber(Math.cos(-2*k*Math.PI/data.length),Math.sin(-2*k*Math.PI/data.length))));
	}
	
	public static double[] getFFT(double[] data) { //data length is a power of 2
		double[] fft=new double[data.length];
		for (int i=0;i<fft.length;i++) {  //evaluate all of the points
			System.out.println(data.length-i);
			fft[i]=findFFT(data,i).abs();
		}
		return fft;
	}
}

public class ComplexNumber {
	
	double real;
	double imaginary;
	
	public ComplexNumber() {
		real=0;
		imaginary=0;
	}
	
	public ComplexNumber(double d) {
		real=d;
		imaginary=0;
	}
	
	public ComplexNumber(double r, double i) {
		real=r;
		imaginary=i;
	}
	
	public ComplexNumber add(ComplexNumber other) {
		return new ComplexNumber(real+other.real,imaginary+other.imaginary);
	}
	
	public ComplexNumber mult(ComplexNumber other) {
		return new ComplexNumber(real*other.real-imaginary*other.imaginary,real*other.imaginary+imaginary*other.real);
	}
	
	public double abs() {
		return Math.sqrt(real*real+imaginary*imaginary);
	}
	
	public String toString() {
		return ""+real+""+imaginary+"i";
	}
	
}

A third class is used to generate a vector of numbers:

double[] time=new double[8192];
double[] signal=new double[8192];
for (int i=0;i<time.length;i++) {
		time[i]=i/4096.0;
		signal[i]=2*Math.sin(2*Math.PI*time[i]);
}

This should in theory, give a bunch of 0's, and one 2, when I call the FFT.getFFT method, but it doesn't. The source for my algorithm is http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#The_radix-2_DIT_case

So I'm not sure if it's a small coding error, or if I completely misread the algorithm. Is there anything wrong with the code?

[Brother Bob]

Your algorithm seems to be fundamentally broken just by looking at your getFFT function and how it relates to the findFFT function. The getFFT function loops through each frequency bin individually to calculate the corresponding frequency bin via findFFT. This is, at best, the naive DFT implementation, which is terribly inefficient. The FFT decomposition calculates all bins at the same time, not one by one, by recursively splitting the FFT. The very link you posted has a pseudo code implementing the Cooley-Tukey algorithm, and you can see how the FFT is split into smaller FFTs, and there is no loop that calculates the bins by their own. All bins are the result of complete recursion, and no bins are completed before any other bins.

Furthermore, the expected result is probably not what you expect it to be. The algorithm looks like an non-normalized FFT, so you should not get one 2 and a bunch of 0. First of all, the amplitude of the sine is 8196 times what you expect, because there is nothing in the algorithm that normalizes the output with regards to the length of the data. As you add more data, the complex exponential correlators will keep multiplying and add with the signal, and longer signals will result in a larger final correlation value. Second, the sine consists of two complex exponentials, one with positive and one with negative frequency, so you will have two bins with half the signal each. Thus, you will have all but two zeros whose magnitude is 2*8196/2.

[HelloSkitty]

Thank you Bob. I did get two peaks at over 8000 when I plotted my fft numbers (I originally took 3 square roots because I didn't know what to do about them), so I added code in the main method to divide the fft numbers by time.length after computation. Now it looks like the picture I attached (a plot of fft result to time). Even if the two pieces have half the signal each, the height of the peaks are 0.85 (I put a red dot where I put my mouse at the screenshot time, coordinates are (1,0.85)), twice that is 1.7. But that doesn't make sense since my function was 2*sin(2*PI*t). And there are still many bars of non-zero magnitude which I know shouldn't be there.

Also, to your first point, I was planning on later evaluating all of the numbers at once, rather than one-by-one, because right now my program runs at O(N^2logN), which I know is even slower than the literal summation computation of a DFT, but I wanted to first make sure I could get each point as an individual correct, before taking them as a group.

Attached Thumbnails

• 

[Brother Bob]

It just doesn't make sense to implement the FFT and calculate only one bin. You need to calculate all of them to get just one value, and once you have that one value, you already have all of them. So to get that one value you're asking for just for testing, you need to fully implement the FFT algorithm.

[HelloSkitty]

Okay, I changed my FFT class to make it use arrays instead of individual numbers, so it evaluates everything at the same time instead of each number individually (ComplexNumber class is unchanged, and main was edited to call this one now).

import java.lang.*;
import java.util.*;

public class FFT2 { //Uses Cooley-Tukey Algorithm
	
	public static ComplexNumber[] findFFT2(double[] data,int[] k) { //recursive function
		if (data.length==1) {  //base case
			return new ComplexNumber[]{new ComplexNumber(data[0])};
		}
		double[] dataE=new double[data.length/2];
		double[] dataO=new double[data.length/2];
		for (int i=0;i<data.length;i++) {  //radix-2
			if (i%2==0) {
                 dataE[i/2]=data[i];
			} else dataO[i/2]=data[i];
		}
		int[] kH=new int[k.length/2];
		for (int i=0;i<kH.length;i++) {
			kH[i]=i;
		}
		double[] f=mult(k,-2*Math.PI/data.length);
		double[] cos=cos(f);
		double[] sin=sin(f);
		ComplexNumber[] cis=cis(cos,sin);
		ComplexNumber[] half=add(findFFT2(dataE,kH),mult(findFFT2(dataO,k),cis));
		ComplexNumber[] result=new ComplexNumber[half.length*2];
		for (int i=0;i<result.length;i++) {
			result[i]=half[i%half.length];
		}
		return result;
	}
	
	public static double[] getFFT2(double[] data) {
		int[] k=new int[data.length];
		for (int i=0;i<k.length;i++) {
			k[i]=i;
		}
		ComplexNumber[] fft=findFFT2(data,k);
		double[] rfft=new double[fft.length];
		for (int i=0;i<fft.length;i++) {
			rfft[i]=fft[i].abs();
		}
		return rfft;
	}
	
	public static double[] mult(int[] a,double b) {
		double[] c=new double[a.length];
		for (int i=0;i<c.length;i++) {
			c[i]=b*a[i];
		}
		return c;
	}
	
	public static ComplexNumber[] add(ComplexNumber[] a,ComplexNumber[] b) {
		ComplexNumber[] c=new ComplexNumber[a.length];
		for (int i=0;i<c.length;i++) {
			c[i]=a[i].add(b[i]);
		}
		return c;
	}
	
	public static ComplexNumber[] mult(ComplexNumber[] a,ComplexNumber[] b) {
		ComplexNumber[] c=new ComplexNumber[a.length];
		for (int i=0;i<c.length;i++) {
			c[i]=a[i].mult(b[i]);
		}
		return c;
	}
	
	public static double[] cos(double[] a) {
		double[] b=new double[a.length];
		for (int i=0;i<b.length;i++) {
			b[i]=Math.cos(a[i]);
		}
		return b;
	}
	
	public static double[] sin(double[] a) {
		double[] b=new double[a.length];
		for (int i=0;i<b.length;i++) {
			b[i]=Math.sin(a[i]);
		}
		return b;
	}
	
	public static ComplexNumber[] cis(double[] a,double[] b) {
		ComplexNumber[] c=new ComplexNumber[a.length];
		for (int i=0;i<c.length;i++) {
			c[i]=new ComplexNumber(a[i],b[i]);
		}
		return c;
	}

}

It's much faster now (almost instant), but it still gives me the picture in my second post when I plot the fft vector to my time vector, which means something is still wrong with my interpretation of the wikipedian algorithm.

[Brother Bob]

In the pseudo code, you seem to calculate the row
Y_k ← t + exp(−2πi k/N) Y_k+N/2
but not
Y_k+N/2 ← t − exp(−2πi k/N) Y_k+N/2
so you're missing half the bins. Instead you replicate the first half into the second half. Ignoring syntax issues, I would expect it to looks something like this.

ComplexNumber[] fftEven = findFFT2(dataE,kH);
ComplexNumber[] fftOdd = findFFT2(dataO,k);
ComplexNumber[] halfLow = add(fftEven, mult(fftOdd, cis));
ComplexNumber[] halfHigh = sub(fftEven, mult(fftOdd, cis));

Then concatenate halfLow and halfHigh for the recursion result. I may be missing some detail though, but either way, your code appears to replicates the lower half into the upper half, contrary to the pseudo code, and that immediately looks suspicious to me.

[HelloSkitty]

Yes, the problem was that I added both times... I fixed it and it seems to work well now. Thank you again for your help and the fast replies!