DIT-FFT C实现

//代码部分：
//
//函数功能：序列  A(n)的FFT
//入口参数：结构体数组
//			N点傅里叶变换
//返回值：  无
/
void FFT_DEAL(complex *A, int N_point)
{
	complex T;
	FFT_deal FFT_D = {0};
	complex WN_p = {0};
	point_T.inver(A, N_point);
	FFT_D.M =(int) (log(N_point) / log(2));
	for (FFT_D.L=1;  FFT_D.L<=FFT_D.M;  FFT_D.L++)
	{
		FFT_D.B = (int)pow( 2, FFT_D.L-1 );
		for (FFT_D.J=0;  FFT_D.J <= FFT_D.B-1;  FFT_D.J++)
		{
			FFT_D.P = ( (int)pow(2,FFT_D.M- FFT_D.L) )*FFT_D.J;
			
			WN_p.Real = (double)    cos( (double)3.1415926 * 2 * FFT_D.P / N_point);
			WN_p.Img  = (double)(-1*sin( (double)3.1415926 * 2 * FFT_D.P / N_point));
			
			for (FFT_D.K = FFT_D.J; FFT_D.K <= (N_point - 1); FFT_D.K += (int)pow(2, FFT_D.L))
			{
				T = point_T.com_ADD( &A[FFT_D.K], &point_T.com_MUL(&A[FFT_D.K+ FFT_D.B], &WN_p) );
				A[FFT_D.K + FFT_D.B] = point_T.com_del( &A[FFT_D.K], &point_T.com_MUL( &A[FFT_D.K + FFT_D.B], &WN_p ) );
				A[FFT_D.K] = T;
			}
		}
	}
}


//
//函数功能：序列  A(n)的FFT
//入口参数：double 数组
//			N点傅里叶变换
//			mode: 
//				1:相频特性
//				0:幅频特性
//结果：
//			原数组存放FFT之后的数据
/
void FFT_angle_range(double* Xn, int N_point, char mode)
{
	int i = -1;
	double* exechange;
	complex* L = (complex*)malloc(sizeof(complex) * N_point);
	if (NULL==L)exit(1);
	while (i++ < (N_point - 1))
	{
		L[i].Real = (double)Xn[i];
		L[i].Img = 0;
	}
	point_T.FF_D(L, N_point);
	if (mode) exechange = point_T.get_agle(L, N_point);
	else     exechange = point_T.get_rage(L, N_point);
	i = -1;
	while (i++ < (N_point - 1)) Xn[i] = exechange[i];
	free(exechange);
	free(L);

}
///
///
//函数功能：复数乘
///

complex complex_MUL(complex* A, complex* B)
{
	complex L;
	L.Real = A->Real * B->Real - A->Img * B->Img;
	L.Img  = A->Real * B->Img  + A->Img * B->Real;
	return L;
}
///
//函数功能：复数加
//
///

complex complex_ADD(complex* A, complex* B)
{
	complex L;
	L.Real = A->Real+ B->Real;
	L.Img =  A->Img + B->Img;
	return L;
}

//
//函数功能：复数减
//
//
complex complex_del(complex* A, complex* B)
{
	complex L;
	L.Real = A->Real - B->Real;
	L.Img  = A->Img  - B->Img ;
	return L;
}

//
//函数功能：数据的倒序
//入口参数：结构体数组；N点傅里叶变换
//返回值：  无
//

void inverted_order(complex *A, int N_point)
{
	int LH = N_point / 2;
	int J = LH;
	int N1 = N_point - 2;
	int I = 0, K = 0;
	complex T;
	for (I = 1; I  <= N1; I++)
	{
		if (I < J)
		{
			T = A[I];
			A[I] = A[J];
			A[J] = T;
		}
		K = LH;
		while (1)
		{
			if (J < K)break;
			J = J - K;
			K = K / 2;
		}
		J = J + K;
	}

}
//
//函数功能：得到FFT之后的  X(K)的幅频特性
//入口参数：结构体数组 ：x(K)序列
//			N点傅里叶变换
//返回值：  X(K)的幅频特性
//          结构体数组
//
double* get_range(complex* p, int length)
{
	int i = -1;
	double* range = (double*)malloc(sizeof(double) * length);
	if (NULL==range)exit(1);
	while (i++ < length-1)
		range[i] = sqrt(p[i].Real * p[i].Real + p[i].Img * p[i].Img);
	return range;
}

//
//函数功能：得到FFT之后的  X(K)的相频特性
//入口参数：结构体数组 ：x(K)序列
//			N点傅里叶变换
//返回值：  X(K)的幅频特性
//          结构体数组
/
double* get_angle(complex* p, int length)
{
	int i = -1;
	double* angle = (double*)malloc(sizeof(double) * length);
	if (NULL==angle)exit(1);
	while (i++ < length-1)
		angle[i] = atan(p[i].Real * p[i].Real + p[i].Img * p[i].Img);
	return angle;
}
//"fft.h"
#ifndef fft
#define fft
typedef struct {
	double Real;
	double Img;
}complex;

typedef struct {
	int L;
	int M;
	int B;
	int J;
	int P;
	int K;

}FFT_deal;

typedef struct {
	void (*inver) (complex* A, int N_point);
	void (*FF_D)  (complex* A, int N_point);
	complex (*com_MUL)(complex* A, complex* B) ;
	complex (*com_ADD)(complex* A, complex* B) ;
	complex (*com_del)(complex* A, complex* B) ;
	double* (*get_rage) (complex* p, int length);
	double* (*get_agle) (complex* p, int length);
	void (*FFT_agle_rage)(double* Xn, int N_point, char mode);
}point_to;


void inverted_order(complex* A, int N_point);
void FFT_DEAL(complex* A, int N_point);
complex complex_MUL(complex* A, complex* B);
complex complex_ADD(complex* A, complex* B);
complex complex_del(complex* A, complex* B);
double* get_range(complex* p, int length);
double* get_angle(complex* p, int length);
void FFT_angle_range(double* Xn, int N_point, char mode);
#endif

//fft.c
#include <stdio.h>
#include "fft.h"
#include <math.h>
#include <malloc.h>
#include <stdlib.h>
//
point_to point_T = {        //
inverted_order,				//
FFT_DEAL,					//
complex_MUL,				//
complex_ADD,				//
complex_del,				//
get_range,					//
get_angle,					//
FFT_angle_range				//
};							//
//


///
例：
矩形序列R4(n)的256点DFT
//
#include <stdio.h>
#include "fft.h"
extern point_to point_T;
double a[256] = {1,1,1,1};

int main()
{
	int i = 0;
	point_T.FFT_agle_rage(a, 256,0);
	while (i < 256 )printf("%f\r\n",a[i++]);
	return 0;
}

结果：
4.000000
3.998494
3.993979
3.986459
3.975944
3.962447
3.945985
3.926578
3.904250
3.879029
3.850945
3.820033
3.786332
3.749882
3.710729
3.668921
3.624510
3.577549
3.528097
3.476214
3.421965
3.365415
3.306633
3.245693
3.182667
3.117634
3.050673
2.981865
2.911294
2.839046
2.765208
2.689871
2.613126
2.535065
2.455783
2.375375
2.293939
2.211572
2.128372
2.044441
1.959877
1.874782
1.789257
1.703403
1.617322
1.531116
1.444886
1.358733
1.272759
1.187062
1.101744
1.016903
0.932635
0.849039
0.766210
0.684242
0.603228
0.523258
0.444423
0.366809
0.290504
0.215590
0.142149
0.070260
0.000000
0.068556
0.135337
0.200272
0.263297
0.324347
0.383361
0.440283
0.495056
0.547630
0.597956
0.645990
0.691690
0.735016
0.775936
0.814416
0.850430
0.883953
0.914964
0.943446
0.969385
0.992773
1.013601
1.031869
1.047576
1.060728
1.071333
1.079402
1.084952
1.088002
1.088574
1.086694
1.082392
1.075702
1.066659
1.055303
1.041678
1.025829
1.007806
0.987662
0.965452
0.941233
0.915069
0.887021
0.857158
0.825547
0.792261
0.757373
0.720960
0.683100
0.643873
0.603361
0.561649
0.518823
0.474969
0.430176
0.384535
0.338137
0.291073
0.243438
0.195326
0.146830
0.098047
0.049071
0.000000
0.049071
0.098047
0.146830
0.195326
0.243438
0.291073
0.338137
0.384535
0.430176
0.474969
0.518823
0.561649
0.603361
0.643873
0.683100
0.720960
0.757373
0.792261
0.825547
0.857158
0.887021
0.915069
0.941233
0.965452
0.987662
1.007806
1.025829
1.041678
1.055303
1.066659
1.075702
1.082392
1.086694
1.088574
1.088002
1.084952
1.079402
1.071333
1.060728
1.047576
1.031869
1.013601
0.992773
0.969385
0.943446
0.914964
0.883953
0.850430
0.814416
0.775936
0.735016
0.691690
0.645990
0.597956
0.547630
0.495056
0.440283
0.383362
0.324347
0.263297
0.200273
0.135337
0.068556
0.000000
0.070260
0.142148
0.215590
0.290504
0.366809
0.444422
0.523258
0.603227
0.684242
0.766210
0.849039
0.932635
1.016902
1.101744
1.187062
1.272759
1.358733
1.444886
1.531116
1.617322
1.703403
1.789256
1.874782
1.959877
2.044441
2.128372
2.211572
2.293939
2.375375
2.455783
2.535065
2.613126
2.689871
2.765208
2.839046
2.911294
2.981865
3.050673
3.117634
3.182667
3.245693
3.306633
3.365415
3.421965
3.476214
3.528097
3.577549
3.624510
3.668921
3.710729
3.749882
3.786332
3.820033
3.850945
3.879029
3.904250
3.926578
3.945985
3.962447
3.975944
3.986459
3.993979
3.998494