FFT 的C 语言
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 : EOF
code file : fft.h
code date : 2014.09.17
e-mail : jasonleaster@gmail.com
****************************************************************/
#ifndef _FFT_IN_C_H
#define _FFT_IN_C_H
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.1415926
#define DEBUG
struct 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);
#endifcomplex_add.c
/***************************************************************
code writer : EOF
code file : complex_add.c
code date : 2014.09.17
e-mail : jasonleaster@gmail.com
code purpose :
We need add operation between complex number, so here is
my method.
If you find something wrong with my code, please touch
me 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 : EOF
code file : complex_mul.c
code date : 2014.09.17
e-mail : jasonleaster@gmail.com
code purpose :
We need multiple(*) operation between complex number, so here is
my method.
If you find something wrong with my code, please touch
me 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 : EOF
code file : get_r_in_Wn.c
code date : 2014.09.17
e-mail : jasonleaster@gmail.com
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);
}
/***************************************************************
code writer : EOF
code file : reverse_bits.c
code date : 2014.09.17
e-mail : jasonleaster@gmail.com
code purpose :
Reverse the bits of input number
If you find something wrong with my code, please touch
me 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 : EOF
code file : show_complex_number.c
code date : 2014.09.17
e-mail : jasonleaster@gmail.com
If you find something wrong with my code, please touch
me 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);
}
/***************************************************************
code writer : EOF
code file : show_signal.c
code date : 2014.09.17
e-mail : jasonleaster@gmail.com
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 touch
me 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 : EOF
code file : fft.c
code date : 2014.09.17
e-mail : jasonleaster@gmail.com
code purpose :
This code core-part of fft in my implementation.
If you find something wrong with my code, please touch
me 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 DEBUG
printf(" tmp:%5d index:%5d ",tmp,index);
show_complex_number(&p_signal->points[index]);
#endif
}
for(layyer = 1;layyer <= bits;layyer++)
{
#ifdef DEBUG
printf("layyer %d input\n",layyer);
show_signal(p_input_signal);
#endif
for(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 DEBUG
printf("r: %d\n",r);
#endif
x = -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 DEBUG
show_complex_number(&complex_tmp);
#endif
p_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 DEBUG
printf("layyer%d output\n",layyer);
show_signal(p_out_put_signal);
#endif
for(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;
}
版权声明：本文博客原创文章，博客，未经同意，不得转载。
