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		:	[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 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);
#endif

complex_add.c

/***************************************************************
code writer	:	EOF
code file	:	complex_add.c
code date	:	2014.09.17
e-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 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		:	[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 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		:	[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	:	EOF
code file	:	reverse_bits.c
code date	:	2014.09.17
e-mail		:	[email protected]

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		:	[email protected]

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

show_signal.c

/***************************************************************
code writer	:	EOF
code file	:	show_signal.c
code date	:	2014.09.17
e-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 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		:	[email protected]

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