DIT-FFT[C语言实现]

复数运算
DIT-FFT基本原理
- 旋转因子的周期性
- DIT-FFT四步骤
- - step1：选定长度
  - step2：奇偶分解成子序列
  - step3：利用可约性进行转化
  - step4：利用周期性和对称性将X(k)分段表示
C语言实现DIT-FFT

复数运算

首先，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的基本原理：

旋转因子的周期性

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

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

对称性：[WNN−m]∗=Wnm[W_N^{N-m}]^*=W_n^m[WNN−m]∗=Wnm,

可约性：WNmk=WNkmW_N^{mk}=W_{\frac{N}{k}}^mWNmk=WkNm.

DIT-FFT四步骤

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

step1：选定长度

N=2MN = 2^M N=2M

step2：奇偶分解成子序列

x1(r)=x(2r),r=0≤r≤N/2−1x2(r)=x(2r+1),r=0≤r≤N/2−1x_1(r)=x(2r)\quad ,r =0\le r \le N/2 -1\\ x_2(r)=x(2r+1)\quad ,r =0\le r \le N/2 -1 x1(r)=x(2r),r=0≤r≤N/2−1x2(r)=x(2r+1),r=0≤r≤N/2−1

step3：利用可约性进行转化

X(k)=∑nevenx(n)WNkn+∑noddx(n)WNkn=∑r=0N/2−1x(2r)WN2kr+WNk∑r=0N/2−1x(2r+1)WN2kr=∑r=0N/2−1x1(r)WN/2kr+WNk∑r=0N/2−1x2(r)WN/2krX(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} X(k)=neven∑x(n)WNkn+nodd∑x(n)WNkn=r=0∑N/2−1x(2r)WN2kr+WNkr=0∑N/2−1x(2r+1)WN2kr=r=0∑N/2−1x1(r)WN/2kr+WNkr=0∑N/2−1x2(r)WN/2kr

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

令
X1(k)=∑r=0N/2−1x1(r)WN/2krX2(k)=∑r=0N/2−1x2(r)WN/2kr0≤k≤N−1X_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 X1(k)=r=0∑N/2−1x1(r)WN/2krX2(k)=r=0∑N/2−1x2(r)WN/2kr0≤k≤N−1
则
X(k)=∑r=0N/2−1x1(r)WN/2kr+WNk∑r=0N/2−1x2(r)WN/2kr=X1(k)+WNkX2(k),0≤k≤N−1X(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) \quad ,0\le k \le N-1 X(k)=r=0∑N/2−1x1(r)WN/2kr+WNkr=0∑N/2−1x2(r)WN/2kr=X1(k)+WNkX2(k),0≤k≤N−1
由于
X1(k+N/2)=X1(k),0≤k≤N/2−1X2(k+N/2)=X2(k),0≤k≤N/2−1WNk+N/2=−WNkX_1(k+N/2)=X_1(k)\quad, 0\le k \le N/2-1\\ X_2(k+N/2)=X_2(k)\quad, 0\le k \le N/2-1\\ W_N^{k+N/2}=-W_N^k X1(k+N/2)=X1(k),0≤k≤N/2−1X2(k+N/2)=X2(k),0≤k≤N/2−1WNk+N/2=−WNk
可得
X(k)=X1(k)+WNkX2(k),0≤k≤N/2−1X(k+N/2)=X1(k)−WNkX2(k),0≤k≤N/2−1X(k)=X_1(k)+W_N^k X_2(k) \quad ,0\le k \le N/2-1\\ X(k+N/2)=X_1(k)-W_N^k X_2(k) \quad ,0\le k \le N/2-1 X(k)=X1(k)+WNkX2(k),0≤k≤N/2−1X(k+N/2)=X1(k)−WNkX2(k),0≤k≤N/2−1

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.9497i0.0000 + 5.0000i-4.7071 + 6.9497i1.0000 + 0.0000i-4.7071 - 6.9497i0.0000 - 5.0000i-3.2929 - 2.9497i

结果相吻合。