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

1、傅里叶变换（FT）

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

傅里叶变换：

$F(\omega ) = \int_{-\infty }^{+\infty}f(t)e^{-j\omega t}dt$

傅里叶逆变换：

$f(t ) = \frac{1}{2\pi }\int_{-\infty }^{+\infty}F(\omega )e^{j\omega t}d\omega$

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

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

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

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

有信号 $f(t)$ ，现对其采样。若采样频率为 $f_{s}$ ，则采样间隔 $T_{s} = \frac{1}{f_{s}}$ ，用以采样的冲激序列串定义为：

$\delta_{s}(t) = \sum_{n=-\infty}^{+\infty}\delta(t-nT_{s})$

因此采样后的信号为：

$f_{s}(t)=\sum_{n=-\infty}^{+\infty}f(t)\delta(t-nT_{s})$

此时将采样后信号 $f_{s}(t)$ 代入傅里叶变换公式：

$F_s(\omega ) = \int_{-\infty }^{+\infty}\left [ \sum_{n=-\infty}^{+\infty}f(t)\delta(t-nT_{s}) \right ]e^{-j\omega t}dt$

交换积分与求和次序：

$F_s(\omega ) = \sum_{n=-\infty}^{+\infty}\int_{-\infty }^{+\infty}f(t)\delta(t-nT_{s})e^{-j\omega t}dt$

由 $\delta (t)$ 函数的筛选性 $\int_{-\infty}^{+\infty}x(t)\delta (t-t_0{})dt=x(t_0)$ ，将上式中 $f(t)e^{-j\omega t}$ 看成 $x(t)$ ，得到离散时间傅里叶变换（DTFT）：

$F_s(\omega ) = \sum_{n=-\infty}^{+\infty}f(nT_{s})e^{-j\omega nT_s}$

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

$F_s(\omega ) = \sum_{k=-\infty}^{+\infty}f(k)e^{-j\omega k}$

2.2 离散傅里叶变换（DFT）

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

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

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

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

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

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

在上述周期延拓的基础上，假设采样间隔为 $T_{s}$ ，则 $N$ 个采样点的采样周期 $T_{0}=NT_{s}$ ，从连续信号 $f(t)$ 中截取的 $N$ 个采样点的信号可表示为：

$x(t) = \sum_{n=0}^{N-1}f(t)\delta (t-nT_s)$

因为周期延拓， $x(t)$ 为是周期信号，周期函数的傅里叶级数为：

$F(k\omega)=\frac{1}{T}\int_{0}^{T}f(t)e^{-jk\omega t}dt$

将 $x(t)$ 代入其中，得：

$X(k\omega)=\frac{1}{T_0}\int_{0}^{T_0}\left [ \sum_{n=0}^{N-1}f(t)\delta (t-nT_s) \right ]e^{-jk\omega_{}t}dt$

交换积分与求和次序：

$X(k\omega)=\frac{1}{T_0}\sum_{n=0}^{N-1}\int_{0}^{T_0} f(t)\delta (t-nT_s) e^{-jk\omega_{}t}dt$

同理，由 $\delta (t)$ 函数的筛选性，上式变为：

$X(k\omega)=\frac{1}{T_0}\sum_{n=0}^{N-1}f(nT_s)e^{-jk\omega nT_s}$

因为 $T_{0}=NT_{s}$ ， $\omega =\frac{2\pi }{T0}=\frac{2\pi }{NT_{s}}$ ，代入上式：

$X(k\omega)=\frac{1}{NT_{s}}\sum_{n=0}^{N-1}f(nT_s)e^{-j\frac{2\pi }{N} kn}$

令 $x[n]=f(nT_s)$ ， $X(k\omega)T_s=X[k]$ ，得离散傅里叶变换（DFT）的表达式为：

$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi }{N} kn}$ ， $0\leq k< N$ 且 $k\in Z$

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

$x[n]=\frac{1}{N}\sum_{n=0}^{N-1}X[k]e^{j\frac{2\pi }{N} kn}$

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

补充：令 $W_{N}^{kn}=e^{-j\frac{2\pi }{N} kn}$ ，则：

$X[k]=\sum_{n=0}^{N-1}x[n]W_{N}^{kn}$

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

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

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

#define PI 3.1415927

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

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

void 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);
	}

	free(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;
}

int main() {
	double xn[9] = {1,2,3,4,5,0,0,0,85};
	dft(xn, sizeof(xn) / sizeof(double));
	system("pause");
	return 0;
}

运行结果：