FFT算法和DFT算法C语言实现（赋详解）

声明：

本人在校期间主修过《数字信号处理》这门课程
对离散傅里叶变换（DFT）和快速傅里叶变换（FFT）深有了解
现编写了基于C语言的FFT算法，已完成对抽样序列的FFT变换并通过窗口输出。

编写思路：

由于FFT变换里面含有对虚数的运算，现将输出序列定义成一个结构体函数
结构体函数里面包含输出序列的实部和虚部，并定义处理复数运算的几种函数
快速傅里叶变换的核心在于倒序运算和蝶形运算
现将倒序运算和蝶形运算的核心思想用代码的形式表现出来

代码部分

## 倒序代码：

for (int I = 1; I < N1; I++)						//定义倒序序列函数
	{
		if (I < J)
		{
			T = x[I];
			x[I] = x[J];
			x[J] = T;
		}
		K = LH;
		if (J >= K)
		{
			do
			{
				J = J - K;
				K = K / 2;
			} while (J >= K);
		}
		J = J + K;
	}

## 倒序代码解释：

用J表示当前倒序数的十进制数值。
对于N=2^M，M位二进制数最高位的十进制权值为N/2，且从左向右二进制位的权值依次为N/4，N/8，……，2，1。
因此，最高位加一相当于十进制运算J+N/2。如果最高位是0（J<N/2）,则直接由J+N/2得下一个倒序值；
如果最高位是1（J>=N/2），则先将最高位变成0（J=J-N/2），然后次高位加1（J+N/4）。
但次高位加1时，同样需要判断0、1值，如果为0（J<N/4），则直接加1（J=J+N/4），否则将次高位变成0（J=J-N/4），再判断下一位；
依次类推，直到完成最高位加1，逢2向右进位的运算

## FFT代码：

for (L = 1; L <= M; L++)							//FFT运算
	{
		B = (int)(pow(2, L - 1));
		for (J = 0; J < B; J++)
		{
			P = (int)(J*pow(2, M - L));
			for (k = J; k < N; k = (int)(k + pow(2, L)))
			{
				K1 = k + B;
				complex wn, t;
				Wn_i(N, P, &wn);
				c_mul(f[K1], wn, &t);					//。。。。。。。。。。。。
				c_sub(f[k], t, &(f[K1]));				//蝶形运算
				c_plus(f[k], t, &(f[k]));				//。。。。。。。。。。。。
			}
		}
	}

## FFT代码解释：

第L级中，每个蝶形的两个输入数据相距B=2^(L-1)个点
每级有B个不同的旋转因子；同一旋转因子对应着间隔为2^L点的2^(M-1)个蝶形
先从输入端（第1级）开始，逐级进行，共进行M级运算
在进行第L级运算时，依次求出B个不同的旋转因子，每求出一个旋转因子，就计算完它对应的所有2^(M-L)个蝶形

## 复数函数定义与运算代码部分：

void c_plus(complex a, complex b, complex *c)			//复数加法
{
	c->real = a.real + b.real;
	c->imag = a.imag + b.imag;
}
void c_sub(complex a, complex b, complex *c)			//复数减法
{
	c->real = a.real - b.real;
	c->imag = a.imag - b.imag;
}
void c_mul(complex a, complex b, complex *c)			//复数乘法
{
	c->real = a.real*b.real - a.imag*b.imag;
	c->imag = a.real*b.imag + a.imag*b.real;
}

void Wn_i(int n1, int i, complex *Wn)				//定义FFT旋转因子
{
	Wn->real = cos(2 * PI*i / n1);
	Wn->imag = -sin(2 * PI*i / n1);
}

查看完整代码请点击这里 ←

## 完整代码实现的功能：

1.基础功能：
 - 可自定义序列的抽样点数 
 - 对序列进行FFT和DFT运算输出
2.拓展功能：
 - 在定义抽样点的基础上比较DFT运算和FFT运算的时间