FFT 的C 語言實現

標籤：fft   c language   

FFT 的C 語言實現

說好的C 語言實現，必須搞定它！

理論介紹:

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

這裡有之前matlab & Octave 的實現

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

先介紹一下總體的實現檔案

最主要的是fft.c這個檔案是演算法實現的核心

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;}

FFT 的C 語言實現