专栏地址：『youcans 的图像处理学习课』
文章目录：『youcans 的图像处理学习课 - 总目录』
【youcans 的图像处理学习课】8. 频率域图像滤波（上）
图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制，是常用的图像预处理操作。
频率域图像处理先将图像进行傅里叶变换，然后在变换域进行处理，最后进行傅里叶反变换转换回空间域。
频率域滤波是滤波器传递函数与输入图像的傅里叶变换的对应像素相乘。频率域中的滤波概念更加直观，滤波器设计也更容易。使用快速傅里叶变换算法，比空间卷积运算速度快很多。
本文提供上述各种算法的完整例程和运行结果，并通过一个应用案例示范综合使用多种图像增强方法。
1. 傅里叶变换
滤波通常是指对图像中特定频率的分量进行过滤或抑制。图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制，是常用的图像预处理操作。
图像滤波不仅可以在空间域进行还可以在频率域进行。空间滤波是图像与各种空间滤波器（模板、核）的卷积，而空间卷积的傅里叶变换是频率域中相应变换的乘积，因此频率域滤波是频率域滤波器（传递函数）与图像的傅里叶变换相乘。
频率域图像滤波，先将图像进行傅里叶变换，然后在变换域进行处理，最后进行傅里叶反变换转换回空间域。
空间域滤波器和频率域滤波器形成一个傅里叶变换对：
f
(
x
,
y
)
⊗
h
(
x
,
y
)
⇔
F
(
u
,
v
)
H
(
u
,
v
)
f
(
x
,
y
)
h
(
x
,
y
)
⇔
F
(
u
,
v
)
⊗
H
(
u
,
v
)
f(x,y) \otimes h(x,y) \Leftrightarrow F(u,v)H(u,v) \\f(x,y) h(x,y) \Leftrightarrow F(u,v) \otimes H(u,v)
f(x,y)⊗h(x,y)⇔F(u,v)H(u,v)f(x,y)h(x,y)⇔F(u,v)⊗H(u,v)
也就是说，空间域滤波器和频率域滤波器实际上是相互对应的，有些空间域滤波器在频率域通过傅里叶变换实现会更方便、更快速。例如，空间域的拉普拉斯滤波器就是频率域的高通滤波器。
1.1 傅里叶级数与傅里叶变换
傅里叶级数（Fourier series）在数论、组合数学、信号处理、概率论、统计学、密码学、声学、光学等领域都有着广泛的应用。
傅里叶级数公式指出，任何周期函数都可以表示为不同频率的正弦函数和/或余弦函数的加权之和：
f
(
t
)
=
A
0
+
∑
n
=
1
∞
A
n
s
i
n
(
n
ω
t
+
ψ
n
)
=
A
0
+
∑
n
=
1
∞
[
a
n
c
o
s
(
n
ω
t
)
+
b
n
s
i
n
(
n
ω
t
)
]
\begin{aligned} f(t) &= A_0 + \sum^{\infty}_{n=1} A_n sin(n \omega t + \psi _n)\\ &= A_0 + \sum^{\infty}_{n=1} [a_n cos(n \omega t) + b_n sin(n \omega t)] \end{aligned}
f(t)=A0+n=1∑∞Ansin(nωt+ψn)=A0+n=1∑∞[ancos(nωt)+bnsin(nωt)]
这个和就是傅里叶级数。
进一步地，任何非周期函数也可以表示为不同频率的正弦函数和/或余弦函数乘以加权函数的积分：
F
(
ω
)
=
∫
−
∞
+
∞
f
(
t
)
e
−
j
ω
t
d
t
f
(
t
)
=
1
2
π
∫
−
∞
+
∞
F
(
ω
)
e
j
ω
t
d
ω
\begin{aligned} F(\omega) &= \int_{-\infty}^{+\infty} f(t) e^{-j\omega t} dt\\ f(t) &= \frac{1}{2 \pi} \int_{-\infty}^{+\infty} F(\omega) e^{j\omega t} d \omega \end{aligned}
F(ω)f(t)=∫−∞+∞f(t)e−jωtdt=2π1∫−∞+∞F(ω)ejωtdω
这个公式就是傅里叶变换（Fourier transform ）和逆变换。
*傅里叶变换存在的充分条件是：f(t) 的绝对值的积分是有限的，在信号处理、图像处理领域这一条件都能满足。
例程 8.1：连续周期信号的傅立叶级数
# 8.1：连续周期信号的傅立叶级数
from scipy import integrate
nf = 30
T = 10
tao = 1.0
y = 1
k = np.arange(0, nf)
an = np.zeros(nf)
bn = np.zeros(nf)
amp = np.zeros(nf)
pha = np.zeros(nf)
half0, err0 = integrate.quad(lambda t: y, -tao/2, tao/2)
an[0] = 2 * half0 / T
for n in range(1, nf):
half1, err1 = integrate.quad(lambda t: 2*y * np.cos(2.0/T * np.pi * n * t), -tao/2, tao/2)
an[n] = half1 / 10
half2, err2 = integrate.quad(lambda t: 2*y * np.sin(2.0/T * np.pi * n * t), -tao/2, tao/2)
bn[n] = half2 / 10
amp[n] = np.sqrt(an[n]**2 + bn[n]**2)
for i in range(0, nf):
pha[i] = 0.0 if an[i]>=0 else np.pi
plt.figure(figsize=(9, 6))
plt.subplot(211), plt.title("Amplitude spectrum"), plt.stem(k, amp)
plt.subplot(212), plt.title("Phase spectrum"), plt.stem(k, pha)
plt.show()
例程 8.2：连续非周期信号的傅立叶系数
# 8.2：连续非周期信号的傅立叶系数
dx = 0.001
x = np.pi * np.arange(-1+dx, 1+dx, dx)
n = len(x)
nquart = int(np.floor(n/4))
f = np.zeros_like(x)
f[nquart: 2*nquart] = (4/n) * np.arange(1, nquart+1)
f[2*nquart: 3*nquart] = np.ones(nquart) - (4/n) * np.arange(0, nquart)
plt.figure(figsize=(9, 6))
plt.title("Fourier series of hat function")
plt.plot(x, f, '-', color='k',label="Original")
# Compute Fourier series
A0 = np.sum(f * np.ones_like(x)) * dx
fFS = A0 / 2
orders = 3
A = np.zeros(orders)
B = np.zeros(orders)
for k in range(orders):
A[k] = np.sum(f * np.cos((k+1) * x)) * dx # Inner product
B[k] = np.sum(f * np.sin((k+1) * x)) * dx
fFS = fFS + A[k] * np.cos((k + 1) * x) + B[k] * np.sin((k + 1) * x)
plt.plot(x, fFS, '-', label="{} order".format(k))
print('Line ', k, ': A = ', A[k], ' B = ', B[k], fFS.shape)
plt.legend(loc="upper right")
plt.show()
例程 8.3：一维连续函数的傅里叶变换（盒式函数）
以一维盒式函数为例，其傅里叶变换为：
F ( ω ) = ∫ − ∞ + ∞ f ( t ) e − j ω t d t = ∫ − π π f ( t ) e − j ω t d t = s i n ( ω ) ω \begin{aligned} F(\omega) &= \int_{-\infty}^{+\infty} f(t) e^{-j\omega t} dt\\ &= \int_{-\pi}^{\pi} f(t) e^{-j\omega t} dt\\ &= \frac{sin(\omega)}{\omega} \end{aligned} F(ω)=∫−∞+∞f(t)e−jωtdt=∫−ππf(t)e−jωtdt=ωsin(ω)
# 8.3：一维连续函数的傅里叶变换（盒式函数）
# 盒式函数 (Box_function)
x = np.arange(-5, 5, 0.1)
w = 2 * np.pi
halfW = w / 2
y = np.where(x, x > -halfW, 0)
y = np.where(x < halfW, y, 0)
plt.figure(figsize=(9, 4))
plt.subplot(131), plt.title("Box_function")
plt.plot(x, y, '-b')
plt.subplot(132), plt.title("Fourier transform")
fu = np.sinc(x)
plt.plot(x, fu, '-b')
plt.subplot(133), plt.title("Spectrum of box_func")
fu = abs(fu)
plt.plot(x, fu, '-b')
plt.show()
1.2 连续函数的取样
连续函数必须经过取样和量化转换为离散函数，才能用计算机进行处理。
考虑一个连续函数 f ( t ) f(t) f(t)，以自变量 t 的均匀间隔 Δ T \Delta T ΔT 对函数取样，取样序列中的任意取样值为：
f k = ∫ − ∞ + ∞ f ( t ) δ ( t − k Δ T ) d t = f ( k Δ T ) f_k = \int_{-\infty}^{+\infty} f(t) \delta (t-k \Delta T) dt = f(k \Delta T) fk=∫−∞+∞f(t)δ(t−kΔT)dt=f(kΔT)
取样后的函数的傅里叶变换为：
F ~ ( μ ) = ( F ⋆ S ) ( μ ) = ∫ − ∞ + ∞ F ( τ ) S ( μ − τ ) d τ = 1 Δ T ∑ n = − ∞ ∞ F ( μ − n Δ T ) \begin{aligned} \tilde{F}(\mu) &= (F \star S) (\mu) \\ &= \int_{-\infty}^{+\infty} F(\tau) S(\mu-\tau) d \tau\\ &= \frac{1}{\Delta T} \sum_{n=-\infty}^{\infty}F(\mu-\frac{n}{\Delta T}) \end{aligned} F~(μ)=(F⋆S)(μ)=∫−∞+∞F(τ)S(μ−τ)dτ=ΔT1n=−∞∑∞F(μ−ΔTn)
香农（Shannon）取样定理指出，对于一个连续信号，用大于信号最高频率 2倍的取样率来取样，则不会丢失信号的有效信息。或者说，以 1 / Δ T 1/\Delta T 1/ΔT 的取样率对信号取样点的的最大频率是 μ m a x = 1 / 2 Δ T \mu_{max}=1/2\Delta T μmax=1/2ΔT 。
例程 8.4：连续函数的取样
# 8.4：连续函数的取样
# 定义函数，用于计算所有基矢的频率
def gen_freq(N, fs):
k = np.arange(0, np.floor(N/2) + 1, 1)
return (k * fs) / N
T = 100
# 定义多个不同频率的基频信号
fk = [2/T, 5/T, 12/T] # 频率
A = [7, 3, 2] # 振幅
phi = [np.pi, 2, 2*np.pi] # 初始相位
n = np.arange(T)
s0 = A[0] * np.sin(2 * np.pi * fk[0] * n + phi[0])
s1 = A[1] * np.sin(2 * np.pi * fk[1] * n + phi[1])
s2 = A[2] * np.sin(2 * np.pi * fk[2] * n + phi[2])
s = s0 + s1 + s2 # 叠加生成混合信号
g = np.fft.rfft(s) # 傅里叶变换
plt.figure(figsize=(8, 6))
plt.subplot(311)
plt.plot(n, s0, n, s1, n, s2, ':', marker='+', alpha=0.5)
plt.plot(n, s, 'r-', lw=2)
plt.title("Sampling of continuous functions")
plt.subplot(312)
fs = 1 # 采样间隔为 1
freq = gen_freq(T, fs=fs) # 计算频率序列
ck = np.abs(g) / T # 计算振幅
plt.plot(freq, ck, '.') # 频率-振幅图
for f in fk:
ck0 = round(ck[np.where(freq==f*fs)][0], 1)
plt.annotate('$({},{})$'.format(f*fs, ck0), xy=(f*fs, ck0), xytext=(5, -10), textcoords='offset points')
plt.subplot(313)
fs = 10 # 采样间隔为 10
freq = gen_freq(T, fs=fs) # 计算频率序列
ck = np.abs(g) / T # 计算振幅
plt.plot(freq, ck, '.') # 频率-振幅图
for f in fk:
ck0 = round(ck[np.where(freq==f*fs)][0], 1)
plt.annotate('$({},{})$'.format(f*fs, ck0), xy=(f*fs, ck0), xytext=(5, -10), textcoords='offset points')
plt.show()
1.3 一维离散傅里叶变换
数字信号和数字图像都是采样得到的离散变量。
对原函数的变换取样后的数据
F ~ ( μ ) = ∫ − ∞ + ∞ f ~ ( t ) e − j 2 π μ t d t = ∫ − ∞ + ∞ ∑ n = − ∞ ∞ f ( t ) δ ( t − n Δ T ) e − j 2 π μ t d t = ∑ n = − ∞ ∞ ∫ − ∞ + ∞ f ( t ) δ ( t − n Δ T ) e − j 2 π μ t d t = ∑ n = − ∞ ∞ f n e − j 2 π μ n Δ T \begin{aligned} \tilde{F}(\mu) &= \int_{-\infty}^{+\infty} \tilde{f}(t) e^{-j 2 \pi \mu t} dt\\ &= \int_{-\infty}^{+\infty} \sum_{n=-\infty}^{\infty} f(t) \delta (t-n{\Delta T}) e^{-j 2 \pi \mu t} dt\\ &= \sum_{n=-\infty}^{\infty} \int_{-\infty}^{+\infty} f(t) \delta (t-n{\Delta T}) e^{-j 2 \pi \mu t} dt\\ &= \sum_{n=-\infty}^{\infty} f_n e^{-j 2 \pi \mu n{\Delta T}}\\ \end{aligned} F~(μ)=∫−∞+∞f~(t)e−j2πμtdt=∫−∞+∞n=−∞∑∞f(t)δ(t−nΔT)e−j2πμtdt=n=−∞∑∞∫−∞+∞f(t)δ(t−nΔT)e−j2πμtdt=n=−∞∑∞fne−j2πμnΔT
当取样频率为 μ = m / M Δ T \mu = m/M \Delta T μ=m/MΔT 时，可以得到离散傅里叶变换（DFT）和逆变换公式为：
F m = ∑ n = 0 M − 1 f n e − j 2 π μ m n / M , m = 0 , . . . M − 1 f n = 1 M ∑ m = 0 M − 1 F m e j 2 π μ m n / M , n = 0 , . . . M − 1 \begin{aligned} F_m &= \sum_{n=0}^{M-1} f_n e^{-j\ 2\pi \mu mn/M}, &m=0,...M-1\\ f_n &= \frac{1}{M} \sum_{m=0}^{M-1} F_m e^{j\ 2\pi \mu mn/M}, &n=0,...M-1 \end{aligned} Fmfn=n=0∑M−1fne−j 2πμmn/M,=M1m=0∑M−1Fmej 2πμmn/M,m=0,...M−1n=0,...M−1
由于任何周期或非周期函数都可以表示为不同频率的正弦函数和余弦函数之和的形式，因此用傅里叶变换表示的函数特征完全可以通过傅里叶反变换来重建，而且不会丢失任何信息。这是频率域图像处理的数学基础。
离散傅里叶变换 是将离散信号分解为一系列离散三角函数分量，每个三角函数分量都有对应的幅值 A、频率 f 和相位 φ \varphi φ。通过所有分量叠加可以得到原离散信号。
在图像处理中，通常使用 x , y x, y x,y 表示离散的空间域坐标变量，用 u , v u,v u,v 表示离散的频率域变量。于是将一维离散傅里叶变换表示为：
F ( u ) = ∑ x = 0 M − 1 f ( x ) e − j 2 π u x / M , u = 0 , . . . M − 1 f ( x ) = 1 M ∑ u = 0 M − 1 F ( u ) e j 2 π u x / M , x = 0 , . . . M − 1 \begin{aligned} F(u) &= \sum_{x=0}^{M-1} f(x) e^{-j\ 2\pi u x/M}, &u=0,...M-1\\ f(x) &= \frac{1}{M} \sum_{u=0}^{M-1} F(u) e^{j\ 2\pi u x/M}, &x=0,...M-1 \end{aligned} F(u)f(x)=x=0∑M−1f(x)e−j 2πux/M,=M1u=0∑M−1F(u)ej 2πux/M,u=0,...M−1x=0,...M−1
例程 8.6：一维离散傅里叶变换
# 8.6：一维离散傅里叶变换
# 生成方波信号
N = 200
t = np.linspace(-10, 10, N)
dt = t[1] - t[0]
sint = np.sin(t)
sig = np.sign(sint)
fig = plt.figure(figsize=(10, 4))
plt.subplot(131), plt.title("source"), plt.xticks([]), plt.yticks([])
plt.plot(t, sig)
# 离散傅里叶变换
fft = np.fft.fft(sig, N) # 离散傅里叶变换
fftshift = np.fft.fftshift(fft) # 对称平移
amp = abs(fftshift) / len(fft) # 幅值
pha = np.angle(fftshift) # 相位
fre = np.fft.fftshift(np.fft.fftfreq(d=dt, n=N)) # 频率
plt.subplot(132), plt.title("DFT"), plt.xticks([]), plt.yticks([])
plt.plot(t, fft)
# 信号恢复
recover = np.zeros(N)
for a, p, f in zip(amp, pha, fre):
sigCos = a * np.cos(2 * np.pi * f * np.arange(N) * dt + p) # 根据幅度,相位,频率构造三角函数
recover += sigCos # 把这些三角函数都加起来
plt.subplot(133), plt.title("recover"), plt.xticks([]), plt.yticks([])
plt.plot(t, recover)
plt.show()
2. 二维函数的傅里叶变换
2.1 二维连续傅里叶变换
设 f ( t , z ) f(t,z) f(t,z) 是二维连续变量 t , z t, z t,z 的连续函数，则其二维连续傅里叶变换和逆变换为：
F
(
μ
,
ν
)
=
∫
−
∞
+
∞
∫
−
∞
+
∞
f
(
t
,
z
)
e
−
j
2
π
(
μ
t
+
ν
z
)
d
t
d
z
f
(
t
,
z
)
=
∫
−
∞
+
∞
∫
−
∞
+
∞
F
(
μ
,
ν
)
e
j
2
π
(
μ
t
+
ν
z
)
d
μ
d
ν
\begin{aligned} F(\mu,\nu) &= \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(t,z) e^{-j 2\pi (\mu t + \nu z)} dt \ dz\\ f(t,z) &= \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} F(\mu,\nu) e^{j 2\pi (\mu t + \nu z)} d\mu \ d\nu \end{aligned}
F(μ,ν)f(t,z)=∫−∞+∞∫−∞+∞f(t,z)e−j2π(μt+νz)dt dz=∫−∞+∞∫−∞+∞F(μ,ν)ej2π(μt+νz)dμ dν
对于图像处理，式中的
t
,
z
t, z
t,z 表示连续空间变量，
μ
,
ν
\mu, \nu
μ,ν 表示连续频率变量。
例程 8.7：二维连续函数的傅里叶变换（二维盒式函数）
以盒式函数为例，其傅里叶变换为：
F ( μ , ν ) = ∫ − ∞ + ∞ ∫ − ∞ + ∞ f ( t , z ) e − j 2 π ( μ t + ν z ) d t d z = ∫ − T / 2 + T / 2 ∫ − Z / 2 + Z / 2 A e − j 2 π ( μ t + ν z ) d t d z = A T Z [ s i n ( π μ T ) π μ T ] [ s i n ( π ν Z ) π ν Z ] \begin{aligned} F(\mu,\nu) &= \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(t,z) e^{-j 2\pi (\mu t + \nu z)} dt \ dz\\ &= \int_{-T/2}^{+T/2} \int_{-Z/2}^{+Z/2} A e^{-j 2\pi (\mu t + \nu z)} dt \ dz\\ &= ATZ[\frac{sin(\pi \mu T)}{\pi \mu T}][\frac{sin(\pi \nu Z)}{\pi \nu Z}] \end{aligned} F(μ,ν)=∫−∞+∞∫−∞+∞f(t,z)e−j2π(μt+νz)dt dz=∫−T/2+T/2∫−Z/2+Z/2Ae−j2π(μt+νz)dt dz=ATZ[πμTsin(πμT)][πνZsin(πνZ)]
# 8.7：二维连续函数的傅里叶变换（二维盒式函数）
# 二维盒式函数 (2D-Box_function)
t = np.arange(-5, 5, 0.1) # start, end, step
z = np.arange(-5, 5, 0.1)
height, width = len(t), len(z)
tt, zz = np.meshgrid(t, z) # 将一维数组 xnew, ynew 转换为网格点集（二维数组）
f = np.zeros([len(t), len(z)]) # 初始化，置零
f[30:70, 30:70] = 1 # 二维盒式函数
fu = np.sinc(t)
fv = np.sinc(z)
fuv = np.outer(np.sinc(t), np.sinc(z)) # 由公式计算连续傅里叶变换
print(fu.shape, fv.shape, fuv.shape)
fig = plt.figure(figsize=(10, 6))
ax1 = plt.subplot(121, projection='3d')
ax1.set_title("2D-Box function")
ax1.plot_wireframe(tt, zz, f, rstride=2, cstride=2, linewidth=0.5, color='c')
ax1.set_xticks([]), ax1.set_yticks([]), ax1.set_zticks([])
ax2 = plt.subplot(122, projection='3d')
ax2.set_title("2D-Fourier transform")
ax2.plot_wireframe(tt, zz, fuv, rstride=2, cstride=2, linewidth=0.5, color='c')
ax2.set_xticks([]), ax2.set_yticks([]), ax2.set_zticks([])
plt.show()
例程 8.8：二维连续函数的傅里叶变换（不同参数的盒式函数）
# 8.8：二维连续函数的傅里叶变换（不同参数的盒式函数）
# 二维盒式函数 (2D-Box_function)
height, width = 128, 128
m = int((height - 1) / 2)
n = int((width - 1) / 2)
fig = plt.figure(figsize=(9, 6))
T = [5, 10, 20]
print(len(T))
for i in range(len(T)):
f = np.zeros([height, width]) # 初始化，置零
f[m-T[i]:m+T[i], n-T[i]:n+T[i]] = 1 # 不同尺寸的盒式函数
fft = np.fft.fft2(f)
shift_fft = np.fft.fftshift(fft)
amp = np.log(1 + np.abs(shift_fft))
plt.subplot(2,len(T),i+1), plt.xticks([]), plt.yticks([])
plt.imshow(f, "gray"), plt.title("Box Filter (T={})".format(T[i]))
#
plt.subplot(2,len(T),len(T)+i+1), plt.xticks([]), plt.yticks([])
plt.imshow(amp, "gray"), plt.title("FFT Spectrum")
plt.tight_layout()
plt.show()
2.2 图像的混叠和重取样
由于无法对一个函数无限地取样，因此在数字图像中总是会出现混叠。
混叠分为空间混叠和时间混叠。空间混叠是由欠取样造成的，在陌生重复的图像中较为明显。时间混叠与动态图像序列中图像的时间间隔相关，如辐条倒转的“车轮效应”。
空间混叠的主要问题是会引入伪影，即原始图像中未出现的线条锯齿、虚假高光和频率模式。
对图像进行缩放都会引入混叠。图像放大可以视为过取样，使用双线性、双三次内插可以降低图像放大中的混叠。图像缩小可以视为欠取样，混叠通常更为严重。
为了降低混叠，在重取样前要使用低通滤波器来平滑，以衰减图像的高频分量，可以有效地抑制混叠，但图像的清晰度也有所下降。
例程 8.9：图像的抗混叠
# 8.9：图像的抗混叠
imgGray = cv2.imread("../images/Fig0417a.tif", flags=0) # flags=0 读取为灰度图像
height, width = imgGray.shape[:2] # 图片的高度和宽度
# 先缩小图像，再放大图像
imgZoom1 = cv2.resize(imgGray, (int(0.33*width), int(0.33*height)))
imgZoom2 = cv2.resize(imgZoom1, None, fx=3, fy=3, interpolation=cv2.INTER_AREA)
# 先对原图像进行5x5的低通滤波，再缩小图像，再放大图像
kSize = (5, 5)
imgBoxFilter = cv2.boxFilter(imgGray, -1, kSize) # cv2.boxFilter 方法
imgZoom3 = cv2.resize(imgBoxFilter, (int(0.33*width), int(0.33*height)))
imgZoom4 = cv2.resize(imgZoom3, None, fx=3, fy=3, interpolation=cv2.INTER_AREA)
fig = plt.figure(figsize=(9,5))
plt.subplot(131), plt.axis('off'), plt.title("Origin")
plt.imshow(imgGray, cmap='gray')
plt.subplot(132), plt.axis('off'), plt.title("Resample")
plt.imshow(imgZoom2, cmap='gray')
plt.subplot(133), plt.axis('off'), plt.title("Box Resample")
plt.imshow(imgZoom4, cmap='gray')
plt.tight_layout()
plt.show()
2.3 二维离散傅里叶变换（DFT）
对于二维图像处理，通常使用 x , y x, y x,y 表示离散的空间域坐标变量，用 u , v u,v u,v 表示离散的频率域变量。二维离散傅里叶变换（DFT）和反变换（IDFT）为：
F
(
u
,
v
)
=
∑
x
=
0
M
−
1
∑
y
=
0
N
−
1
f
(
x
,
y
)
e
−
j
2
π
(
u
x
/
M
+
v
y
/
N
)
f
(
x
,
y
)
=
1
M
N
∑
u
=
0
M
−
1
∑
v
=
0
N
−
1
F
(
u
,
v
)
e
j
2
π
(
u
x
/
M
+
v
y
/
N
)
\begin{aligned} F(u,v) &= \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y) e^{-j 2\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^{j 2\pi (ux/M+vy/N)} \end{aligned}
F(u,v)f(x,y)=x=0∑M−1y=0∑N−1f(x,y)e−j2π(ux/M+vy/N)=MN1u=0∑M−1v=0∑N−1F(u,v)ej2π(ux/M+vy/N)
二维离散傅里叶变换也可以用极坐标表示：
F
(
u
,
v
)
=
R
(
u
,
v
)
+
j
I
(
u
,
v
)
=
∣
F
(
u
,
v
)
∣
e
j
ϕ
(
u
,
v
)
F(u,v) = R(u,v) + j I(u,v) = |F(u,v)| e^{j \phi (u,v)}
F(u,v)=R(u,v)+jI(u,v)=∣F(u,v)∣ejϕ(u,v)
傅里叶频谱（Fourier spectrum）为：
∣
F
(
u
,
v
)
∣
=
[
R
2
(
u
,
v
)
+
I
2
(
u
,
v
)
]
1
/
2
|F(u,v)| = [R^2(u,v) + I^2(u,v)]^{1/2}
∣F(u,v)∣=[R2(u,v)+I2(u,v)]1/2
傅里叶相位谱（Fourier phase spectrum）为：
ϕ
(
u
,
v
)
=
a
r
c
t
a
n
[
I
(
u
,
v
)
/
R
(
u
,
v
)
]
\phi (u,v) = arctan[I(u,v)/R(u,v)]
ϕ(u,v)=arctan[I(u,v)/R(u,v)]
傅里叶功率谱（Fourier power spectrum）为：
P
(
u
,
v
)
=
∣
F
(
u
,
v
)
∣
2
=
R
2
(
u
,
v
)
+
I
2
(
u
,
v
)
P(u,v) = |F(u,v)|^2 = R^2(u,v) + I^2(u,v)
P(u,v)=∣F(u,v)∣2=R2(u,v)+I2(u,v)
空间取样和频率间隔是相互对应的，频率域所对应的离散变量间的间隔为： Δ u = 1 / M Δ T ， Δ v = 1 / N Δ Z \Delta u = 1/M \Delta T，\Delta v = 1/N \Delta Z Δu=1/MΔT，Δv=1/NΔZ。即：频域中样本之间的间隔，与空间样本之间的间隔及样本数量的乘积成反比。
空间域滤波器和频率域滤波器也是相互对应的，二维卷积定理是在空间域和频率域滤波之间建立等价关系的纽带：
(
f
⋆
h
)
(
x
,
y
)
⇔
(
F
⋅
H
)
(
u
,
v
)
(f \star h)(x,y) \Leftrightarrow (F \cdot H)(u,v)
(f⋆h)(x,y)⇔(F⋅H)(u,v)
这表明 F 和 H 分别是 f 和 h 的傅里叶变换；f 和 h 的空间卷积的傅里叶变换，是它们的变换的乘积。
因此，计算两个函数的空间卷积，可以直接在空间域计算，也可以在频率域计算：先计算每个函数的傅里叶变换，再对两个变换相乘，最后进行傅里叶反变换转换回空间域。
也就是说，空间域滤波器和频率域滤波器实际上是相互对应的，有些空间域滤波器在频率域通过傅里叶变换实现会更方便、更快速。
2.4 Numpy 实现图像傅里叶变换
Numpy 中的 np.fft.fft2() 函数可以实现图像的傅里叶变换 。
函数说明：
numpy.fft.fft2(a, s=None, axes=(- 2, - 1), norm=None) → out
参数说明：
- a：输入数组，一维或多维数组，可以是复数形式
- s：输出数组的形状（每个变换轴的长度），对应于 fft(x,n) 中的 n，可选项
- out：输出数组，复数形式的一维或多维数组（complex ndarray）
注意事项：
- 用于二维图像傅里叶变换时，输入数组 a 是 numpy 二维数组（灰度图像）或三维数组（彩色图像）。
- 输出结果是复数形式 ( R e a l + j ∗ I m a g ) (Real + j * Imag) (Real+j∗Imag) 的数组，不能直接用于显示图像。为了显示傅里叶变换结果的图像，需要将数组的值调整到 [0,255] 的灰度空间内。
经过 np.fft.fft2() 函数实现傅里叶变换得到的图片频谱信息，幅度谱的最大值（低频分量）在左上角 (0,0) 处。为了便于观察，通常用 np.fft.fftshift() 函数将低频分量移动到频域图像的中心位置。
函数说明：
numpy.fft.fftshift(x, axes=None) → y
参数说明：
- x：输入数组，一维或多维数组
- axes：整数，或输入数组形状的元组，用于指定移动的轴，可选项
- y：输出数组
例程 8.10：二维图像的离散傅里叶变换（Numpy）
# 8.10：Numpy 实现二维离散傅里叶变换
normalize = lambda x: (x - x.min()) / (x.max() - x.min() + 1e-6)
imgGray = cv2.imread("../images/Fig0424a.tif", flags=0) # flags=0 读取为灰度图像
# imgGray = cv2.imread("../images/imgBall.png", flags=1) # flags=0 读取为灰度图像
# 傅里叶变换
# fft = np.fft.fft2(imgGray.astype(np.float32))
fft = np.fft.fft2(imgGray) # np.fft.fft2 实现傅里叶变换
# 非中心化，计算幅度谱和相位谱
ampSpectrum = np.sqrt(np.power(fft.real, 2) + np.power(fft.imag, 2)) # 幅度谱
print("ampSpectrum max={}, min={}".format(ampSpectrum.max(), ampSpectrum.min()))
# phase = np.arctan2(fft.imag, fft.real) # 计算相位角(弧度制)
# phiSpectrum = phase / np.pi*180 # 将相位角转换为 [-180, 180]
phiSpectrum = np.angle(fft)
# 中心化，将低频分量移动到频域图像的中心
fftShift = np.fft.fftshift(fft) # 将低频分量移动到频域图像的中心
# 中心化后的幅度谱
ampSpeShift = np.sqrt(np.power(fftShift.real, 2) + np.power(fftShift.imag, 2))
ampShiftNorm = np.uint8(normalize(ampSpeShift)*255) # 归一化为 [0,255]
# 幅度谱做对数变换
ampSpeLog = np.log(1 + ampSpeShift) # 幅度谱做对数变换以便于显示
ampSpeLog = np.uint8(normalize(ampSpeLog)*255) # 归一化为 [0,255]
# np.fft.ifft2 实现图像的逆傅里叶变换
invShift = np.fft.ifftshift(fftShift) # 将低频逆转换回图像四角
imgIfft = np.fft.ifft2(invShift) # 逆傅里叶变换，返回值是复数数组
imgRebuild = np.abs(imgIfft) # 将复数数组调整至灰度空间
plt.figure(figsize=(9, 6))
plt.subplot(231), plt.title("Original image"), plt.axis('off')
plt.imshow(imgGray, cmap='gray')
plt.subplot(232), plt.title("FFT phase spectrum"), plt.axis('off')
plt.imshow(phiSpectrum, cmap='gray')
plt.subplot(233), plt.title("Rebuild image with IFFT"), plt.axis('off')
plt.imshow(imgRebuild, cmap='gray')
plt.subplot(234), plt.title("FFT amplitude spectrum"), plt.axis('off')
plt.imshow(ampSpectrum, cmap='gray')
plt.subplot(235), plt.title("FFT-shift amplitude"), plt.axis('off')
plt.imshow(ampSpeShift, cmap='gray')
plt.subplot(236), plt.title("Log-trans of FFT amp"), plt.axis('off')
plt.imshow(ampSpeLog, cmap='gray')
plt.tight_layout()
plt.show()
程序说明：
图中未中心化的幅度谱（FFT amp spe）也并不是完全黑色，在图像的四角位置都有微小的亮区域，但是很难观察到，这也是对幅度谱进行中心化处理（fftShift）的原因。
2.5 OpenCV 实现图像傅里叶变换（cv.dft）
使用 OpenCV 中的 cv.dft() 函数也可以实现图像的傅里叶变换，cv.idft() 函数实现图像傅里叶逆变换。
函数说明：
cv.dft(src[, dst[, flags[, nonzeroRows]]]) → dst
cv.idft(src[, dst[, flags[, nonzeroRows]]]) → dst
参数说明：
- src：输入图像，单通道灰度图像，使用 np.float32 格式
- dst：输出图像，图像大小与 src 相同，数据类型由 flag 决定
- flag：转换标识符
- cv.DFT_INVERSE：用一维或二维逆变换取代默认的正向变换
- cv.DFT_SCALE：缩放比例标识，根据元素数量求出缩放结果，常与DFT_INVERSE搭配使用
- cv.DFT_ROWS: 对输入矩阵的每行进行正向或反向的傅里叶变换，常用于三维或高维变换等复杂操作
- cv.DFT_COMPLEX_OUTPUT：对一维或二维实数数组进行正向变换，默认方法，结果是由 2个通道表示的复数阵列，第一通道是实数部分，第二通道是虚数部分
- cv.DFT_REAL_OUTPUT：对一维或二维复数数组进行逆变换，结果通常是一个尺寸相同的复数矩阵
注意事项：
- 输入图像 src 是 np.float32 格式，如图像使用 np.uint8 格式则必须先转换 np.float32 格式。
- 默认方法 cv.DFT_COMPLEX_OUTPUT 时，输入 src 是 np.float32 格式的单通道二维数组，输出 dst 是 2个通道的二维数组，第一通道 dft[:,:,0] 是实数部分，第二通道 dft[:,:,1] 是虚数部分。
- 不能直接用于显示图像。可以使用 cv.magnitude() 函数将傅里叶变换的结果转换到灰度 [0,255]。
- idft(src, dst, flags) 等价于 dft(src, dst, flags=DFT_INVERSE)。
- OpenCV 实现傅里叶变换，计算速度比 Numpy 更快。
转换标识符为 cv.DFT_COMPLEX_OUTPUT 时，cv.dft() 函数的输出是 2个通道的二维数组，使用 cv.magnitude() 函数可以实现计算二维矢量的幅值 。
函数说明：
cv.magnitude(x, y[, magnitude]) → dst
参数说明：
- x：一维或多维数组，也表示复数的实部，浮点型
- y：一维或多维数组，也表示复数的虚部，浮点型，数组大小必须与 x 相同
- dst：输出数组，数组大小和数据类型与 x 相同，运算公式为：
d s t ( I ) = x ( I ) 2 + y ( I ) 2 dst(I) = \sqrt{x(I)^2 + y(I)^2} dst(I)=x(I)2+y(I)2
傅里叶变换及相关操作的取值范围可能不适于图像显示，需要进行归一化处理。 OpenCV 中的 cv.normalize() 函数可以实现图像的归一化。
函数说明：
cv.normalize(src, dst[, alpha[, beta[, norm_type[, dtype[, mask]]]]]) → dst
参数说明：
- src：输入图像
- dst：输出结果，与输入图像同尺寸同类型
- alpha：归一化后的最小值，可选项，默认值为0
- beta：归一化后的最大值，可选项，默认值为1
- norm_type：归一化类型
- NORM_INF：Linf 范数（绝对值的最大值）
- NORM_L1：L1 范数（绝对值的和）
- NORM_L2：L2 范数（欧几里德距离），默认类型
- NORM_MINMAX：线性缩放，常用类型
- dtype：可选项，默认值 -1，表示输出矩阵与输入图像类型相同
- mask：掩模遮罩，可选项，默认无遮罩
傅里叶变换在理论上需要 O ( M N ) 2 O(MN)^2 O(MN)2 次运算，非常耗时；快速傅里叶变换只需要 O ( M N l o g ( M N ) ) O(MN log (MN)) O(MNlog(MN)) 次运算就可以完成。
OpenCV 中的傅里叶变换函数 cv.dft() 对于行数和列数都可以分解为 2 p ∗ 3 q ∗ 5 r 2^p * 3^q * 5^r 2p∗3q∗5r 的矩阵的计算性能最好。为了提高运算性能，可以对原矩阵的右侧和下方补 0，以满足该分解条件。OpenCV 中的 cv.getOptimalDFTSize() 函数可以实现图像的最优 DFT 尺寸扩充，适用于 cv.dft() 和 np.fft.fft2()。
函数说明：
cv.getOptimalDFTSize(versize) → retval
参数说明：
- versize：数组大小
- retval：DFT 扩充的最优数组大小
例程 8.11：二维图像的离散傅里叶变换（OpenCV）
# 8.11：OpenCV 实现二维图像的离散傅里叶变换
imgGray = cv2.imread("../images/Fig0424a.tif", flags=0) # flags=0 读取为灰度图像
# cv2.dft 实现图像的傅里叶变换
imgFloat32 = np.float32(imgGray) # 将图像转换成 float32
dft = cv2.dft(imgFloat32, flags=cv2.DFT_COMPLEX_OUTPUT) # 傅里叶变换
dftShift = np.fft.fftshift(dft) # 将低频分量移动到频域图像的中心
# 幅度谱
# ampSpe = np.sqrt(np.power(dft[:,:,0], 2) + np.power(dftShift[:,:,1], 2))
dftAmp = cv2.magnitude(dft[:,:,0], dft[:,:,1]) # 幅度谱，未中心化
dftShiftAmp = cv2.magnitude(dftShift[:,:,0], dftShift[:,:,1]) # 幅度谱，中心化
dftAmpLog = np.log(1 + dftShiftAmp) # 幅度谱对数变换，以便于显示
# 相位谱
phase = np.arctan2(dftShift[:,:,1], dftShift[:,:,0]) # 计算相位角(弧度制)
dftPhi = phase / np.pi*180 # 将相位角转换为 [-180, 180]
print("dftMag max={}, min={}".format(dftAmp.max(), dftAmp.min()))
print("dftPhi max={}, min={}".format(dftPhi.max(), dftPhi.min()))
print("dftAmpLog max={}, min={}".format(dftAmpLog.max(), dftAmpLog.min()))
# cv2.idft 实现图像的逆傅里叶变换
invShift = np.fft.ifftshift(dftShift) # 将低频逆转换回图像四角
imgIdft = cv2.idft(invShift) # 逆傅里叶变换
imgRebuild = cv2.magnitude(imgIdft[:,:,0], imgIdft[:,:,1]) # 重建图像
plt.figure(figsize=(9, 6))
plt.subplot(231), plt.title("Original image"), plt.axis('off')
plt.imshow(imgGray, cmap='gray')
plt.subplot(232), plt.title("DFT Phase"), plt.axis('off')
plt.imshow(dftPhi, cmap='gray')
plt.subplot(233), plt.title("Rebuild image with IDFT"), plt.axis('off')
plt.imshow(imgRebuild, cmap='gray')
plt.subplot(234), plt.title("DFT amplitude spectrum"), plt.axis('off')
plt.imshow(dftAmp, cmap='gray')
plt.subplot(235), plt.title("DFT-shift amplitude"), plt.axis('off')
plt.imshow(dftShiftAmp, cmap='gray')
plt.subplot(236), plt.title("Log-trans of DFT amp"), plt.axis('off')
plt.imshow(dftAmpLog, cmap='gray')
plt.tight_layout()
plt.show()
例程 8.12：OpenCV 快速傅里叶变换
# 8.12：OpenCV 快速傅里叶变换
imgGray = cv2.imread("../images/Fig0429a.tif", flags=0) # flags=0 读取为灰度图像
rows,cols = imgGray.shape[:2] # 图像的行(高度)/列(宽度)
# 快速傅里叶变换(要对原始图像进行矩阵扩充)
rPad = cv2.getOptimalDFTSize(rows) # 最优 DFT 扩充尺寸
cPad = cv2.getOptimalDFTSize(cols) # 用于快速傅里叶变换
imgEx = np.zeros((rPad,cPad,2),np.float32) # 对原始图像进行边缘扩充
imgEx[:rows,:cols,0] = imgGray # 边缘扩充，下侧和右侧补0
dftImgEx = cv2.dft(imgEx, cv2.DFT_COMPLEX_OUTPUT) # 快速傅里叶变换
# 傅里叶逆变换
idftImg = cv2.idft(dftImgEx) # 逆傅里叶变换
idftMag = cv2.magnitude(idftImg[:,:,0], idftImg[:,:,1]) # 逆傅里叶变换幅值
# 矩阵裁剪，得到恢复图像
idftMagNorm = np.uint8(cv2.normalize(idftMag, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]
imgRebuild = np.copy(idftMagNorm[:rows, :cols])
print("max(imgGray-imgRebuild) = ", np.max(imgGray-imgRebuild))
print("imgGray:{}, dftImg:{}, idftImg:{}, imgRebuild:{}".
format(imgGray.shape, dftImgEx.shape, idftImg.shape, imgRebuild.shape))
plt.figure(figsize=(9, 6))
plt.subplot(131), plt.title("Original image"), plt.axis('off')
plt.imshow(imgGray, cmap='gray')
plt.subplot(132), plt.title("Log-trans of DFT amp"), plt.axis('off')
dftAmp = cv2.magnitude(dftImgEx[:,:,0], dftImgEx[:,:,1]) # 幅度谱，中心化
dftAmpLog = np.log(1 + dftAmp) # 幅度谱对数变换，以便于显示
plt.imshow(cv2.normalize(dftAmpLog, None, 0, 255, cv2.NORM_MINMAX), cmap='gray')
plt.subplot(133), plt.title("Rebuild image with IDFT"), plt.axis('off')
plt.imshow(imgRebuild, cmap='gray')
plt.tight_layout()
plt.show()
版权声明：
本文中部分原始图片来自 Rafael C. Gonzalez “Digital Image Processing, 4th.Ed.”，特此致谢。
youcans@xupt 原创作品，转载必须标注原文链接：(https://blog.csdn.net/youcans/article/details/120995650)
Copyright 2022 youcans, XUPT
欢迎关注 『youcans 的 OpenCV 学习课』 系列，持续更新
【youcans 的 OpenCV 学习课】1. 安装与环境配置
【youcans 的 OpenCV 学习课】2. 图像读取与显示
【youcans 的 OpenCV 学习课】3. 图像的创建与修改
【youcans 的 OpenCV 学习课】4. 图像的叠加与混合
【youcans 的 OpenCV 学习课】5. 图像的几何变换
【youcans 的 OpenCV 学习课】6. 灰度变换与直方图处理
【youcans 的 OpenCV 学习课】7. 空间域图像滤波
【youcans 的 OpenCV 学习课】8. 频率域图像滤波（上）
【youcans 的 OpenCV 学习课】9. 频率域图像滤波（下）
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
