Lode 的计算机图形学教程

傅里叶变换

目录

返回索引

简介

傅里叶变换是图像处理中的重要工具，与滤波理论直接相关——滤波器在空间域（即图像本身）中是卷积运算，在频谱域（即图像的 FT）中则只是简单的乘法运算！

大多数关于图像傅里叶变换的教程都使用枯燥的灰度图像。甚至许多可用于计算图像 FT 的程序和插件也只支持灰度图像。在 24 位和 32 位彩色已经日渐普及的今天，我们不能再接受这种局限！本教程以 24 位彩色图像呈现，通过对每个颜色通道分别进行 FT 来实现。以下是一个利用图像 FT 改善可平铺草地纹理效果的示例：



另外，还有一种计算彩色图像 FT 的方法——使用四元数代替复数，可以将一幅 4 通道图像转换为另一幅 4 通道图像。

在开始介绍图像的傅里叶变换（即二维 FT）之前，我们先从更简单的一维 FT 入手，它常用于音频和电磁信号处理。

为了充分理解本教程，了解复数、其模（振幅）和辐角（相位）会很有帮助。如果您希望深入了解这些内容，可以参阅关于复数的附录（其他 html 文件之一），或通过搜索引擎查找相关教程。

以下两个源文件可供下载，包含本文所描述的所有代码：

fourier2d.cpp
fourier1d.cpp

库

本教程包含用于计算傅里叶变换的演示代码，但如果您希望在实际应用中使用傅里叶变换，市面上有更好、更快的库，我强烈建议使用这些库，本教程中的代码仅供学习参考。

特别需要指出的是，本教程后续的快速傅里叶变换代码只支持 2 的幂次大小的输入，而优质的库可以对任意大小执行 FFT。

以下是几个优质的免费库（请仔细查看其许可证，以确认是否可用于您的项目）：

信号与频谱

在开始介绍傅里叶变换之前，有必要先了解信号与频谱的基础知识。就目前而言，信号指的是时域信号，例如音频和电磁信号。如果将这样的信号绘制成关于时间的函数图像，一个简单的信号可能如下所示：



这是一个正弦波，如果它是音频信号，您会听到一个只有 1 个频率的纯音（就像 PC 扬声器的蜂鸣声，或非常纯净的长笛声）。红色曲线背后的灰色曲线表示信号的振幅，始终为正值。

信号的频谱描述了每个频率在信号中所占的比重。由于上述信号是正弦波，只有一个频率，因此频谱非常简单：只有一个值为正（即正弦曲线的频率），其余均为零。所以频谱只有一个单峰。但频谱有两侧：左侧为负频率侧，右侧为正频率侧。负侧包含负频率。对于实数信号（没有虚部），例如音频信号，频谱的负侧始终是正侧的镜像。因此，对于上述正弦信号，正侧有一个单峰，该峰会在负侧被镜像。其频谱如下所示：



白色曲线表示信号的振幅。红色和绿色曲线分别表示频谱的实部和虚部，但频谱很少以这种方式表示。通常，频谱同时给出振幅和相位，但由于现在相位对我们来说没有太多有趣的信息，这里不包含相位。

现在让我们研究另一个例子：同样的正弦波，但加上了直流分量：sin(x)+C。直流分量这个名称来自电子学，代表直流电（与交流电相对）。任何信号中的直流分量就是其时间平均值。正弦信号的时间平均值为 0，但如果在正弦波上加一个常数，信号的时间平均值就会变成该常数。添加这个常数后，正弦信号会整体上移，高于水平轴：



这样的常数（即直流分量）的频率为 0。因此我们可以预期频谱在原点处会出现一个额外的峰值，实际上确实如此：



因此，如果信号的频谱在原点处不为零，就说明其时间平均值不为零，存在直流分量！峰值距原点越远（向左或向右），代表的频率越高，因此远在右侧（及左侧）的峰值意味着信号包含高频分量。需要注意的是，音频信号从不含直流分量，因为没有人能产生或听到它。但电气信号可以有直流分量。

现在让我们看一个由两个正弦函数之和构成的信号：第二个正弦函数的频率是第一个的两倍，因此曲线形状为 sin(x)+sin(2*x)：



由于现在有两个不同频率的正弦波，我们可以预期频谱的正侧有两个峰值（负侧也有两个，因为它是镜像）：



如果两个正弦波的相位不同（即其中一个正弦波水平移动了），频谱的振幅仍然相同，只有相位会不同，但这里没有显示相位。

任意声音信号，例如一个说话的词语，由无限多个正弦函数的无穷级数（或积分）组成，每个正弦函数频率不同，且各有其振幅和相位，而频谱是观察信号中每个频率成分多少的完美工具。这样的频谱不只有几个峰值，而是一个连续函数。

有几种特殊函数具有特殊的频谱：

狄拉克冲激是一个除原点处为无穷大外其他地方均为零的信号。这是连续函数的理想版本；对于计算机上的离散函数，它可以表示为原点处具有有限高度的单个峰值。



这样的峰值在音频信号中听起来像爆裂声，并且包含所有频率。这就是为什么它的频谱看起来是这样的（水平黑线）：



频谱在各处均为正值，因此信号包含所有频率。这也意味着，要将这样的峰值表示为正弦函数之和，需要叠加无穷多个振幅相同、相位各异的基本正弦函数，它们在原点以外相互抵消，从而形成峰值。

注意，可以将频谱解释为时域信号（一个恒定的直流信号），此时红色函数可以看作其频谱（在零点处有一个峰值，因为直流信号的频率为 0）。这种对偶性是傅里叶变换有趣性质之一，我们稍后会看到。

最后，另一个特殊信号是 sinc(x) 函数；sinc(x) = sin(x)/x：



其频谱是一个矩形脉冲（这里只是近似，因为频谱仅用 128 个点计算）：



同样由于信号与其频谱之间的对偶性，矩形时域信号的频谱是 sinc 函数。

频谱也用于研究光的颜色。这完全是同一回事，因为光也是一维信号，包含不同频率，其中一些频率恰好是眼睛敏感的。频谱显示每个频率在光中的含量，就像音频频谱对音频信号所做的那样。

像 Winamp 这样的 MP3 播放器也会显示正在播放的音频信号的频谱（但它随时间变化，因为每次都对最近几毫秒重新计算频谱），音频信号的频谱可用于检测音乐中的低音、高音和人声，用于语音识别、音频滤波等。由于本教程是关于计算机图形学的，这里就不深入探讨了。

通常，频谱绘制在对数刻度而非线性刻度上，此时振幅以分贝（dB）表示。

傅里叶变换

那么如何计算给定信号的频谱呢？答案是使用傅里叶变换！

用于连续信号的连续傅里叶变换定义如下：



逆连续傅里叶变换（可以从频谱还原回信号）定义为：



F(w) 是频谱，其中 w 表示频率；f(x) 是时域信号，其中 x 表示时间。i 是 sqrt(-1)，参见复数理论。

注意两个变换之间的相似性，这正好解释了信号与其频谱之间的对偶性。

计算机无法处理连续信号，只能处理有限的离散信号。这些信号在时间上是有限的，并且只有一组离散的点（即信号是采样的）。傅里叶变换及其逆变换的一个性质是，离散信号的 FT 是周期性的。由于计算机要求信号和频谱都是离散的，信号和频谱都将是周期性的。但通过只取其中一个周期，我们就得到了有限信号。因此，如果您在计算机上对信号或图像进行 DFT，从数学上讲，该信号是无限重复的，或者图像是无限平铺的，频谱也是如此。一个很好的性质是，信号和频谱的离散点数相同，因此 128×128 像素图像的 DFT 也将有 128×128 个像素。

由于信号在时间上是有限的，积分的无穷边界可以替换为有限边界，积分符号可以替换为求和符号。因此 DFT 定义为：



逆变换定义为：



这在计算机上已经更具可编程性了。要编程实现 FT，需要对每个 n 进行计算，因此需要一个双重 for 循环（一个遍历每个 n，一个遍历每个 k，即求和）。虚数的指数可以用著名公式 e^ix = cosx + isinx 替换为余弦和正弦。本教程后续将进一步介绍编程实现。

DFT 有多种定义方式，例如可以在正向 DFT 而非逆变换中除以 N，或者在两者中都除以 sqrt(N)。在计算机上绘图时，在正向 DFT 中除以 N 可以得到最佳效果。

傅里叶变换的性质

信号通常用小写字母表示，其傅里叶变换或频谱用大写字母表示。信号与其频谱之间的关系通常记为 f(x) <--> F(w)，左侧为信号，右侧为其频谱。

傅里叶变换有一些有趣的性质，其中一些有助于理解为什么某些信号的频谱具有特定形状。

但请注意，这不是一个完整的列表，一些缩放因子（如 2*pi）被省略了，这些因子并不那么重要，可以省略，也取决于您使用频率还是角频率。信号处理手册以及维基百科和 Mathworld 等网站包含更完整、数学上更严格的 FT 性质表。我们最感兴趣的是函数的形状而非比例。

线性性

f(x)+g(x) <--> F(w)+G(w)

a*f(x) <--> a*F(w)

这意味着，如果您将两个信号相加/相减，它们的频谱也会相加/相减；如果增大/减小信号的振幅，其频谱的振幅也会以相同的比例增大/减小。

尺度性

f(a*x) <--> (1/a) * F(w/a)

这意味着，如果您在 x 方向上拉宽函数，其频谱在 x 方向上会变窄，反之亦然。振幅也会相应改变。

时移性

f(x-x0) <--> (exp(-i*w*x0)) * F(w)

由于时移时唯一发生的变化是傅里叶变换与一个虚数指数的乘法，因此时移不会在频谱的振幅中体现出来，只会影响相位。

频移性

(exp(-i*w0*x))*f(x)<-->F(w-w0)

这是时移性的对偶。

对偶性或对称性

若 f(x) <--> F(w)
则 F(x) <--> f(-w)

至少在忽略某些缩放因子的情况下成立。正是由于这个性质，例如，矩形脉冲的频谱是 sinc 函数，同时 sinc 函数的频谱也是矩形脉冲。

对称规则

以下只是部分对称规则：

卷积定理

卷积是两个函数之间的一种运算，其定义为积分，但也可以如下解释：

取 2 个函数，例如两个矩形函数 f1(x) 和 f2(x)。将其中一个矩形固定在某个位置。将第二个矩形关于 y 轴镜像（对矩形影响不大）。然后将第二个函数沿 u 值平移，u 从 -∞ 到 +∞。

结果函数 g(u)，即 f1(x) 和 f2(x) 的卷积，以 u 为参数。对于某个特定的 u，将 f1 与平移镜像后的 f2 相乘，然后取结果下方的面积（通过积分）。这个面积就是该 u 处 g(u) 的值！需要对每个 u 执行此操作才能得到完整结果。例如，对于两个矩形函数，当从 u=-∞ 开始时，g(u) 为 0，因为两个矩形不重叠，两个函数相乘的结果为零函数，面积为零。随着 u 向右移动，g(u) 始终保持为 0，直到两个矩形开始重叠。此后越向右，面积越大，直到达到最大值，g(u) 持续上升。然后，随着第二个矩形继续向右移动，重叠越来越少，g(u) 的值再次减小。最终 g(u) 再次变为零并保持到 +∞。因此结果 g(u) 是一个三角函数。

这就是不使用 FT 的绘图程序在二维中对滤波器所做的事情。图像是保持固定位置的函数，而滤波器矩阵在整幅图像上移动，您通过将矩阵的每个元素与矩阵所覆盖的图像中对应像素相乘来计算每个像素的值。对于大型滤波器矩阵，这需要大量计算。您可以在 Paint Shop Pro 和 Photoshop 等绘图程序中使用"用户自定义"滤波器来创建自己的滤波器矩阵。

相关运算与卷积几乎相同，只是不需要对第二个函数进行镜像。

卷积定理指出，在时域中对两个函数进行卷积，等价于在傅里叶域中将它们的频谱相乘，反之亦然。这种乘法运算比卷积简单得多。

由于滤波器也可以在时域和傅里叶域中描述（分别称为滤波器的冲激响应和传递函数），线性滤波器在时域中执行卷积，在频域中对应于将滤波器的传递函数与信号的频谱相乘。时域中滤波器的冲激响应是滤波器对狄拉克冲激（见前文）的响应，而滤波器的传递函数是冲激响应的傅里叶变换。

对于我们而言，卷积定理在实验信号和图像的傅里叶变换时将非常有用。

编程实现一维傅里叶变换

以下是一个简单程序的代码，它将计算给定信号的 FT，并绘制信号和 FT。将其放入 main.cpp 文件中。

首先声明变量和函数。N 将是信号中离散点的数量。对于信号及其频谱，分别使用小写字母 f 或 g 以及大写字母 F 或 G，这与数学中的惯例一致。

// working with fixed sizes for simplicity in this tutorial
const int N = 128;

double fRe[N]; //the function's real part, imaginary part, and amplitude
double fIm[N];
double fAmp[N];
double FRe[N]; //the FT's real part, imaginary part and amplitude
double FIm[N];
double FAmp[N];
const double pi = 3.1415926535897932384626433832795;

void DFT(int n, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm); //Calculates the DFT
void plot(int yPos, const double *g, double scale, bool trans, ColorRGB color); //plots a function
void calculateAmp(int n, double *ga, const double *gRe, const double *gIm); //calculates the amplitude of a complex function

主函数创建一个正弦信号，然后使用其他函数计算 FT 并绘图。

int main(int /*argc*/, char */*argv*/[])
{
  screen(640, 480, 0, "1D DFT");
  cls(RGB_White);

  for(int x = 0; x < N; x++)
  {
    fRe[x] = 25 * sin(x / 2.0); //generate a sine signal.  Put anything else you like here.
    fIm[x] = 0;
  }

  //Plot the signal
  calculateAmp(N, fAmp, fRe, fIm);
  plot(100, fAmp, 1.0, 0, ColorRGB(160, 160, 160));
  plot(100, fIm, 1.0, 0, ColorRGB(128, 255, 128));
  plot(100, fRe, 1.0, 0, ColorRGB(255, 0, 0));

  DFT(N, 0, fRe, fIm, FRe, FIm);

  //Plot the FT of the signal
  calculateAmp(N, FAmp, FRe, FIm);
  plot(350, FRe, 12.0, 1, ColorRGB(255, 128, 128));
  plot(350, FIm, 12.0, 1, ColorRGB(128, 255, 128));
  plot(350, FAmp, 12.0, 1, ColorRGB(0, 0, 0));

  drawLine(w / 2, 0, w / 2, h - 1, ColorRGB(128, 128, 255));
  redraw();
  sleep();

  return 0;
}

接下来是执行实际计算的函数（尽管这里使用的算法并不是最优的）：DFT 函数

void DFT(int n, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm)
{
  for(int w = 0; w < n; w++)
  {
    GRe[w] = GIm[w] = 0;
    for(int x = 0; x < n; x++)
    {
      double a = -2 * pi * w * x / float(n);
      if(inverse) a = -a;
      double ca = cos(a);
      double sa = sin(a);
      GRe[w] += gRe[x] * ca - gIm[x] * sa;
      GIm[w] += gRe[x] * sa + gIm[x] * ca;
    }
    if(!inverse)
    {
      GRe[w] /= n;
      GIm[w] /= n;
    }
  }
}

该函数接受 4 个数组作为参数：*gRe 和 *gIm 用于读取输入（这是信号的实部和虚部），*GRe 和 *GIm 用于写入结果（这是信号的 FT）。输入和输出都是复数，因此需要两个数组（一个实部，一个虚部）。n 是信号（和数组）的大小，即信号中离散点的数量。

第一个 for 循环，w 从 0 到 n-1，表示需要对结果的每个点进行计算。第二个 for 循环，x 从 0 到 n，表示输入信号的每个点都需要被考虑，这就是 DFT 数学定义中的求和。指数函数通过公式 e^ix = cosx + isinx 替换为余弦和正弦，实部和虚部分别计算。之后将值除以 n，使值足够小以便在屏幕上绘图。

将"inverse"设置为 true，还可以计算 IDFT（逆离散傅里叶变换），但本例中未使用。唯一的区别是指数函数的符号不同，以及不除以 n。

接下来，plot 函数接受一个数组作为参数，使用线段和实心圆将该数组中包含的信号或频谱以您指定的颜色绘制到屏幕上。您还可以使用"shift"参数将边缘移到中心，反之亦然。这对于更好地表示傅里叶变换很有用：通常 DFT 的结果在两侧有低频分量，中心有最高频率。由于结果实际上是周期性的，将两侧移到中心等同于将图移动半个周期。

void plot(int yPos, const double *g, double scale, bool shift, ColorRGB color)
{
  drawLine(0, yPos, w - 1, yPos, ColorRGB(128, 128, 255));
  for(int x = 1; x < N; x++)
  {
    int x1, x2, y1, y2;
    int x3, x4, y3, y4;
    x1 = (x - 1) * 5;
    //Get first endpoint, use the one half a period away if shift is used
    if(!shift) y1 = -int(scale * g[x - 1]) + yPos;
    else y1 = -int(scale * g[(x - 1 + N / 2) % N]) + yPos;
    x2 = x * 5;
    //Get next endpoint, use the one half a period away if shift is used
    if(!shift) y2 = -int(scale * g[x]) + yPos;
    else y2 = -int(scale * g[(x + N / 2) % N]) + yPos;
    //Clip the line to the screen so we won't be drawing outside the screen
    clipLine(x1, y1, x2, y2, x3, y3, x4, y4);
    drawLine(x3, y3, x4, y4, color);
    //Draw circles to show that our function is made out of discrete points
    drawDisk(x4, y4, 2, color);
  }
}

最后，calculateAmp 函数根据实部和虚部计算振幅，以便您也可以绘制振幅。

//Calculates the amplitude of *gRe and *gIm and puts the result in *gAmp
void calculateAmp(int n, double *gAmp, const double *gRe, const double *gIm)
{
  for(int x = 0; x < n; x++)
  {
    gAmp[x] = sqrt(gRe[x] * gRe[x] + gIm[x] * gIm[x]);
  }
}

该程序的输出如下：



红色函数是输入函数，即正弦波；底部的黑色函数是其频谱的振幅，由 DFT 函数计算得出。绿色曲线是虚部。

快速傅里叶变换

上述 DFT 函数可以正确计算离散傅里叶变换，但使用了两个 n 次的 for 循环，因此需要 O(n²) 次算术运算。更快的算法是快速傅里叶变换（FFT），它只需要 O(n*logn) 次运算。这对于非常大的 n 来说差别很大：如果 n 为 1024，DFT 函数需要 1048576（约 100 万）次循环，而 FFT 只需 10240 次。

FFT 的工作原理是将信号分成两半：一半包含所有偶数索引的值，另一半包含所有奇数索引的值。然后递归地再次分割这两半，依此类推。因此，要求信号的大小 n 必须是 2 的幂。这里描述的算法是 Radix-2 Cooley-Tukey FFT 算法（由 Cooley 和 Tukey 于 1965 年开发）。

还存在其他算法，可以快速处理非 2 的幂大小的信号，甚至是大质数。例如混合基 Cooley-Tukey 算法以及其他算法。但这超出了本教程的范围。上面库部分链接的代码支持非 2 的幂，因此对于实际应用我建议使用类似的库。但在本教程中，我们继续使用更简单的代码进行演示。

以下是 DFT 函数的新版本，称为 FFT，使用上述 radix-2 算法。您可以将其添加到之前的程序中，并将函数调用改为 FFT 来测试它。

首先，该函数将计算 n 的对数，这是算法所需要的。

void FFT(int n, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm)
{
  //Calculate m=log_2(n)
  int m = 0;
  int p = 1;
  while(p < n)
  {
    p *= 2;
    m++;
  }

接下来执行位逆序（Bit Reversal）。由于 FFT 每次将信号分成偶数和奇数分量，算法要求输入信号已经按照奇数分量在前、偶数分量在后的形式排列，每个半部分依此类推。这等同于将每个值的索引的二进制位逆序（例如，索引为 0001 的第二个值现在将获得索引 1000，移至中间位置），因此得名"位逆序"。

  //Bit reversal
  GRe[n - 1] = gRe[n - 1];
  GIm[n - 1] = gIm[n - 1];
  int j = 0;
  for(int i = 0; i < n - 1; i++)
  {
    GRe[i] = gRe[j];
    GIm[i] = gIm[j];
    int k = n / 2;
    while(k <= j)
    {
      j -= k;
      k /= 2;
    }
    j += k;
  }

接下来执行实际的 FFT。从数学上讲，它是一个递归函数，但这是一个修改版本，以非递归方式工作——每次将完整数组传递给递归函数实在太麻烦了。ca 和 sa 仍然表示余弦和正弦，但使用半角公式和初始值 -1.0 和 0.0（pi 的余弦和正弦）计算：每次循环，ca 和 sa 表示前一次循环角度一半的余弦和正弦。

  //Calculate the FFT
  double ca = -1.0;
  double sa = 0.0;
  int l1 = 1, l2 = 1;
  for(int l = 0; l < m; l++)
  {
    l1 = l2;
    l2 *= 2;
    double u1 = 1.0;
    double u2 = 0.0;
    for(int j = 0; j < l1; j++)
    {
      for(int i = j; i < n; i += l2)
      {
        int i1 = i + l1;
        double t1 = u1 * GRe[i1] - u2 * GIm[i1];
        double t2 = u1 * GIm[i1] + u2 * GRe[i1];
        GRe[i1] = GRe[i] - t1;
        GIm[i1] = GIm[i] - t2;
        GRe[i] += t1;
        GIm[i] += t2;
      }
      double z =  u1 * ca - u2 * sa;
      u2 = u1 * sa + u2 * ca;
      u1 = z;
    }
    sa = sqrt((1.0 - ca) / 2.0);
    if(!inverse) sa =- sa;
    ca = sqrt((1.0 + ca) / 2.0);
  }

最后，如果不是逆 DFT，则将值除以 n。

  //Divide through n if it isn't the IDFT
  if(!inverse)
  for(int i = 0; i < n; i++)
  {
    GRe[i] /= n;
    GIm[i] /= n;
  }
}

该函数再次从 *gRe 和 *gIm 读取，并将结果存储在 *GRe 和 *GIm 中。

FFT 部分现在有 3 个嵌套循环，但只有第一个从 0 到 n-1，其他两个合起来只循环 log_2(n) 次。对于这个只有 128 个值的小信号，您可能不会注意到速度的提升，因为 DFT 和 FFT 都是瞬间计算完成的（除非您在 40 年代的计算机上工作），但对于更大的信号、二维信号，以及需要实时反复计算 FT 的应用，这将产生重要差异。

滤波器

在傅里叶域中，应用滤波器非常简单。您只需将频谱与滤波器的传递函数相乘即可。

低通滤波器只允许低频分量通过（例如音乐中的低音），而高通滤波器只允许高频分量通过。您还可以制作带通和带阻滤波器，分别允许或阻止特定频率的分量。放大器将整个频谱乘以一个常数，当然您也可以制作放大特定频率的滤波器，或者使高频越来越弱的三角滤波器等。

例如，要制作低通滤波器，将频谱与原点附近的矩形函数相乘即可。所有低频分量将被保留，而高频分量则乘以零，从而被滤除。

以下是一个程序的部分代码，该程序首先允许您选择不同的信号并查看其 FT，按下空格键后，可以选择对频谱应用不同的滤波器，并计算逆 FT 以查看原始信号受滤波器影响的效果。

完整代码可在此处下载：fourier1d.cpp。

首先声明所有变量和函数。

const int N = 128; //yeah ok, we use old fashioned fixed size arrays here

double fRe[N]; //the function's real part, imaginary part, and amplitude
double fIm[N];
double fAmp[N];
double FRe[N]; //the FT's real part, imaginary part and amplitude
double FIm[N];
double FAmp[N];
const double pi = 3.1415926535897932384626433832795;

void FFT(int n, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm); //Calculates the DFT
void plot(int yPos, const double *g, double scale, bool trans, ColorRGB color); //plots a function
void calculateAmp(int n, double *ga, const double *gRe, const double *gIm); //calculates the amplitude of a complex function

主函数从这里开始。它进入第一个循环，该循环允许您选择一个函数。参数 p1 和 p2 可以修改以调整某些函数的振幅和/或宽度。

int main(int /*argc*/, char */*argv*/[])
{
  screen(640, 480, 0, "Fourier Transform and Filters");
  for(int x = 0; x < N; x++) fRe[x] = fIm[x] = 0;
  bool endloop = 0, changed = 1;
  int n = N;
  double p1 = 25.0, p2 = 2.0;
  while(!endloop)
  {
    if(done()) end();
    readKeys();

这是在循环内读取按键后的处理。每个按键都关联一个动作来创建特定函数，例如按下 g 键（SDLK_g）会将两个正弦波之和放入 fRe[x]。方向键也有一些动作，用于将虚部设置为 0、等于实部，或实部的修改版本。最后，空格键将"endloop"设置为 true，从而结束第一个循环。

这里没有包含所有按键，因为代码非常繁杂冗长。这里只展示几个作为示例，其余的在可下载的源文件中。

    if(keyPressed(SDLK_z)) {for(int x = 0; x < n; x++) fRe[x] = 0; changed = 1;} //zero
    if(keyPressed(SDLK_a)) {for(int x = 0; x < n; x++) fRe[x] = p1; changed = 1;} //constant
    if(keyPressed(SDLK_b)) {for(int x = 0; x < n; x++) fRe[x] = p1 * sin(x / p2); changed = 1;} //sine


    //ETCETERA... *snip*

    if(keyPressed(SDLK_SPACE)) endloop = 1;

这是第一个循环的最后部分，在这里绘制所选函数，计算其 FT，并绘制 FT。

    if(changed)
    {
      cls(RGB_White);
      calculateAmp(N, fAmp, fRe, fIm);
      plot(100, fAmp, 1.0, 0, ColorRGB(160, 160, 160));
      plot(100, fIm, 1.0, 0, ColorRGB(128, 255, 128));
      plot(100, fRe, 1.0, 0, ColorRGB(255, 0, 0));

      FFT(N, 0, fRe, fIm, FRe, FIm);
      calculateAmp(N, FAmp, FRe, FIm);
      plot(350, FRe, 12.0, 1, ColorRGB(255, 128, 128));
      plot(350, FIm, 12.0, 1, ColorRGB(128, 255, 128));
      plot(350, FAmp, 12.0, 1, ColorRGB(0, 0, 0));

      drawLine(w / 2, 0, w / 2, h - 1, ColorRGB(128, 128, 255));
      print("Press a-z to choose a function, press the arrows to change the imaginary part, press space to accept.", 0, 0, RGB_Black);
      redraw();
    }
    changed = 0;
  }

第一个循环结束后，主函数的第二部分开始。首先创建一个额外的频谱缓冲区，以便同时保存原始频谱和其滤波版本。然后开始第二个循环。

  //The original FRe2, FIm2 and FAmp2 will become multiplied by the filter
  double FRe2[N];
  double FIm2[N];
  double FAmp2[N];
  for(int x = 0; x < N; x++)
  {
    FRe2[x] = FRe[x];
    FIm2[x] = FIm[x];
  }
  cls();
  changed = 1; endloop = 0;
  while(!endloop)
  {
    if(done()) end();
    readKeys();

现在所有按键都以某种方式改变频谱：它们将应用滤波器。低通滤波器（LP）将频谱设置为 0，原点附近除外。高通滤波器（HP）则相反。

这里没有显示所有按键，因为代码非常繁杂冗长。只展示几个作为示例，其余的以更紧凑的形式包含在可下载的源文件中。


    if(keyPressed(SDLK_z))
    {
      for(int x = 0; x < n; x++)
      {
        FRe2[x] = FRe[x];
        FIm2[x] = FIm[x];
      }
      changed = 1;
    } //no filter
    if(keyPressed(SDLK_a))
    {
      for(int x = 0; x < n; x++)
      {
        if(x < 44 || x > N - 44)
        {
          FRe2[x] = FRe[x];
          FIm2[x] = FIm[x];
        }
        else FRe2[x] = FIm2[x] = 0;
      }
      changed = 1;
    } //LP

    //ETCETERA... *snip*

同样，如果按下某个键，将绘制新的频谱，并执行逆 FFT 以绘制原始函数的滤波版本。按下 Escape 键程序将退出。


    if(changed)
    {
      cls(RGB_White);
      FFT(N, 1, FRe2, FIm2, fRe, fIm);
      calculateAmp(N, fAmp, fRe, fIm);
      plot(100, fAmp, 1.0, 0, ColorRGB(160, 160, 160));
      plot(100, fIm, 1.0, 0, ColorRGB(128, 255, 128));
      plot(100, fRe, 1.0, 0, ColorRGB(255, 0, 0));

      calculateAmp(N, FAmp2, FRe2, FIm2);

      plot(350, FRe2, 12.0, 1, ColorRGB(255, 128, 128));
      plot(350, FIm2, 12.0, 1, ColorRGB(128, 255, 128));
      plot(350, FAmp2, 12.0, 1, ColorRGB(0, 0, 0));

      drawLine(w/2, 0,w/2,h-1, ColorRGB(128, 128, 255));
      print("Press a-z to choose a filter. Press esc to quit.", 0, 0, RGB_Black);
      redraw();
    }
    changed=0;
  }

  return 0;
}

FFT、plot 和 calculateAmp 函数与前面解释的相同，这里不再重复。

这段代码没有太多新内容，FFT、plot 和 calculateAmp 函数仍然相同，但通过运行这个程序，您可以看到低通（LP）、高通（HP）以及其他几种滤波器的效果。只需按照屏幕上的说明操作即可。

例如，这里展示了一个由 5 个正弦波和直流分量之和构成的信号（按"p"键获得），以及其频谱，然后分别对其应用 LP 和 HP：

原始信号：


LP 滤波器：


HP 滤波器

这里对阶跃函数（按"e"键获得）执行相同操作：


对其应用 LP 滤波器：



对其应用 HP 滤波器。该滤波器几乎只滤除直流分量，因此信号中只剩下两条几乎平坦的直线：



图像的二维傅里叶变换

将傅里叶变换扩展到二维实际上非常简单。首先对图像的每一行进行一维 FT，然后对结果的每一列再进行一维 FT。

实现方法是：使用一维 FT 函数逐行对图像进行 FT，然后转置图像（使列变为行），再对每行进行一维 FT，最后再次转置图像以恢复原始方向。这种方法的优点是操作的内存更加连续，有利于 CPU 缓存（速度更快），但需要中间的转置操作。另一种方法是为一维 FT 函数添加"步幅（stride）"参数，即当前处理的相邻值在内存中的间距。处理列时将步幅设置为行大小，处理行时设置为 1。此外，由于我们有 3 个颜色通道，步幅还要乘以 3。这种方法稍微简单一些，因此本教程将使用此方法。


以下是一个实现上述功能的程序代码，其中提供了两个实现相同功能的替代函数：一个使用慢速 DFT 的二维版本，另一个使用快速傅里叶变换（FFT）的二维版本。程序首先计算图像的 FT，然后计算结果的逆 FT，以验证公式是否正确工作：如果能还原出原始图像则说明正确。您可以在主函数中自由切换 DFT2D() 和 FFT2D() 的调用，您会注意到 DFT2D() 函数非常慢，而 FFT2D() 函数运行非常快。

由于 RGB 图像有 3 个颜色通道，FT 分别对每个颜色通道计算，实际上计算了 3 个灰度 FT。

程序需要一个 24 位彩色的 128×128 PNG 图像 pics/test.png（路径相对于程序）。可在 photos.zip 中找到一张。
.
再次首先声明所有变量和函数。
W 和 H 是图像的宽度和高度（以像素为单位）。
信号及其频谱的数组现在是三维的：2 个维度用于大小，1 个维度用于 RGB 颜色分量。这里不能使用 ColorRGB 类，因为 FT 需要更高的精度，因此颜色分量存储在 double 数组中；索引 0 表示红色，索引 1 表示绿色，索引 2 表示蓝色。

这里有两个版本的 FT 函数：慢速的 DFT2D 和快速的 FFT2D。

/// working with fixed sizes for simplicity in this tutorial
const int W = 128; //the width of the image
const int H = 128; //the height of the image

double fRe[H][W][3], fIm[H][W][3], fAmp[H][W][3]; //the signal's real part, imaginary part, and amplitude
double FRe[H][W][3], FIm[H][W][3], FAmp[H][W][3]; //the FT's real part, imaginary part and amplitude
double fRe2[H][W][3], fIm2[H][W][3], fAmp2[H][W][3]; //will become the signal again after IDFT of the spectrum
double FRe2[H][W][3], FIm2[H][W][3], FAmp2[H][W][3]; //filtered spectrum

double pi = 3.1415926535897932384626433832795;

void draw(int xpos, int yPos, int w, int h, const double *g, bool shift, bool neg128);
void DFT2D(int w, int h, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm);
void FFT2D(int w, int h, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm);
void calculateAmp(int w, int h, double *ga, const double *gRe, const double *gIm);

主函数首先加载一幅图像，然后将信号设置为该图像。

接下来绘制图像。绘制实部、虚部（值为 0！）以及振幅（看起来与实部相同）。

然后使用 DFT2D 函数计算图像的 FT，您也可以将其改为 FFT2D：速度会更快。然后绘制刚刚计算的频谱。

最后，再次计算该频谱的逆 FT 并绘制结果，这仅用于验证函数是否正确工作：如果能得到原始图像，则说明正确。

int main(int /*argc*/, char */*argv*/[])
{
  screen(N * 3, M * 3, 0, "2D DFT and FFT");

  std::vector<ColorRGB> img;
  unsigned long dummyw, dummyh;
  if(loadImage(img, dummyw, dummyh, "pics/test.png"))
  {
    print("image pics/test.png not found");
    redraw();
    sleep();
    cls();
  }

  //set signal to the image
  for(int y = 0; y < H; y++)
  for(int x = 0; x < W; x++)
  {
    fRe[y][x][0] = img[W * y + x].r;
    fRe[y][x][1] = img[W * y + x].g;
    fRe[y][x][2] = img[W * y + x].b;
  }

  //draw the image (real, imaginary and amplitude)
  calculateAmp(W, H, fAmp[0][0], fRe[0][0], fIm[0][0]);
  draw(0, 0, W, H, fRe[0][0], 0, 0);
  draw(W, 0, W, H, fIm[0][0], 0, 0);
  draw(2 * W, 0, W, H, fAmp[0][0], 0, 0);

  //calculate and draw the FT of the image
  DFT2D(W, H, 0, fRe[0][0], fIm[0][0], FRe[0][0], FIm[0][0]); //Change this to FFT2D to make it faster
  calculateAmp(W, H, FAmp[0][0], FRe[0][0], FIm[0][0]);
  draw(0, H, W, H, FRe[0][0], 1, 1);
  draw(W, H, W, H, FIm[0][0], 1, 1);
  draw(2 * W, H, W, H, FAmp[0][0], 1, 0);

  //calculate and draw the image again
  DFT2D(W, H, 1, FRe[0][0], FIm[0][0], fRe2[0][0], fIm2[0][0]); //Change this to FFT2D to make it faster
  calculateAmp(W, H, fAmp2[0][0], fRe2[0][0], fIm2[0][0]);
  draw(0, 2 * H, W, H, fRe2[0][0], 0, 0);
  draw(W, 2 * H, W, H, fIm2[0][0], 0, 0);
  draw(2 * W, 2 * H, W, H, fAmp2[0][0], 0, 0);

  redraw();
  sleep();
  return 0;
}

这里再次出现一维 DFT 函数，但这次添加了步幅（stride）和因子（factor）参数。如上所述，步幅参数允许更大的步进，以便在我们的三维数组中处理列和颜色通道。因子参数的存在是因为我们需要控制用什么值（而不仅仅是 n）来除结果。

void DFT(int n, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm, int stride, double factor)
{
  for(int n = 0; n < n; n++)
  {
    GRe[n * stride] = GIm[n * stride] = 0;
    for(int x = 0; x < n; x++)
    {
      double a = -2 * pi * n * x / float(n);
      if(inverse) a = -a;
      double ca = cos(a);
      double sa = sin(a);
      GRe[n * stride] += gRe[x * stride] * ca - gIm[x * stride] * sa;
      GIm[n * stride] += gRe[x * stride] * sa + gIm[x * stride] * ca;
    }
    //Divide through the factor, e.g. n
    GRe[n * stride] /= factor;
    GIm[n * stride] /= factor;
  }
}

然后，DFT2D 函数对每一行和每一列分别调用一维 DFT 函数，并对每个颜色分量分别处理。该函数现在有两个维度参数：w 和 h（有时也称为 n 和 m），因为图像是二维的。参数 inverse 可以设置为 true 以执行逆变换。同样，*gRe 和 *gIm 是输入数组，*GRe 和 *GIm 是输出数组。在这个二维版本中，正向 DFT 的值除以 h，逆 DFT 的值除以 w。也有定义要求在一个方向上除以 h*w，另一个方向上除以 1，或者两个方向都除以 sqrt(h*w)，但只要正向 FT 和逆变换合起来能还原原始图像，选哪种都没有太大关系。不过，在每个方向上使用大约 h*w 的平方根，可以给出与输入图像类似的值范围，绘图效果最好。

void DFT2D(int w, int h, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm)
{
  //temporary buffers
  std::vector<double> Gr2(m * n * 3);
  std::vector<double> Gi2(m * n * 3);

  for(int y = 0; y < h; y++) // for each row
  for(int c = 0; c < 3; c++) // for each color channel
  {
    DFT(w, inverse, &gRe[y * w * 3 + c], &gIm[y * w * 3 + c], &Gr2[y * w * 3 + c], &Gi2[y * w * 3 + c], 3, 1);
  }
  for(int x = 0; x < w; x++) // for each column
  for(int c = 0; c < 3; c++) // for each color channel
  {
    DFT(h, inverse, &Gr2[x * 3 + c], &Gi2[x * 3 + c], &GRe[x * 3 + c], &GIm[x * 3 + c], w * 3, inverse ? w : h);
  }
}

这里再次给出快速傅里叶变换，现在加入了步幅和因子参数。

void FFT(int n, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm, int stride, double factor)
{
  //Calculate h=log_2(n)
  int h = 0;
  int p = 1;
  while(p < n)
  {
    p *= 2;
    h++;
  }
  //Bit reversal
  GRe[(n - 1) * stride] = gRe[(n - 1) * stride];
  GIm[(n - 1) * stride] = gIm[(n - 1) * stride];
  int j = 0;
  for(int i = 0; i < n - 1; i++)
  {
    GRe[i * stride] = gRe[j * stride];
    GIm[i * stride] = gIm[j * stride];
    int k = n / 2;
    while(k <= j)
    {
      j -= k;
      k /= 2;
    }
    j += k;
  }
  //Calculate the FFT
  double ca = -1.0;
  double sa = 0.0;
  int l1 = 1, l2 = 1;
  for(int l = 0; l < h; l++)
  {
    l1 = l2;
    l2 *= 2;
    double u1 = 1.0;
    double u2 = 0.0;
    for(int j = 0; j < l1; j++)
    {
      for(int i = j; i < n; i += l2)
      {
        int i1 = i + l1;
        double t1 = u1 * GRe[i1 * stride] - u2 * GIm[i1 * stride];
        double t2 = u1 * GIm[i1 * stride] + u2 * GRe[i1 * stride];
        GRe[i1 * stride] = GRe[i * stride] - t1;
        GIm[i1 * stride] = GIm[i * stride] - t2;
        GRe[i * stride] += t1;
        GIm[i * stride] += t2;
      }
      double z =  u1 * ca - u2 * sa;
      u2 = u1 * sa + u2 * ca;
      u1 = z;
    }
    sa = sqrt((1.0 - ca) / 2.0);
    if(!inverse) sa =- sa;
    ca = sqrt((1.0 + ca) / 2.0);
  }
  //Divide through the factor, e.g. n
  for(int i = 0; i < n; i++)
  {
    GRe[i * stride] /= factor;
    GIm[i * stride] /= factor;
  }
}

扩展到二维与 DFT2D 完全相同，只是使用 FFT 函数

void FFT2D(int w, int h, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm)
{
  //temporary buffers
  std::vector<double> Gr2(m * n * 3);
  std::vector<double> Gi2(m * n * 3);

  for(int y = 0; y < h; y++) // for each row
  for(int c = 0; c < 3; c++) // for each color channel
  {
    FFT(w, inverse, &gRe[y * w * 3 + c], &gIm[y * w * 3 + c], &Gr2[y * w * 3 + c], &Gi2[y * w * 3 + c], 3, 1);
  }
  for(int x = 0; x < w; x++) // for each column
  for(int c = 0; c < 3; c++) // for each color channel
  {
    FFT(h, inverse, &Gr2[x * 3 + c], &Gi2[x * 3 + c], &GRe[x * 3 + c], &GIm[x * 3 + c], w * 3, inverse ? w : h);
  }
}

draw 函数的作用与一维版本中的 plot 函数相同：在屏幕上显示信号。现在该函数将把您传入的信号或频谱数组绘制为 24 位 RGB 图像，左上角位于 xpos,yPos。shift 参数可用于将四角移到中心、将中心移到四角，这对于 FT 将直流分量移到中心非常有用。neg128 参数可用于绘制可能有负像素值的图像：此时灰色（RGB 颜色 128, 128, 128）表示 0，较暗的颜色表示负值，较亮的颜色表示正值。

//Draws an image
void draw(int xPos, int yPos, int n, int m, const double *g, bool shift, bool neg128) //g is the image to be drawn
{
  for(int y = 0; y < m; y++)
  for(int x = 0; x < n; x++)
  {
    int x2 = x, y2 = y;
    if(shift) {x2 = (x + n / 2) % n; y2 = (y + m / 2) % m;} //Shift corners to center
    ColorRGB color;
    //calculate color values out of the floating point buffer
    color.r = int(g[3 * w * y2 + 3 * x2 + 0]);
    color.g = int(g[3 * w * y2 + 3 * x2 + 1]);
    color.b = int(g[3 * w * y2 + 3 * x2 + 2]);

    if(neg128) color = color + RGB_Gray;
    //negative colors give confusing effects so set them to 0
    if(color.r < 0) color.r = 0;
    if(color.g < 0) color.g = 0;
    if(color.b < 0) color.b = 0;
    //set color components higher than 255 to 255
    if(color.r > 255) color.r = 255;
    if(color.g > 255) color.g = 255;
    if(color.b > 255) color.b = 255;
    //plot the pixel
    pset(x + xPos, y + yPos, color);
  }
}

calculateAmp 函数也与一维版本相同，只是增加了一个用于第二维度的循环和一个用于颜色分量的循环。该函数用于生成信号和频谱的振幅，以便绘制它们。

//Calculates the amplitude of *gRe and *gIm and puts the result in *gAmp
void calculateAmp(int n, int m, double *gAmp, const double *gRe, const double *gIm)
{
  for(int y = 0; y < m; y++)
  for(int x = 0; x < n; x++)
  for(int c = 0; c < 3; c++)
  {
    gAmp[w * 3 * y + 3 * x + c] = sqrt(gRe[w * 3 * y + 3 * x + c] * gRe[w * 3 * y + 3 * x + c] + gIm[w * 3 * y + 3 * x + c] * gIm[w * 3 * y + 3 * x + c]);
  }
}

注意，由于 C++ 的工作方式，数组以特殊方式传递给函数，在函数内部必须将其视为一维数组。


请注意，上述代码并不是最优的，仅为演示而保持简单。有更快的库可以实现这一功能，例如 FFTW、kissfft 等。这些库速度更快、功能更完整，并且也能快速支持非 2 的幂大小。对于实际问题，建议使用这类经过充分测试的快速实现。

程序的输出如下：

2D DFT

顶行是图像的实部、虚部和振幅。图像的虚部为 0。

中间行是图像频谱的实部、虚部和振幅。实部和虚部可以有负值，因此启用了"neg128"选项来绘制它们，灰色表示 0。由于振幅始终为正，黑色表示 0。

底行是图像，通过对频谱进行逆 DFT 计算得到。看起来与原始图像几乎相同，由于舍入误差，少数像素的颜色可能略有不同。

文本"99% done"只在使用 DFT2D 函数时才会出现，因为它非常慢（在 Athlon 1700+ 处理器上处理 128×128 图像需要 30 秒），进度计数器很有用。FFT2D 立即完成相同的计算，因此不需要这样的计数器。

如果研究图像的频谱，您会看到穿过中心的水平和垂直线。这是由图像的边缘引起的：从数学上讲，FT 是对无限平铺的图像计算的，从图像一侧到另一侧时会有突变，从而在频谱中产生这些线。

一些说明：

图像的频谱

您可以将 test.png 替换为任何 24 位 128×128 彩色 PNG。您甚至可以使用其他尺寸，但需要修改代码顶部的参数 N 和 M。注意，FFT2D 函数只适用于大小为 2 的幂的图像（64×64、128×256、256×256 等均可），否则您将不得不使用较慢的 DFT2D 函数（或支持非 2 的幂的其他优质 FFT 库）。

图像的频谱本身不是最有趣的（对其应用滤波器然后再次进行逆 FT 更有趣），但通过观察其频谱，您仍然可以了解一些关于图像的信息：以下是几个示例，每张图片显示图像的实部和频谱的振幅。要查看频谱的实部和虚部，请自行将图像通过程序处理。

这幅图像是加了直流分量的正弦函数（因为图像只能有正像素值）：每个像素的颜色为 128 + 128 * sin(pi*y/16.0)。如果您还记得一维信号的频谱，您会记得正弦函数在负侧和正侧各有一个峰值，直流分量在零处有一个峰值。在频谱中可以看到同样的情况：中心点是直流分量，另外两个点表示正弦函数的频率，一个点只是另一个点的镜像。在 x 方向上没有像素，因为图像在该方向上处处相同。



这里使用更高频率的正弦函数生成图像，因此两个点离原点更远，表示更高的频率。还记得傅里叶变换的尺度性吗？这里很好地说明了这一点：图像收缩，其频谱变宽。



在下一幅图像中，正弦函数不能平铺：如果将同一图像叠放两次，在 y 方向上将看不到正弦波了（而在前两个例子中可以看到）。频谱也不再将其视为真正的正弦波，它不只包含两个点，而是一条线。



二维 FT 的一个性质是，如果旋转图像，频谱也会沿同一方向旋转：



下图显示了 2 个正弦函数之和，每个正弦函数方向不同（等离子体效果通常就是这样制作的）：



下图显示了图像中的线条如何在频谱中生成垂直线：



频谱中的斜线显然是由天空到山峰的剧烈过渡引起的。



矩形函数的 FT 是二维 sinc 函数：



这是一个可平铺纹理的 FT，由于它是可平铺的，水平和垂直边缘没有突变，因此频谱中没有水平和垂直线：



二维滤波器

现在来介绍 FT 在图像上最有趣的应用：您可以轻松地对频谱应用各种滤波器，同样通过将其与二维传递函数相乘，或将频谱的某些部分置零，或对频谱应用某些奇特公式。之后，对修改后的频谱进行逆 FT 即可得到滤波后的图像。

以下是一个程序的部分代码，该程序允许您选择不同的图像（一个 PNG、其修改版本以及一些生成的函数），然后按下空格键后，可以应用不同的滤波器并查看其效果。完整代码可在此处下载：fourier2d.cpp。

这段代码没有太多新内容，运行它并研究滤波器的结果更有趣。

再次，首先是所有函数和变量的声明。f2 和 F2 将成为原始图像和频谱的滤波版本的缓冲区。

// working with fixed sizes for simplicity in this tutorial
const int W = 128; //the width of the image
const int H = 128; //the height of the image

double fRe[H][W][3], fIm[H][W][3], fAmp[H][W][3]; //the signal's real part, imaginary part, and amplitude
double FRe[H][W][3], FIm[H][W][3], FAmp[H][W][3]; //the FT's real part, imaginary part and amplitude
double fRe2[H][W][3], fIm2[H][W][3], fAmp2[H][W][3]; //will become the signal again after IDFT of the spectrum
double FRe2[H][W][3], FIm2[H][W][3], FAmp2[H][W][3]; //filtered spectrum

double pi = 3.1415926535897932384626433832795;

void draw(int xpos, int yPos, int w, int h, const double *g, bool shift, bool neg128);
void DFT2D(int w, int h, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm);
void FFT2D(int w, int h, bool inverse, const double *gRe, const double *gIm, double *GRe, double *GIm);
void calculateAmp(int w, int h, double *ga, const double *gRe, const double *gIm);

主函数首先加载一幅 PNG 图像，然后开始第一个循环，该循环让您尝试 PNG 图像的不同修改版本以及其他一些生成的信号：

int main(int /*argc*/, char */*argv*/[])
{
  screen(3 * W,4 * H + 8, 0,"2D FFT and Filters");
  std::vector img;
  unsigned long dummyw, dummyh;
  if(loadImage(img, dummyw, dummyh, "pics/test.png"))
  {
    print("image pics/test.png not found");
    redraw();
    sleep();
    cls();
  }
  //set signal to the image
  for(int y = 0; y < H; y++)
  for(int x = 0; x < W; x++)
  {
    fRe[y][x][0] = img[W * y + x].r;
    fRe[y][x][1] = img[W * y + x].g;
    fRe[y][x][2] = img[W * y + x].b;
  }

  int ytrans=8; //translate everything a bit down to put the text on top

  //set new FT buffers
  for(int y = 0; y < H; y++)
  for(int x = 0; x < W; x++)
  for(int c = 0; c < 3; c++)
  {
    FRe2[y][x][c] = FRe[y][x][c];
    FIm2[y][x][c] = FIm[y][x][c];
  }

  bool changed = 1, endloop = 0;
  while(!endloop)
  {
    if(done()) end();
    readKeys();

以下是第一个循环输入部分的一小段代码，这里看起来相当杂乱，只显示了几个按键，所有按键的代码都包含在可下载文件中。

    if(keyPressed(SDLK_z) || changed)
    {
      for(int y = 0; y < H; y++)
      for(int x = 0; x < W; x++)
      for(int c = 0; c < 3; c++)
      {
        fRe2[y][x][c] = fRe[y][x][c];
      } changed = 1;
    } //no effect
    if(keyPressed(SDLK_a))
    {
      for(int y = 0; y < H; y++)
      for(int x = 0; x < W; x++)
      for(int c = 0; c < 3; c++)
      {
        fRe2[y][x][c] = 255 - fRe[y][x][c];
      }
    changed = 1;} //negative
    if(keyPressed(SDLK_b))
    {
      for(int y = 0; y < H; y++)
      for(int x = 0; x < W; x++)
      for(int c = 0; c < 3; c++)
      {
        fRe2[y][x][c] = fRe[y][x][c]/2;
      }
      changed = 1;
    } //half amplitude


    //ETCETERA... *snip*

    if(keyPressed(SDLK_SPACE)) endloop = 1;

如果按下某个键，将绘制图像的新版本及其频谱。如果第一个循环完成（通过按空格键），则初始化第二个循环：

    if(changed)
    {
      //Draw the image and its FT
      calculateAmp(W, H, fAmp2[0][0], fRe2[0][0], fIm2[0][0]);
      draw(0, ytrans, W, H, fRe2[0][0], 0, 0);  draw(W, 0+ytrans, W, H, fIm2[0][0], 0, 0);  draw(2 * W, 0+ytrans, W, H, fAmp2[0][0], 0, 0); //draw real, imag and amplitude
      FFT2D(W, H, 0, fRe2[0][0], fIm2[0][0], FRe[0][0], FIm[0][0]);
      calculateAmp(W, H, FAmp[0][0], FRe[0][0], FIm[0][0]);
      draw(0,H + ytrans, W, H, FRe[0][0], 1, 1);  draw(W, H + ytrans, W, H, FIm[0][0], 1, 1);  draw(2 * W, H + ytrans, W, H, FAmp[0][0], 1, 0); //draw real, imag and amplitude
      print("Press a-z for effects and space to accept", 0, 0, RGB_White, 1, ColorRGB(128, 0, 0));
      redraw();
    }
    changed = 0;
  }
  changed = 1;
  while(!done())
  {
    readKeys();

以下是应用滤波器的输入代码的一小部分：给出了低通、高通和带通滤波器。其余按键的代码在可下载文件中。

    if(keyPressed(SDLK_g))
    {
      for(int y = 0; y < H; y++)
      for(int x = 0; x < W; x++)
      for(int c = 0; c < 3; c++)
      {
        if(x * x + y * y < 64 || (W - x) * (W - x) + (H - y) * (H - y) < 64 || x * x + (H - y) * (H - y) < 64 || (W - x) * (W - x) + y * y < 64)
        {
          FRe2[y][x][c] = FRe[y][x][c];
          FIm2[y][x][c] = FIm[y][x][c];
        }
        else
        {
          FRe2[y][x][c] = 0;
          FIm2[y][x][c] = 0;
        }
      }
      changed = 1;
    } //LP filter

    if(keyPressed(SDLK_h))
    {
      for(int y = 0; y < H; y++)
      for(int x = 0; x < W; x++)
      for(int c = 0; c < 3; c++)
      {
        if(x * x + y * y > 16 && (W - x) * (W - x) + (H - y) * (H - y) > 16 && x * x + (H - y) * (H - y) > 16 && (W - x) * (W - x) + y * y > 16)
        {
          FRe2[y][x][c] = FRe[y][x][c];
          FIm2[y][x][c] = FIm[y][x][c];
        }
        else
        {
          FRe2[y][x][c] = 0;
          FIm2[y][x][c] = 0;
        }
      }
      changed = 1;
    } //HP filter

    if(keyPressed(SDLK_t))
    {
      for(int y = 0; y < H; y++)
      for(int x = 0; x < W; x++)
      for(int c = 0; c < 3; c++)
      {
        if((x * x + y * y < 256 || (W - x) * (W - x) + (H - y) * (H - y) < 256 || x * x + (H - y) * (H - y) < 256 || (W - x) * (W - x) + y * y < 256) && (x * x + y * y > 128 && (W - x) * (W - x) + (H - y) * (H - y) > 128 && x * x + (H - y) * (H - y) > 128 && (W - x) * (W - x) + y * y > 128))
        {
          FRe2[y][x][c] = FRe[y][x][c];
          FIm2[y][x][c] = FIm[y][x][c];
        }
        else
        {
          FRe2[y][x][c] = 0;
          FIm2[y][x][c] = 0;
        }
      }
      changed = 1;
    } //BP filter

    //ETCETERA... *snip*

这是第二个循环中绘制滤波后频谱和信号的部分，也是主函数的最后部分：

    if(changed)
    {
      //after pressing a key: calculate the inverse!
      FFT2D(W, H, 1, FRe2[0][0], FIm2[0][0], fRe[0][0], fIm[0][0]);
      calculateAmp(W, H, fAmp[0][0], fRe[0][0], fIm[0][0]);
      draw(0, 3 * H + ytrans, W, H, fRe[0][0], 0, 0);  draw(W, 3 * H + ytrans, W, H, fIm[0][0], 0, 0);  draw(2 * W, 3 * H + ytrans, W, H, fAmp[0][0], 0, 0);
      calculateAmp(W, H, FAmp2[0][0], FRe2[0][0], FIm2[0][0]);
      draw(0, 2 * H + ytrans, W, H, FRe2[0][0], 1, 1);  draw(W, 2 * H + ytrans, W, H, FIm2[0][0], 1, 1);  draw(2 * W, 2 * H + ytrans, W, H, FAmp2[0][0], 1, 0);
      print("Press a-z for filters and arrows to change DC", 0, 0, RGB_White, 1, ColorRGB(128, 0, 0));
      redraw();
    }
    changed=0;
  }
  return 0;
}

FFT2D、draw 和 calculateAmp 函数与之前相同，这里不再显示。它们也包含在可下载文件中。

以下是程序可能的输出：



顶行显示原始图像，下方是其频谱（右侧为振幅），第三行显示应用低通滤波器后的频谱：只保留低频分量。底行显示对新频谱进行逆 FT 的结果：图像的所有高频分量被去除，图像变得模糊。

以下是另一个示例，显示将频谱的虚部设置为 0 时会发生什么。要实现这一效果，先按空格键，然后在程序中按 b 键。



后续章节将介绍更多不同类型的滤波器。

直流分量

频谱的中心是图像的直流分量，代表图像的平均颜色：如果对只保留直流分量的频谱计算 IDFT，得到的是一幅只有一种颜色的平坦图像，即原始图像的平均颜色。在上面给出代码的程序中，先按空格键接受图像，然后可以使用方向键添加或去除直流分量，按 a 键（或 azerty 键盘上的 q 键）恢复完整频谱，按 m 键（azerty 键盘上的 , 键）将频谱置零。要只获得直流分量，先按 m 键去除频谱，然后按向上方向键只恢复直流分量。

这是原始图像（频谱中的亮蓝线可能来自散热片的锐利边缘）：



这是去除直流分量后的同一图像（注意频谱中心的黑点）。实际上您看不到完整情况，图像中看起来是黑色的部分实际上是负值，因为这幅图像表示从原始图像中减去平均颜色后的结果。右侧是图像的振幅，揭示了负值部分。

去除直流分量后的图像去除直流分量后图像的振幅

下面只保留了直流分量。频谱中只剩一个点。在计算机屏幕上它看起来是白色的，但实际上是 3 个包含直流颜色信息的大浮点数。图像频谱的直流分量通常非常大，包含最多的"能量"：如前面的图像所示，从频谱中去除这个单点会使图像变暗很多。这里的"能量"应理解为信息量。

图像中只保留直流分量

这是图像的平均颜色，您可以称之为终极模糊——当您站在离图像很远的地方所看到的，或者将所有像素的颜色求和后除以像素总数所得到的结果。

低通滤波器：模糊

音频中的低频代表音乐的低音，而在图像中，低频代表最模糊的部分和渐变。锐利的边缘和线条具有高频。

图像可以通过对每个像素计算该像素及其邻域的加权平均值来进行模糊，使用卷积矩阵：这就是 Photoshop 和 Paint Shop Pro 的做法，例如，您可以在这些程序中使用"用户自定义"滤波器来创建这样的卷积矩阵。但是，由于空间域中的卷积等同于频域中的乘法，我们也可以通过将频谱与低通滤波器的传递函数相乘来模糊图像。低通滤波器的传递函数使所有高频分量为零或至少削弱它们，同时保留或放大低频分量。

对于非常强烈的模糊，需要一个非常大的卷积矩阵，这需要大量计算。进行快速傅里叶变换，然后乘以频谱，再进行逆 FFT，对于强模糊来说速度快得多。使用频谱还可以更简单、更直观地创建更好的模糊效果。

从频谱中去除的高频越多，图像越模糊，在空间域中所需的计算也越多：








在上述图像中，保留了频谱的一个圆形区域。如果不使用这样"锐利"的圆形，而是使用从中心开始指数衰减的高斯核，可以创建更好的模糊效果，即高斯模糊。

您也可以保留一个正方形区域，但圆形区域应该能产生更好的模糊效果：

附注：对于某些类型的模糊，例如高斯模糊，还有比 FFT 更快的算法，不需要 FFT 也不需要大卷积核。对于高斯模糊，可以通过连续进行三次对任意大小的盒子都能在线性时间内完成的盒式模糊来很好地近似，更多内容请参见本系列的滤波教程。

高通滤波器

图像上的高通滤波器只保留高频分量。由于应用高通滤波器时，频谱的直流分量也被去除，结果图像可能有无法绘制的负像素值。为了解决这个问题，Photoshop 的高通滤波器会在结果像素的每个颜色分量上加 128（最大值 255 的一半），使灰色表示零。这里，我们将查看图像的正值实部、振幅（揭示负值），以及重新加入直流分量后的图像。

高通滤波器与低通滤波器完全相反，因此现在去除频谱中心的一个圆形区域。左侧是图像的正值实部，中间是去除圆形后的频谱，右侧是图像的振幅：



为了更好地看清发生了什么，再次加入直流分量：



可以看到只剩下图像的高频部分：细微的渐变和颜色变化消失了，只有物体的边缘，尤其是散热片仍然可以辨认。

Photoshop 不是加入直流分量，而是向图像添加灰色以使所有像素为正，因此 Photoshop 中高通滤波器的结果如下所示：



如果我们对去除的圆形使用更大的半径，只有最高频率保留下来，只有散热片仍然看起来相当清晰：



高通滤波器也可以通过先对图像进行模糊，然后从原始图像中减去模糊图像来实现。

看了这些图片后，您可能会想知道高通滤波器是否真的有什么实际用途？确实有，例如它可以用来制作更好的可平铺纹理，以这个草地纹理为例：



如果平铺它，看起来是这样的：



看起来相当难看，因为图像中有明显的亮点和暗点，而且这些看起来非常重复。这些大面积的亮色和暗色区域是非常低频的分量，因此可以用一个细微的高通滤波器去除！

让我们将图像导入程序，对其应用一个半径非常小的高通滤波器，然后恢复直流分量：直流分量是所有频率中最低的，只是一个平坦的颜色，因此不会显得重复，而且非常重要，因为它包含草地的绿色。如截图右侧第三行所示，频谱中只去除了一个非常小的圆形区域，圆心是一个白色点，即直流分量。



底行出现了草地纹理的新版本，在保持一定距离观看时，平铺效果已经好多了：



您仍然可以看出边缘在哪里，因为源图像实际上并不是一个完美可平铺的图像。如果源图像是完美可平铺的，结果也会如此。对于一个美丽的岩石墙壁写实纹理，如果其中有一个大的暗斑使岩石墙壁在平铺时显得重复，这个高通滤波器的结果将非常有用。

程序的源代码经过了少量修改，以在"SDLK_j"按钮下获得半径更小的高通滤波器。

附注：对于锐化，也有比 FFT 更快的算法，例如使用快速模糊然后计算以下图像之和：original + (original - blur) * amount。这就是"反锐化蒙版"（unsharp masking）。

带通滤波器

带通滤波器只允许介于某个最低频率和最高频率之间的频率通过。您的收音机就使用其一维版本来收听特定频率的电台。同样，这种滤波器也会去除直流分量，因此我们再次将其加回，因为没有它图像太暗看不清。以下是原始图像，以及带直流分量的带通滤波器示例：





许多这样的图像之和（当然直流分量只加一次）就能还原出原始图像。以下是几个允许更高频率通过的带通滤波器：







以下是一个更窄的带通滤波器的结果，带通滤波器越窄，越能看出图像是如何由正弦函数之和构成的：



注意：程序默认不包含所有带通示例，但您可以通过修改代码中带通滤波器的数值（第二个输入循环中的 SDLK_r 和 SDLK_s 下）轻松获得它们。

保持实数性

到目前为止，我们总是分别计算每个通道的 FT，并将虚部留为黑色。这是一种很大的资源浪费，因为有 3 个未使用的虚部通道。

由于实函数的 FT 结果在复数域中是对称的，一半的信息实际上是冗余的。因此可以将 FT 的结果存储在一个实数数组中：将直流分量放在中心，一侧放实部的一半，另一侧放虚部的一半。对于偶数大小的输入，虚部还会有一个值为零的情况，其实部有某个非零值。该值可以放在剩余的索引中。这种变换是完全可逆的，但需要正确处理索引。在这个已经很长的教程中，提供相应代码目前超出了范围。

另一种避免复数的方法是使用实数变换，例如 DCT（离散余弦变换）。存在几种类似的变换，也有类似的"快速"变体，但这也超出了本文目前的讨论范围。

版权所有 (c) 2004-2018 Lode Vandevenne，保留所有权利。