(十一)OpenCV实现图像频率域滤波
1.基础
见《数字图像处理第四版》P137-P209
1.1傅里叶变换Fourier Transform
Fourier Transform由法国的一位数学家和物理学家Jean-Baptiste Joseph Fourier (1768-1830)提出，主要包括适用于周期函数的傅里叶级数(Fourier Serier)和应用于非周期函数(曲线下面积有限)的傅里叶变换。
- 傅里叶级数
傅里叶级数是指，周期为TTT的连续变量ttt的周期信号f(t)f(t)f(t),可表示为乘以适当系数的正弦函数和余弦函数之和。
借助欧拉公式和欧拉恒等变换，
能将傅里叶级数表示成复数指数的形式,其形式为，
其中，系数cnc_ncn可表示为:cn=1T∫−T/2T/2f(t)e−2πntTdt,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=T1∫−T/2T/2f(t)e−T2πntdt,n=0,±1,±2,...，j=−1j=\sqrt-1j=−1
从傅里叶级数的的公式可以看出，周期函数被分解成了使用不同频率，幅度，相角的正弦和/或余弦函数表示的形式。能够将信号从时域在频域上展开来进行分析。这正是傅里叶变换的意义，关于其形象化的描述可参考
- 连续函数的傅里叶变换
非周期函数可看成周期是无穷大的函数，参考傅里叶级数中cnc_ncn，其频率和系数幅度之间的关系，即傅里叶变换，可表示为：
F(μ)=∫−∞∞f(t)e−j2πμtdtF(\mu) = \int_{-\infty}^{\infty}f(t)e^{-j2\pi\mu t}dtF(μ)=∫−∞∞f(t)e−j2πμ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(f⋆h)(t)=∫−∞∞f(τ)h(t−τ)dτ
其中，−τ-\tau−τ表示对卷积核函数旋转180度，ttt是一个函数滑过另一个函数的位移，τ\tauτ是积分虚拟变量。上式的傅里叶变换表示为：
式中方括号内是h(t−τ)h(t-\tau)h(t−τ)的傅里叶变换，F(h(t−τ))=H(μ)e−j2πμτF(h(t-\tau))=H(\mu)e^{-j2\pi\mu\tau}F(h(t−τ))=H(μ)e−j2πμτ,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)(f⋆h)(t)⇔(F∙H)(μ)
(f∙h)(t)⇔(F⋆H)(μ)(f\bullet h)(t) \Leftrightarrow (F\star H)(\mu)(f∙h)(t)⇔(F⋆H)(μ)
1.3冲激函数及其取样性质
连续变量ttt在t=0t=0t=0处的单位冲激表示为δt\delta{t}δt,表示为：
冲激函数满足如下性质：
- 满足恒等式：∫−∞∞δ(t)dt=1\int_{-\infty}^{\infty}\delta(t)dt = 1∫−∞∞δ(t)dt=1
- 取样性质：∫−∞∞f(t)δ(t−t0)dt=f(t0)\int_{-\infty}^{\infty}f(t)\delta(t-t_0)dt = f(t_0)∫−∞∞f(t)δ(t−t0)dt=f(t0)，取样性质只是得到f(t)f(t)f(t)在冲激位置t0t_0t0处的值。
实际对f(t)f(t)f(t)取样时，需要使用一系列冲激信号，以ΔT\Delta TΔT为间隔进行取样，由此可定义冲击串：
上述，将ttt解释为连续时间变量,对于离散变量，冲激函数定义为：
满足的性质等效于：
- 恒等式：- 取样性质：
1.4冲激和冲激串的傅里叶变换
根据连续单变量函数傅里叶变换的定义，δ(t)\delta(t)δ(t)函数的傅里叶变换可表示为：
ζ{δ(t)}=F(μ)=∫−∞∞δ(t)e−j2πμtdt=∫−∞∞e−j2πμtδ(t)dt=e−j2πμ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)e−j2πμtdt=∫−∞∞e−j2πμtδ(t)dt=e−j2πμ0=1
在t=t0t=t_0t=t0处，一个冲激的傅里叶变换为：
ζ{δ(t−t0)}=F(μ)=∫−∞∞δ(t−t0)e−j2πμtdt=∫−∞∞e−j2πμtδ(t−t0)dt=e−j2πμt0\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}ζ{δ(t−t0)}=F(μ)=∫−∞∞δ(t−t0)e−j2πμtdt=∫−∞∞e−j2πμtδ(t−t0)dt=e−j2πμt0
为了得到离散信号，通常需要使用周期冲激串来进行采样，因此推导离散信号的傅立叶变换时需要先求得周期冲激串的傅立叶变换。
首先，由上面的定义的傅立叶变换和反变换可知，若函数f(t)f(t)f(t)有傅立叶变换F(μ)F(\mu)F(μ),则函数F(t)F(t)F(t)的傅里叶变换为f(−μ)f(-\mu)f(−μ),由此性质,冲激δ(t−t0)\delta(t-t_0)δ(t−t0)的傅立叶变换为e−j2πμt0e^{-j2\pi \mu t_0}e−j2πμt0，则e−j2πμt0e^{-j2\pi \mu t_0}e−j2πμt0的傅里叶变换为δ(−μ−t0)\delta(-\mu-t_0)δ(−μ−t0)。记a=−t0a=-t_0a=−t0,则ej2πμae^{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)，也即ζ[ej2πμ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]上的积分仅包含位于原点的冲激，故cnc_ncn可以表示为：
cn=1ΔT∫−ΔT/2ΔT/2δ(t)e−j2πnΔTdt=1ΔTe0=1ΔTc_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)e−jΔT2πndt=ΔT1e0=ΔT1
那么，傅立叶级数可表示为：
结合前述ζ[ej2πμa]=δ(μ−a)\zeta[e^{j2\pi \mu a}]=\delta(\mu - a)ζ[ej2πμa]=δ(μ−a),则
ζ[ej2πnΔTt]=δ(μ−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)δ(t−nΔT)
f~(t)\tilde f(t)f~(t)表示取样后的函数，取样后，取样序列中任意一个取样的值fkf_kfk由下公式给出：
fk=∫−∞∞f(t)δ(t−kΔT)dt=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(\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)}=(F⋆S)(μ)
由上面的推导可以知道冲激串的傅里叶变换为：
由一维卷积的定义，则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/ΔT1/\Delta T1/ΔT是取样函数的取样率，从上图可以看到取样率要高到足以在各个周期之间提供有效的间隔，以保持F(μ)F(\mu)F(μ)的完整性。图c和图d分别对应的是临界取样和欠取样。
对以原点为中心的有限区间[−μmax,μmax][-\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>μmax1/2\Delta T\gt\mu_{max}1/2ΔT>μmax时可以保证足够的间隔：1ΔT>2μmax\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(μ)ej2πμtf(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/ΔT1/\Delta T1/ΔT的无限周期连续函数，因此表征F~(μ）\tilde F(\mu）F~(μ）只需要一个周期。假设从μ=0\mu=0μ=0到μ=1/ΔT\mu=1/\Delta Tμ=1/ΔT的一个周期内等间隔的取F~(μ）\tilde F(\mu）F~(μ）的M个样本，在如下频率处取样可实现：μ=mMΔ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,...,M−1把μ\muμ代入上面的公式中可得F~(μ)\tilde F(\mu)F~(μ)的表达式为:
已知一个由f(t)的M个样本组成的集合{fm}\{f_m\}{fm}时,上式给出了一个与输入样本集合的离散傅里叶变换相对应的M个复值集合{Fm}\{F_m\}{Fm}。反之，求得{Fm}\{F_m\}{Fm}时也可由傅里叶反变换求得{fm}\{f_m\}{fm}
在图像中通常使用xxx,yyy表示图像坐标，并使用μ\muμ,vvv表示频率变量，这里的这些变量都理解成整数，则上述傅里叶变换对可表示为：
2.2二维函数的傅里叶变换：
将前述概念扩充到二维上。
二维冲激函数： δ(t,z){=1,t=z=00,others,且∫−∞∞∫−∞∞δ(t,z)dtdz=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−t0,z−z0)dtdz=f(t0,z0)\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)δ(t−t0,z−z0)dtdz=f(t0,z0)
二维连续函数的傅里叶变换对：
F(μ,v)=∫−∞∞∫−∞∞f(t,z)e−j2π(μt+vz)dtdzF(\mu,v) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(t,z)e^{-j2\pi(\mu t+vz)}dtdzF(μ,v)=∫−∞∞∫−∞∞f(t,z)e−j2π(μt+vz)dtdz和f(t,z)=∫−∞∞∫−∞∞F(μ,v)ej2π(μt+vz)dμdvf(t,z) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(\mu,v)e^{j2\pi(\mu t+vz)}d\mu d vf(t,z)=∫−∞∞∫−∞∞F(μ,v)ej2π(μt+vz)dμdv
二维冲激串(二维取样函数)：
二维取样定理：
{ΔT<12μmaxΔZ<12vmax\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)ej2π(μ0x/M+v0y/M)⇔F(μ−μ0,v−v0)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,v−v0)
f(x−x0,y−y0)⇔F(μ−μ0,v−v0)e−j2π(x0μ/M+y0v/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(x−x0,y−y0)⇔F(μ−μ0,v−v0)e−j2π(x0μ/M+y0v/M)
用所示的指数项乘以f(x,y)f(x,y)f(x,y)将使DFT的原点移动到(μ0,v0)(\mu_0,v_0)(μ0,v0)处，用负指数乘以F(μ,v)F(\mu, v)F(μ,v)将使f(x,y)f(x,y)f(x,y)的原点移动到点f(x0,y0)f(x_0,y_0)f(x0,y0)处。平移不影响F(μ,v)F(\mu,v)F(μ,v)的幅度
旋转:
使用极坐标：x=rcosθ,y=rsinθ,μ=ωcosφ,v=sinφx=rcos\theta,y=rsin \theta, \mu=\omega cos \varphi, v = sin \varphix=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μ和vvv方向是无限周期的，即：
F(μ,v)=F(u+k1M,v)=F(μ,v+k2N)=F(u+k1M,v+k2N)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+k1M,y)=f(x,y+k2N)=f(x+k1M,y+k2N)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/2M/2M/2处相遇的两个半周期组成，可以发现该周期内较第部分的F(μ)F(\mu)F(μ)出现在较高的频率处，不便于分析和滤波。根据上述的平移性质，
f(x)ej2π(μ0x/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,则指数项变为ejπxe^{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)x⇔F(μ−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+y⇔F(μ−M/2,v−N/2)
频谱，相谱和功率谱：
从前面的分析可知，二维DFT通常是复函数，用极坐标形式表示为：
F(μ,v)=R(μ,v)+jI(μ,v)=∣F(μ,v)∣ejϕ(μ,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)∣=[R2(μ,v)+I2(μ,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)=arctan[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=R2(μ,v)+I2(μ,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)(f⋆h)(x,y)⇔(F∙H)(μ,v)
使用卷积定理时的问题：
如果使用DFT和卷积定理来求得如下图左列相同的卷积结果时，需考虑DFT表达式中固有的周期性。
上图中左列，是根据定义使用尺寸为400的卷积核对尺寸为400的原信号进行卷积，先对卷积核函数旋转180度，再依次滑过原信号，得到图左列最后所示的卷积结果。
若使用卷积定理，则f,h都具有周期性，若根据卷积定理，先计算f和h的DFT，然后让这个变换相乘，再计算傅里叶反变换，将得到如上图右侧所示的结果，显然如此计算的卷积结果是不对的。
交叠错误很容易解决。考虑f(x)和h(x),他们分别由A个样本和B个样本组成，可证明，若在这两个函数中填充零，使他们长度P相同，则选择
P≥A+B−1P\ge A+B-1P≥A+B−1
可避免交叠错误。
令f(x,y)f(x,y)f(x,y)和h(x,y)h(x,y)h(x,y)分别表示大小为AXBAXBAXB像素和CXDCXDCXD像素的图像阵列.它们循环卷积中的交叠错误可通过对这两个函数填充零来避免，
fp(x,y)={f(x,y),0≤x≤A−1和0≤y≤B−10,A≤x≤P或B≤y≤Qf_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),0≤x≤A−1和0≤y≤B−10,A≤x≤P或B≤y≤Q
hp(x,y)={h(x,y),0≤C−1和0≤y≤D−10,C≤x≤P或D≤y≤Qh_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),0≤C−1和0≤y≤D−10,C≤x≤P或D≤y≤Q
其中，P≥A+C−1P \ge A+C-1P≥A+C−1, Q≥B+D−1Q \ge B+D-1Q≥B+D−1
填充零后图像大小为PXQPXQPXQ,若两个阵列大小相同，都为MXNMXNMXN,则要求P≥2M−1P\ge2M-1P≥2M−1和Q≥2N−1Q\ge2N-1Q≥2N−1.DFT算法执行偶数大小的阵列时速度较快，因此通常选择满足上述公式的最小偶整数作为P和Q的值，若两个阵列相同，则P=2MP=2MP=2M和Q=2NQ=2NQ=2N
3.频率域滤波步骤
1）一幅大小为MXNMXNMXN的输入图像f(x,y)f(x,y)f(x,y),由前述公式得到填充零后的尺寸P=2M和Q=2N
2)使用零填充，镜像填充或者复制填充，形成大小为PXQPXQPXQ的填充后的图像fP(x,y)f_P(x,y)fP(x,y)
3)根据前述二维傅里叶变换的平移性质，将fP(x,y)f_P(x,y)fP(x,y)乘以(−1)x+y(-1)^{x+y}(−1)x+y，使傅里叶变换位于PXQPXQPXQ大小的频率矩形的中心
4）计算步骤3得到的图像的DFT，即F(μ,v)F(\mu,v)F(μ,v)
5）建一个实对称滤波器传递函数H(μ,v)H(\mu, v)H(μ,v),其大小为PXQPXQPXQ,中心在(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−1G(i,k)=H(i,k)F(i.k),i=0,1,2,...,M-1G(i,k)=H(i,k)F(i.k),i=0,1,2,...,M−1和k=0,1,2,...,N−1k=0,1,2,...,N-1k=0,1,2,...,N−1
7）计算G(μ,v)G(\mu, v)G(μ,v)的IDFT得到滤波后的图像(大小为PXQ)：gp(x,y)={real[ζ−1[G(μ,v)]]}(−1)x+yg_p(x,y)=\{real[\zeta^{-1}[G(\mu,v)]]\}(-1)^{x+y}gp(x,y)={real[ζ−1[G(μ,v)]]}(−1)x+y
8）最后，从gP(x,y)g_P(x,y)gP(x,y)的左上象限提取一个大小为MXNMXNMXN的区域，得到与输入图像大小相同的滤波结果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 )将MatMatMat数组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 )计算xxx和yyy的模，dst(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 sizeint m = getOptimalDFTSize( I.rows );int n = getOptimalDFTSize( I.cols ); // on the border add zero valuescopyMakeBorder(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] = magnitudeMat magI = planes[0];//! [magnitude]//! [log]magI += Scalar::all(1); // switch to logarithmic scalelog(magI, magI);//! [log]//! [crop_rearrange]// crop the spectrum, if it has an odd number of rows or columnsmagI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));// rearrange the quadrants of Fourier image so that the origin is at the image centerint cx = magI.cols/2;int cy = magI.rows/2;Mat q0(magI, Rect(0, 0, cx, cy)); // Top-Left - Create a ROI per quadrantMat q1(magI, Rect(cx, 0, cx, cy)); // Top-RightMat q2(magI, Rect(0, cy, cx, cy)); // Bottom-LeftMat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-RightMat 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 resultimshow("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
(十一)OpenCV实现图像频率域滤波相关推荐
- OpenCV —— 频率域滤波（傅里叶变换，低通和高通滤波，带通和带阻滤波，同态滤波）
频率域滤波 基本概念 傅里叶变换 二维离散的傅里叶变换 快速傅里叶变换 傅里叶幅度谱与相位谱 谱残差显著性检测 卷积与傅里叶变换的 频率域滤波 低通滤波和高通滤波 带通和带阻滤波 同态滤波 基本概念 ...
- 【OpenCV 例程200篇】86. 频率域滤波应用：指纹图像处理
[OpenCV 例程200篇]86. 频率域滤波应用:指纹图像处理 欢迎关注 『OpenCV 例程200篇』 系列,持续更新中 欢迎关注 『Python小白的OpenCV学习课』 系列,持续更新中 4 ...
- OpenCV实现频率域滤波——以高斯低通滤波去噪为例
最近由于作业原因,试着用OpenCV实现频率域滤波,但是OpenCV中并没有像MATLAB中fftshift这样的中心化操作,所以我写了一个频率域滤波的函数,以后用频率域滤波的时候在主函数中调用即可. ...
- 第4章 Python 数字图像处理(DIP) - 频率域滤波11 - 使用高通滤波器锐化图像
目录 使用高通滤波器锐化图像 由低通滤波器得到理想.高斯和巴特沃斯高通滤波器 指纹增强 频域中的拉普拉斯 钝化掩蔽.高提升滤波和高频强调滤波 同态滤波 使用高通滤波器锐化图像 由低通滤波器得到理想.高 ...
- 第4章 Python 数字图像处理(DIP) - 频率域滤波10 - 使用低通频率域滤波器平滑图像 - 理想、高斯、巴特沃斯低通滤波器
目录 使用低通频率域滤波器平滑图像 理想低通滤波器(ILPF) 高斯低通滤波器(GLPF) 巴特沃斯低通滤波器 低通滤波的例子 使用低通频率域滤波器平滑图像 理想低通滤波器(ILPF) 在以原点为中心 ...
- MATLAB中实现图像的空间域滤波和频率域滤波
1. 空间域滤波 空间域滤波是指在图像空间中借助模板对图像领域进行操作,处理图像每一个像素值.主要分为线性滤波和非线性滤波两类,根据功能可分为平滑滤波器和锐化滤波器.平滑可通过低通来实现,平滑的目的有 ...
- 传统基本图像处理方法：图像增强（灰度变换、直方图增强、空间域滤波、频率域滤波）、图像分割、图像配准等
图像处理设计主要有以下几种处理:图像增强(灰度变换.直方图增强.空间域滤波.频率域滤波).图像分割.图像配准等等. 图像增强: 图像增强作为基本的图像处理技术,目的在于通过对图像进行加工使其比原始图像 ...
- 图像 快速傅里叶变换 及 频率域滤波 java 实现
首先感谢中山大学12级软件学院计算机应用方向和我同班的乔勃大God,以及软件学院副院长.数图seisei朝老师的帮助!马屁还是要拍的o(*￣▽￣*)ブ //---------------------- ...
- 数字图像处理第四章频率域滤波（低通滤波器、高通滤波器、拉普拉斯滤波、同态滤波器）
本章节的主要内容具体包括:傅里叶变换的概念及处理的相关知识.频率域卷积概念.三种低通滤波器的原理及代码实现.三种高通滤波器的原理及代码实现.频率域拉普拉斯算法原理及实现.同态滤波器原理及代码实现. 4 ...
最新文章
- 三层架构和MVC一样吗？（区别）
- 基于 Jenkins 快速搭建持续集成环境
- c++对象长度之静态数据成员（3）
- 【IOS学习基础】OC类的相关
- MS SQL 2008 发布订阅配置错误总结
- Android开发之API29以上Environment.getExternalStoragePublicDirectory废弃的问题
- 本周ASP.NET英文技术文章推荐[03/25 - 03/31]
- svn 地址中文_iGit自助迁移SVN工程解密
- Amoeba+Mysql实现读写分离+java连接amoeba
- C++ MFC WebBrowser 探索（二）
- 腾讯云首发智能网关流控，公有云进入网络精细管控时代
- （娱乐项目）Python图片转换成矩阵数据，矩阵数据转换成图片
- word添加目录和更新
- 期权、期货及其他衍生产品 第三章读书笔记 利用期货的对冲策略
- 腾讯云 配置短信验证
- 周杰 清华大学计算机学院,清华大学自动化系主任周杰教授访问我院并做学术报告...
- 【友盟】 微博分享缺少C8998文件
- ion-slides 图片只能滑动一半、图片索引错误导致图片显示错误（缓存问题导致的）
- 免费网络硬盘、FTP、大容量邮箱、电子相册合集
- Java简单案例练习