C-language implementation of FFT

C-language implementation of FFT

A good C language implementation must be done!

Theoretical introduction:

Http://blog.csdn.net/cinmyheart/article/details/39052739

Here is the implementation of MATLAB & Ave ave.

Http://blog.csdn.net/cinmyheart/article/details/39042623

First, we will introduce the overall implementation file.

The most important thing is that the FFT. c file is the core of algorithm implementation.

FFT. h

/***************************************************************code writer:EOFcode file:fft.hcode date:2014.09.17e-mail:[email protected]****************************************************************/#ifndef _FFT_IN_C_H#define _FFT_IN_C_H#include <stdio.h>#include <stdlib.h>#include <math.h>#define PI 3.1415926#define DEBUGstruct complex_number{double real;double imagine;};struct signal{int size;//how many points in this domain.struct complex_number points[0];};int reverse_bits(int num,int bits);int get_r_in_Wn(int k, int m, int bits);void init_signal(struct signal* p_signal,double* array,int size);struct signal* fft(struct signal* p_signal);struct complex_number complex_mul(struct complex_number* x,struct complex_number* y);struct complex_number complex_add(struct complex_number* x, struct complex_number *y);struct complex_number complex_sub(struct complex_number* x, struct complex_number *y);void show_signal(struct signal* const signal);void show_complex_number(struct complex_number * x);#endif

Complex_add.c

/***************************************************************code writer:EOFcode file:complex_add.ccode date:2014.09.17e-mail:[email protected]code purpose:We need add operation between complex number, so here is my method.If you find something wrong with my code, please touchme by e-mail. Thank you :)****************************************************************/#include "fft.h"#include "fft.h"struct complex_number complex_add(struct complex_number* x, struct complex_number *y){struct complex_number ret;if(!x || !y){printf("You passed NULL into %s()\n",__FUNCTION__);goto out ;}ret.real    = x->real    + y->real;ret.imagine = x->imagine + y->imagine;out:return ret;}

Complex_mul.c

/***************************************************************code writer:EOFcode file:complex_mul.ccode date:2014.09.17e-mail:[email protected]code purpose:We need multiple(*) operation between complex number, so here is my method.If you find something wrong with my code, please touchme by e-mail. Thank you :)****************************************************************/#include "fft.h"struct complex_number complex_mul(struct complex_number* x,struct complex_number *y){struct complex_number ret;if(!x || !y){printf("You passed NULL into %s()\n",__FUNCTION__);goto out ;}ret.real    = (x->real) * (y->real)    - (x->imagine)*(y->imagine);ret.imagine = (x->real) * (y->imagine) + (x->imagine)* (y->real);out:return ret;}

Complex_sub.c

#include "fft.h"struct complex_number complex_sub(struct complex_number* x, struct complex_number *y){struct complex_number ret;if(!x || !y){printf("You passed NULL into %s()\n",__FUNCTION__);goto out ;}ret.real    = x->real    - y->real;ret.imagine = x->imagine - y->imagine;out:return ret;}

Get_r_in_wn.c

/******************************************************************************code writer:EOFcode file:get_r_in_Wn.ccode date:2014.09.17e-mail:[email protected]Input Parameter : @k, the location of input signal                          @m, the current layyer                          @N, the total number of inputed signal                          @bits, how many bits should be used to represent 'N'Output Parameter: @ret , the value of 'r'*******************************************************************************/int get_r_in_Wn(int k, int m, int bits){int tmp = k<<(bits-m);return tmp&((1<<m) -1);}

Reverse_bits.c

/***************************************************************code writer:EOFcode file:reverse_bits.ccode date:2014.09.17e-mail:[email protected]code purpose: Reverse the bits of input numberIf you find something wrong with my code, please touchme by e-mail. Thank you :)****************************************************************/int reverse_bits(int num,int bits){int ret  = 0;int copy_num = 0;for(ret = 0,copy_num = num; bits > 0; bits--){ret  += (copy_num % 2) * (1<<(bits-1));copy_num >>= 1;}return ret;}

Show_complex_number.c

/***************************************************************code writer:EOFcode file:show_complex_number.ccode date:2014.09.17e-mail:[email protected]If you find something wrong with my code, please touchme by e-mail. Thank you :)****************************************************************/#include "fft.h"void show_complex_number(struct complex_number * x){printf("real:%lf  imagine:%lf\n",x->real,x->imagine);}

Show_signal.c

/***************************************************************code writer:EOFcode file:show_signal.ccode date:2014.09.17e-mail:[email protected]code purpose:If you want to see the detail about signal that @p_signal point to, just call this API.If you find something wrong with my code, please touchme by e-mail. Thank you :)****************************************************************/#include "fft.h"void show_signal(struct signal* const p_signal){if(!p_signal){printf("You passed NULL into function %s  line:%d\n",__FUNCTION__,__LINE__);return ;}int tmp = 0;for(tmp = 0; tmp < p_signal->size;tmp++){printf("point %d real : %lf imagine %lf\n",tmp,p_signal->points[tmp].real,p_signal->points[tmp].imagine);}printf("\n\n");}

FFT. c

/***************************************************************code writer:EOFcode file:fft.ccode date:2014.09.17e-mail:[email protected]code purpose:This code core-part of fft in my implementation.If you find something wrong with my code, please touchme by e-mail. Thank you :)****************************************************************/#include "fft.h"struct signal* fft(struct signal* p_signal){struct signal*  p_input_signal = (struct signal*) malloc(sizeof(struct complex_number)*(p_signal->size) + sizeof(p_signal->size));struct signal*  p_out_put_signal =  (struct signal*)malloc(sizeof(struct complex_number)*(p_signal->size) + sizeof(p_signal->size));*p_input_signal   = *p_signal;*p_out_put_signal = *p_signal;int tmp   = 0;int index = 0;int bits  = 0;int layyer= 0;int selected_point = 0;int pre_half   = 0;int    r  = 0;double x  = 0;struct complex_number W_rN ;struct complex_number complex_tmp ;/***We caculate how many bits should be used to ** represent the size of the number of input signal.*/for(tmp = p_signal->size-1;tmp > 0;tmp>>=1){bits++;}/***We should re-sequence the input signal** by bit-reverse.*/for(tmp = 0;tmp < p_signal->size;tmp++){index = reverse_bits(tmp,bits);p_input_signal->points[tmp] = p_signal->points[index];#ifdef DEBUGprintf(" tmp:%5d index:%5d  ",tmp,index);show_complex_number(&p_signal->points[index]);#endif}for(layyer = 1;layyer <= bits;layyer++){#ifdef DEBUGprintf("layyer %d input\n",layyer);show_signal(p_input_signal);#endiffor(selected_point = 0;selected_point < p_signal->size;selected_point += 1<<(layyer)){for(pre_half = selected_point,r = 0,x = 0;    pre_half < (selected_point + (1<<(layyer-1)));    pre_half++){r = get_r_in_Wn(pre_half,layyer,bits);#ifdef DEBUGprintf("r: %d\n",r);#endifx = -2*PI*r/((double)(p_input_signal->size));W_rN.real    = cos(x);W_rN.imagine = sin(x);complex_tmp = complex_mul(&W_rN , &(p_input_signal->points[pre_half + (1<<(layyer-1))])  );#ifdef DEBUGshow_complex_number(&complex_tmp);#endifp_out_put_signal->points[pre_half] = complex_add(&p_input_signal->points[pre_half],&complex_tmp);p_out_put_signal->points[pre_half + (1<<(layyer-1))] = complex_sub(&p_input_signal->points[pre_half],&complex_tmp);}}#ifdef DEBUGprintf("layyer%d output\n",layyer);show_signal(p_out_put_signal);#endiffor(tmp = 0;tmp < p_out_put_signal->size;tmp++){p_input_signal->points[tmp] = p_out_put_signal->points[tmp];}}free(p_input_signal);return p_out_put_signal;}

C-language implementation of FFT