1 信號處理中的卷積

無論是信號處理、圖像處理還是其他一些領(lǐng)域，我們經(jīng)常會在一些相互關(guān)聯(lián)的數(shù)據(jù)處理中使用卷積。卷積可以說是算法中一個非常重要的概念。這個概念最早起源于信號處理之中。
假設(shè)對于一個線性系統(tǒng)其在單位脈沖δ(t)的響應(yīng)下，輸出為h(t).那么在Aδ(t)的響應(yīng)下輸出為Ah(t).而所有的信號都可以用δ(t)乘以一個系數(shù)的和來表示。即

1. 信號響應(yīng)的計算
2. 和多項式的關(guān)系
3. 濾波器的應(yīng)用
4. 圖像處理中的應(yīng)用

2 matlab中卷積函數(shù)

matlab中提供了兩個常用的卷積函數(shù)conv和conv2.分別用來計算一維卷積和二維卷積。這兩個函數(shù)都是對離散的數(shù)據(jù)進(jìn)行卷積，接下來我們就分別討論這兩種卷積。

2.1一維卷積conv函數(shù)

通過對之前信號處理中卷積的介紹，我們將其進(jìn)行離散化處理可以得到如下離散的一維卷積表達(dá)式：

卷積和多項式

其實兩個多項式之間的乘法我們可以看作是它們系數(shù)矩陣之間的卷積。例如

clearclcu = [1 0 1];v = [2 7];w = conv(u,v)

執(zhí)行結(jié)果為w=[2 7 2 7],表示u,v得到的對應(yīng)多項式為：

卷積和濾波器

從另一方面我們也可以將卷積看作是一個濾波的過程（FIR），將u看成被濾波的信號，將v看成一個濾波器。此時濾波后得到的結(jié)果序列長度和u相同。在v里面有一個濾波中心點，在卷積乘加的過程中，中心點對應(yīng)“時刻”的濾波值始終和對應(yīng)卷積結(jié)果時刻的值相乘，通過平移濾波器序列改變中心點的位置，后面圖示會直觀的解釋這一概念。對于信號處理由于只能知道過去時刻的信號，所以中心點應(yīng)該選擇v的第一個元素上（即v[0]為中心點）。對于已知整個序列的信號（如圖像處理）則中心點的位置一般選擇在size(v)/2+1(如果不能整除則選在這個值的左邊或右邊)上。特殊的如果選擇的濾波器各個值都相等且和為1則這是一個均值濾波器，下面是一個中心點在v[0]，size(v)=3的一個均值濾波器濾波示意圖：

上圖中為了方便觀測卷積濾波器的物理意義，這里直接求w[2]和w[3].圖中的連線表示被連接的兩個數(shù)需要相乘?？梢钥闯龃藭rv就相當(dāng)于一個濾波因子對u的對應(yīng)值進(jìn)行濾波,我們將v看出一個單位脈沖響應(yīng)，因此也稱為FIR(有限脈沖響應(yīng))濾波器。
對于圖像處理或者已知整個信號序列值的場合，濾波器的中心點往往選在在v的中心（或左右），matlab中conv函數(shù)在調(diào)用的過程中如果加入了same參數(shù)則表示把該卷積看成一個濾波器，且濾波器的中心點在size(v)/2+1或者round(size(v)/2)(v長度為奇數(shù)時)。下面是matlab中將卷積看成濾波器且濾波器中心點在中間的計算代碼和計算過程的示意圖：

clearclcu = [1,2,3,4,3,1];v = [1/3 1/3 1/3];w = conv(u,v,'same')

執(zhí)行結(jié)果：

w =    1.0000    2.0000    3.0000    3.3333    2.6667    1.3333

我們發(fā)現(xiàn)將卷積看成是一個濾波器（濾波算子）后，在乘加的過程中我們總需要將v進(jìn)行交叉。這也是為什么在圖像處理中為什么各種算子和對應(yīng)的卷積核總是倒轉(zhuǎn)了180°的原因。針對這個問題有兩種解決方法，一個是改變在濾波過程中的卷積公式，將其中的減改成加即定義濾波卷積公式為：

在信號濾波中濾波器中心點一般選擇在最左邊，而在圖像處理中濾波器中心一般選擇在中心位置。其他位置的使用較少。事實上，可以通過對上節(jié)通用的卷積結(jié)果進(jìn)行移位和截斷能得到任意中心位置的卷積結(jié)果。
除了same參數(shù)以外，matlab還提供了valid參數(shù)用來計算有效的濾波點數(shù)。所謂的有效點數(shù)就是中心點在濾波器中心的條件下，能全部使用濾波器數(shù)據(jù)進(jìn)行計算的數(shù)據(jù)點。比如上面中的w[0]就不是一個有效的數(shù)據(jù)點。用圖像處理的角度來說就是，算子在計算邊界時沒有足夠的數(shù)據(jù)進(jìn)行計算然后舍棄這些點。計算后的w長度=max(length(u)-length(v)+1,0).下面是一個matlab示例，可以體會一下：

clearclcu = [1,2,3,4,3,1];v = [1/3 1/3 1/3];w = conv(u,v,'valid')

執(zhí)行結(jié)果如下：

w =    2.0000    3.0000    3.3333    2.6667

2.2二維卷積conv2函數(shù)

和前面一維卷積一樣，我們可以將卷積的概念推廣到二維卷積上。二維卷積最常用的地方就是圖像的濾波處理。這里我們將卷積所用的濾波器稱為卷積核，一些常用的卷積核和我們又稱之為某某算子。比如常用于邊沿檢測的Sobel算子。事實上，由于卷積的概念來自于數(shù)字信號處理，所以在圖像處理中算子和卷積核剛好關(guān)于中心點成中心對稱（前面我們已經(jīng)解釋了）。但是不引起混淆的情況下，我默認(rèn)它們是一樣的不再進(jìn)行區(qū)分了。
首先我們我們先看一下二維卷積的定義，可以使用doc conv2參考matlab中的說明文檔，但需要注意matlab中角標(biāo)是從1開始而不是0開始，公式會有些不一致。
兩個二維矩陣a,b的卷積結(jié)果c可以用如下的公式表示：

在matlab在調(diào)用二維卷積的時候同樣也提供了三個可選參數(shù):

full:     計算兩個二維矩陣的常規(guī)卷積（默認(rèn)參數(shù)），和一維卷積一樣same:     和一維卷積一樣valid:    和一維卷積一樣

上面三個參數(shù)中當(dāng)參數(shù)為full時表示的意義和信號處理中卷積的意義相同。當(dāng)參數(shù)為same時，對應(yīng)的卷積公式為：

same參數(shù)對應(yīng)公式:

可以看出full參數(shù)計算卷積求c(0,0)時，從圖像上來看首先是將b矩陣旋轉(zhuǎn)180°(注意b矩陣左上角為(2,2),右下角為(0,0)),然后讓a(0,0)和b(0,0)對齊，最后對應(yīng)乘加計算c(0,0)的值，這里有可能會涉及到a矩陣邊界以外的點。
對應(yīng)same參數(shù)來說大體計算流程并沒有改變多少，只是這里不是讓a(0,0)和b(0,0)對齊，而是讓a(0,0)和b矩陣的中心點b(mb/2,nb/2)對齊(行列為奇數(shù)時向下取整)，同樣這里有可能會涉及到a矩陣邊界以外的點。
由于邊界以外的點的情況我們往往并不知道，所以matlab中提供了一個valid參數(shù)，該參數(shù)下計算卷積時，如果計算某個卷積點c(i,j)用到了a矩陣邊界外的點，那么就會去掉這個卷積結(jié)果。其對應(yīng)的公式為：

除了上面和一維卷積相同的調(diào)用形式以外，matlab中還提供了另外一種調(diào)用形式，即C=conv2(h1,h2,A).它首先讓A的每一列和h1進(jìn)行卷積，然后再讓結(jié)果的每一行和h2進(jìn)行卷積，如果n1=length(h1),n2=length(h2),那么mc=max([ma+n1-1,ma,n1]),nc=max([na+n2-1],na,n2).一般情況下這種調(diào)用形式我們用的比較少。

2.3 邊界的擴(kuò)充

matlab的卷積計算過程中使用的是信號處理中求系統(tǒng)響應(yīng)的標(biāo)準(zhǔn)計算公式（圖像處理等實際上就是在標(biāo)準(zhǔn)卷積結(jié)果中進(jìn)行截取），這時候針對圖像處理中的“邊界”以外的點都是認(rèn)為為0的。然而在計算機(jī)圖形學(xué)中有時候我們并不認(rèn)為這是一個合理的假設(shè)。比如這幅圖是在一個大的圖像中截取的我們認(rèn)為邊界點的以外的點約等于邊界上的點更合適。或者這幅圖是一個周期循環(huán)的圖，我們在這里截取的是一個周期上的圖像，這時候認(rèn)為邊界以外的點等于另一邊界上的點更合適。
此外，當(dāng)卷積核一部分跑到邊界以外的時候，我們想舍去這些點。這時候我們就可以用valid參數(shù)進(jìn)行舍棄。

3 卷積的C++實現(xiàn)

為了使我們的代碼在計算系統(tǒng)響應(yīng)、多項式乘法、濾波器設(shè)計、圖像處理中有更廣泛的通用性，除了實現(xiàn)上面matlab中提供的卷積實現(xiàn)以外，我們還要增加matlabconv2函數(shù)中沒有的邊界處理條件選擇。
可選的邊界條件有：邊界以外的點為0（默認(rèn)），邊界以外的點和邊界值相等，邊界以外的點和另一邊界上的點組成周期信號。
這里實現(xiàn)C++卷積算法的時候使用的是C++矩陣計算庫Eigen.
需要注意的是這里卷積實現(xiàn)的時候?qū)onv和conv2合在一起了。事實上，在matlab中的大部分場合使用conv的地方都可以使用conv2代替，因此在C++實現(xiàn)里不再區(qū)分conv和conv2.但是需要注意的是，使用conv計算一個行向量和一個列向量的卷積得到的是一個向量，而使用conv2得到的卻是一個矩陣，在這一點上conv和conv2是有區(qū)別的。因此在使用時如果你確保你求的是一維卷積，請保證兩個向量要么都是行向量要么都是列向量。下面是matlab的一個示例：

clearclca = rand(1,3);b = rand(3,1);c = conv(a,b)d = conv2(a,b)e = conv2(a,b')執(zhí)行結(jié)果：c =    0.1409    0.4958    0.2941    0.1608    0.0380d =    0.1409    0.0357    0.0384    0.4601    0.1165    0.1255    0.1391    0.0352    0.0380e =    0.1409    0.4958    0.2941    0.1608    0.0380

3.1 C++中卷積的簡單實現(xiàn)

有了上面的分析之后，首先我們先編寫一個較為簡單的卷積函數(shù)來驗證一下功能。之后再將這個函數(shù)封裝成一個類。簡單的C++卷積代碼如下：

#include <iostream>  using namespace std;enum BoundaryCondition{    zero    ,bound    ,period};enum Method{    full    ,same    ,valid};//參數(shù)設(shè)置const int A_Row = 4;const int A_Col = 5;const int B_Row = 2;const int B_Col = 3;BoundaryCondition Bc = zero;Method method = full;float A[A_Row][A_Col] = {{1,2,3,4,5},                         {2,3,4,5,6},                         {3,4,5,6,7},                         {4,5,6,7,8}};float B[B_Row][B_Col] = {{1,2,1},                         {2,3,2}};float** C = 0;int cR;int cC;template <typename T>T GetA_Ele(int row,int col);template <typename T>void conv(const T a[][A_Col],const T b[][B_Col],T**& c);int main(){    conv(A,B,C);    for(int i = 0;i < cR;i++)    {        for(int j = 0;j < cC;j++)            cout<<C[i][j]<<"\t";        cout<<endl;    }    delete[] C;}//計算A和B的卷積  template<typename T>void conv(const T a[][A_Col],const T b[][B_Col],T**& c){    int offsetR = 0;    int offsetC = 0;        switch(method)    {        case full:            cR = A_Row + B_Row - 1;            cC = A_Col + B_Col - 1;            break;        case same:            cR = A_Row;            cC = A_Col;            offsetR = B_Row/2;            offsetC = B_Col/2;            break;        case valid:            cR = A_Row - B_Row + 1;            cC = A_Col - B_Col + 1;            if((cR < 1)|(cC < 1))                return;            offsetR = B_Row - 1;            offsetC = B_Col - 1;            break;        default:            return;    }    c = new T*[cR];             //給二維數(shù)組分配空間    for(int i = 0;i < cR;i++)        c[i] = new T[cC];    for(int i = 0;i < cR;i++)    {        for(int j = 0;j < cC;j++)        {            c[i][j] = 0;            for(int k1 = 0;k1 < B_Row;k1++)            {                for(int k2 = 0;k2 < B_Col;k2++)                    c[i][j] += b[k1][k2]*GetA_Ele<float>(i - k1 + offsetR,j - k2 + offsetC);            }        }    }}//根據(jù)邊界條件獲取A矩陣的元素template <typename T>T GetA_Ele(int row,int col){    switch(Bc)    {        case zero:      //索引超出界限認(rèn)為0            if((row < 0)|(row > A_Row)|(col < 0)|(col > A_Col))                return 0;        case bound:     //超出索引部分和邊界值相等            if(row < 0)                row = 0;            else if(row >= A_Row)                row = A_Row - 1;            if(col < 0)                col = 0;            else if(col >= A_Col)                col = A_Col - 1;            return A[row][col];        case period:            while((row < 0)|(row >= A_Row))            {                if(row < 0)                    row += A_Row;                else                    row -= A_Row;            }            while((col < 0)|(col >= A_Col))            {                if(col < 0)                    col += A_Col;                else                    col -= A_Col;            }            return A[row][col];        default:            return T(0);    }}

上面函數(shù)執(zhí)行結(jié)果為:

1       4       8       12      16      19      154       14      26      37      48      56      437       22      37      48      59      67      5110      30      48      59      70      78      5912      35      55      66      77      85      64

由于沒有使用矩陣運行庫且算法較為簡單，所以參數(shù)設(shè)置起來比較麻煩。需要手動給出矩陣的行和列。這里方法中的三個參數(shù)full,same,valid和MATLAB中相同。此外為了更方便的進(jìn)行圖像處理，有提供了3種可選的卷積邊界條件。由于在卷積計算的過程中，對A矩陣的取值有可能超出索引，因此必須使用一個Get方法進(jìn)行封裝。
為了使用matlab對我們的算法進(jìn)行驗證，這里取Bc = zero;，然后判斷三個不同的方法計算結(jié)果是否和matlab中相同。matlab對應(yīng)代碼如下：

clearclca = [1,2,3,4,5;...    2,3,4,5,6;...    3,4,5,6,7;...    4 5 6 7 8];b = [1 2 1;2 3 2];c = conv2(a,b,'full')

有了上面的簡單移植代碼做鋪墊，我們就可以編寫一個較為正式的conv模板來供我們以后使用了。對應(yīng)的C++類模板conv使用方法如下（源碼見文章的最后面）：

#include "Eigen\core"#include <iostream>#include "Conv.h"using namespace std;using namespace Eigen;using namespace ConvSpace;int main(){    MatrixXd a(4,5);    MatrixXd b(2,3);    a << 1, 2, 3, 4, 5,        2, 3, 4, 5, 6,        3, 4, 5, 6, 7,        4, 5, 6, 7, 8;    b << 1, 2, 1,        2, 3, 2;    try    {        Conv<MatrixXd> myConv(a, b);        //第一種調(diào)用形式        //Conv<MatrixXd> myConv(a, b,ConvType::Full,BoundaryConditon::Period);  //不使用默認(rèn)參數(shù)        cout << myConv.Eval() << endl;      //直接輸出結(jié)果        //const MatrixXd &r = myConv.Eval();    //第二種獲取結(jié)果的形式        //cout << r << endl;        //Conv<MatrixXd> myConv;            //第二種調(diào)用形式        //myConv.SetA(a);        //myConv.SetB(b);        //cout << myConv.Eval() << endl;    }    catch(ConvError ce)    {        cout<<ce.GetMessage()<<endl;    }    int c;    cin >> c;    //程序暫停}

執(zhí)行結(jié)果如下:

該算法需要使用到Eigen矩陣運行庫，因此需要包含Eigen相關(guān)文件和命名空間。此外這里將該算法在ConvSpace命名空間中實現(xiàn)，因此在使用函數(shù)之前必須要包含該命名空間和相關(guān)文件。
該函數(shù)大體上有兩種調(diào)用形式，一種是在構(gòu)造函數(shù)中給出被卷積矩陣然后進(jìn)行計算，另外一種是使用空構(gòu)造函數(shù)，然后通過Set方法設(shè)置相關(guān)矩陣，最后計算。第一種方式調(diào)用比較簡單如下：

Conv<MatrixXd> myConv(a, b);        //第一種調(diào)用形式cout << myConv.Eval() << endl;      //直接輸出結(jié)果

在使用上面代碼時應(yīng)該將所有的代碼放到一個try語句塊中，并且捕獲ConvError異常。當(dāng)出現(xiàn)參數(shù)配置錯誤或其他錯誤時，會拋出該異常通過查看該異常的說明可以判斷異常出現(xiàn)的原因。
上面使用了默認(rèn)的邊界條件zero和默認(rèn)的卷積方法full。當(dāng)然也可以自己選擇邊界條件和卷積方法（如上面注釋）。
前面在介紹原理的時候提到過，matlab中conv2還有另外一種調(diào)用形式，是計算兩個向量和一個矩陣的卷積。這里并沒有實現(xiàn)這個算法。實際上向量也是一種特殊的矩陣，有興趣的同學(xué)可以在代碼中補(bǔ)全這一部分內(nèi)容，最好將補(bǔ)全后的代碼再分享出來。
源代碼

4 關(guān)于圖像卷積的思考

由前面卷積的定義，我們可以發(fā)現(xiàn)把卷積的概念引入圖像處理中并不能直接的看出其物理意義，主要由以下三個方面的移植問題：

1. 它將其中的一個輸入序列“向前”進(jìn)行了平移，這就要求計算某個點k的卷積結(jié)果時，需要知道k+i“時刻”的輸入序列值。而這在時間序列作為輸入時往往是不可能做到的。
   
2. 有圖像卷積的圖示可以知道，卷積實際上是將卷積核先進(jìn)行中心對稱旋轉(zhuǎn)后，再和另外一個矩陣對應(yīng)乘加。由于并不是直接乘加因此沒有濾波算子的所表達(dá)的物理意義明確。
   
3. 在圖像處理中卷積往往有一個“邊界的概念”。而在數(shù)字信號處理中并不會有這種概念，因為在數(shù)字信號零狀態(tài)下，線性系統(tǒng)沒有輸入就沒有輸出。
   由上面幾點可以看出，卷積的概念推廣到圖像處理中，已經(jīng)有何很大的改變，也賦予了卷積更多的意義。由于卷積的概念在圖像處理中并實現(xiàn)完全的“本地化”給其賦予一個直觀的物理意義。因此我們可以對卷積的定義稍作修改，定義一個在圖像處理中意義更加明確的圖像卷積如下：
   
   
   
   其對應(yīng)的示意圖如下：

上面定義的圖像卷積和數(shù)組信號處理中的卷積略有不同，這樣將卷積的概念移植到圖像處理中并不會將卷積核進(jìn)行顛倒，注意b矩陣左上角為b(0,0).而且可以通過兩個可調(diào)的偏移量l和m調(diào)整卷積核在乘加過程中的中心位置。上面的圖像卷積公式有著更加直觀的物理意義。
當(dāng)然有利有弊，在數(shù)字信號處理中我們關(guān)于卷積的討論已經(jīng)有了很多非常成熟的結(jié)論。因此如果我們使用的是普通的卷積概念，這些討論和結(jié)論在圖像處理中是可以直接使用的。如卷積等于頻域中的乘積、卷積的交換律和分配率等。雖然上面卷積通過討論也能得到相似的性質(zhì)，但是畢竟增加了許多的工作量。
總之，有利有弊，如果我們需要直觀的在圖像處理中使用卷積的話，上面的公式在某種方面也有助于我們對基本卷積公式的理解。