离散傅里叶变换（DFT）的推导及C语言实现

1、傅里叶变换（FT）

傅里叶变换（连续时间傅里叶变换）是该部分内容的理论基础，回顾一下：

傅里叶变换：

        

傅里叶逆变换：

        

以上是连续时间傅里叶变换，但计算机只能处理离散的数据。因此有了离散傅里叶变换（DFT），下面进行详细推导。

---

2、离散傅里叶变换（DFT）

2.1 采样和离散时间傅里叶变换（DTFT）

使用采样可将连续域上的数据离散化。信号学科里面使用冲激序列串来对连续信号采样。

有信号，现对其采样。若采样频率为 ，则采样间隔  ，用以采样的冲激序列串定义为：

        

因此采样后的信号为：

        

此时将采样后信号  代入傅里叶变换公式：

        

交换积分与求和次序：

        

由  函数的筛选性 ，将上式中  看成， 得到离散时间傅里叶变换（DTFT）：

        

因为  表示采样的时间点，所以可令 ，因此离散时间傅里叶变换的更一般形式为：

        

---

2.2 离散傅里叶变换（DFT）

2.1中通过采样使得信号在时域离散化，但并不能保证其频域也离散化，同样不利于计算机处理。

由性质：时域离散，频域周期化；频域离散，时域周期化。因此若想让信号在频域也离散，则需要该信号在时域上为周期信号。

但原信号不一定为周期信号。解决方式是周期延拓：截取原无限长信号的个采样点，假设：

① 个采样点为原信号的一个周期；

  

② N个采样点外为该N点的周期延拓。

 

这样原先的离散非周期信号就变成离散周期信号，因此频域得以离散化。

在上述周期延拓的基础上，假设采样间隔为  ，则个采样点的采样周期 ，从连续信号中截取的个采样点的信号可表示为：

        

因为周期延拓，为是周期信号，周期函数的傅里叶级数为：​ 

         

将代入其中，得：

        

交换积分与求和次序：

        

同理，由  函数的筛选性，上式变为：

        

因为，，代入上式：

        

令，，得离散傅里叶变换（DFT）的表达式为：

        ，   且 

离散傅里叶逆变换（IDFT）的表达式为：

        

注意：为了满足正逆变换的自洽，放在正逆变换其中之一前就可以，通常放在逆变换前。

补充：令，则：

        

---

2.3 离散傅里叶变换（DFT）的C语言实现

根据上述推导出的公式，很容易编码实现。如下：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PI 3.141593

double realComput(double xn[], int ndft, int k);
double imageComput(double xn[], int ndft, int k);

typedef struct {
	double real;
	double image;
} Complex;

Complex* dft(double x[], int ndft) {
	Complex* dftRes = (Complex*)malloc(ndft * sizeof(Complex));
	if (dftRes == NULL) {
		return;
	}

	for (int i = 0; i < ndft; ++i) {
		dftRes[i].real = realComput(x, ndft, i);
		dftRes[i].image = imageComput(x, ndft, i);
		// printf("%lf + %lfi\n", dftRes[i].real, dftRes[i].image);
	}
	return dftRes;
}

double realComput(double xn[], int ndft, int k) {
	double realPart = 0;
	for (int i = 0; i < ndft; ++i) {
		realPart += xn[i] * cos(2 * PI / ndft * k * i);
	}
	return realPart;
}

double imageComput(double xn[], int ndft, int k) {
	double imagePart = 0;
	for (int i = 0; i < ndft; ++i) {
		imagePart -= xn[i] * sin(2 * PI / ndft * k * i);
	}
	return imagePart;
}

double* ampSpectrum(Complex* dftRes, int ndft) {
	double* amp = (double*)malloc(sizeof(double) * ndft);
	if (amp == NULL) {
		return;
	}
	for (int i = 0; i < ndft; ++i) {
		amp[i] = sqrt(dftRes[i].real * dftRes[i].real+ dftRes[i].image* dftRes[i].image);
	}
	return amp;
}

double* phaseSpectrum(Complex* dftRes, int ndft) {
	double* phase = (double*)malloc(sizeof(double) * ndft);
	if (phase == NULL) {
		return;
	}
	for (int i = 0; i < ndft; ++i) {
		phase[i] = atan2(dftRes[i].image, dftRes[i].real);
	}
	return phase;
}

void testDFT() {
	double xn[] = { 1, 2, 3, 4, 6, 41, 0, 855, 954, -1 };
	int ndft = sizeof(xn) / sizeof(double);
	Complex* dftRes = dft(xn, ndft);

	double* ampRes = ampSpectrum(dftRes, ndft);
	double* phaRes = phaseSpectrum(dftRes, ndft);

	for (int i = 0; i < ndft; ++i) {
		printf("%lf + %lfi\n", dftRes[i].real, dftRes[i].image);
	}
	printf("\n");

	for (int i = 0; i < ndft; ++i) {
		printf("%lf, ", ampRes[i]);
	}
	printf("\n");

	for (int i = 0; i < ndft; ++i) {
		printf("%lf, ", phaRes[i]);
	}
	printf("\n");

	free(dftRes);
	free(ampRes);
	free(phaRes);
}

int main() {
	testDFT();
	return 0;
}

运行结果：

暂未发现问题。