歡迎關注 『OpenCV 例程200篇』 系列，持續更新中
歡迎關注 『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.頻率域圖像濾波（下）

---

youcans 的 OpenCV 學習課—8.頻率域圖像濾波（上）

本系列面向 Python 小白，從零開始實戰解說 OpenCV 項目實戰。
圖像濾波是在盡可能保留圖像細節特征的條件下對目標圖像的噪聲進行抑制，是常用的圖像預處理操作。
頻率域圖像處理先將圖像進行傅里葉變換，然后在變換域進行處理，最后進行傅里葉反變換轉換回空間域。
頻率域濾波是濾波器傳遞函數與輸入圖像的傅里葉變換的對應像素相乘。頻率域中的濾波概念更加直觀，濾波器設計也更容易。使用快速傅里葉變換算法，比空間卷積運算速度快很多。
本文提供上述各種算法的完整例程和運行結果，并通過一個應用案例示范綜合使用多種圖像增強方法。

---

1. 傅里葉變換

濾波通常是指對圖像中特定頻率的分量進行過濾或抑制。圖像濾波是在盡可能保留圖像細節特征的條件下對目標圖像的噪聲進行抑制，是常用的圖像預處理操作。

圖像濾波不僅可以在空間域進行還可以在頻率域進行。空間濾波是圖像與各種空間濾波器（模板、核）的卷積，而空間卷積的傅里葉變換是頻率域中相應變換的乘積，因此頻率域濾波是頻率域濾波器（傳遞函數）與圖像的傅里葉變換相乘。

頻率域圖像濾波，先將圖像進行傅里葉變換，然后在變換域進行處理，最后進行傅里葉反變換轉換回空間域。

空間域濾波器和頻率域濾波器形成一個傅里葉變換對：
$$\begin{gathered} f(x,y)?h(x,y)?F(u,v)H(u,v) \\ f(x,y)h(x,y)?F(u,v)?H(u,v) \end{gathered}$$
也就是說，空間域濾波器和頻率域濾波器實際上是相互對應的，有些空間域濾波器在頻率域通過傅里葉變換實現會更方便、更快速。例如，空間域的拉普拉斯濾波器就是頻率域的高通濾波器。

1.1 傅里葉級數與傅里葉變換

傅里葉級數（Fourier series）在數論、組合數學、信號處理、概率論、統計學、密碼學、聲學、光學等領域都有著廣泛的應用。

傅里葉級數公式指出，任何周期函數都可以表示為不同頻率的正弦函數和/或余弦函數的加權之和：
$$\begin{array}{rl}f(t)& =A_0+\sum _{n=1}^{\infty }{A}_{n}sin(n\omega t+{\psi }_n)\\ & =A_0+\sum _{n=1}^{\infty }[a_{n}cos(n\omega t)+b_{n}sin(n\omega t)]\end{array}$$
這個和就是傅里葉級數。

進一步地，任何非周期函數也可以表示為不同頻率的正弦函數和/或余弦函數乘以加權函數的積分：
$$\begin{array}{rl}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{array}$$
這個公式就是傅里葉變換（Fourier transform ）和逆變換。

*傅里葉變換存在的充分條件是：f(t) 的絕對值的積分是有限的，在信號處理、圖像處理領域這一條件都能滿足。

例程 8.1：連續周期信號的傅立葉級數

# 8.1：連續周期信號的傅立葉級數from scipy import integratenf = 30T = 10tao = 1.0y = 1k = 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 / Tfor 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 / 10half2, err2 = integrate.quad(lambda t: 2*y * np.sin(2.0/T * np.pi * n * t), -tao/2, tao/2)bn[n] = half2 / 10amp[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.piplt.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.001x = 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 seriesA0 = np.sum(f * np.ones_like(x)) * dxfFS = A0 / 2orders = 3A = np.zeros(orders)B = np.zeros(orders)for k in range(orders):A[k] = np.sum(f * np.cos(np.pi * (k+1) * x / np.pi)) * dx # Inner productB[k] = np.sum(f * np.sin(np.pi * (k+1) * x / np.pi)) * dxfFS = fFS + A[k] * np.cos((k + 1) * np.pi * x / np.pi) + B[k] * np.sin((k + 1) * np.pi * x / np.pi)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：一維連續函數的傅里葉變換（盒式函數）

以一維盒式函數為例，其傅里葉變換為：

$$\begin{array}{rl}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{array}$$

# 8.3：一維連續函數的傅里葉變換（盒式函數）# 盒式函數 (Box_function)x = np.arange(-5, 5, 0.1)w = 2 * np.pihalfW = w / 2y = 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)$，以自變量 t 的均勻間隔 $\Delta T$ 對函數取樣，取樣序列中的任意取樣值為：

$$f_k={\int }_{?\infty }^{+\infty }f(t)\delta (t?k\Delta T)dt=f(k\Delta T)$$

取樣后的函數的傅里葉變換為：

$$\begin{array}{rl}\tilde{F}(\mu )& =(F?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{array}$$

香農（Shannon）取樣定理指出，對于一個連續信號，用大于信號最高頻率 2倍的取樣率來取樣，則不會丟失信號的有效信息。或者說，以 $1/\Delta T$ 的取樣率對信號取樣點的的最大頻率是 ${\mu }_{max}=1/2\Delta T$ 。

例程 8.4：連續函數的取樣

# 8.4：連續函數的取樣# 定義函數，用于計算所有基矢的頻率def gen_freq(N, fs):k = np.arange(0, np.floor(N/2) + 1, 1)return (k * fs) / NT = 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 # 采樣間隔為 1freq = 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 # 采樣間隔為 10freq = 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 一維離散傅里葉變換

數字信號和數字圖像都是采樣得到的離散變量。

對原函數的變換取樣后的數據

$$\begin{array}{rl}\tilde{F}(\mu )& ={\int }_{?\infty }^{+\infty }\tilde{f}(t)e^{?j2\pi \mu t}dt\\ & ={\int }_{?\infty }^{+\infty }\sum _{n=?\infty }^{\infty }f(t)\delta (t?n\Delta T)e^{?j2\pi \mu t}dt\\ & =\sum _{n=?\infty }^{\infty }{\int }_{?\infty }^{+\infty }f(t)\delta (t?n\Delta T)e^{?j2\pi \mu t}dt\\ & =\sum _{n=?\infty }^{\infty }{f}_n{e}^{?j2\pi \mu n\Delta T}\end{array}$$

當取樣頻率為 $\mu =m/M\Delta T$ 時，可以得到離散傅里葉變換（DFT）和逆變換公式為：

$$\begin{array}{rlr}{F}_m& =\sum _{n=0}^{M?1}{f}_n{e}^{?j2\pi \mu mn/M},& m=0,...M?1\\ f_n& =\frac{1}{M}\sum _{m=0}^{M?1}{F}_m{e}^{j2\pi \mu mn/M},& n=0,...M?1\end{array}$$

由于任何周期或非周期函數都可以表示為不同頻率的正弦函數和余弦函數之和的形式，因此用傅里葉變換表示的函數特征完全可以通過傅里葉反變換來重建，而且不會丟失任何信息。這是頻率域圖像處理的數學基礎。

離散傅里葉變換 是將離散信號分解為一系列離散三角函數分量，每個三角函數分量都有對應的幅值 A、頻率 f 和相位 $\phi$。通過所有分量疊加可以得到原離散信號。

在圖像處理中，通常使用 $x,y$ 表示離散的空間域坐標變量，用 $u,v$ 表示離散的頻率域變量。于是將一維離散傅里葉變換表示為：

$$\begin{array}{rlr}F(u)& =\sum _{x=0}^{M?1}f(x)e^{?j2\pi ux/M},& u=0,...M?1\\ f(x)& =\frac{1}{M}\sum _{u=0}^{M?1}F(u)e^{j2\pi ux/M},& x=0,...M?1\end{array}$$

例程 8.6：一維離散傅里葉變換

# 8.6：一維離散傅里葉變換# 生成方波信號N = 200t = 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)$ 是二維連續變量 $t,z$ 的連續函數，則其二維連續傅里葉變換和逆變換為：

$$\begin{array}{rl}F(\mu ,\nu )& ={\int }_{?\infty }^{+\infty }{\int }_{?\infty }^{+\infty }f(t,z)e^{?j2\pi (\mu t+\nu z)}dtdz\\ f(t,z)& ={\int }_{?\infty }^{+\infty }{\int }_{?\infty }^{+\infty }F(\mu ,\nu )e^{j2\pi (\mu t+\nu z)}d\mu d\nu \end{array}$$
對于圖像處理，式中的 $t,z$ 表示連續空間變量， $\mu ,\nu$ 表示連續頻率變量。

例程 8.7：二維連續函數的傅里葉變換（二維盒式函數）

以盒式函數為例，其傅里葉變換為：

$$\begin{array}{rl}F(\mu ,\nu )& ={\int }_{?\infty }^{+\infty }{\int }_{?\infty }^{+\infty }f(t,z)e^{?j2\pi (\mu t+\nu z)}dtdz\\ & ={\int }_{?T/2}^{+T/2}{\int }_{?Z/2}^{+Z/2}A{e}^{?j2\pi (\mu t+\nu z)}dtdz\\ & =ATZ[\frac{sin(\pi \mu T)}{\pi \mu T}][\frac{sin(\pi \nu Z)}{\pi \nu Z}]\end{array}$$

# 8.7：二維連續函數的傅里葉變換（二維盒式函數）# 二維盒式函數 (2D-Box_function)t = np.arange(-5, 5, 0.1) # start, end, stepz = 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, 128m = 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$ 表示離散的空間域坐標變量，用 $u,v$ 表示離散的頻率域變量。二維離散傅里葉變換（DFT）和反變換（IDFT）為：

$$\begin{array}{rl}F(u,v)& =\sum _{x=0}^{M?1}\sum _{y=0}^{N?1}f(x,y)e^{?j2\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^{j2\pi (ux/M+vy/N)}\end{array}$$
二維離散傅里葉變換也可以用極坐標表示：
$$F(u,v)=R(u,v)+jI(u,v)=|F(u,v)|e^{j?(u,v)}$$
傅里葉頻譜（Fourier spectrum）為：
$$|F(u,v)|=[R^2(u,v)+I^2(u,v){]}^{1/2}$$
傅里葉相位譜（Fourier phase spectrum）為：
$$?(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)$$

空間取樣和頻率間隔是相互對應的，頻率域所對應的離散變量間的間隔為：$\Delta u=1/M\Delta T，\Delta v=1/N\Delta Z$。即：頻域中樣本之間的間隔，與空間樣本之間的間隔及樣本數量的乘積成反比。

空間域濾波器和頻率域濾波器也是相互對應的，二維卷積定理是在空間域和頻率域濾波之間建立等價關系的紐帶：
$$(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 二維數組（灰度圖像）或三維數組（彩色圖像）。

輸出結果是復數形式 $(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]]]) → dstcv.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 相同，運算公式為：

$$dst(I)=\sqrt{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(MN{)}^2$ 次運算，非常耗時；快速傅里葉變換只需要 $O(MNlog(MN))$ 次運算就可以完成。

OpenCV 中的傅里葉變換函數 cv.dft() 對于行數和列數都可以分解為 $2^p?3^q?5^r$ 的矩陣的計算性能最好。為了提高運算性能，可以對原矩陣的右側和下方補 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) # 將圖像轉換成 float32dft = 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 # 邊緣擴充，下側和右側補0dftImgEx = 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 的 OpenCV 學習課 原創作品
轉載必須標注原文鏈接：https://blog.csdn.net/youcans/article/details/122917951
Copyright 2021 youcans, XUPT
Crated：2022-02-14

歡迎關注 『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.頻率域圖像濾波（下）

總結

以上是生活随笔為你收集整理的youcans 的 OpenCV 学习课—8.频率域图像滤波（上）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。