图像的频域滤波（基于严格的傅里叶分析理论）

本文将基于严格的傅里叶分析理论基础给出严谨的图像频域滤波算法，并利用代码演示图像频域滤波的正确结果。

一、离散傅里叶变换
1.1 离散傅里叶变换的产生
离散信号分析和处理的主要手段是利用计算机去实现，然而序列f(n)的离散时间傅里叶变换是频率的连续函数，并且其逆变换为积分运算无法直接用计算机实现。显然，要在计算机上完成这些变换必须要将连续函数转换为离散数据，同时将求和范围从无限宽收缩至一个有限区间。对于离散傅里叶级数，无论是在时域还是在频域，只对有限项求和，故而可以利用计算机计算。因此借助于离散傅里叶级数的概念，把有限长的序列作为周期性离散序列的一个周期来处理，从而定义了离散傅里叶变换。这样，在允许一定程度近似的条件下，有限长序列的离散时间傅里叶变换可由计算机实现。亦即，可由计算机计算有限长序列的频谱。
1.2 离散傅里叶变换的物理意义
由离散傅里叶变换可知，有限长序列f(n)的N点离散傅里叶变换F(k)正好是f(n)的周期延拓序列的离散傅里叶级数系数的主值序列；又因为周期序列的频谱完全取决于其离散傅里叶级数系数，因此，F(k)实质上是f(n)的周期延拓序列的频谱特性。
1.3 小结
只有时域周期信号才有离散傅里叶级数（DFS），相应的也就有离散傅里叶级数系数。其中离散傅里叶级数系数是频域上的周期序列。而所谓离散傅里叶变换对实质上就是时域或者频域上周期序列的主值序列。例如对于时域上某有限长序列f(n)求其N点离散傅里叶变换（DFT），那么就可以看做先将f(n)以时域上的N点为周期进行周期延拓，然后求该时域周期信号的离散傅里叶级数系数（离散傅里叶级数系数是频域上的周期序列），最后取离散傅里叶级数系数的主值序列。反过来，对于频域上某有限长序列F(k)求其N点离散傅里叶反变换（IDFT），可视为先将F(k)以频域上的N点为周期进行周期延拓，然后求该频域周期信号的IDFS，得到时域周期信号，最后取该周期信号的主值序列。

二、关于循环卷积与线性卷积
基于以上小结关于离散傅里叶变换的论述（周期延拓—>DFS或IDFS—>取主值序列），此时计算时域中两个有限序列的线性卷积的L点DFT等价于计算两个有限序列循环卷积的L点DFT，同时这也解释清楚了为什么离散傅里叶变换的性质全部都是基于循环卷积，而不是线性卷积；而这也正是我们探讨线性卷积与循环卷积的等价条件的原因。
Note：在这里严格推导本小节字体加粗部分
要证两个有限序列的线性卷积的L点DFT等价于两个有限序列循环卷积的L点DFT，首先就要证明下述等式

证明过程如下：

紧接着等式两边同时取主值序列，得：
等式右边代表的是x(n)与f(n)的L点循环卷积，上式的意义为：x(n)与f(n)的L点循环卷积等于x(n)与f(n)的线性卷积以L为周期进行周期延拓后的主值序列。
设序列x(n)的长度为M，序列f(n)的长度为N，那么线性卷积x(n)*f(n)的长度为N+M-1。显然，若L<N+M-1，则上式左侧会出现时域混叠现象，此时x(n)与y(n)的L点循环卷积无法正确的反映x(n)与f(n)的线性卷积结果。因此容易得到x(n)与f(n)的L点循环卷积等于x(n)与f(n)的线性卷积的条件是L≥N+M-1。

更进一步的，对式
两边同时计算DFS，然后取主值序列得：
根据DFT与DFS的关系，于是有
这就是说：
时域序列的线性卷积的L点离散傅里叶变换等价于时域序列的循环卷积的L点离散傅里叶变换。
频域与时域也是同样的道理，设序列x(n)的长度为M，序列f(n)的长度为N，那么线性卷积x(n)*f(n)为有限长序列，其长度为N+M-1。那么对线性卷积x(n)*f(n)进行离散傅里叶变换，其有效的频谱长度也为N+M-1。
对于上式，根据DFT与DFS关系得，线性卷积x(n)*f(n)的L点离散傅里叶变换等于将线性卷积x(n)*f(n)的结果以L为周期进行周期延拓，计算其离散傅里叶级数系数，然后取其长为L的主值序列，得到的频谱长度为L，而线性卷积x(n)*f(n)的有效频谱长度为N+M-1，若L<N+M-1，则会出现频谱混叠的现象，那么x(n)与f(n)循环卷积的DFT无法正确计算线性卷积x(n)*f(n)的频谱。因此x(n)与f(n)循环卷积的DFT与线性卷积x(n)*f(n)等价的条件为L≥N+M-1。
综上，x(n)与f(n)循环卷积与线性卷积x(n)*f(n)等价的条件为L≥N+M-1。

三、图像的频域滤波算法步骤
图像进行频域滤波时同样需要考虑循环卷积与线性卷积的等价条件。设图像的尺寸为M·N，频域滤波器的尺寸为m·n，那么为避免时域混叠，需要将图像以及频域滤波器扩充至(M+m-1)·(N+n-1)的尺寸。
因此频域滤波一般基于以下步骤：
1、给定一幅大小为M·N的输入图像f(x,y)，需要对其添加必要数量的0，扩充成尺寸为(2M)·(2N)的图像fp(x,y)。
2、用(-1)^(x+y)乘以fp(x,y)，将其频谱移至中心。
3、计算上一步的DFT，得到F(u,v)。
4、构建一个实对称的频域滤波函数H(u，v)，尺寸为(2M)·(2N)，中心位于(M，N)处。
5、阵列相乘得到G(u，v)=F(u,v)H(u，v)
6、进行IDFT，得到处理后的图像：Re{IDFT[G(u，v)]}·(-1)^(x+y)；随后再进行图像的裁剪。
其中，为忽略由于计算不准确导致的寄生复成分，选择了实部；另外，由于在算法最开始时对图像进行了中心化，因此需要额外乘以(-1)^(x+y)将其补偿。

四、图像在频域上的高斯高通滤波
基于以上算法步骤，将图像进行高斯高通滤波，结果如下所示：
本例使用泰勒斯威夫特的图像作为案例，将输入图像灰度化，随后将图像扩充，其结果如下所示：

将扩充后的图像进行DFT，其幅度谱如下所示：

随后构建高斯高通滤波器，其频率响应如下所示：

将输入图像的频谱与该滤波器的频谱对应相乘，结果如下所示：

显然低频成分被衰减。
将上述频谱进行IDFT，取其实部并乘以(-1)^(x+y)，得到的结果如下所示：

显然图像的边缘很清楚的提取了出来。
下面展示在使用同一个高斯高通滤波函数的情况下，未严格按照上述算法流程（图像没有扩充足够多的零像素点）得到的高通滤波图像，结果如下所示：

很显然，图像存在严重的空域混叠现象，无法得到精确的高斯高通滤波图像。当然如果最后对图像的精度要求不高，只是大致的对图像做滤波处理，那么便无需严格按照上述算法流程进行滤波。
最后，完整代码如下所示：

#include<iostream>
#include