快速傅里叶变换
一些使离散傅里叶变换更快的想法，并实现了一个精简的DFT算法版本。
带有复数的基本DFT算法
根据傅里叶变换的主要原则，一个在N内周期性的信号f(n)，并且我们有它的N个采样，可以像下面这样转换成它的傅里叶分量
用
复数傅里叶分量 ck 的另一个等式看起来像
和
这意味着复数傅里叶分量 ck 是 f(n) 和围绕单位圆旋转的复数向量的复数和。对于一个完整的傅里叶变换，我们需要 N 个傅里叶分量 (ck)，并且对于每个分量，我们需要 N 个正弦和 N 个余弦值。对于 1000 个样本，这将需要计算 100 万个正弦和余弦值。这消耗了大量的 CPU 能力和时间。但它确实有效 J
public void CalcDFT(int N)
{
int k, n;
TKomplex w;
if (N > 0)
{
for (k = 0; k < N; k++)
{
c[k].real = 0;
c[k].imag = 0;
for (n = 0; n < N; n++)
{
w.real = Math.Cos((double)(2.0 * Math.PI * (double)(k * n) / (double)(N)));
w.imag = -Math.Sin((double)(2.0 * Math.PI * (double)(k * n) / (double)(N)));
c[k] = ksum(c[k], kprod(w, y[n]));
}
c[k].real = c[k].real / (double)(N) * 2.0;
c[k].imag = -c[k].imag / (double)(N) * 2.0;
}
}
c[0].real = c[0].real / 2;
c[0].imag = c[0].imag / 2;
}
通过这个标准实现计算，1000 个样本的傅里叶变换在我的电脑上花费了 140 毫秒（当然世界上有更快的电脑 J）。
可以在 Visual C# 项目 DFT_std 中找到此标准实现的源代码。
ksum() 和 kprod() 是用于复数加法和乘法的简单函数。
public TKomplex ksum(TKomplex a, TKomplex b)
{
TKomplex res;
res.real = a.real + b.real;
res.imag = a.imag + b.imag;
return (res);
}
public TKomplex kprod(TKomplex a, TKomplex b)
{
TKomplex res;
res.real = a.real * b.real - a.imag * b.imag;
res.imag = a.real * b.imag + a.imag * b.real;
return (res);
}
关于周期性的思考
现在，正弦和余弦函数都是 2π 的周期函数，对于 k = 1，我们必须计算正弦和余弦的值，从 0 到 2π，步长为 2π /N。 对于 k = 2，我们有 0 到 4 π，步长为 2/N，依此类推。 所以对于 k = 2，我们得到的值与 k = 1 时已经得到的值相同。 只是，我们得到每个第二个值，当 n 大于 N/2 时，我们完成了一圈单位圆并再次从与 n = 0 相同的值开始。 k 越大，步长越大，我们经过的单位圆圈数越多。 但我们得到的正弦和余弦值仍然与 k = 1 时得到的相同。 这意味着基本上只需要计算 0 到 2π 之间的 N 个正弦和 N 个余弦值，步长为 2π/N，并将它们保存在查找表中以供整个变换使用。 对于 1000 个样本，我们只需要计算 1000 个正弦值和 1000 个余弦值，而不是 999000 个，这节省了相当多的时间。
我们使用以下方法构建查找表：
for (k = 1; k < N; k++)
{
we[k].real = Math.Cos((2.0 * Math.PI * (double)(k) / (double)(N)));
we[k].imag = -Math.Sin((2.0 * Math.PI * (double)(k) / (double)(N)));
}
我们可以通过 k * n 除以 N 的简单取模运算从查找表中获得相应的正弦和余弦值。 这样，离散傅里叶变换变得不是“快”，而是相当“快”和精简 J
public void CalcDFT() // Quick and lean Fourier transformation
{
int k, n;
if (N > 0)
{
for (k = 0; k < N; k++)
{
c[k].real = 0;
c[k].imag = 0;
for (n = 0; n < (N - 1); n++)
{
c[k] = ksum(c[k], kprod(we[(k * n) % N], y[n]));
}
c[k].real = c[k].real / N * 2;
c[k].imag = -c[k].imag / N * 2;
}
c[0].real = c[0].real / 2;
c[0].imag = c[0].imag / 2;
}
}
对于与上述相同的 1000 个样本，此实现花费 53 毫秒。 快了很多。
此实现的 code 在 DFT_lean 中。
复数？
现在，我的脑海中又出现了一些想法。 在主要条款中没有复数。 输入是实数序列（至少在我的应用程序中来自电信号分析 J），我们最终得到另一个实数序列。 那么，为什么我们在中间使用复数呢？
在我们的算法中，我们有一个复数乘法和一个复数加法。 所以让我们看看它们。 首先是乘法，we[x] * y[n]。 看起来像
(we[x].re * y[n].re - we[x].im * y[n].im) + j(we[x].re * y[n].im + we[x].im * y[n].re)
y[n].im 总是 0，因此第二部分和第三部分也是 0，我们的复数乘法简化为
(we[x].re * y[n].re) + j(we[x].im * y[n].re)
这意味着
y[n].re * cosine(x) + j(y[n].re * sine(x))
这看起来很像主要条款中的 ak 和 bk。 在下面的加法中，我们只需要对实部和虚部求和。 因此，在整个计算中，实部和虚部没有混合，最后我们在 ck 的实部得到 ak 部分，在 ck 的虚部得到 bk 部分，我们可以忘记所有复杂的東西 J
有了输入序列不包含虚数的限制，我们可以实现没有复数的傅里叶变换
public void CalcDFT()
{
int k, n;
if (N > 0)
{
for (k = 0; k < N; k++)
{
a[k] = 0;
b[k] = 0;
for (n = 0; n < (N - 1); n++)
{
a[k] = a[k] + ((cosine[(k * n) % N] * y[n]));
b[k] = b[k] + ((sine[(k * n) % N] * y[n]));
}
a[k] = a[k] / N * 2;
b[k] = b[k] / N * 2;
}
a[0] = a[0] / 2;
b[0] = b[0] / 2;
}
}
初始化变成这样
public TFftAlgorithm(int order)
{
int k;
N = order;
y = new double[N + 1];
a = new double[N + 1];
b = new double[N + 1];
xw = new double[N + 1];
sine = new double[N + 1];
cosine = new double[N + 1];
cosine[0] = 1.0; // we don't have to calculate cos(0) = 1
sine[0] = 0; // and sin(0) = 0
for (k = 1; k < N; k++) // init vectors of unit circle
{
cosine[k] = Math.Cos((2.0 * Math.PI * (double)(k) / (double)(N)));
sine[k] = Math.Sin((2.0 * Math.PI * (double)(k) / (double)(N)));
}
}
通过这种实现，1000 个样本在 20 毫秒内被转换。 这几乎比本文顶部的第一种实现快 5 倍。 它不是快速傅里叶变换。 但至少这是一个快速傅里叶变换，它看起来非常精简 J
这个实现位于项目 DFT_lean_quick 中。
对于矩形信号，结果应该看起来像这样
蓝色形状是输入序列：一个有 1000 个样本的矩形信号。幅度为 20。红色形状是用前 30 个傅里叶分量进行的反傅里叶变换。
欲了解更多数学知识，请访问 www.mosismath.com :-)