C-language implementation of FFT

Source: Internet
Author: User
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

Related Article

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.