(十一)OpenCV实现图像频率域滤波

1.基础

见《数字图像处理第四版》P137-P209

1.1傅里叶变换Fourier Transform

Fourier Transform法国的一位数学家和物理学家Jean-Baptiste Joseph Fourier (1768-1830)提出,主要包括适用于周期函数的傅里叶级数(Fourier Serier)和应用于非周期函数(曲线下面积有限)的傅里叶变换。

  • 傅里叶级数
    傅里叶级数是指,周期为 T T T的连续变量 t t t的周期信号 f ( t ) f(t) f(t),可表示为乘以适当系数的正弦函数和余弦函数之和
    借助欧拉公式和欧拉恒等变换

在这里插入图片描述

能将傅里叶级数表示成复数指数的形式,其形式为,

在这里插入图片描述
其中,系数 c n c_n cn可表示为: c n = 1 T ∫ − T / 2 T / 2 f ( t ) e − 2 π n t T d t , n = 0 , ± 1 , ± 2 , . . . c_n = \frac{1}{T}\int_{-T/2}^{T/2} f(t)e^{-\frac{2\pi nt}{T}}dt, n=0,\pm 1,\pm 2,... cn=T1T/2T/2f(t)eT2πntdt,n=0,±1,±2,... j = − 1 j=\sqrt-1 j= 1
从傅里叶级数的的公式可以看出,周期函数被分解成了使用不同频率,幅度,相角的正弦和/或余弦函数表示的形式。能够将信号从时域在频域上展开来进行分析。这正是傅里叶变换的意义,关于其形象化的描述可参考

  • 连续函数的傅里叶变换

    非周期函数可看成周期是无穷大的函数,参考傅里叶级数中 c n c_n cn,其频率和系数幅度之间的关系,即傅里叶变换,可表示为:

    F ( μ ) = ∫ − ∞ ∞ f ( t ) e − j 2 π μ t d t F(\mu) = \int_{-\infty}^{\infty}f(t)e^{-j2\pi\mu t}dt F(μ)=f(t)ej2πμtdt

1.2卷积定理

连续变量t的两个连续函数 f ( t ) f(t) f(t) h ( t ) h(t) h(t)的卷积的定义为:

( f ⋆ h ) ( t ) = ∫ − ∞ ∞ f ( τ ) h ( t − τ ) d τ (f\star h)(t) = \int _{-\infty}^{\infty}f(\tau)h(t-\tau)d\tau (fh)(t)=f(τ)h(tτ)dτ


其中, − τ -\tau τ表示对卷积核函数旋转180度, t t t是一个函数滑过另一个函数的位移, τ \tau τ是积分虚拟变量。上式的傅里叶变换表示为:

在这里插入图片描述
式中方括号内是 h ( t − τ ) h(t-\tau) h(tτ)的傅里叶变换, F ( h ( t − τ ) ) = H ( μ ) e − j 2 π μ τ F(h(t-\tau))=H(\mu)e^{-j2\pi\mu\tau} F(h(tτ))=H(μ)ej2πμτ, H ( μ ) H(\mu) H(μ) h ( t ) h(t) h(t)的傅里叶变换,式(5)可表示为:
在这里插入图片描述式中 ∙ \bullet 表示相乘, ζ \zeta ζ表示傅里叶变换。

由上式可知,空间域中两个函数的卷积的傅里叶变换,等于频域中两个函数的傅里叶变换的乘积。反之,如果有频域中两个傅里叶变换的乘积,可以通过计算傅里叶反变换得到空间域的卷积。综上,卷积定理可以表示为:

( f ⋆ h ) ( t ) ⇔ ( F ∙ H ) ( μ ) (f\star h)(t) \Leftrightarrow (F\bullet H)(\mu) (fh)(t)(FH)(μ)

( f ∙ h ) ( t ) ⇔ ( F ⋆ H ) ( μ ) (f\bullet h)(t) \Leftrightarrow (F\star H)(\mu) (fh)(t)(FH)(μ)

1.3冲激函数及其取样性质

连续变量 t t t t = 0 t=0 t=0处的单位冲激表示为 δ t \delta{t} δt,表示为:

在这里插入图片描述
冲激函数满足如下性质:

  • 满足恒等式: ∫ − ∞ ∞ δ ( t ) d t = 1 \int_{-\infty}^{\infty}\delta(t)dt = 1 δ(t)dt=1
  • 取样性质: ∫ − ∞ ∞ f ( t ) δ ( t − t 0 ) d t = f ( t 0 ) \int_{-\infty}^{\infty}f(t)\delta(t-t_0)dt = f(t_0) f(t)δ(tt0)dt=f(t0),取样性质只是得到 f ( t ) f(t) f(t)在冲激位置 t 0 t_0 t0处的值。

实际对 f ( t ) f(t) f(t)取样时,需要使用一系列冲激信号,以 Δ T \Delta T ΔT为间隔进行取样,由此可定义冲击串:

在这里插入图片描述
上述,将 t t t解释为连续时间变量,对于离散变量,冲激函数定义为:

在这里插入图片描述

满足的性质等效于:

  • 恒等式:在这里插入图片描述- 取样性质:在这里插入图片描述

1.4冲激和冲激串的傅里叶变换

根据连续单变量函数傅里叶变换的定义, δ ( t ) \delta(t) δ(t)函数的傅里叶变换可表示为:

ζ { δ ( t ) } = F ( μ ) = ∫ − ∞ ∞ δ ( t ) e − j 2 π μ t d t = ∫ − ∞ ∞ e − j 2 π μ t δ ( t ) d t = e − j 2 π μ 0 = 1 \zeta \{\delta(t)\}=F(\mu)=\int_{-\infty}^{\infty}\delta(t)e^{-j2\pi \mu t}dt=\int_{-\infty}^{\infty}e^{-j2\pi \mu t}\delta(t)dt=e^{-j2\pi \mu 0}=1 ζ{δ(t)}=F(μ)=δ(t)ej2πμtdt=ej2πμtδ(t)dt=ej2πμ0=1

t = t 0 t=t_0 t=t0处,一个冲激的傅里叶变换为:

ζ { δ ( t − t 0 ) } = F ( μ ) = ∫ − ∞ ∞ δ ( t − t 0 ) e − j 2 π μ t d t = ∫ − ∞ ∞ e − j 2 π μ t δ ( t − t 0 ) d t = e − j 2 π μ t 0 \zeta \{\delta(t-t_0)\}=F(\mu)=\int_{-\infty}^{\infty}\delta(t-t_0)e^{-j2\pi \mu t}dt=\int_{-\infty}^{\infty}e^{-j2\pi \mu t}\delta(t-t_0)dt=e^{-j2\pi \mu t_0} ζ{δ(tt0)}=F(μ)=δ(tt0)ej2πμtdt=ej2πμtδ(tt0)dt=ej2πμt0

为了得到离散信号,通常需要使用周期冲激串来进行采样,因此推导离散信号的傅立叶变换时需要先求得周期冲激串的傅立叶变换。

首先,由上面的定义的傅立叶变换和反变换可知,若函数 f ( t ) f(t) f(t)有傅立叶变换 F ( μ ) F(\mu) F(μ),则函数 F ( t ) F(t) F(t)的傅里叶变换为 f ( − μ ) f(-\mu) f(μ),由此性质,冲激 δ ( t − t 0 ) \delta(t-t_0) δ(tt0)的傅立叶变换为 e − j 2 π μ t 0 e^{-j2\pi \mu t_0} ej2πμt0,则 e − j 2 π μ t 0 e^{-j2\pi \mu t_0} ej2πμt0的傅里叶变换为 δ ( − μ − t 0 ) \delta(-\mu-t_0) δ(μt0)。记 a = − t 0 a=-t_0 a=t0,则 e j 2 π μ a e^{j2\pi \mu a} ej2πμa的傅立叶变换为 δ ( − μ + a ) \delta(-\mu + a) δ(μ+a),对于冲激函数,除 μ = a \mu = a μ=a外, δ \delta δ都为0,因此 δ ( − μ + a ) = δ ( μ − a ) \delta(-\mu + a)=\delta(\mu - a) δ(μ+a)=δ(μa),也即 ζ [ e j 2 π μ a ] = δ ( μ − a ) \zeta[e^{j2\pi \mu a}]=\delta(\mu - a) ζ[ej2πμa]=δ(μa)

由1.3知 S Δ T ( t ) S_{\Delta T}(t) SΔT(t)是周期为 Δ T \Delta T ΔT的周期函数,其傅立叶级数表示为:

在这里插入图片描述
在这里插入图片描述
周期在 [ − Δ T / 2 , Δ T / 2 ] [-\Delta T/2, \Delta T/2] [ΔT/2ΔT/2]上的积分仅包含位于原点的冲激,故 c n c_n cn可以表示为:

c n = 1 Δ T ∫ − Δ T / 2 Δ T / 2 δ ( t ) e − j 2 π n Δ T d t = 1 Δ T e 0 = 1 Δ T c_n=\frac{1}{\Delta T}\int_{-\Delta T/2}^{\Delta T/2}\delta(t)e^{-j\frac{2\pi n}{\Delta T}}dt=\frac{1}{\Delta T}e^0=\frac{1}{\Delta T} cn=ΔT1ΔT/2ΔT/2δ(t)ejΔT2πndt=ΔT1e0=ΔT1

那么,傅立叶级数可表示为:
在这里插入图片描述结合前述 ζ [ e j 2 π μ a ] = δ ( μ − a ) \zeta[e^{j2\pi \mu a}]=\delta(\mu - a) ζ[ej2πμa]=δ(μa),则

ζ [ e j 2 π n Δ T t ] = δ ( μ − n Δ T ) \zeta[e^{j\frac{2\pi n}{\Delta T}t}]=\delta(\mu - \frac{n}{\Delta T}) ζ[ejΔT2πnt]=δ(μΔTn)
结合和的傅里叶变换等于各分量的傅里叶变换之和,则周期冲激串的傅里叶变换为

在这里插入图片描述由上式可知,周期为 Δ T \Delta T ΔT的冲激串的傅里叶变换仍为冲激串,其周期为 1 Δ T \frac{1}{\Delta T} ΔT1,其周期互为反比。

1.5取样后函数的傅里叶变换和取样定理

在这里插入图片描述
使用冲激串取样的过程如上图所示,其数学表示为:

f ~ ( t ) = f ( t ) S Δ T ( t ) = ∑ n = − ∞ ∞ f ( t ) δ ( t − n Δ T ) \tilde f(t)=f(t)S_{\Delta T}(t)=\sum_{n=-\infty}^{\infty}f(t)\delta (t-n\Delta T) f~(t)=f(t)SΔT(t)=n=f(t)δ(tnΔT)
f ~ ( t ) \tilde f(t) f~(t)表示取样后的函数,取样后,取样序列中任意一个取样的值 f k f_k fk由下公式给出:

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)δ(tkΔT)dt=f(kΔT)

取样后函数的傅里叶变换表示:
F ( μ ) F(\mu) F(μ)表示连续函数 f ( t ) f(t) f(t)函数的傅里叶变换,取样后的函数 f ~ ( t ) \tilde f(t) f~(t) f ( t ) f(t) f(t)与冲激串的乘积。由卷积定理,空间域中两个函数乘积的傅里叶变换,是频率域中两个变换的卷积。因此 f ~ ( t ) \tilde f(t) f~(t)的傅里叶变换 F ~ ( μ ) \tilde F(\mu) F~(μ)可表示为:

F ~ ( μ ) = ζ { f ~ ( t ) } = ζ { f ( t ) S Δ T ( t ) } = ( F ⋆ S ) ( μ ) \tilde F(\mu) = \zeta\{\tilde f(t)\}=\zeta\{f(t)S_{\Delta T}(t)\}=(F\star S)(\mu) F~(μ)=ζ{f~(t)}=ζ{f(t)SΔT(t)}=(FS)(μ)
由上面的推导可以知道冲激串的傅里叶变换为:
在这里插入图片描述
由一维卷积的定义,则 F ( μ ) F(\mu) F(μ) S ( μ ) S(\mu) S(μ)的卷积可以表示为:

在这里插入图片描述
由上述推导可以知,取样后函数 f ~ ( t ) \tilde f(t) f~(t)的傅里叶变换 F ~ ( μ ) \tilde F(\mu) F~(μ)是原函数傅里叶变换的一个无限的周期的副本序列。考虑一个带限函数的傅里叶变换示意图:

在这里插入图片描述

上图中,图a是函数 f ( t ) f(t) f(t)的傅里叶变换 F ( μ ) F(\mu) F(μ),图b是函数 f ~ ( t ) \tilde f(t) f~(t)的傅里叶变换 F ~ ( μ ) \tilde F(\mu) F~(μ) 1 / Δ T 1/\Delta T 1/ΔT是取样函数的取样率,从上图可以看到取样率要高到足以在各个周期之间提供有效的间隔,以保持 F ( μ ) F(\mu) F(μ)的完整性。图c和图d分别对应的是临界取样和欠取样。

对以原点为中心的有限区间 [ − μ m a x , μ m a x ] [-\mu_{max},\mu_{max}] [μmax,μmax]外的频率值,傅里叶变换为零的函数 f ( t ) f(t) f(t)称为带限函数。从上图中可知,当 Δ T \Delta T ΔT较大时, F ~ ( μ ) \tilde F(\mu) F~(μ)会发生交叠,将不能从 F ~ ( μ ) \tilde F(\mu) F~(μ)中分离完整的 F ( μ ) F(\mu) F(μ),即无法通过傅里叶反变换复原原信号。因此为了复原原信号,必须提高采样频率,如下图:
在这里插入图片描述
1 / 2 Δ T > μ m a x 1/2\Delta T\gt\mu_{max} 1/2ΔT>μmax时可以保证足够的间隔: 1 Δ T > 2 μ m a x \frac{1}{\Delta T}\gt2\mu_{max} ΔT1>2μmax,该式表明,如果以超过函数最高频率2倍的取样率来得到样本,那么连续带限函数就能由其样本集合复原,该结论即取样定理

1.6由取样后的数据复原函数

定义函数 H ( μ ) H(\mu) H(μ)为:
在这里插入图片描述 F ( μ ) F(\mu) F(μ)为: F ( μ ) = H ( μ ) F ~ ( μ ) F(\mu)=H(\mu)\tilde F(\mu) F(μ)=H(μ)F~(μ),得到 F ( μ ) F(\mu) F(μ)后可用傅里叶反变换复原 f ( t ) f(t) f(t):
f ( t ) = ∫ − ∞ ∞ F ( μ ) e j 2 π μ t f(t)=\int_{-\infty}^{\infty}F(\mu)e^{j2\pi\mu t} f(t)=F(μ)ej2πμt

F ( μ ) = H ( μ ) F ~ ( μ ) F(\mu)=H(\mu)\tilde F(\mu) F(μ)=H(μ)F~(μ)可得:
f ( t ) = ζ − 1 { F ( μ ) } = ζ − 1 { H ( μ ) F ~ ( μ ) } = h ( t ) ⋆ f ~ ( t ) f(t)=\zeta^{-1}\{F(\mu)\}=\zeta^{-1}\{H(\mu)\tilde F(\mu)\}=h(t)\star\tilde f(t) f(t)=ζ1{F(μ)}=ζ1{H(μ)F~(μ)}=h(t)f~(t)
即频率域对函数的处理等同于空间域中的卷积。

2.离散傅里叶变换(Discrete Fourier Transform)

2.1一维离散傅里叶变换

前面讨论采样定理时,求出了采样后函数 f ~ ( t ) \tilde f(t) f~(t)的傅里叶变换 F ~ ( μ ) \tilde F(\mu) F~(μ)与原函数 f ( t ) f(t) f(t)的傅里叶变换 F ( μ ) F(\mu) F(μ)之间的关系,但是并没有给出 F ~ ( μ ) \tilde F(\mu) F~(μ)的表达式,这里利用冲激串的定义来求:
在这里插入图片描述
F ~ ( μ ) \tilde F(\mu) F~(μ F ( μ ) F(\mu) F(μ的关系,可知, F ~ ( μ ) \tilde F(\mu) F~(μ是周期为 1 / Δ T 1/\Delta T 1/ΔT的无限周期连续函数,因此表征 F ~ ( μ ) \tilde F(\mu) F~(μ只需要一个周期。假设从 μ = 0 \mu=0 μ=0 μ = 1 / Δ T \mu=1/\Delta T μ=1/ΔT的一个周期内等间隔的取 F ~ ( μ ) \tilde F(\mu) F~(μ的M个样本,在如下频率处取样可实现: μ = m M Δ T , m = 0 , 1 , 2 , . . . , M − 1 \mu=\frac{m}{M\Delta T}, m=0,1,2,...,M-1 μ=MΔTm,m=0,1,2,...,M1 μ \mu μ代入上面的公式中可得 F ~ ( μ ) \tilde F(\mu) F~(μ)的表达式为:

在这里插入图片描述
已知一个由f(t)的M个样本组成的集合 { f m } \{f_m\} {fm}时,上式给出了一个与输入样本集合的离散傅里叶变换相对应的M个复值集合 { F m } \{F_m\} {Fm}。反之,求得 { F m } \{F_m\} {Fm}时也可由傅里叶反变换求得 { f m } \{f_m\} {fm}

在这里插入图片描述在图像中通常使用 x x x, y y y表示图像坐标,并使用 μ \mu μ, v v v表示频率变量,这里的这些变量都理解成整数,则上述傅里叶变换对可表示为:

在这里插入图片描述

2.2二维函数的傅里叶变换:

将前述概念扩充到二维上。

二维冲激函数: δ ( t , z ) { = 1 , t = z = 0 0 , o t h e r s , 且 ∫ − ∞ ∞ ∫ − ∞ ∞ δ ( t , z ) d t d z = 1 \delta(t,z)\begin{cases} =1,t=z=0\\ 0,others \end{cases}, 且\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\delta(t,z)dtdz = 1 δ(t,z){=1,t=z=00,others,δ(t,z)dtdz=1

二维冲激函数的取样性质:

∫ − ∞ ∞ ∫ − ∞ ∞ f ( t , z ) δ ( t − t 0 , z − z 0 ) d t d z = f ( t 0 , z 0 ) \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(t,z)\delta(t-t_0,z-z_0)dtdz = f(t_0,z_0) f(t,z)δ(tt0,zz0)dtdz=f(t0,z0)

二维连续函数的傅里叶变换对:

F ( μ , v ) = ∫ − ∞ ∞ ∫ − ∞ ∞ f ( t , z ) e − j 2 π ( μ t + v z ) d t d z F(\mu,v) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(t,z)e^{-j2\pi(\mu t+vz)}dtdz F(μ,v)=f(t,z)ej2π(μt+vz)dtdz f ( t , z ) = ∫ − ∞ ∞ ∫ − ∞ ∞ F ( μ , v ) e j 2 π ( μ t + v z ) d μ d v f(t,z) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(\mu,v)e^{j2\pi(\mu t+vz)}d\mu d v f(t,z)=F(μ,v)ej2π(μt+vz)dμdv

二维冲激串(二维取样函数):
在这里插入图片描述
二维取样定理:

{ Δ T < 1 2 μ m a x Δ Z < 1 2 v m a x \begin{cases} \Delta T \lt \frac{1}{2\mu_{max}} \\ \Delta Z \lt \frac{1}{2 v_{max}} \end{cases} {ΔT<2μmax1ΔZ<2vmax1

二维信号的混叠:

在这里插入图片描述二维离散傅里叶变换对:

在这里插入图片描述
二维离散傅里叶变换的平移旋转和周期性:

由定义可以求得:

  • 平移:

    f ( x , y ) e j 2 π ( μ 0 x / M + v 0 y / M ) ⇔ F ( μ − μ 0 , v − v 0 ) f(x,y)e^{j2\pi(\mu_0x/M+v_0 y/M)}\Leftrightarrow F(\mu -\mu_0, v- v_0) f(x,y)ej2π(μ0x/M+v0y/M)F(μμ0,vv0)

    f ( x − x 0 , y − y 0 ) ⇔ F ( μ − μ 0 , v − v 0 ) e − j 2 π ( x 0 μ / M + y 0 v / M ) f(x-x_0,y-y_0)\Leftrightarrow F(\mu -\mu_0, v- v_0)e^{-j2\pi(x_0\mu/M+y_0 v/M)} f(xx0,yy0)F(μμ0,vv0)ej2π(x0μ/M+y0v/M)

    用所示的指数项乘以 f ( x , y ) f(x,y) f(x,y)将使DFT的原点移动到 ( μ 0 , v 0 ) (\mu_0,v_0) (μ0,v0)处,用负指数乘以 F ( μ , v ) F(\mu, v) F(μ,v)将使 f ( x , y ) f(x,y) f(x,y)的原点移动到点 f ( x 0 , y 0 ) f(x_0,y_0) f(x0,y0)处。平移不影响 F ( μ , v ) F(\mu,v) F(μ,v)的幅度

  • 旋转:

    使用极坐标: x = r c o s θ , y = r s i n θ , μ = ω c o s φ , v = s i n φ x=rcos\theta,y=rsin \theta, \mu=\omega cos \varphi, v = sin \varphi x=rcosθ,y=rsinθ,μ=ωcosφ,v=sinφ,则有如下变换对:$f(r,\theta+\theta_0)\Leftrightarrow F(\omega, \varphi + \theta_0) 即 , 若 f ( x , y ) 旋 转 即,若f(x,y)旋转 f(x,y)\theta_0 角 度 , 则 角度,则 F(\mu,v) 也 将 旋 转 相 同 的 角 度 。 反 之 , 若 也将旋转相同的角度。反之,若 F(\mu,v) 旋 转 某 个 角 度 , 旋转某个角度, f(x,y)$也旋转相同的角度。

  • 周期性:

    同一维情况相同,二维离散傅里叶变换及其反变换在 μ \mu μ v v v方向是无限周期的,即:

    F ( μ , v ) = F ( u + k 1 M , v ) = F ( μ , v + k 2 N ) = F ( u + k 1 M , v + k 2 N ) F(\mu,v)=F(u+k_1M,v)=F(\mu,v+k_2N)=F(u+k_1M,v+k_2N) F(μ,v)=F(u+k1M,v)=F(μ,v+k2N)=F(u+k1M,v+k2N)


    f ( x , y ) = f ( x + k 1 M , y ) = f ( x , y + k 2 N ) = f ( x + k 1 M , y + k 2 N ) f(x,y)=f(x+k_1M,y)=f(x,y+k_2N)=f(x+k_1M,y+k_2N) f(x,y)=f(x+k1M,y)=f(x,y+k2N)=f(x+k1M,y+k2N)

    从第一部分的分析,可发现,区间0到M-1上的变换数据由在点 M / 2 M/2 M/2处相遇的两个半周期组成,可以发现该周期内较第部分的 F ( μ ) F(\mu) F(μ)出现在较高的频率处,不便于分析和滤波。根据上述的平移性质,

    f ( x ) e j 2 π ( μ 0 x / M ) ⇔ F ( μ − μ 0 ) f(x)e^{j2\pi (\mu_0x/M)}\Leftrightarrow F(\mu -\mu_0) f(x)ej2π(μ0x/M)F(μμ0)
    该式将把数据的原点 F ( 0 ) F(0) F(0)平移到 μ 0 \mu_0 μ0,若另 μ 0 = M / 2 \mu_0=M/2 μ0=M/2,则指数项变为 e j π x e^{j\pi x} ejπx,因为x是整数,故根据欧拉公式,其可表示为 ( − 1 ) x (-1)^x (1)x,代入上式可得:

    f ( x ) ( − 1 ) x ⇔ F ( μ − M / 2 ) f(x)(-1)^x\Leftrightarrow F(\mu - M/2) f(x)(1)xF(μM/2)

    如下图所示。
    在这里插入图片描述 对于二维情况可以表示为:

    f ( x , y ) ( − 1 ) x + y ⇔ F ( μ − M / 2 , v − N / 2 ) f(x,y)(-1)^{x+y}\Leftrightarrow F(\mu - M/2, v-N/2) f(x,y)(1)x+yF(μM/2,vN/2)

频谱,相谱和功率谱:

从前面的分析可知,二维DFT通常是复函数,用极坐标形式表示为:

F ( μ , v ) = R ( μ , v ) + j I ( μ , v ) = ∣ F ( μ , v ) ∣ e j ϕ ( μ , v ) F(\mu,v)=R(\mu,v)+jI(\mu,v)=|F(\mu,v)|e^{j\phi(\mu,v)} F(μ,v)=R(μ,v)+jI(μ,v)=F(μ,v)ejϕ(μ,v)

傅里叶频谱:
∣ F ( μ , v ) ∣ = [ R 2 ( μ , v ) + I 2 ( μ , v ) ] 1 / 2 |F(\mu,v)|=[R^2(\mu,v)+I^2(\mu,v)]^{1/2} F(μ,v)=[R2(μ,v)+I2(μ,v)]1/2

相谱: ϕ ( μ , v ) = a r c t a n [ I ( μ , v ) R ( μ , v ) ] \phi(\mu,v)=arctan[\frac{I(\mu,v)}{R(\mu,v)}] ϕ(μ,v)=arctan[R(μ,v)I(μ,v)]

功率谱: P ( μ , v ) = ∣ F ( μ , v ) ∣ 2 = R 2 ( μ , v ) + I 2 ( μ , v ) P(\mu,v)=|F(\mu,v)|^2 = R^2(\mu,v) + I^2(\mu,v) P(μ,v)=F(μ,v)2=R2(μ,v)+I2(μ,v)

一个图像的频谱示意图如下图:
在这里插入图片描述

图A是原图像,图B是直接对原图像进行傅里叶变换得到的频谱,可以看到四个角较亮,但因为其值分散在四个角落,因此很难观察,使用前述的变换,对原函数 f ( x , y ) f(x,y) f(x,y)乘以 ( − 1 ) x + y (-1)^{x+y} (1)x+y将其变换的中心移动到图像的中心可得图C,更便于观察分析。图D是对频谱值取对数后的结果。

2.2二维离散卷积定理

将一维离散卷积定理扩展到两个变量,可得:

在这里插入图片描述二维卷积定理可表示为:

( f ⋆ h ) ( x , y ) ⇔ ( F ∙ H ) ( μ , v ) (f \star h)(x,y) \Leftrightarrow (F\bullet H)(\mu,v) (fh)(x,y)(FH)(μ,v)
使用卷积定理时的问题:

如果使用DFT和卷积定理来求得如下图左列相同的卷积结果时,需考虑DFT表达式中固有的周期性。
在这里插入图片描述
上图中左列,是根据定义使用尺寸为400的卷积核对尺寸为400的原信号进行卷积,先对卷积核函数旋转180度,再依次滑过原信号,得到图左列最后所示的卷积结果。

若使用卷积定理,则f,h都具有周期性,若根据卷积定理,先计算f和h的DFT,然后让这个变换相乘,再计算傅里叶反变换,将得到如上图右侧所示的结果,显然如此计算的卷积结果是不对的。

交叠错误很容易解决。考虑f(x)和h(x),他们分别由A个样本和B个样本组成,可证明,若在这两个函数中填充零,使他们长度P相同,则选择

P ≥ A + B − 1 P\ge A+B-1 PA+B1

可避免交叠错误。

f ( x , y ) f(x,y) f(x,y) h ( x , y ) h(x,y) h(x,y)分别表示大小为 A X B AXB AXB像素和 C X D CXD CXD像素的图像阵列.它们循环卷积中的交叠错误可通过对这两个函数填充零来避免,
f p ( x , y ) = { f ( x , y ) , 0 ≤ x ≤ A − 1 和 0 ≤ y ≤ B − 1 0 , A ≤ x ≤ P 或 B ≤ y ≤ Q f_p(x,y) = \begin{cases} 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{cases} fp(x,y)={f(x,y),0xA10yB10,AxPByQ

h p ( x , y ) = { h ( x , y ) , 0 ≤ C − 1 和 0 ≤ y ≤ D − 1 0 , C ≤ x ≤ P 或 D ≤ y ≤ Q h_p(x,y)=\begin{cases} h(x,y), 0 \le C-1和 0 \le y \le D-1\\ 0, C \le x \le P 或 D \le y \le Q \end{cases} hp(x,y)={h(x,y),0C10yD10,CxPDyQ
其中, P ≥ A + C − 1 P \ge A+C-1 PA+C1, Q ≥ B + D − 1 Q \ge B+D-1 QB+D1

填充零后图像大小为 P X Q PXQ PXQ,若两个阵列大小相同,都为 M X N MXN MXN,则要求 P ≥ 2 M − 1 P\ge2M-1 P2M1 Q ≥ 2 N − 1 Q\ge2N-1 Q2N1.DFT算法执行偶数大小的阵列时速度较快,因此通常选择满足上述公式的最小偶整数作为P和Q的值,若两个阵列相同,则 P = 2 M P=2M P=2M Q = 2 N Q=2N Q=2N

3.频率域滤波步骤

  • 1)一幅大小为 M X N MXN MXN的输入图像 f ( x , y ) f(x,y) f(x,y),由前述公式得到填充零后的尺寸P=2M和Q=2N

  • 2)使用零填充,镜像填充或者复制填充,形成大小为 P X Q PXQ PXQ的填充后的图像 f P ( x , y ) f_P(x,y) fP(x,y)

  • 3)根据前述二维傅里叶变换的平移性质,将 f P ( x , y ) f_P(x,y) fP(x,y)乘以 ( − 1 ) x + y (-1)^{x+y} (1)x+y,使傅里叶变换位于 P X Q PXQ PXQ大小的频率矩形的中心

  • 4)计算步骤3得到的图像的DFT,即 F ( μ , v ) F(\mu,v) F(μ,v)

  • 5)建一个实对称滤波器传递函数 H ( μ , v ) H(\mu, v) H(μ,v),其大小为 P X Q PXQ PXQ,中心在 ( P / 2 , Q / 2 ) (P/2,Q/2) (P/2,Q/2)

  • 6)采用对应像素相乘得到 G ( μ , v ) = H ( μ , v ) F ( μ , v ) G(\mu,v)=H(\mu,v)F(\mu,v) G(μ,v)=H(μ,v)F(μ,v),即 G ( i , k ) = H ( i , k ) F ( i . k ) , i = 0 , 1 , 2 , . . . , M − 1 G(i,k)=H(i,k)F(i.k),i=0,1,2,...,M-1 G(i,k)=H(i,k)F(i.k),i=0,1,2,...,M1 k = 0 , 1 , 2 , . . . , N − 1 k=0,1,2,...,N-1 k=0,1,2,...,N1

  • 7)计算 G ( μ , v ) G(\mu, v) G(μ,v)的IDFT得到滤波后的图像(大小为PXQ): g p ( x , y ) = { r e a l [ ζ − 1 [ G ( μ , v ) ] ] } ( − 1 ) x + y g_p(x,y)=\{real[\zeta^{-1}[G(\mu,v)]]\}(-1)^{x+y} gp(x,y)={real[ζ1[G(μ,v)]]}(1)x+y

  • 8)最后,从 g P ( x , y ) g_P(x,y) gP(x,y)的左上象限提取一个大小为 M X N MXN MXN的区域,得到与输入图像大小相同的滤波结果 g ( x , y ) g(x,y) g(x,y)

4.OpenCV API

  • int cv::getOptimalDFTSize ( int vecsize )
    在进行傅里叶变换时,需要对原图像进行扩充以提高处理速度,该方法正是为了获取满足最优处理速度的最小扩充
  • void cv::copyMakeBorder ( InputArray src, OutputArray dst, int top, int bottom, int left, int right, int borderType, const Scalar & value = Scalar() )对图像进行扩充,可设置上下左右的扩充大小和类型
  • void cv::merge ( const Mat * mv, size_t count, OutputArray dst ) M a t Mat Mat数组mv合并成count channel的dst,与split刚好相反
  • void cv::dft ( InputArray src, OutputArray dst, int flags = 0, int nonzeroRows = 0 )执行1d或2d离散傅里叶变换
  • void cv::magnitude ( InputArray x, InputArray y, OutputArray magnitude )计算 x x x y y y的模, 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)

示例:

// https://docs.opencv.org/4.x/d8/d01/tutorial_discrete_fourier_transform.html
#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"

#include <iostream>

using namespace cv;
using namespace std;

static void help(char ** argv)
{
    cout << endl
        <<  "This program demonstrated the use of the discrete Fourier transform (DFT). " << endl
        <<  "The dft of an image is taken and it's power spectrum is displayed."  << endl << endl
        <<  "Usage:"                                                                      << endl
        << argv[0] << " [image_name -- default lena.jpg]" << endl << endl;
}

int main(int argc, char ** argv)
{
    help(argv);

    const char* filename = argc >=2 ? argv[1] : "lena.jpg";

    Mat I = imread( samples::findFile( filename ), IMREAD_GRAYSCALE);
    if( I.empty()){
        cout << "Error opening image" << endl;
        return EXIT_FAILURE;
    }

    //! [expand]
    Mat padded;                            //expand input image to optimal size
    int m = getOptimalDFTSize( I.rows );
    int n = getOptimalDFTSize( I.cols ); // on the border add zero values
    copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
    //! [expand]

    //! [complex_and_real]
    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
    merge(planes, 2, complexI);         // Add to the expanded another plane with zeros
    //! [complex_and_real]

    //! [dft]
    dft(complexI, complexI);            // this way the result may fit in the source matrix
    //! [dft]

    // compute the magnitude and switch to logarithmic scale
    // => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
    //! [magnitude]
    split(complexI, planes);                   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
    magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
    Mat magI = planes[0];
    //! [magnitude]

    //! [log]
    magI += Scalar::all(1);                    // switch to logarithmic scale
    log(magI, magI);
    //! [log]

    //! [crop_rearrange]
    // crop the spectrum, if it has an odd number of rows or columns
    magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

    // rearrange the quadrants of Fourier image  so that the origin is at the image center
    int cx = magI.cols/2;
    int cy = magI.rows/2;

    Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
    Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
    Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
    Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right

    Mat tmp;                           // swap quadrants (Top-Left with Bottom-Right)
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)
    q2.copyTo(q1);
    tmp.copyTo(q2);
    //! [crop_rearrange]

    //! [normalize]
    normalize(magI, magI, 0, 1, NORM_MINMAX); // Transform the matrix with float values into a
                                            // viewable image form (float between values 0 and 1).
    //! [normalize]
    imshow("Input Image"       , I   );    // Show the result
    imshow("spectrum magnitude", magI);
    waitKey();

    return EXIT_SUCCESS;
}

结果:
在这里插入图片描述
REF:

1.https://zhuanlan.zhihu.com/p/351173878
2.https://zhuanlan.zhihu.com/p/19763358
3.https://zhuanlan.zhihu.com/p/31371519
4.https://docs.opencv.org/4.x/d8/d01/tutorial_discrete_fourier_transform.html
5.https://www.cnblogs.com/HL-space/p/10546602.html

  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
opencv中的频域滤波使用傅里叶变换(FFT)进行操作。以下是一些常见的频域滤波操作: 1. 频域低通滤波:通过去除高频部分来平滑图像,常用的方法有理想低通滤波、巴特沃斯低通滤波和高斯低通滤波。 2. 频域高通滤波:通过去除低频部分来增强图像的边缘和细节,常用的方法有理想高通滤波、巴特沃斯高通滤波和高斯高通滤波。 3. 频域带通滤波:通过保留一定频率范围内的信息来增强图像的某些频率成分。 在OpenCV中,可以使用dft函数进行傅里叶变换,使用idft函数进行逆变换。可以使用getOptimalDFTSize函数来获取最佳尺寸。下面是一个示例代码,演示如何进行频域滤波: ```python import cv2 import numpy as np # 读取图像 image = cv2.imread('image.jpg', 0) # 将图像转换为频域 dft = cv2.dft(np.float32(image), flags=cv2.DFT_COMPLEX_OUTPUT) dft_shift = np.fft.fftshift(dft) # 创建一个掩模,中心为1,其余为0 rows, cols = image.shape mask = np.zeros((rows, cols, 2), np.uint8) mask[int(rows/2)-30:int(rows/2)+30, int(cols/2)-30:int(cols/2)+30] = 1 # 应用滤波器 dft_shift = dft_shift * mask # 将频域图像转换回空域图像 f_ishift = np.fft.ifftshift(dft_shift) image_filtered = cv2.idft(f_ishift) image_filtered = cv2.magnitude(image_filtered[:, :, 0], image_filtered[:, :, 1]) # 显示原始图像滤波后的图像 cv2.imshow('Original', image) cv2.imshow('Filtered', image_filtered) cv2.waitKey(0) cv2.destroyAllWindows() ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值