Click here to Skip to main content
15,077,377 members
Please Sign up or sign in to vote.
0.00/5 (No votes)
See more:
Hello,everyone,recently I had met a problem about FFT algorithm:I tried to program the algorithm in C with reference to the textbook,but it didn't work as the matlab' fft function.
Firstly,I defined a data type for complex operation in "Complex.h"
C++


CSS
typedef struct
{
    double x_real;
    double x_imag;
}Complex;
// the multiplication for complex x1 and x2
Complex CMul(Complex x1,Complex x2)
{

    Complex result;
    result.x_real=x1.x_real*x2.x_real-x1.x_imag*x2.x_imag;
    result.x_imag=x1.x_real*x2.x_imag+x1.x_imag*x2.x_real;
    return  result;
}
//the addition for the complex x1 and x2
Complex CAdd(Complex x1,Complex x2)
{

    Complex result;
    result.x_real=x1.x_real+x2.x_real;
    result.x_imag=x1.x_imag+x2.x_imag;
    return  result;
}
//the subtract for the complex x1 and x2
Complex CSub(Complex x1,Complex x2)
{

    Complex result;
    result.x_real=x1.x_real-x2.x_real;
    result.x_imag=x1.x_imag-x2.x_imag;
    return  result;
}

Then I began to program the FFT in"fft.cpp"
C++


#include <stdio.h>
#include <math.h>
#include "Complex.h"
#define pi 3.1415926
//bit-reverse operation for the input 
void RSort(Complex x[],int length )
{
    int LH,J,N1,I,K;
	Complex T;
	LH=length/2;
	J=LH;
	N1=length-2;
	for (I=1;I<=N1;I++)
	{
		if (I<J)
		{
           T=x[I];
		   x[I]=x[J];
		   x[J]=T;
		}
		K=LH;
		while(J>=K)
		{
			J=J-K;
			K=K/2;
		}
		J=J+K;
	}
}
/*
void FFT(Complex x[],int N ,int M)---compute the FFT of the origin signal
x[]--the input of origin and output of its FFT
N--the length of x[];
M--the power of N;
*/
void FFT(Complex x[],int N ,int M)
{
/*
    L--the level of the Butterfly operation
    B--the distance between two signal in every Butterfly operation
    wn[]---the rotation factor
    p--the exponent of rotation factor
*/
    int L,B,J,p,k,i;
    Complex wn[4];
   //to get the rotation factor
    for (i=0;i<N/2;i++)
    {
	wn[i].x_real=cos(2*pi*i/N);
	wn[i].x_imag=-sin(2*pi*i/N);
    }
   //bit-reversal for the input signal
   RSort(x,N);
  //the main part of the FFT algorithm
   for (L=1;L<=M;L++)
	{
		B=int(pow(2,L-1));
		for (J=0;J<=B-1;J++)
		{
			p=int(J*pow(2,M-L));
			for (k=J;k<=N-1;)
			{
				x[k]=CAdd(x[k],CMul(wn[p],x[k+B]));
				x[k+B]=CSub(x[k],CMul(wn[p],x[k+B]));
				k=k+int(pow(2,L));
			}
		}
	}

}
void main()
{
    Complex *x;
	int i;
	for (i=0;i<=15;i++)
	{
		x[i].x_real=i;
		x[i].x_imag=0;
	}
    printf("the input signal:");
	for (i=0;i<=15;i++)
	{
		printf("%lf  ",x[i].x_real);
	}
	printf("\n");
	FFT(x,16,4);
    printf("the result after FFT:");
	for (i=0;i<=15;i++)
	{
		printf("%lf i%lf  ",x[i].x_real,x[i].x_imag);
	}
    printf("\n");
}

For convenience,I defined the input signal to be a signal x[]=[0,1,2,3,4,5,6,7].The output of the result after executing program is :

Quote:

the input signal:0.000000 1.000000 2.000000 3.000000 4.000000 5.000000 6.
000000 7.000000
the result after FFT:28.000000 i0.000000 -1.414213 i-4.828427 4.000000 i-6.00
0000 -0.707107 i-0.707107 12.000000 i0.000000 0.000000 i-2.000000 4.000000 i
0.000000 0.000000 i0.000000
Press any key to continue

And the result using FFT function in matlab is:
Quote:

>> x=[0 1 2 3 4 5 6 7];
>> y=fft(x)

y =

28.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 -4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i


I'm looking forward to your opinions,Thanks very much!^_^^_^...
Posted

Just a note: your question is not really a request for expert advice or something like that; this is really a request to do considerable amount of job for you: testing, debugging, detecting bugs and trying to fix them. Maybe someone will find something for you, but for many like myself it would be much easier to implement the algorithm from scratch.

All of this really goes beyond to format of CodeProject Quick Questions & Answers.

So I don't know how else to help…

At least use C++ complex library instead you rudimentary implementation. FFT has many implementations you could use. Read some articles, such as:
How to implement the FFT algorithm[^],
http://www.librow.com/articles/article-10[^].

You can find a lot more.

—SA
   
Comments
Rslboy 10-Nov-11 4:32am
   
Thanks for your advice!Yeah,It does have many jobs for me to do,not only just ask for answers,but aslo do it on my own! Thanks for the reference you give again!^_^...
Sergey Alexandrovich Kryukov 10-Nov-11 11:13am
   
Yes, I surely respect your desire to make things on your own; this is what I always try to do. I consider "reinventing the wheel" as a productive business, it already gave me some interesting results.
If you agree that my advice makes sense, please consider accepting it formally (green button) -- thanks.
--SA
Rslboy 12-Nov-11 7:20am
   
Thanks..^_^。。。
The errors which are logical errors occur in the function of FFT():
x[k]=CAdd(x[k],CMul(wn[p],x[k+B]));
x[k+B]=CSub(x[k],CMul(wn[p],x[k+B]));

These codes should be modified to(before these variable T1 and T2 of Complex should be defined):
C++
T2=CMul(wn[p],x[k+B]);
T1=x[k];
x[k]=CAdd(T1,T2);
x[k+B]=CSub(T1,T2);

And then,it works correctly!I felt excited about that.Now I wanted to give thanks to SAKryukov for his help firstly!:-):-)I had began to gain strong confidence in myself! Thanks again.
This program didn't work efficiently and it still needed to be improved well!Let's go for it!
   

This content, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)




CodeProject, 20 Bay Street, 11th Floor Toronto, Ontario, Canada M5J 2N8 +1 (416) 849-8900