卷積和傅立葉變換

本文版權屬于重慶大學計算機學院劉骥，禁止轉載

卷積和傅立葉變換
- 什麼是卷積
  - 卷積的數學定義
  - 卷積的計算
  - 卷積的作用
- 傅立葉變換
  - 傅立葉變換的用途
  - 時域頻域
- 卷積與傅立葉變換
- 後記

什麼是卷積

卷積的數學定義

f(x) 和 g(x) 是兩個連續函數，那麼卷積定義為：

h(x)=(f∗g)(x)=∫∞−∞f(u)g(x−u)du

若 f[n] 和 g[n] 是兩個數列，那麼卷積定義為：

h[n]=(f∗g)[n]=∑u=−∞∞f[u]g[n−u]

h(x) 稱為 g(x) 的卷積。 f 稱為卷積核，g是被卷函數（或者數列）。在圖像進行中，通常在推導公式時使用連續卷積，因為它有更好的數學性質。而在實際程式設計中，通常用離散卷積，因為計算機隻能處理離散值。卷積并非圖像處理的專有概念，它是從泛函分析、信号進行中引入的一個概念，是其他領域早就充分研究的一個概念。

卷積的計算

先用離散卷積公式舉個例子，假設 f 和g的定義如下：

f={1,1,1}g={2,3,2,6}

用離散卷積計算如下：

h[n]=(f∗g)[n]=f[1]g[n−1]+f[2]g[x−2]+f[3]g[n−3]h[1]=(f∗g)[1]=f[1]g[1−1]+f[2]g[1−2]+f[3]g[1−3]=unkownh[2]=(f∗g)[2]=f[1]g[2−1]+f[2]g[2−2]+f[3]g[2−3]=2h[3]=(f∗g)[3]=f[1]g[3−1]+f[2]g[3−2]+f[3]g[3−3]=5h[4]=(f∗g)[4]=f[1]g[4−1]+f[2]g[4−2]+f[3]g[4−3]=7h[5]=(f∗g)[5]=f[1]g[5−1]+f[2]g[5−2]+f[3]g[5−3]=11h[6]=(f∗g)[6]=f[1]g[6−1]+f[2]g[6−2]+f[3]g[6−3]=8h[7]=(f∗g)[7]=f[1]g[7−1]+f[2]g[7−2]+f[3]g[7−3]=6h[8]=(f∗g)[8]=f[1]g[8−1]+f[2]g[8−2]+f[3]g[8−3]=unkown

在計算中，如果下标超過數列的範圍(例如出現 g[−1] 或者 g[5] )，設定該值為0。按照該規則， h[2]−h[7] 能夠計算出值。取 h[2]−h[7] 構成新的數列：

h={2,5,7,11,8,6}

這就是 f∗g 的結果。

卷積不難計算，但為什麼書上或者其他資料上看起來那麼難呢？因為大多數教材和資料，都是用連續卷積來舉例。例如 f(x)=x ， g(x)=ax ，求 (f∗g)(x) 就需要用到積分運算，而大多數人在學完高等數學之後，微積分就還給老師啦！由于計算機隻能處理離散值，是以實際項目中面對的大多數是離散卷積。在函數庫的幫助下，這項工作甚至不用自己動手。如果使用python，那麼：

import numpy as np
f=np.array([,,])
g=np.array([,,,])
fg=np.convolve(f,g)
print(fg)

輸出結果如下：

卷積和傅立葉變換卷積和傅立葉變換

與手工計算沒有差别。是以

卷積的計算不會讓人困惑，真正讓人困惑的是卷積的作用。

卷積的作用

初學者都希望知道卷積的實體意義，知乎上有不錯的主題可以參考。但你真的需要了解的那麼清楚嗎？事實上，你隻需要了解卷積對于程式員的價值：

在程式中，卷積往往作為權重求和的替代物，因為程式員懶得再寫一個權重求和函數。

重要的的事情說三遍。姑且不論卷積有什麼重大的實體意義，在程式設計過程中，大多數程式員用卷積函數，隻是為了替代權重和。考慮這樣一個應用場景，有一個數列 g ，程式希望按照如下方式對數列中的資料進行處理：

g[n]′=2∗g[n−1]+3∗g[n]+4∗g[n−1]

具體來說這個程式如下：

# coding=utf-8
import numpy as np
g=np.array([,,,])
h=np.zeros(f.shape)
for i  in range(,f.shape[-]):
    h[i]=*g[i-]+*g[i]+*g[i+] #注意這裡有錯誤
print(h)

上述程式沒有考慮邊界條件，是以運作的時候會出現索引超出範圍的異常：

卷積和傅立葉變換卷積和傅立葉變換

當然可以修改這個程式，加入邊界條件的判斷（或者修改疊代範圍）：

# coding=utf-8
import numpy as np
g=np.array([,,,])
h=np.zeros(f.shape)
for i  in range(,f.shape[-]):
    if i-<:
        a=
    else:
        a=g[i-]
    if i+>f.shape[-]-:
        b=
    else:
        b=g[i+]
    h[i]=*a+*g[i]+*b 
print(h)

輸出結果為：

卷積和傅立葉變換卷積和傅立葉變換

這麼寫也太麻煩了吧？有沒有辦法簡化一下？有的，方法如下：

# coding=utf-8
import numpy as np
f=np.array([,,]) #注意f中數值的順序
g=np.array([,,,])
h=np.convolve(f,g) #用卷積替代
print(h)

結果如下：

卷積和傅立葉變換卷積和傅立葉變換

好像多了兩個值。很正常，因為卷積預設計算了兩個邊界。是以隻要去掉邊界就是想要的結果，代碼需要改為：

h=np.convolve(f,g,mode='same')

現在留意一下這段代碼：

f=np.array([,,])

請注意，之前定義的權重順序為 (2,3,4) ，這裡是 (4,3,2) 正好反向。如果對照卷積的定義，再略微思考，不難發現其中的秘密：

卷積公式：h[n]=(f∗g)[n]=∑u=−∞∞f[u]g[n−u]權重和公式：h[n]=∑i=1uf[i]g[n+(i−u−12−1)]

換言之，當把 f 反向之後，再做卷積，那麼就與權重和等價。并且如果要用卷積來做權重和，那麼數列f的長度最好是奇數（否則 u−12 除不盡），這也是為什麼圖像卷積使用的卷積核通常是3x3、5x5的原因。

再來做一次練習，如果要做以下操作：

g[n]′=(g[n+1]−g[n−1])/2

那麼對應的權重相當于 (−0.5,0,0.5) ，于是程式如下：

# coding=utf-8
import numpy as np
f=np.array([,,-]) #注意f中數值的順序
g=np.array([,,,])
h=np.convolve(f,g,mode='same') #用卷積替代
print(h)

結果就是：

卷積和傅立葉變換卷積和傅立葉變換

最後一項由于超出邊界，是以是 0.5∗0−0.5∗3=−1.5 。

至此，我們可以得到如下結論：

當程式員需要權重和時，他們會用卷積替代。計算機中所謂卷積，有兩種含義：真正意義上的卷積，或者權重和的另一種說法。

重要的事情說三遍。最後我們來談談圖像進行中的這樣卷積，那樣卷積。比如高斯卷積。現在你知道了，程式員的意思是

高斯函數為權重系數的權重和

。很多matlab代碼，真的是用權重和而不是卷積公式來實作高斯卷積的！雖然這麼寫沒錯，但每當看到這些代碼，很懷疑作者是否真的了解卷積。

好吧，大多數情況下，你不用了解卷積，你隻需要記住，

将卷積核反過來就是權重，卷積就變成了權重和。

傅立葉變換

傅立葉變換的用途

傅立葉變換是圖像進行中最複雜的一個概念。它之是以複雜是因為它是信号處理的知識。就像軟體架構，架構是一個建築學術語。另一個學科的常識，在其他學科眼中就變成了高深的學問。例如時間複雜度O(n)是計算機學科的常識，試試給學曆史的朋友解釋一下吧。那麼要搞懂傅立葉需要學信号處理嗎？Yes！不學是以搞不懂，沒有其他原因。那讀完這篇文章可以搞懂傅立葉變換嗎？做夢！這篇文章隻是一個科普，目的僅僅在于科普傅立葉變換的應用，而具體的原理肯定需要系統的學習。

首先我們來看一個問題，如圖：

卷積和傅立葉變換卷積和傅立葉變換

圖上藍色的線是一條直線 g1 ，桔紅色的是一條曲線 g2 。現在我們希望找到一個函數，将桔紅色的曲線變換為藍色的曲線：

g1(x)=f(g2(x))

能夠做到這一點嗎？當然能夠做到這一點。如果站在上帝的視角，可以知道：

g1(x)=2x+1

g2(x)=2x+1+sin(1111x)

是以函數 f 就是：

f(x)=2∗g−12(x)+1

問題在于，首先沒有上帝視角，其次 g−12 怎麼求呢？

下面再舉一個例子：

g1={1,2,3,4,5}

g2={1,2,7,4,5}

問如何将 g2 轉化為 g1 。從上帝視角可以推出，如果有一個數列：

f={0,0,−4,0,0}

那麼

g1[n]=g2[n]+f[n]

問題在于，現實中，如果隻能觀察到 g2 ，那麼如何才能轉化為 g1 呢？

再來看一個例子，圖像是這樣的：

卷積和傅立葉變換卷積和傅立葉變換

怎麼知道她其實是這樣的？

卷積和傅立葉變換卷積和傅立葉變換

這些問題都是信号處理的典型問題。而結論是：傅立葉變換可以解決。

時域、頻域

了解傅立葉變換所碰到的第一個問題就是何謂時域、頻域，而

傅立葉變換的實質就是将時域的函數變換為頻域的函數

。好吧，無論怎麼說，反正就是不懂。是以我準備給你一點感性認識。假設有這樣兩個數列：

g1={1,1,1,1,1,1,1,1}

g2={1,2,1,2,1,2,1,2}

如果把它們的圖像畫出來，就是這樣的：

卷積和傅立葉變換卷積和傅立葉變換

現在對它們做傅裡葉變換，python代碼如下：

# coding=utf-8
import numpy as np
g1=np.array([,,,,,,,]) #時域值
g2=np.array([,,,,,,,]) #時域值
fftg1=np.fft.rfft(g1)   #頻域值
fftg2=np.fft.rfft(g2)   #頻域值
print(fftg1)
print(fftg2)

得到結果如下：

卷積和傅立葉變換卷積和傅立葉變換

是複數啊？！當然是複數，請參考傅立葉變換的公式。這些值就是時域對應的頻域值。

它們是原始資料的另一種表現形式，好比把C語言程式翻譯成了彙編碼

。沒錯你看到的就是彙編碼。那麼它們有什麼用呢？如果修改彙編碼，是不是可以改變程式的運作呢？傅立葉變換的用途也在于此。

應用傅裡葉變換，重點不在傅立葉變換本身，而在于變換之後的處理。

變換之後的處理，則是一門專門的學問，不是一兩篇文章可以說得清楚的。

是以為什麼看不懂傅立葉變換呢？其原因在于：

傅立葉變換是加減乘除，懂加減乘除并沒有卵用。

重要的事情說三遍。下面我們來看看傅立葉變換之後的一種用途（實際的用途有成百上千種，這是一門專門的學問）。下面來繪制頻譜。什麼是頻譜呢？通俗的講是每個複數的模值，不通俗的講是信号中每個頻率的占比（我也不是很懂）。代碼如下：

# coding=utf-8
import numpy as np
from pylab import *
g1=np.array([,,,,,,,]) #時域值
g2=np.array([,,,,,,,]) #時域值
fftg1=np.fft.rfft(g1)   #頻域值
fftg2=np.fft.rfft(g2)   #頻域值
figure()
subplot()
plot(abs(fftg1)) #一直沒搞懂求模值為什麼是abs而不是norm
subplot()
plot(abs(fftg2))
show()

顯示結果如下：

卷積和傅立葉變換卷積和傅立葉變換

不一樣對吧？肯定不一樣啦！資料不一樣，肯定不一樣啦！簡直是廢話！我們看看下面的程式呢？

# coding=utf-8
import numpy as np
from pylab import *
g1=np.array([,,,,,,,]) #時域值
g2=np.array([,,,,,,,]) #時域值
fftg1=np.fft.rfft(g1)   #頻域值
fftg2=np.fft.rfft(g2)   #頻域值
figure()
subplot()
plot(abs(fftg1)) 
subplot()
plot(abs(fftg2))
show()

卷積和傅立葉變換卷積和傅立葉變換

可以發現，盡管資料不同，曲線的形式是相同的。再來看一個更一般的程式：

# coding=utf-8
import numpy as np
from pylab import *
x=np.linspace(,*pi,)
g1=sin(x) #時域值
g2=*sin(x) #時域值
fftg1=np.fft.rfft(g1)   #頻域值
fftg2=np.fft.rfft(g2)   #頻域值
figure()
subplot()
plot(abs(fftg1))
subplot()
plot(abs(fftg2))
show()

結果如下：

卷積和傅立葉變換卷積和傅立葉變換

再試試這個：

# coding=utf-8
import numpy as np
from pylab import *
x=np.linspace(,*pi,)
g1=sin(x) #時域值
g2=sin(*x) #時域值，頻率不同
fftg1=np.fft.rfft(g1)   #頻域值
fftg2=np.fft.rfft(g2)   #頻域值
figure()
subplot()
plot(abs(fftg1)) 
subplot()
plot(abs(fftg2))
show()

結果如下：

卷積和傅立葉變換卷積和傅立葉變換

是不是隐約有點感覺呢？不同數列的頻率分布是不同的，那麼它們對應的頻譜就不同。那是否可以通過調整頻率分布來對信号進行處理呢？例如：

# coding=utf-8
import numpy as np
from pylab import *
x=np.linspace(,*pi,)
g1=sin(x) #時域值
fftg1=np.fft.rfft(g1)   #頻域值
s=abs(fftg1)    #頻譜
fftg1=fftg1/s   
rs=roll(s,) #頻譜向右移動100
figure()
subplot()
plot(s)
subplot()
plot(rs)
fftg1=fftg1*rs  #恢複信号
figure()
subplot()
plot(x,g1) 
subplot()
plot(x,np.fft.irfft(fftg1)) #逆傅立葉變換
show()

結果如下：

卷積和傅立葉變換卷積和傅立葉變換

程式利用roll函數将頻譜右移，相當于把頻率變高。将信号還原，就可以看到變換後的結果：

卷積和傅立葉變換卷積和傅立葉變換

最後來看一個更複雜的操作：

# coding=utf-8
import numpy as np
from pylab import *
x=np.linspace(,*pi,)
g1=*x++sin(*x) #時域值
fftg1=np.fft.rfft(g1)   #頻域值
s=abs(fftg1)
mask=np.ones(s.shape,dtype=np.float32)
fftg1=fftg1/s
mask[:]=
rs=s*mask #用mask處理頻域
figure()
subplot()
plot(s)
subplot()
plot(rs)
fftg1=fftg1*rs
figure()
subplot()
plot(x,g1) 
subplot()
plot(x,np.fft.irfft(fftg1))
show()

頻域的變換結果為：

卷積和傅立葉變換卷積和傅立葉變換

轉換後的信号變為：

卷積和傅立葉變換卷積和傅立葉變換

從上面的例子可以看出：

傅立葉變換的研究重點不是如何進行變換，而是如何處理變換後的結果。

這個問題太複雜了，可以研究一輩子。數字圖像處理書籍中通常隻會介紹幾種比較簡單的、常用的處理。而什麼時候需要把時域函數轉換為頻域函數呢？當然是頻域處理起來更友善的時候！

卷積與傅立葉變換

現在讨論本文的最後一個問題。卷積的公式忘記了嗎？重新貼一下：

h(x)=(f∗g)(x)=∫∞−∞f(u)g(x−u)du

h[n]=(f∗g)[n]=∑u=−∞∞f[u]g[n−u]

有沒有想過為什麼公式這麼奇葩，這麼不好了解呢？有一個原因：

時域上的卷積對應着頻域上的乘積（傅立葉變換）

什麼意思呢？看下面這段代碼：

# coding=utf-8
import numpy as np
f=np.array([,,])
g=np.array([,,,])
fg=np.convolve(f,g) #直接卷積
n = len(f)+len(g)-
N = **(int(np.log2(n))+)
a=np.fft.rfft(f,N) #因為f和g不等長，是以擴充為相同的N位
b=np.fft.rfft(g,N)
c=a*b   #f和g傅立葉變換後的乘積
fft_fg=np.fft.irfft(c)[:n] #逆傅立葉變換，前n位有效
print(fg) #直接卷積
print(fft_fg) #傅立葉卷積
print(fg-fft_fg) #直接卷積和傅立葉卷積的差異

執行結果就是：

卷積和傅立葉變換卷積和傅立葉變換

如果寫成數學公式，那就是：

(f∗g)(x)=F−1(F(f(x))∗F(g(x)))

這條性質意味着什麼？

權重和 ⟶ 卷積 ⟶ 傅立葉變換後的乘積

權重和 ⟵ 卷積 ⟵ 傅立葉變換後的乘積

簡化之後就是：

權重和 ⟷ 傅立葉變換後的乘積

這有什麼用嗎？至少有以下兩個用途：

1. 加速算法運作。因為權重和，又要加又要乘，計算量是比較大的，而頻域的乘法是按位相乘，計算量要小得多。不過實際應用時要權衡傅立葉變換和逆變的時間。

2. 某些處理在頻域很容易（比如屏蔽特定頻率）。正常的做法是先做傅立葉變換，處理之後再逆變。如果能夠找到與之對應的權重和（卷積），那麼就可以直接用權重和（卷積）來計算。例如高斯濾波的實質是過濾高頻、低頻，保留中頻。實際應用中用卷積來實作高斯濾波，因為高斯函數是一種非常特殊的函數，它的傅裡葉變換是高斯函數，逆變仍然是高斯函數。

後記

本文所述内容隻能作為科普，若想真正掌握，還需要閱讀專業書籍，做大量練習，期待某一天突然開竅。如果你覺得本文通俗易懂，那麼說明我是搞懂了，但你未必。推薦一本書：

卷積和傅立葉變換卷積和傅立葉變換

本書對傅立葉變換的描述是非常清楚的，隻需要很少的數學知識也能讀懂。當然某些部分需要較長時間的思考，

有價值的知識，都是有門檻的

。

卷積和傅立葉變換卷積和傅立葉變換

卷積和傅立葉變換

什麼是卷積

卷積的數學定義

卷積的計算

卷積的作用

傅立葉變換

傅立葉變換的用途

時域、頻域

卷積與傅立葉變換

後記

繼續閱讀

（三十六：2021.01.13）Pre-MICCAI 2019學習（一）《從植入粒子後的CT圖像自動定位植入的粒子》講在前面摘要論文内容

Jetson TX2 opencv3.4的編譯安裝

二維圖像三維化錄影機标定視覺計算理論（Computational Theory of Vision）

Photoshop最常見的問答集

Visual C++多媒體設計及圖形、圖像處理

虛拟機下 linux 連接配接網絡并設定固定ip

鐵三角無線返聽系統ATW-3255。适合在演講、背景通訊聯絡和表演場所使用優化的内部數字信号處理提供了良好的立體聲分離度

3d相機機器視覺檢測技術

使用自相關函數的紋理圖像周期計算

2014ECCV視覺跟蹤錄用文章

opencv視覺跟蹤——消除背景模組化

JPEG壓縮技術JPEG壓縮的步驟

JPEG壓縮技術的原理

Python3 色情圖檔識别 Python 3 色情圖檔識别

BMP檔案結構及圖像每行位元組計算方法

如何用Opencv求圖像的灰階投影曲線