DIT-FFT[C语言实现]

复数运算

首先，FFT运算涉及复数运算，如果不想用C语言的复数库，可以选择自定义复数结构。先定义如下结构体以及相关运算：

typedef struct{
	float real;
	float imag;
}complex;

complex complex_add(complex a, complex b)
{
	complex out;
	out.real = a.real + b.real;
	out.imag = a.imag + b.imag;
	return out;
}

complex complex_sub(complex a, complex b)
{
	complex out;
	out.real = a.real - b.real;
	out.imag = a.imag - b.imag;
	return out;
}

complex complex_multiply(complex a, complex b)
{
	complex out;
	out.real = a.real * b.real - a.imag * b.imag;
	out.imag = a.real * b.imag + a.imag * b.real;
	return out;
}

complex rotate_factor(int id)
{
	complex out;
	out.real = (float) cos(-2*PI*id/N);
	out.imag = (float) sin(-2*PI*id/N);
	return out;
}

DIT-FFT基本原理

下面讲一讲DIT-FFT的基本原理：

旋转因子的周期性

旋转因子，具有周期性、对称性、可约性。

周期性：$W_N^{m+lN}=W_N^m$,

对称性：$[W_N^{N-m}{]}^{\ast }=W_n^m$,

可约性：$W_N^{mk}=W_{\frac{N}{k}}^m$.

DIT-FFT四步骤

DIT-FFT，即Decimation In Time - Fast Fourier Transform, 时域基2FFT算法，主要可分为4个步骤。

step1：选定长度

$$N=2^M$$

step2：奇偶分解成子序列

$$\begin{gathered} x_1(r)=x(2r),r=0\le r\le N/2-1 \\ {x}_2(r)=x(2r+1),r=0\le r\le N/2-1 \end{gathered}$$

step3：利用可约性进行转化

$$\begin{gathered} X(k)=\sum _{n_{even}}x(n)W_N^{kn}+\sum _{n_{odd}}x(n)W_N^{kn} \\ =\sum _{r=0}^{N/2-1}x(2r)W_N^{2kr}+W_N^k\sum _{r=0}^{N/2-1}x(2r+1)W_N^{2kr} \\ =\sum _{r=0}^{N/2-1}{x}_1(r)W_{N/2}^{kr}+W_N^k\sum _{r=0}^{N/2-1}{x}_2(r)W_{N/2}^{kr} \end{gathered}$$

step4：利用周期性和对称性将X(k)分段表示

令
$$\begin{gathered} X_1(k)=\sum _{r=0}^{N/2-1}{x}_1(r)W_{N/2}^{kr} \\ {X}_2(k)=\sum _{r=0}^{N/2-1}{x}_2(r)W_{N/2}^{kr} \\ 0\le k\le N-1 \end{gathered}$$
则
$$\begin{gathered} X(k)=\sum _{r=0}^{N/2-1}{x}_1(r)W_{N/2}^{kr}+W_N^k\sum _{r=0}^{N/2-1}{x}_2(r)W_{N/2}^{kr} \\ =X_1(k)+W_N^k{X}_2(k),0\le k\le N-1 \end{gathered}$$
由于
$$\begin{gathered} X_1(k+N/2)=X_1(k),0\le k\le N/2-1 \\ {X}_2(k+N/2)=X_2(k),0\le k\le N/2-1 \\ {W}_N^{k+N/2}=-W_N^k \end{gathered}$$
可得
$$\begin{gathered} X(k)=X_1(k)+W_N^k{X}_2(k),0\le k\le N/2-1 \\ X(k+N/2)=X_1(k)-W_N^k{X}_2(k),0\le k\le N/2-1 \end{gathered}$$

C语言实现DIT-FFT

在编写算法之前，现在文件开头，引入相关的库，并做一些声明：

#include "stdio.h"
#include "math.h"
#define N 8
#define M (int)log2(N)
#define PI 3.1415926

为保证输出自然序，需要对输入进行位反序。此处采用变换索引的方式进行：

int vector[N];
	int i, j, k;
	for(i=0; i<N;i++)
	{
		vector[i] = i;
	}
	int binlist[N][M];
	int (*p)[M] = binlist;
	int transvector[N];
	for(p=binlist,i=0; p<binlist + N; p++,i++)
	{
		int Devidend = vector[i];
		int Devisor, Remainder;
		for(j=0; j<M; j++)
		{
			Devisor = Devidend / 2;
			Remainder = Devidend % 2;
			*(*p+M-1-j) = Remainder;
			Devidend = Devisor;
		}
		int tempnum;
		for(k=0; k<=M/2-1; k++)
		{
			tempnum = *(*p+k);
			*(*p+k) = *(*p+M-1-k);
			*(*p+M-1-k) = tempnum;
		}
		transvector[i] = 0;
		for(k=0; k<M; k++)
		{
			transvector[i] += *(*p+k) * (int)pow(2, M-1-k);
		}	
	}

主要思路为：生成0到N-1索引、转换为二进制、进行位反转、转换为十进制。

下面提供一组测试数据：

	complex x[N];
	x[0].real = 1;	x[0].imag = 0;
	x[1].real = 0;	x[1].imag = 0;
	x[2].real = 4;	x[2].imag = 0;
	x[3].real = 2;	x[3].imag = 0;
	x[4].real = 5;	x[4].imag = 0;
	x[5].real = 3;	x[5].imag = 0;
	x[6].real = 2;	x[6].imag = 0;
	x[7].real = 6;	x[7].imag = 0;

当然，也可以通过交互方式输入：

	for(i=0; i<N; i++)
	{
		printf("请输入第%d个数的实部:",i+1);
		scanf("%f", &x[i].real);
		printf("请输入第%d个数的虚部:",i+1);
		scanf("%f", &x[i].imag);
	}

于是，可以获得重排序列如下：

	complex X[N];
	for(i=0; i<N; i++)
	{
		X[vector[i]].real = x[transvector[i]].real;
		X[vector[i]].imag = x[transvector[i]].imag;
	}

而后便是算法的核心的部分了，如下：

	int L, B, level_id, J, P, base_id, K;
	complex temp_num;
	for(L=1; L<=M; L++)
	{
		B = (int) pow(2, L-1);
		level_id = 0;
		for(J=0; J<B;J++)
		{
			P = J * (int) pow(2, M-L);
			base_id = 0;
			for(K=1; K<=( (N/2) / ( (int) pow(2, L-1) )); K++)
			{
				temp_num = X[base_id+level_id];
				X[base_id+level_id] = complex_add(X[base_id+level_id], complex_multiply(X[base_id + B+level_id], rotate_factor(P)));
				X[base_id + B+level_id] = complex_sub(temp_num, complex_multiply(X[base_id + B+level_id], rotate_factor(P)));
				base_id += (int) pow(2, L);
			}
			level_id += 1;
		}
	} 

这里根据FFT变换的性质，采用三组for循环实现。第一层循环为级数循环，第二层循环根据单级内不同类型的蝶形算子进行计算，第三层循环根据每一级相同类型的旋转因子的个数进行计算。

结果显示部分如下：

	printf("The result is:\n");
	for(i=0; i<N;i++)
	{
		printf("X[%d] = %8.4f  +  j %8.4f\n", i, X[i].real, X[i].imag);
		
	}

运算结果为：

The result is:
X[0] =  23.0000  +  j   0.0000
X[1] =  -3.2929  +  j   2.9497
X[2] =  -0.0000  +  j   5.0000
X[3] =  -4.7071  +  j   6.9497
X[4] =   1.0000  +  j   0.0000
X[5] =  -4.7071  +  j  -6.9497
X[6] =   0.0000  +  j  -5.0000
X[7] =  -3.2929  +  j  -2.9497

在matlab中进行验证：

>> fft([1,0,4,2,5,3,2,6]).'
ans =

  23.0000 + 0.0000i
  -3.2929 + 2.9497i
   0.0000 + 5.0000i
  -4.7071 + 6.9497i
   1.0000 + 0.0000i
  -4.7071 - 6.9497i
   0.0000 - 5.0000i
  -3.2929 - 2.9497i

结果相吻合。