数字图像处理第四章：频率域滤波

1.背景

傅立叶的基本思想就是：任何简单或复杂的周期函数都可表示为不同频率的正弦/余弦加权和的形式。因此，在图像处理中，图像从时域角度处理往往不容易，通过傅立叶变换形成频率域再做处理就方便很多。由上章可知，有可分离的和不可分离的两种核，分别考虑大小为MxM和mxm的图像与核，在可分离核中，采用快速傅立叶变化算法(FFT)的滤波优势定义为:$$C_n(m)=\frac{M^2{m}^2}{2M^2{\log }_2{M}^2}=\frac{m^2}{4{\log }_{2}M}$$可分离的优势定义为：$$C_s(m)=\frac{2M^{2}m}{2M^2{\log }_2{M}^2}=\frac{m}{2{\log }_{2}M}$$当C(m)>1时，FFT优势更大，且FFT优势随m增大迅增，反之则是空间滤波优势大

2.基本概念

复数

就是高中学过的那个负数没错，公式定义为（R表实部，I表虚部）：$$C=R+jI$$

傅立叶级数

可由正弦余弦组合形成的周期为T的函数f(t)，其形式为：
$$f(t)=\sum _{n=-\infty }^{\infty }{c}_n{e}^{j\frac{2\pi n}{T}t}$$$$c_n=\frac{1}{T}{\int }_{-T/2}^{T/2}f(t)e^{-j\frac{2\pi n}{T}t}dt,n=0,\pm 1,\pm \mathrm{2...}$$

冲激函数及其取样性质

冲激函数是一个理想化的信号，通常用符号 $\delta (t)$ 或 $\delta [n]$ 表示。它在 $t=0$ 或 $n=0$ 时取值为无穷大，而在其他时间或时间点处均为零。冲激函数的积分等于1或求和等于1，即：$${\int }_{-\infty }^{\infty }\delta (t)dt=1$$$$\sum _{n=-\infty }^{\infty }\delta [n]=1$$具体来说，将一个信号与冲激函数相乘，并对其进行积分或求和，可以得到信号在冲激函数的时间点上的采样值。这个过程被称为取样，冲激具有关于积分的所谓取样性质：$${\int }_{-\infty }^{\infty }f(t)\delta (t)dt=f(0)$$取样性质只是得到$f(t)$在冲激位置的值，冲激串$S_{\Delta t}(t)$定义为无穷多个冲激$\Delta T$单位的和:$$S_{\Delta T}(t)=\sum _{k=-\infty }^{\infty }\delta (t-k\Delta T)$$
冲激函数的取样性质可以描述为：

1. 取样定理：一个信号如果没有高于采样频率一半的频率分量，那么它可以通过在等间隔的时间点上采样，并使用插值函数对采样值进行重构，得到与原信号完全相同的信号。
2. 周期性：冲激函数的离散版本是周期性的，即 $\delta [n]=\delta [n+kN]$，其中 $N$ 是一个正整数，$k$ 是任意整数。
3. 平移性：如果将冲激函数向右平移 $n_0$ 个单位，则其离散版本的表达式为 $\delta [n-n_0]$。
4. 放大性：如果将冲激函数的幅度放大 $a$ 倍，则其离散版本的表达式为 $a\delta [n]$。
   这些性质使冲激函数成为数字信号处理中非常有用的工具。

连续单变量函数的傅立叶变换

傅立叶变换可定义为 ξ { f ( t ) } = ∫ − ∞ ∞ f ( t ) e − j 2 π μ t d t \xi\{f(t)\}=\int_{-\infty}^{\infty}f(t)e^{-j2\pi\mu t}dt ξ{f(t)}=∫−∞∞​f(t)e−j2πμtdt由于$t$和$\mu$都是连续变量，上式可以改写为以下两个式子：$$F(\mu )={\int }_{-\infty }^{\infty }f(t)e^{-j2\pi \mu t}dt,f(t)={\int }_{-\infty }^{\infty }F(\mu )e^{j2\pi \mu t}d\mu$$上面两个式子共同构成傅立叶变换对，使用欧拉公式，有如下公式：$$F(\mu )={\int }_{-\infty }^{\infty }f(t)\left[\cos (2\pi \mu t)-j\sin (2\pi \mu t)\right]$$该方程积分变量是t，剩下的唯一变量$\mu$是正余弦函数的频率，因此傅立叶变换域是频率域，频率变量$\mu$取决于t的单位。或者说频率域的单位是输入函数自变量每秒的周期数

以下是一个“盒式”函数及对应的傅立叶变换和频谱图：

卷积

将两个函数（或信号）通过一定的规则组合在一起，产生一个新的函数（或信号）。通俗一点，可以将卷积理解为一种“混合”或“融合”过程主要有以下步骤：

1. 取两个函数（或信号），其中一个通常被称为输入函数（或信号），另一个被称为卷积核（或滤波器）。
2. 对于输入函数中的每个时间点（或位置），将卷积核与输入函数的对应部分进行相乘。
3. 将乘积的结果在时间（或空间）上进行平移，并对相邻的乘积结果进行求和。
4. 将所有求和的结果组合起来，形成一个新的函数（或信号），称为卷积结果。

卷积定理：空间域中两个函数乘积的傅立叶变换，是俩函数变换在频率域中的卷积

3.取样和取样函数的傅立叶变换

连续函数转换为离散值需要取样和量化

取样

就是通过周期为$\Delta T$时间的冲激串将原函数进行点的划分，形成离散的数值量，以下是一个函数和取样的结果图

取样后的傅立叶变换

令原函数为$f(t)$，经过傅立叶变换后为$F(\mu )$，则取样后的函数$\tilde{f}(t)$的傅立叶变换$\tilde{F}(\mu )$是$$\tilde{F}(\mu )=\frac{1}{\Delta T}\sum _{n=-\infty }^{\infty }F\left(\mu -\frac{n}{\Delta T}\right)$$上式表明$\tilde{F}(\mu )$是$F(\mu )$的一个无线的、周期的副本序列，间隔由1/$\Delta T$的值决定，且$\tilde{F}(\mu )$依然是连续函数

取样定理与混叠

原点为中心的有限区间$[-{\mu }_{max},{\mu }_{max}]$外的频率值，其FFT为零的函数$f(t)$叫带限函数
取样定理指出，为了准确地重构一个连续时间信号，我们需要以至少两倍的采样率对该信号进行采样。以下左边的子图显示原始信号的图形，右边的子图显示采样信号的图形：

具体来说，采样频率应该大于等于信号中最高频率的两倍，等于时的取样率称为奈奎斯特律（Nyquist定理）。这是因为信号中的频率成分可以通过采样进行捕捉，但是如果采样频率低于信号中最高频率的两倍，就会出现混叠现象，即高于采样频率一半的频率成分被误认为低于采样频率一半的频率导致取样后的信号不能与其他信号区分开来，而混叠对是指具有相同数字化结果的两个不同连续函数。

低通滤波器指函数$H(\mu )$，该函数定义为：$$H(\mu )=\{\begin{array}{ll}\Delta T,& -{\mu }_{max}\le \mu \le {\mu }_{max}\\ 0,& 其他\end{array}$$

以下左边显示原始信号的图形，右边的子图显示有混叠问题的采样信号图形：

上述两幅图的代码如下：

import numpy as np
import matplotlib.pyplot as plt

# 定义原始信号
t = np.linspace(0, 1, 1000)  # 连续时间
f = 10  # 原始信号频率
signal = np.sin(2 * np.pi * f * t)  # 原始信号

# 进行采样
Fs = 40  # 正常的采样频率
# Fs = 2 * f  # 采样频率设置为原始信号频率的一半，有混叠问题
Ts = 1 / Fs  # 采样周期
samples = np.arange(0, 1, Ts)  # 采样时间点
sampled_signal = np.sin(2 * np.pi * f * samples)  # 采样信号

# 绘制原始信号和采样信号
plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.plot(t, signal)
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Original Signal')

plt.subplot(122)
plt.stem(samples, sampled_signal)
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Sampled Signal')

plt.tight_layout()
plt.show()

重建函数

就是由一组样本重建取样过后的函数，由$F(\mu )=H(\mu )\tilde{F}(\mu )$可知 f ( t ) = ξ { F ( μ ) } = ξ { H ( μ ) F ~ ( μ ) } = h ( t ) ∗ f ~ ( t ) f(t)=\xi\{F(\mu)\}=\xi\{H(\mu)\tilde{F}(\mu)\}=h(t)*\tilde{f}(t) f(t)=ξ{F(μ)}=ξ{H(μ)F~(μ)}=h(t)∗f~​(t)最后可以得到如下表达式：$$f(t)=\sum _{n=-\infty }^{\infty }f(n\Delta T)sinc[(t-n\Delta T)/\Delta T]$$由于$sinc$函数是$H(\mu )$的傅立叶反变换函数，就有了以上的公式

4.单变量的离散傅立叶变换

由取样函数变换的DFT

离散傅立叶变换（DFT）：对$\tilde{F}(\mu )$的一个周期进行取样是DFT的基础

若一个函数f是另一个函数g的DFT，则函数g是函数f的IDFT，这两函数就是傅立叶变换对

循环卷积公式如下所示：$$F(x)\ast h(x)=\sum _{m=0}^{M-1}f(m)h(x-m),x=0,1,\mathrm{2....},M-1$$之所以叫这个名这是由于函数是周期的，它们的卷积也是周期的，这也是由DFT及其反变换的周期性导致的结论

取样和频率间隔的关系

若以$\Delta T$各单位间隔对$f(t)$取样后的$f(x)$由M个样本组成，则记录的长度是$$T=M\Delta T$$对应间隔$\Delta u$为$$\Delta u=\frac{1}{M\Delta T}=\frac{1}{T}$$DFT的$M$个分量跨越的整个频率范围是$$R=M\Delta u=\frac{1}{\Delta T}$$可知DFT的频率分辨率$\Delta u$与记录的长度T成反比，DFT的频率范围取决于取样间隔$\Delta T$，已知连续函数$f(t)$经$\Delta T$间隔取样后可以计算出其傅立叶变换$F(u)$，反之亦然

5.二变量函数的傅立叶变换

二维冲激及其取样性质

同一维的冲击类似，二维冲激函数比一维多了一次积分，定义为$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\delta (t,z)dtdz=1$$取样性质为：$$\sum _{x=-\infty }^{\infty }\sum _{y=-\infty }^{\infty }f(x,y)\delta (x,y)=f(0,0)$$对$(x_0,y_0)$出的一个冲激，取样性质为：$$\sum _{x=-\infty }^{\infty }\sum _{y=-\infty }^{\infty }f(x,y)\delta (x-x_0,y-y_0)=f(x_0,y_0)$$

二维连续傅立叶变换对

以下是一个二位函数和其变换后的频率谱：傅立叶变换的算法公式如下：$$F(u,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(t,z)e^{-2iw(ut+vz)}dtdz$$当 u,v=0 时，由$\sin c(m)=\sin (\pi \ast m)/\pi \ast m$,因此sinc(0)=1，所以在u,v=0时取得最大幅值。直流成分，越均匀无变化的部分，值就越大

二维傅立叶的一般公式如下，其中M和N分别是上图中的长和宽：$$F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-2iw(ux/M+vy/N)}$$

二维取样及取样定理

类似于一维的取样，二维取样也可以用一个取样函数$S_{\Delta T\Delta Z}(t,z)$建模，并且将其乘以原二维图像函数$f(t,z)$就得到了取样后的函数；当函数$f(t,z)$的傅立叶变换为零时，该函数称为带限函数。
二维取样定理：若一个带限连续二维函数在$\mu$和$v$两方向上可由大于该函数最高频率2倍的取样率样本表示，则信息无丢失

图像中的混叠

同一维的混叠类似，在对使函数存在于有限的时间范围内的处理时，会引入破坏图像频谱的频率分量，故无法无限取样，也就会出现混叠。
导致图像的混叠主要有两种形式：空间混叠和时间混叠

• 时间混叠：如“车轮”效应，在电影中车在高速行驶时车轮辐条的倒转现象
• 空间混叠：由欠取样造成，使产生的信号频率远低于原信号频率，会产生伪影现象，例如出现原图没有的锯齿（边缘失真）、高光等。

下图是一个由二维信号频率与数字化该信号的取样率相互作用产生的混叠效应：

可以看到点块和条形块重叠的部分产生了与原图像无关的图像，产生了混叠效应

如何可以减少混叠效应呢？
答：通过稍微散焦待数字化的图像来衰减高频分量，即在图像数字化前进行抗混叠滤波
例如，图像的缩放总是会引入混叠，因此在取样前通过合适的滤波器平滑，再进行放大缩小就可以减少混叠的影响

莫尔效应：由近似间隔的两个光栅叠加产生的视觉现象，如网格纱窗、条纹衣服等

DFT及IDFT

二维离散傅立叶变换及其反变换公式如下：$$F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-2i\pi (ux/M+vy/N)}$$$$f(x,y)=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{2i\pi (ux/M+vy/N)}$$其中$f(x,y)$是大小为M*N的图像，对离散变量$u,v$进行间隔点求值得到$F(u,v)$，反之同理

6.二维离散傅立叶变换和反变换的性质

主要有平移旋转、周期、对称性

空间间隔和频率间隔的关系

一维的取样和频率间隔解释同4中的相同，连续函数$f(t,z)$的数字图像$f(x,y)$由$t,z$方向的M*N个样本组成，$\Delta T$和KaTeX parse error: Undefined control sequence: \Delat at position 1: \̲D̲e̲l̲a̲t̲ ̲Z表示样本间隔，则频率间隔分别为$$\Delta u=\frac{1}{M\Delta T}，\Delta v=\frac{1}{N\Delta Z}$$频率域样本间隔与空间样本间隔成反比

平移和旋转

即DFT和IDFT的平移与旋转都是相互对应的，有以下两个公式$$f(x,y)e^{2i\pi (u_{0}x/M+v_{0}y/N)}\iff F(u-u_0,v-v_0)$$$$f(x-x_0,y-y_0)\iff F(u,v)e^{-2i\pi (x_{0}u/M+y_{0}v/N)}$$说明将变换对的其中一个函数乘以对应的指数项就可将另一函数移至某点，平移不影响$F(u,v)$幅度，极坐标也类似，只是将平移换成了旋转一定的角度

周期性

由一维的知识可知二维的傅立叶变换和反变换在$u,v$方向是无限周期的。一维情况下，通过周期性，我们能够让傅立叶变换中心化，方便对DFT的操作，如$$f(x)e^{i2\pi (u_{0}x/M)}\iff F(u-u_0)$$将该指数项乘以$f(x)$把变换的原点$F(0)$移到了$u_0$，若$u_0=M/2$则指数项就变成了$e^{i\pi x}$此时有$$f(x)(-1{)}^x\iff F(u-M/2)$$这也是使得$F(u)$在区间内居中，二维时原理相同，式子如下$$f(x,y)(-1{)}^{x+y}\iff F(u-M/2,v-N/2)$$这个式子也就将原点$F(0,0)$移到了$(floor(M/2),floor(N/2)$处

对称性

偶函数是对称的，奇函数是反对称的，同初等数学知识一样，对于任意两个离散的偶函数$w_e$和奇函数$w_o$有$$\sum _{x=0}^{M-1}\sum _{y=0}^{M-1}{w}_e(x,y)w_o(x,y)=0$$
通过上面的概念，可以产生DFT及其反变换的许多性质，其中实函数$f(x,y)$的傅立叶变换是共轭对称的这一性质很重要，即$$F^{\ast }(u,v)=F(-u,-v)$$

傅立叶频谱和相角

二维DFT通常是复函数，其相关的公式如下表：

公式名称	公式
极坐标	$$F(u,v)=R(u,v)+iI(u,v)=\mid F(u,v)\mid e^{i\varphi (u,v)}$$
幅度(傅立叶频谱)	$$\mid F(u,x)\mid =[R^2(u,v)+I^2(u,v){]}^{1/2}$$
相角(相位谱)	$$\varphi (u,v)=\arctan \left[\frac{I(u,v)}{R(u,v)}\right]$$
功率谱	$$P(u,v)=\mid F(u,v){\mid }^2=R^2(u,v)+I^2(u,v)$$

$R(u,v)$和$I(u,v)$代表实部和虚部，以上都是针对离散变量而言的，由DFT的共轭对称可知谱关于原点偶对称，相角关于原点奇对称，由公式可知当(u,v)=0时，(0,0)处的指数项为1，有$$F(0,0)=MN\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)\Rightarrow \mid F(0,0)\mid =MN\mid \overline{f}\mid$$由于$MN$通常很大，因此$\mid F(0,0)\mid$通常是频谱的最大分量，有时也被称为变换的直流分量，而直流项是一个信号或波形中的恒定分量或平均值部分，它表示信号的偏移或平均能量级。

由书中对于矩形函数的频谱分析可知，平移对于频谱的影响几乎没有，但是频谱的旋转与原图的旋转一致

二维离散卷积定理

二维循环卷积表达式为$$(f\ast h)(x,y)=\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}f(m,n)h(x-m,y-n)$$$x,y$为离散变量，而二维卷积定理为$$(f\ast h)(x,y)\iff (F\cdot H)(u,v)$$或者反写为$$(f\cdot h)(x,y)\iff \frac{1}{MN}(F\ast H)(u,v)$$
计算两个函数的空间卷积有两种：

1. 直接在空间域计算
2. 计算每个函数的傅立叶变换，让两个变换相乘并计算反变换

一维中，在DFT表达式中固有的周期性$F(u)=F(u+kM)$和$f(x)=f(x+kM)$会使两个函数的处理过程中各个周期相互干扰，产生交叠错误。交叠错误可以通过填充零来纠正，如函数$f(x)$和$h(x)$，分别由A和B个样本组成，在这两函数中填充零，使它们的长度P相同，则选择$$P\ge A+B-1$$就可避免交叠错误

二维也是同理，填充零的方法如下：$$f_p(x,y)=\{\begin{array}{ll}f(x,y),& 0\le x\le A-1和0\le y\le B-1\\ 0,& A\le x\le P或B\le y\le Q\end{array}$$$$h_p(x,y)=\{\begin{array}{ll}h(x,y),& 0\le x\le C-1和0\le y\le D-1\\ 0,& C\le x\le P或D\le y\le Q\end{array}$$其中，$$P\ge A+C-1,Q\ge B+D-1$$
但是在二维两个函数零填充时，有时会导致函数不连续，类似于将函数乘个盒式函数，在频率域中这个会产生频率泄漏，会使图像块效应，这个通常用取样函数乘一个平滑过渡到零的函数就可以解决

由前面的知识可知，在两个信号之间进行二维离散卷积运算时，等效于将这两个信号的二维离散傅里叶变换（DFT）进行逐元素相乘，然后再进行二维离散傅里叶反变换（IDFT），而这种特性有下面两种应用优势：

1. 加速计算：在时域进行二维离散卷积运算可能会非常耗时，而利用二维离散卷积定理可以将计算转移到频域中，通过对频域信号进行逐元素相乘后再进行逆变换，可以显著加速卷积计算的过程。
2. 滤波效果分析：通过观察频域中的频谱信息，可以更好地理解滤波器对信号的影响。频域中的乘法可以看作是在不同频率分量上进行加权，从而实现对不同频率的滤波效果。

7.频率域滤波基础

频率域的其他性质

频率与空间变化率直接相关，其中变化最慢的频率分量$(u=v=0)$与图像的平均灰度成正比，远离变换的原点且两者成反比

频率域滤波基础知识

滤波的步骤：修改一幅图像的傅立叶变化，计算反变换，得到空间域表示的结果
已知图像$f(x,y)$则基本滤波公式为 g ( x , y ) = R e a l { ℑ − 1 [ H ( u , v ) F ( u , v ) ] } g(x,y)=Real\lbrace\Im^{-1}\lbrack H(u,v)F(u,v)\rbrack\rbrace g(x,y)=Real{ℑ−1[H(u,v)F(u,v)]}其中${\Im }^{-1}$表示IDFT，$F(u,v)$是$f(x,y)$的DFT，$H(u,v)$是滤波器传递参数（滤波器或滤波函数），$g(x,y)$是滤波后的图像，Real表示该值的实部

通常我们把衰减高频通过低频模糊图像的函数$H(u,v)$称作低通滤波器，而清晰图像但降低对比度的叫高通滤波器

以下是一些滤波器（上）和相应的滤波图像（下），分别是低通、高通、偏移高通滤波器

同时，对图像在卷积滤波前进行需要进行正确的零填充，否则可能会出现交叠错误或者导致出现不希望的结果

$H(u,v)$是滤波器传递函数填充的方法一般先对图像零填充，再在频率域中创造期望的滤波传递函数，函数大小与填充后的图像大小相同，之后可以平滑传递函数，减少交叠错误和振铃效应

频率域滤波小结

频率域滤波过程如下：

1. 从文件中读取原始图像$f(x,y)$，计算得到填充零后的尺寸
2. 进行合适的填充零操作，得到填充后的图像$f_p(x,y)$
3. 将$f_p(x,y)$乘以$(-1{)}^{x+y}$，使傅立叶变化位于矩阵中心
4. 计算图像的DFT，即$F(u,v)$，并构一个实对称滤波器传递函数$H(u,v)$
5. 用$G(u,v)=H(u,v)F(u,v)$得到$G(u,v)$
6. 计算$G(u,v)$的IDFT得到滤波后的图像$g_p(x,y)$
7. 最后从$g_p(x,y)$左上象限提取MxN的区域，得到结果$g(x,y)$

空间域和频率域滤波间的对应关系

将频率域滤波定义为滤波器传递函数$H(u,v)$与原图像的傅立叶变换$F(u,v)$的对应像素相乘，若已知$H(u,v)$则可求出其空间域的等效核$f(x,y)$，反之亦然，就有了$$h(x,y)\iff H(u,v)$$也就是说空间域滤波和频率域滤波可以通过傅里叶变换和逆变换相互转换

8.使用频率域滤波器平滑图像

主要介绍三个滤波器：理想低通滤波、巴特沃斯低通滤波器、高斯低通滤波器，其中巴特沃斯低通滤波器通过参数调整，可以在另外俩个滤波器之间进行变换。令所有滤波器传递函数大小为$P\ast Q$

理想低通滤波（ILPF）

含义：通过所有圆点为中心的圆内频率，阻止所有圆外的频率的二维低通滤波，传递函数$$H(u,v)=\{\begin{array}{ll}1,& D(u,v)\le D_0\\ 0,& D(u,v)>D_0\end{array}$$$D_0$是常数，$D(u,v)$是频率域中点到频率矩形中心的距离：$$D(u,v)=[(u-P/2{)}^2+(v-Q/2{)}^2{]}^{1/2}$$在$H(u,v)=1$与$H(u,v)=0$间的过渡点叫截止频率，直观图片如下：

相关的代码实现如下：

import numpy as np
import cv2
import matplotlib.pyplot as plt

def ideal_lowpass_filter(image, cutoff_freq):
    # 获取图像尺寸
    height, width = image.shape[:2]

    # 进行二维傅里叶变换
    dft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
    dft_shift = np.fft.fftshift(dft)

    # 构建理想低通滤波器
    filter_mask = np.zeros((height, width, 2), np.float32)
    cx, cy = width // 2, height // 2
    filter_mask[cy - cutoff_freq:cy + cutoff_freq, cx - cutoff_freq:cx + cutoff_freq] = 1

    # 将滤波器应用于频域图像
    dft_shift_filtered = dft_shift * filter_mask

    # 进行逆变换，将图像转换回空域
    dft_filtered_shifted = np.fft.ifftshift(dft_shift_filtered)
    img_filtered = cv2.idft(dft_filtered_shifted)
    img_filtered = cv2.magnitude(img_filtered[:, :, 0], img_filtered[:, :, 1])

    return img_filtered

# 读取图像
image = cv2.imread('./img/Fig0441(a)(characters_test_pattern).tif', 0)  # 以灰度模式读取

# 执行理想低通滤波
cutoff_frequency = 30  # 截止频率
filtered_image = ideal_lowpass_filter(image, cutoff_frequency)

# 显示原始图像和滤波后的图像
plt.subplot(121), plt.imshow(image, cmap='gray')
plt.title('Original Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(filtered_image, cmap='gray')
plt.title('Filtered Image'), plt.xticks([]), plt.yticks([])
plt.show()

高斯低通滤波（GLPF）

传递函数形式如下：$$H(u,v)=e^{-D^2(u,v)/2{\sigma }^2}$$其中的$\sigma$也就是上个滤波中的$D_0$，即截止频率
相关代码和结果图如下：

import numpy as np
import cv2
import matplotlib.pyplot as plt

def gaussian_lowpass_filter(image, cutoff_freq):
    # 获取图像尺寸
    height, width = image.shape[:2]

    # 进行二维傅里叶变换
    dft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
    dft_shift = np.fft.fftshift(dft)

    # 构建高斯低通滤波器
    filter_mask = np.zeros((height, width, 2), np.float32)
    cx, cy = width // 2, height // 2
    for i in range(height):
        for j in range(width):
            dist = np.sqrt((i - cy) ** 2 + (j - cx) ** 2)
            filter_mask[i, j] = np.exp(-0.5 * (dist / cutoff_freq) ** 2)

    # 将滤波器应用于频域图像
    dft_shift_filtered = dft_shift * filter_mask

    # 进行逆变换，将图像转换回空域
    dft_filtered_shifted = np.fft.ifftshift(dft_shift_filtered)
    img_filtered = cv2.idft(dft_filtered_shifted)
    img_filtered = cv2.magnitude(img_filtered[:, :, 0], img_filtered[:, :, 1])

    return img_filtered

# 读取图像
image = cv2.imread('./img/Fig0441(a)(characters_test_pattern).tif', 0)  # 以灰度模式读取

# 执行高斯低通滤波
cutoff_frequency = 30  # 截止频率
filtered_image = gaussian_lowpass_filter(image, cutoff_frequency)

# 显示原始图像和滤波后的图像
plt.subplot(121), plt.imshow(image, cmap='gray')
plt.title('Original Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(filtered_image, cmap='gray')
plt.title('Filtered Image'), plt.xticks([]), plt.yticks([])
plt.show()

虽然图片上可能没有体现出来，但是一般情况下GLPF与ILPF相比，实现的平滑要小一点，平滑过渡随着截止频率的增大而增加

巴特沃斯低通滤波器（BLPF）

其传递函数形式如下：$$H(u,v)=\frac{1}{1+[D(u,v)/D_0{]}^{2n}}$$

代码和结果图如下：

import numpy as np
import cv2
import matplotlib.pyplot as plt

def butterworth_lowpass_filter(image, cutoff_freq, order):
    # 获取图像尺寸
    height, width = image.shape[:2]

    # 进行二维傅里叶变换
    dft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT)
    dft_shift = np.fft.fftshift(dft)

    # 构建巴特沃斯低通滤波器
    filter_mask = np.zeros((height, width, 2), np.float32)
    cx, cy = width // 2, height // 2
    for i in range(height):
        for j in range(width):
            dist = np.sqrt((i - cy) ** 2 + (j - cx) ** 2)
            filter_mask[i, j] = 1 / (1 + (dist / cutoff_freq) ** (2 * order))

    # 将滤波器应用于频域图像
    dft_shift_filtered = dft_shift * filter_mask

    # 进行逆变换，将图像转换回空域
    dft_filtered_shifted = np.fft.ifftshift(dft_shift_filtered)
    img_filtered = cv2.idft(dft_filtered_shifted)
    img_filtered = cv2.magnitude(img_filtered[:, :, 0], img_filtered[:, :, 1])

    return img_filtered

# 读取图像
image = cv2.imread('./img/Fig0441(a)(characters_test_pattern).tif', 0)  # 以灰度模式读取

# 执行巴特沃斯低通滤波
cutoff_frequency = 30  # 截止频率
order = 2  # 阶数
filtered_image = butterworth_lowpass_filter(image, cutoff_frequency, order)

# 显示原始图像和滤波后的图像
plt.subplot(121), plt.imshow(image, cmap='gray')
plt.title('Original Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(filtered_image, cmap='gray')
plt.title('Filtered Image'), plt.xticks([]), plt.yticks([])
plt.show()

与前两个滤波进行比较之后，其实还是能看出来BLPF的平滑效果是介于ILPF和GLPF之间的，空间一阶BLPF没有振铃效应，二三阶的也很难察觉，但是在更高阶中就很明显

9.使用高通滤波器锐化图像

上一节中衰减图像的傅立叶变换中的高频分量可以平滑分量，而要锐化图像，则需要衰减低频分量。同上结类似，令所有滤波器传递函数大小为$P\ast Q$，滤波器传递函数$H(u,v)$视为大小为$P\ast Q$的中心化离散函数

由低通滤波得到对应的理想、高斯、巴特沃斯高通滤波

在空间域中可以由低通核$lp(x,y)$通过$$hp(x,y)=\delta (x,y)-lp(x,y)$$得到高通核，在频率域中则用下列式子也可以得到对应的高通滤波器传递函数：$$H_{HP}(u,v)=1-H_{LP}(u,v)$$由上式就可以得到对应的理想、高斯、巴特沃斯高通滤波传递函数通过高通滤波器实现的锐化效果中，巴特沃斯高通滤波的尖锐过渡效果在理想和高斯之间，阈值处理可突出尖锐特征

频率域中的拉普拉斯

拉普拉斯在频率域和频率矩形中心的传递函数$H(u,v)$如下：$$H(u,v)=-4{\pi }^2(u^2+v^2)$$$$H(u,v)=-4{\pi }^2[(u-P/2{)}^2+(v-Q/2{)}^2]=-4{\pi }^2{D}^2(u,v)$$
$D(u,v)$表示距离，用这个传递函数可以得到图像的拉普拉斯：$${\nabla }^{2}f(x,y)={\Im }^{-1}[H(u,vF(u,v))]$$$F(u,v)$表示原图像$f(x,y)$的$DFT$,增强的式子如下：$$g(x,y)=f(x,y)+c{\nabla }^{2}f(x,y)$$其中c=-1，代码和效果图在上章已有，不再赘述

钝化掩蔽、高提升滤波和高频强调滤波

频率域中定义的模版为：$$g_{mask}(x,y)=f(x,y)-f_{LP}(x,y)$$$$f_{LP}(x,y)={\Im }^{-1}[H_{HP}(u,v)F(u,v)]$$其中$f_{LP}(x,y)$是平滑后的图像$H_{HP}(u,v)$是低通滤波器传递函数。$$g(x,y)=f(x,y)+k{g}_{mask}(x,y)$$式中定义了k=1时的钝化掩蔽以及k>1时的高提升滤波，高频强调滤波的通用公式为$$g(x,y)={\Im }^{-1}[k_1+k_2{H}_{HP}(u,v)]F(u,v)$$

同态滤波

同态滤波的作用是它能够有效地改善低光照、强光照和阴影等情况下的图像质量

已知图像$f(x,y)$可由照射分量$i(x,y)$和反射分量$r(x,y)$的乘积表示，即$$f(x,y)=r(x,y)r(x,y)$$其反变换公式如下$$\Im [z(x,y)]=\Im [\ln f(x,y)]=\Im [\ln i(x,y)]+\Im [\ln r(x,y)]$$或$$Z(u,v)=F_i(u,v)+F_r(u,v)$$其中$F_i(u,v)$和$F_r(u,v)$分别是$\ln i(x,y)$和$\ln r(x,y)$的傅立叶变换，滤波后的图像是$$s(x,y)={\Im }^{-1}[S(u,v)]={\Im }^{-1}[H(u,v)F_i(u,v)]+{\Im }^{-1}[H(u,v)F_r(u,v)]$$令$i^{{}^{\prime }}(x,y)={\Im }^{-1}[H(u,v)F_i(u,v)]$，$r^{{}^{\prime }}(x,y)={\Im }^{-1}[H(u,v)F_r(u,v)]$则有$$s(x,y)=i^{{}^{\prime }}(x,y)+r^{{}^{\prime }}(x,y)$$输出的图像就是$$g(x,y)=e^{s(x,y)}=e^{i^{{}^{\prime }}(x,y)}{e}^{r^{{}^{\prime }}(x,y)}=i_0(x,y)r_0(x,y)$$上面的方法是同态系统的一种特殊情况，整个同态滤波步骤简要概括如下

10.选择性滤波

处理特定频带的滤波器叫频带滤波器，有带阻滤波器和带通滤波器
处理小频率矩形区域的滤波器叫陷波滤波器，有陷波带阻滤波器和陷波带通滤波器

带阻滤波器和带通滤波器

由前面知道，带阻带通滤波器传递函数可由高低通滤波器传递函数组合构建，下面不同带阻滤波器传递函数$H(u,v)$的对应公式

理想带阻滤波器(IBRF)	高斯带阻滤波器(GBRF)	巴特沃斯带阻滤波器(BBRF)
$H(u,v)=\{\begin{array}{ll}0,& C_0-\frac{W}{2}\le D(u,v)\le C_0+\frac{W}{2}\\ 1,& 其他\end{array}$	$H(u,v)=1-e^{-[\frac{D^2(u,v)-C_0^2}{D(u,v)W}{]}^2}$	$H(u,v)=\frac{1}{1+[\frac{D(u,v)W}{D^2(u,v)-C_0^2}{]}^{2n}}$

对应传递函数图像如下：

陷波滤波器

陷波滤波器是最有用的选择性滤波器，是一种可以在某一个频率点迅速衰减输入信号，以达到阻碍此频率信号通过的滤波效果的滤波器，陷波带阻滤波器传递函数一般公式为$$H_{NR}(u,v)=\prod ^Q{H}_k(u,v)H_{-k}(u,v)$$其中$H_k(u,v)$和$H_{-k}(u,v)$分别是中心在$(u_k,v_k)$和$(-u_k,-v_k)$的高通滤波器传递函数，陷波带通滤波器传递函数则是$$H_{NP}(u,v)=1-H_{NR}(u,v)$$陷波滤波能选择性的修改DFT的局部区域，其应用的范围如删除数字化印刷物中的莫尔模式，以下是代码和结果图

import matplotlib.pyplot as plt
import numpy as np
import cv2

def gaussLowPassFilter(shape, radius=10):  # 高斯低通滤波器
    u, v = np.mgrid[-1:1:2.0 / shape[0], -1:1:2.0 / shape[1]]
    D = np.sqrt(u ** 2 + v ** 2)
    D0 = radius / shape[0]
    kernel = np.exp(- (D ** 2) / (2 * D0 ** 2))
    return kernel


def butterworthNRFilter(shape, radius=9, uk=60, vk=80, n=2):  # 巴特沃斯陷波带阻滤波器
    M, N = shape[1], shape[0]
    u, v = np.meshgrid(np.arange(M), np.arange(N))
    Dm = np.sqrt((u - M // 2 - uk) ** 2 + (v - N // 2 - vk) ** 2)
    Dp = np.sqrt((u - M // 2 + uk) ** 2 + (v - N // 2 + vk) ** 2)
    D0 = radius
    n2 = n * 2
    kernel = (1 / (1 + (D0 / (Dm + 1e-6)) ** n2)) * (1 / (1 + (D0 / (Dp + 1e-6)) ** n2))
    return kernel


def imgFrequencyFilter(img, lpTyper="GaussLP", radius=10):
    normalize = lambda x: (x - x.min()) / (x.max() - x.min() + 1e-8)

    # (1) 边缘填充
    imgPad = np.pad(img, ((0, img.shape[0]), (0, img.shape[1])), mode="reflect")
    # (2) 中心化: f(x,y) * (-1)^(x+y)
    mask = np.ones(imgPad.shape)
    mask[1::2, ::2] = -1
    mask[::2, 1::2] = -1
    imgPadCen = imgPad * mask
    # (3) 傅里叶变换
    fft = np.fft.fft2(imgPadCen)

    # (4) 构建 频域滤波器传递函数:
    if lpTyper == "GaussLP":
        print(lpTyper)
        freFilter = gaussLowPassFilter(imgPad.shape, radius=60)
    elif lpTyper == "GaussHP":
        freFilter = gaussLowPassFilter(imgPad.shape, radius=60)
    elif lpTyper == "ButterworthNR":
        print(lpTyper)
        freFilter = butterworthNRFilter(imgPad.shape, radius=9, uk=60, vk=80, n=2)  # 巴特沃斯陷波带阻滤波器
    elif lpTyper == "MButterworthNR":
        print(lpTyper)
        BNRF1 = butterworthNRFilter(imgPad.shape, radius=9, uk=60, vk=80, n=2)  # 巴特沃斯陷波带阻滤波器
        BNRF2 = butterworthNRFilter(imgPad.shape, radius=9, uk=-60, vk=80, n=2)
        BNRF3 = butterworthNRFilter(imgPad.shape, radius=9, uk=60, vk=160, n=2)
        BNRF4 = butterworthNRFilter(imgPad.shape, radius=9, uk=-60, vk=160, n=2)
        freFilter = BNRF1 * BNRF2 * BNRF3 * BNRF4
    else:
        print("Error of unknown filter")
        return -1

    # (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数
    freTrans = fft * freFilter
    # (6) 傅里叶反变换
    ifft = np.fft.ifft2(freTrans)
    # (7) 去中心化反变换的图像
    M, N = img.shape[:2]
    mask2 = np.ones(imgPad.shape)
    mask2[1::2, ::2] = -1
    mask2[::2, 1::2] = -1
    ifftCenPad = ifft.real * mask2
    # (8) 截取左上角，大小和输入图像相等
    imgFilter = ifftCenPad[:M, :N]
    imgFilter = np.clip(imgFilter, 0, imgFilter.max())
    imgFilter = np.uint8(normalize(imgFilter) * 255)
    return imgFilter


# 使用陷波滤波删除数字化印刷图像中的莫尔模式
# (1) 读取原始图像
img = cv2.imread("./img/Fig0464(a)(car_75DPI_Moire).tif", flags=0)  # flags=0 读取为灰度图像
fig = plt.figure(figsize=(10, 5))
plt.subplot(141), plt.title("Original"), plt.axis('off'), plt.imshow(img, cmap='gray')

# (2) 图像高斯低通滤波
imgGLPF = imgFrequencyFilter(img, lpTyper="GaussLP", radius=30)  # 图像高斯低通滤波
plt.subplot(142), plt.title("GaussLP filter"), plt.axis('off'), plt.imshow(imgGLPF, cmap='gray')

# (3) 图像巴特沃斯陷波带阻滤波
imgBNRF = imgFrequencyFilter(img, lpTyper="ButterworthNR", radius=9)
plt.subplot(143), plt.title("ButterworthNR filter"), plt.axis('off'), plt.imshow(imgBNRF, cmap='gray')

# (4) 叠加巴特沃斯陷波带阻滤波
imgSBNRF = imgFrequencyFilter(img, lpTyper="MButterworthNR", radius=9)
plt.subplot(144), plt.title("Superimposed BNRF"), plt.axis('off'), plt.imshow(imgSBNRF, cmap='gray')

plt.tight_layout()
plt.show()

以上代码为youcans@xupt 原创作品，原文链接：https://blog.csdn.net/youcans/article/details/122781106
Copyright 2021 youcans, XUPT
Crated：2022-2-1

此外还能用陷波滤波去除周期干扰，总结一下，陷波滤波器的作用是选择性地抑制或去除图像或信号中的特定频率成分，以消除或减弱噪声、干扰或频率成分，从而改善图像的质量和清晰度

11.快速傅立叶变换

一般而言，图像处理的开销都比较大，下面是一些简化和加速傅立叶变换的方法

二维DFT的可分离性

二维的DFT可以分成多个一维的变换，如$$F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-2i\pi (ux/M+vy/N)}$$可以写成$$F(u,v)=\sum _{x=0}^{M-1}F(x,v)e^{-2i\pi ux/M}$$其中$$F(x,v)=\sum _{y=0}^{N-1}f(x,y)e^{-2i\pi uy/N}$$总结起来就是$f(x,y)$的二维DFT可计算$f(x,y)$的每一行的一维变换，然后沿每一列计算以为变换得到，由此可以由一维到二维的计算

用DFT计算IDFT

由$$f(x,y)=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{2i\pi (ux/M+vy/N)}$$取复共轭可以得到$$MN{f}^{\ast }(x,y)=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}{F}^{\ast }(u,v)e^{-2i\pi (ux/M+vy/N)}$$这个式子表明把$F^{\ast }(u,v)$带入二维傅立叶正变换的结果是$MN{f}^{\ast }(x,y)$，取复共轭再乘$\frac{1}{MN}$后就得到$F(x,y)$的反变换$f(x,y)$

快速傅立叶变换（FFT）

FFT（快速傅里叶变换）是一种高效的算法，用于计算离散傅里叶变换（DFT）和逆离散傅里叶变换（IDFT）。FFT在信号处理、图像处理、数字滤波、频谱分析等领域得到广泛应用。

FFT算法的主要思想是通过将DFT的计算复杂度从O(N^2)降低到O(N log N)，以更快的速度计算信号的频谱信息。它基于蝶形运算的思想，将DFT分解为多个较小规模的子问题，递归地进行计算，然后再进行合并得到最终的结果。

使用FFT算法，可以将一个长度为N的时域序列转换为频域序列，频域序列表示了信号在不同频率下的成分。FFT的输出是一个复数序列，包含了信号的幅度和相位信息。

对于频域序列，可以进行各种频谱分析，如查看频率成分的能量分布、滤波去除不需要的频率成分、频域图像处理等。另外，FFT还可以用于将频域序列转换回时域序列，即进行逆变换IDFT。

FFT算法的应用十分广泛，包括音频信号处理、图像处理、通信系统、数据压缩、图像压缩等领域。它提供了一种高效的方法来分析和处理信号的频谱特性，为各种应用提供了强大的工具。

章总结

本章学习了频率域滤波的相关知识：

1. 其含义和与空间域滤波的区别，在空间域滤波中，滤波器（也称为卷积核）直接应用于图像的每个像素位置，而频率域滤波通过在频率域中修改频谱信息，来实现对图像的滤波操作。
2. 取样、函数重建和混叠的概念也有了认识
3. 频率域卷积和滤波的区别：频率域卷积用傅里叶变换的性质，将两个图像的傅里叶变换相乘，并通过逆傅里叶变换得到卷积结果；而频率域滤波是通过在频率域中修改图像的频谱信息来实现滤波操作，通过将图像的傅里叶变换与滤波器的频谱相乘，并将结果进行逆傅里叶变换，可以得到滤波后的图像
4. 知道如何由空间和得到频率域滤波器函数，以及相反过程：
   • 从空间域到频率域滤波器函数：
     定义一个空间域滤波器函数，对滤波器函数进行傅里叶变换，将其转换为频率域。
     得到的频率域滤波器函数可以表示为复数形式，包括幅度谱和相位谱。
   • 从频率域到空间域滤波器函数：
     定义一个频率域滤波器函数，对其进行逆傅里叶变换，将其转换为空间域。
     得到的空间域滤波器函数可以表示为实数形式。
5. 怎么在频率域中直接构造滤波器传递函数：
   • 确定滤波器在频率域中的大小和形状。通常，滤波器的大小与图像的大小相同。
   • 创建滤波器传递函数：根据所需的滤波效果，在频率域中创建一个复数数组，表示滤波器的传递函数。传递函数可以包含幅度谱和相位谱。
   • 进行逆傅里叶变换：对滤波器传递函数进行逆傅里叶变换，将其转换回空间域。
6. 图像填充的重要性：
   • 保持图像尺寸和边界信息
   • 避免边缘信息丢失
   • 避免图像尺寸不匹配问题：
   • 减少边界影响
7. 频率域中的其他滤波技术，如钝化掩蔽、同态滤波等
8. 快速傅立叶变换的原理