使用 matlab 数字图像处理（七）—— 频率域处理

二维离散傅里叶变换

$$F(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-i2\pi(ux/M+vy/N)}\\ f(x,y)=\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F(u,v)e^{i2\pi(ux/M+vy/N)}$$

则频域原点的傅里叶变换 $F(u,v)$ 为：

$$F(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)$$

显然这是 $f(x,y)$ 各个像素的灰度之和，而如果将系数 $1/MN$ 置于正变换之前，则 $F(0,0)$ 对应于原图像 $f(x,y)$ 的平均灰度。$F(0,0)$ 有时也被称为频率谱的直流分量 DC（direct current）。

FFT（快速傅里叶变换）的必要性

FFT 相对 DFT 不是一种新的变换，而是对 DFT 的快速计算算法。

之所以提出快速傅里叶变换（FFT）方法，是因为在计算离散域上的傅里叶变换时，对于 $N$ 点序列，它的 DFT 变换与反变换对定义为：

$$F(u)=\sum_{x=0}^{N-1}f(x)W_N^{ux},\;u=0,1\ldots,N-1,\\ f(x)=\frac1N\sum_{u=0}^{N-1}F(u)W_N^{-ux},\;x=0,1\ldots,N-1$$
其中 $W_N=e^{-j\frac{2\pi}N}$，

于是不难发现，计算每个 $u$ 值对应的 $F(u)$ 需要 $N$ 次复数乘法和 $N-1$ 复数加法，因此为了计算长度为 $N$ 的序列的傅里叶变换，共需要执行 $N^2$ 次复数加法和 $N(N-1)$ 次复数加法。

事实上，在离散傅里叶变换的运算中有大量重复运算，变量 $W_N$ 是一个复变量，实际上 $W_N$ 只有 $N$ 个独立的值。而这 $N$ 个值也不是完全独立，它们具有一定的对称关系。关于变量 $W_N$ 的周期性和对称性，可以做以下总结：

$$W_N^0=1,\quad W_N^{\frac N2}=-1\\ W_N^{N+r}=W_N^r,\quad W_N^{\frac N2+r}=-W_N^r$$

上面的等式是 $W$ 矩阵中元素的某些特殊值，而下面的等式则说明了 $W$ 矩阵元素的周期性和对称性。利用 $W_N$ 的周期性，DFT 运算中的某些项就可以合并，而利用 $W_N$ 的对称性，则可以仅计算半个 $W_N$ 序列。而根据这两点，就可以将一个长度为 $N$ 的序列分解成两个长度为 $N/2$ 的序列并分别计算 DFT，这样就可以节省大量的运算量。

N = 5;
W_N = exp(-j*2*pi/N);
W_N.^(0:9)

输出为：

1.0000 + 0.0000i   0.3090 - 0.9511i  -0.8090 - 0.5878i  -0.8090 + 0.5878i   0.3090 + 0.9511i

1.0000 + 0.0000i   0.3090 - 0.9511i  -0.8090 - 0.5878i  -0.8090 + 0.5878i   0.3090 + 0.9511i

频域滤波的基本步骤

• （1）计算原始图像 $f(x,y)$ 的 DFT，得到 $F(u,v)$
• （2）将频谱 $F(u,v)$ 的零频点移动到频谱图的中心位置
• （3）计算滤波器函数 $H(u,v)$ 与 $F(u,v)$ 的乘积 $G(u,v)$（时域的卷积等于频域的乘积）
• （4）将频谱 $G(u,v)$ 的零频点移回到频谱图的左上角位置
• （5）计算机步骤（4）的傅里叶反变换 g(x, y)
• （6）取 g(x, y) 的实部作为最终滤波后的结果图像；

由上面的叙述易知，滤波能否取得成功的关键取决于频域滤波函数 $H(u,v)$，常常称之为滤波器，或滤波器传递函数，因为它在滤波中抑制或滤除了频谱中某些频率的分量，而保留其他的一些频率不受影响。

如果滤波器的值只为实数，则进行频域滤波的时候，$F\cdot H$（频域的乘积），$H$ 中的每一个实数元素分别乘以 $F$ 中对应位置的复数元素，从而使 $F$ 中元素的实部和虚部等比例的变化，不会改变 $F$ 的相位谱，这种滤波器因此也被称为“零相移”滤波器。