Jump to content

  • Log In with Google      Sign In   
  • Create Account

We're offering banner ads on our site from just $5!

1. Details HERE. 2. GDNet+ Subscriptions HERE. 3. Ad upload HERE.


FFT


Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

  • You cannot reply to this topic
6 replies to this topic

#1 HelloSkitty   Members   -  Reputation: 152

Like
0Likes
Like

Posted 13 July 2011 - 11:39 AM

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?
A penny for my thoughts? Do you think I have only half a brain?

Not-so-proud owner of blog: http://agathokakologicalartolater.wordpress.com/

Sponsor:

#2 Brother Bob   Moderators   -  Reputation: 8571

Like
2Likes
Like

Posted 13 July 2011 - 12:31 PM

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.

#3 HelloSkitty   Members   -  Reputation: 152

Like
0Likes
Like

Posted 13 July 2011 - 01:08 PM

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

  • FFT.png

A penny for my thoughts? Do you think I have only half a brain?

Not-so-proud owner of blog: http://agathokakologicalartolater.wordpress.com/

#4 Brother Bob   Moderators   -  Reputation: 8571

Like
0Likes
Like

Posted 13 July 2011 - 01:29 PM

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.

#5 HelloSkitty   Members   -  Reputation: 152

Like
0Likes
Like

Posted 13 July 2011 - 02:15 PM

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.
A penny for my thoughts? Do you think I have only half a brain?

Not-so-proud owner of blog: http://agathokakologicalartolater.wordpress.com/

#6 Brother Bob   Moderators   -  Reputation: 8571

Like
1Likes
Like

Posted 13 July 2011 - 02:45 PM

In the pseudo code, you seem to calculate the row
Yk ← t + exp(−2πi k/N) Yk+N/2
but not
Yk+N/2 ← t − exp(−2πi k/N) Yk+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.

#7 HelloSkitty   Members   -  Reputation: 152

Like
0Likes
Like

Posted 14 July 2011 - 08:38 AM

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!
A penny for my thoughts? Do you think I have only half a brain?

Not-so-proud owner of blog: http://agathokakologicalartolater.wordpress.com/




Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.



PARTNERS