FFT 的C 语言

时间：2023-03-09 19:28:40
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.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);

#endif

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

}

reverse_bits.c

/***************************************************************

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

}

show_signal.c

/***************************************************************

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;

}

版权声明：本文博客原创文章，博客，未经同意，不得转载。

相关文章