用快速哈達瑪變換加速濾波反投影算法的濾波過程

喬志偉 韓 焱 魏學業

①(北京交通大學電子信息工程學院 北京 100044)②(中北大學電子測試技術國防科技重點實驗室 太原 030051)

1 引言

CT(Computed Tomography)，即計算機層析技術是當前最好的無損檢測技術之一。工業CT算法主要分為解析法和迭代法兩種，而解析法占據比較優勢的地位。在解析法圖像重建中，濾波反投影算法(FBP)是應用最為廣泛的一種算法。它的主要優點是重建質量好、速度較快、空間分辨率高，缺點是需要完備的投影數據[1]。當要重建的圖像很大的時候，FBP算法的速度仍然是一個問題，所以FBP算法的加速技術研究一直是該算法研究的熱點之一。從算法的角度實現加速可以從兩個方面實現，一個方面是濾波過程的加速；一個方面是反投影過程的加速。本文研究如何進一步加速濾波過程。濾波的過程是一個線性卷積的過程，所以濾波過程的加速問題其實是線性卷積的快速實現問題[2]。

目前，常用的長序列的線性卷積的加速主要有兩種方法。一種方法是用FFT(快速傅里葉變換)算法實現線性卷積；另一種方法是用分段卷積而后重疊相加或重疊保留的方法實現線性卷積[3,4]。FFT算法雖然速度較快，但是需要將信號補零，這就形成了一種冗余并影響了速度。如果信號的點數是一個大素數，則需要補零的點數將很多，其速度優勢將不再明顯。沃爾什變換(Walsh Transform,WT)是正交變換的一種，其變換矩陣的元素取+1或者-1。信號的沃爾什變換只涉及到實數的加法和減法[5,6]，而傅里葉變換涉及到了復數的乘法和加法，并包含三角函數的計算，顯然沃爾什變換是比傅里葉變換更快的一種正交變換。沃爾什變換依據沃爾什函數的序數的排序不同，形成了4種形式的變換：Paley序、反Paley序(Hadamard序)、Walsh序和反Walsh序[5]。基于Hadamard序的沃爾什變換又名哈達瑪變換(Hadamard Transform, HT)。因為HT對應的Hadamard矩陣可以因式分解，所以哈達瑪變換存在快速算法，本文將快速哈達瑪變換簡稱為FHT。既然利用FFT可以加速線性卷積，用更快的FHT則可以進一步加速卷積過程。將FHT用于FBP算法，可以加速FBP算法的濾波過程，進而提高圖像重建的速度。

在本文中，‘*’表示卷積，上角標‘-1’表示傅里葉變換和哈達瑪變換的反變換或者相應的變換矩陣的逆，‘diag’ 表示對角矩陣。

2 濾波反投影算法的濾波過程的計算特性

FBP算法是在傅里葉中心切片定理的基礎上推導出來的。該定理認為：物體f (x,y)在角度θ得到的平行投影的傅里葉變換，等于f (x,y)的2維傅里葉變換的一個切片[7]。投影示意圖如圖1所示。

顯然，如果我們能得到足夠多的各個角度的投影，對每個角度下的投影執行傅里葉變換，就可以得到物體傅里葉變換的各個切片，只要切片多到可以填滿整個傅里葉空間，就相當于得到了物體的傅里葉變換，然后對它執行傅里葉反變換就得到了物體。這種方法被稱為直接傅里葉重建[7]。用直接傅里葉重建法重建圖像時，要進行從極坐標系到直角坐標系的轉換，對于投影的傅里葉變換的離散采樣而言，此變換必然需要插值，而二維頻率域的插值會影響到整個圖像，插值的誤差傳播到整個圖像使圖像精度大大下降，所以這種方法的使用受到了一定的限制。

圖1 平行束投影形成示意圖

在中心切片定理的基礎上可以推導出濾波反投影算法的連續形式：

p(t,θ) 是某角度下投影，g(t,θ)是被單位沖激響應為h(t)的濾波器濾波后的濾波投影。

物體f(x,y)在位置(x,y)的線性衰減系數就是在各個角度下通過該點的所有的濾波投影的累加(嚴格講是對角度的定積分)。對投影進行濾波的濾波器是一個理想的斜變濾波器，由于它不滿足絕對可積條件，所以物理不可實現。考慮到對投影采樣的過程有低通濾波的作用，對該濾波器可以加窗，從而可以形成多種實用濾波器，使用最為廣泛的是RL濾波器和SL濾波器。

FBP算法的過程就是：首先在各個角度下測得相應的投影p(t,θ)；而后對每個投影用斜變濾波器濾波，形成濾波投影g(t,θ)；對于物體的每個點的線性衰減系數，用所有通過這個點的各個角度的濾波投影的值累加得到。顯然，整個過程分為濾波和反投影兩個步驟。本文要研究濾波過程的加速問題。

設探測器的寬度為d，共有M個像元，這樣在每個角度就可以得到一個長度為M的投影信號，可以看成是一個離散時間信號p(n)。不失一般性同時為了論述方便，設M是一個奇數。

采用S-L濾波器做濾波器，則濾波器的單位脈沖響應為

則濾波過程的離散實現公式為

設投影p(n)的取值范圍為(?(M?1)/2,(M ?1)/2)，則濾波以后需要的是同范圍的濾波投影g(n)，而h(n)是一個無限長的單位脈沖響應，根據卷積的定義可知，斜變濾波器的長度應取為2M-1，并且定義域關于0對稱[8]。

如果直接利用線性卷積的公式計算，比較復雜，需要做M2次實數乘法和M(M?1)次實數加法，當信號比較長的時候，花費的時間很多，必須根據濾波過程的特點設計新的加速算法。

3 快速哈達瑪變換(FHT)及其時間復雜度分析

沃爾什變換是一種簡單而又高速的正交變換，本文將利用它的一種特殊形式哈達瑪變換加速濾波的過程，本節將分析HT的定義、特點和FHT的時間復雜度。

3.1 定義

現有一個N點離散時間信號x(n)，且2qN=(q為正整數)，則其哈達瑪變換定義為

其反變換為

其中bi(x)是非負整數的二進制形式的第i位，如5的二進制形式是101，則b0(5)=1;b1(5)=0;b2(5)=1。

可以將如上的定義寫成矩陣的形式。設離散時間信號為一行向量x(n)=[x(0),x(1),x(2),…,x(N?1)]；設該信號的哈達瑪變換為XH(k)=[XH(0),XH(1),XH(2),…,XH(N?1)]，變換核矩陣(哈達瑪矩陣)定義為

3.2 哈達瑪變換的特點

(1)實數變換 通常用到的信號一般都是實信號，其哈達瑪變換仍然是實數的序列，這就避免了FFT所必需的復數運算。

(2)變換只涉及到加減法 從Hadamard矩陣的結構可以看出來，該變換只涉及加法和減法，這就避免了FFT所必需的乘法。

(3)反變換簡單 該變換的正反變換只相差一個系數1/N，這使得反變換可以借助正變換來計算。

(4)存在快速算法 因為Hadamard矩陣可以因式分解，所以該變換存在蝶形快速方法。

3.3 FHT變換的時間復雜度分析

分析HT的矩陣形式可以看到，計算HT的一個點需要做(N-1)次加減法，完成N點的HT共需要做N(N?1)≈N2次加法，當離散時間信號的長度增加時，其加法次數將以二次函數的方式增長，必須尋找快速算法。

利用哈達瑪矩陣可以因式分解的特點，可以得到跟庫利-圖基FFT算法相似的蝶形快速運算結構，顯然，FHT需要Nlog2N次加減法運算，這就使得其效率大大提高。如32點的信號做HT，用普通的方法需要做992次加減法；而用FHT算法，只需要做160次加減法。

4 用哈達瑪變換實現線性卷積的矩陣表達及哈達瑪域濾波器增益矩陣的求解

設有一個信號x(n)，長度為N1點，經過一個濾波器(系統)，其單位脈沖響應為h(n)，長度為N2點，兩個信號的起點都是0；則該信號的響應y(n)=x(n)*h(n),n=0,1,2,…,N1+N2?2，如果直接計算，運算量太大，考慮到用DFT可以間接實現線性卷積，而DFT存在快速算法FFT，可以用FFT快速實現線性卷積?？梢杂脠D2表示。

圖2 FFT實現線性卷積示意圖

在圖2中，L≥N1+N2?1,且L=2m,m 是自然數，這種計算線性卷積的方法就是所謂的FFT卷積法，現在的問題是在哈達瑪域是否存在與傅里葉域類似的規律。

現有理論表明：時域的卷積跟哈達瑪域(也就是沃爾什域)的相乘沒有對應關系，這就必須采用一種間接的方法來尋找這種關系。鑒于傅里葉變換和哈達瑪變換都可以寫成矩陣的形式[9]，而利用矩陣方程可以求解未知的矩陣，下面從矩陣表達的角度討論如何用哈達瑪變換實現線性卷積。

設x(n)和h(n)都已經補零至L點，圖2的矩陣實現如圖3所示。在圖3中，X是x(n)的行向量形式，Y是y(n)的行向量形式，F代表傅里葉變換核矩陣，F?1代表傅里葉反變換核矩陣，而

圖3 DFT實現線性卷積的矩陣形式

是濾波器h(n)的DFTHf(k)構成的一個對角矩陣。顯然可以得到

類比地，我們構造一個哈達瑪域的用HT實現線性卷積的方框圖[10]，見圖4。在圖4中，X是x(n)的行向量形式，Y是y(n)的行向量形式，H代表哈達瑪變換核矩陣(Hadamard矩陣)，?1H代表哈達瑪反變換核矩陣，而Gw是哈達瑪域的濾波器的增益矩陣，這是未知的。

圖4 HT實現線性卷積的矩陣形式

由圖4可知，

顯然，可以通過求解傅里葉域的濾波器增益矩陣，求得哈達瑪域的增益矩陣，從而實現在哈達瑪域的線性卷積，即哈達瑪域的濾波。因為哈達瑪變換有快速算法，所以在哈達瑪域完成濾波可以提高線性卷積的速度。其濾波框圖如圖5所示。

圖5 FHT實現線性卷積的示意圖

5 用FHT實現線性卷積的加速原理及適用性分析

為了可以定量分析用FHT實現線性卷積的時間復雜度，特在第4節的基礎上做出如下假設：

(1)L=N1+N2?1=2m。

(2)一次乘法相當于4次加法。

(3)一次減法等價于一次加法。

(4)將N1和N2都近似為L/2。

5.1 直接利用線性卷積公式計算的計算量

5.2 FFT實現線性卷積的計算量

(1)計算x(n)的FFT需要復數乘法(L/2)log2L次，復數加法Llog2L次，折合11L log2L次實數加法。

(2)計算h(n)的FFT，折合11L log2L次實數加法。

(3)計算Y(k)=X(k) H(k)，L次復數乘法，折合18L次實數加法。

(4)計算IFFT，得到y(n)，折合11L log2L次加法。

總的加法計算量為33L log2L+18L。其復雜度為O(Llog2L)。

5.3 FHT實現線性卷積的計算量

(1)計算L點的x(n)的FHT，需要Llog2L次加法。

(2)用公式(14)計算Gw，共20L3次加法。

(3)計算Xh(k)×Gw，共5L2次加法。

(4)計算IFHT，需要Llog2L次加法。

總的運算量為20L3+5L2+2L log2L 次加法。顯然其復雜度為O(L3)，用這種方法不但比FFT運算量大，甚至比直接計算還要耗費時間。可以看出，第(2)步和第(3)步非常浪費時間，如果能解決這兩個問題，就可以實現加速。

假設要求很多輸入信號通過特定的數字濾波器后的各自的響應。因為濾波器都是一樣的，所以可以將Gw事先計算出來，而后濾波時就無需每次都計算了，以此避免此復雜運算。

文獻已經證明，哈達瑪域濾波器增益矩陣Gw是一個分塊對角矩陣，并且是一個稀疏矩陣[9]，這樣計算Xh(k)?Gw的時候，就沒有必要按照向量乘矩陣的復雜度去計算了。通過統計分析，其計算量大約為3L log2L次乘法和3L log2L次加法。折合15L log2L次加法。

解決了如上兩個問題之后，FHT的總運算量為17Llog2L。顯然復雜度降低為O(Llog2L)，而且與FFT實現線性卷積相比，運算量降低了大約一半。

需要特別注意的是，FHT實現快速線性卷積的使用是有條件的：(1)有很多信號要被同一個濾波器濾波，從而可以事先計算增益濾波器。(2)利用Gw的分塊對角結構和稀疏特性，優化Xh(k)×Gw的計算。而我們要進行的FBP算法的濾波過程恰好符合條件(1)，所以可以將這種方法應用到圖像重建中。

6 數值實驗及結果討論

采用濾波反投影算法對1個由4個細圓柱插到1個粗圓柱的模型的某斷層進行重建。探測器像元共101個，寬度為0.1 cm；采集180o的投影，角度間隔為1o；對投影濾波采用SL濾波器；反投影時用線性插值方法插值；重建出的圖像是以旋轉中心為圖像中心的101×101的圖像，圖像的像素對應于實際物體的大小是0.1 cm×0.1 cm。表1為模型的參數說明。圖6為模型的某斷層零角度采樣示意圖。

表1 模型的參數說明

圖6 模型的某斷層零角度采樣示意圖

PC機的CPU的晶振頻率為1.7 GHz，開發工具使用MATLAB6.5，程序運行時，關閉了除殺毒軟件和必要的監控軟件以外的所有其他程序。

表2給出了直接卷積、使用FFT卷積和使用FHT卷積3種方法的濾波時間對比。表3給出了3種不同濾波方法對應的重建圖像的精度比較，其中dd和rr 分別是歸一化均方距離誤差判據和歸一化平均絕對距離判據[7]。圖7給出了模型原圖和3種方法重建出來的圖像。

表2 3種方法的速度比較

表3 3種方法重建出的圖像的誤差判據的比較

設重建圖像的大小N×N,tuv表示模型圖像的第u行和第v列的線性衰減系數，ruv表示重建圖像的第u行和第v列的線性衰減系數，t表示模型的平均衰減系數，則dd和rr 定義如式(15)。

從表2可以看出，用線性卷積的公式直接計算花費的時間非常長，用了將近13 s的時間。在重建過程中，投影信號共有180個，每個投影信號共101個點，濾波器長度為201個點，共進行了180×101×101=1836180次乘法，所以耗費的時間很長。用FFT算法，大大減少了運算量，花費的時間僅僅為0.063 s；而本文的方法，因為采用了快速哈達瑪變換，進一步加速，耗時僅僅為0.0296 s，比用FFT濾波快了大約一倍的時間。這個結果跟理論分析基本是一致的。

圖7 3種濾波方法重建出來的CT圖像對比

從表3中的誤差判據看出，采用3種方法重建出來的圖像的精度在誤差保留小數點后6位的情況下是一樣的。可以認為采用FHT來實現哈達瑪域的濾波不會影響重建圖像的精度。

從圖7可以看出，3種方法重建出來的圖像基本是一樣的，用肉眼分辨不出來，這也直觀地說明了3種方法只是影響濾波過程的速度，而不會對重建圖像的精度產生影響。

7 結束語

本文提出了一種用快速哈達瑪變換對濾波反投影算法的濾波過程加速的方法。通過引入FHT，使得濾波過程較通常采用的FFT濾波法速度提高了約1倍。需要注意的是，本方法要求事先計算出濾波器增益矩陣，并且要利用該增益矩陣的分塊對角結構減少向量與該矩陣相乘時的運算量。本方法適用于所有的解析法圖像重建的濾波過程。雖然本文是以FBP算法的加速問題為研究對象，但是這種方法同樣可以用到反投影濾波(BPF)算法、直接傅里葉域重建算法以及3維重建的FDK算法中。下一步要研究的工作是哈達瑪變換直接實現線性卷積的方法。

[1] 喬志偉, 魏學業, 韓焱. 解析法圖像重建中的插值技術研究[J].計算機工程與設計, 2009, 30(9): 2213-2216.Qiao Zhi-wei, Wei Xue-ye and Han Yan. Study on interpolation technology in image reconstruction based on analytic method[J]. Computer Engineering and Design, 2009,30(9): 2213-2216.

[2] 種稚萌, 朱世華, 呂剛明. 多徑信道下的異步分布式線性卷積空時編碼[J]. 電子與信息學報, 2009, 31(6): 1415-1419.Zhong Zhi-meng, Zhu Shi-hua, and Lü Gang-ming.Asynchronous distributed linear convolutional space-time code under multipath channels[J]. Journal of Electronics &Information Technology, 2009, 31(6): 1415-1419.

[3] 虞湘賓, 畢光國. 長序列信號快速相關及卷積的算法研究[J].電路與系統學報, 2001, 6(4): 78-82.Yu Xiang-bin and Bi Guang-guo. Algorithms of long sequence fast correlation and convolution[J]. Journal of Circuits and Systems, 2001, 6(4): 78-82.

[4] 呂新華, 武斌. 基于圓周卷積的長序列小波變換快速實現[J].信號處理, 2006, 26(2): 903-905.Lü Xin-hua and Wu Bin. Fast implementation of long sequence wavelet transform based on cyclic convolution[J].Signal Proccessing, 2006, 26(2): 903-905.

[5] 黃曉萍，桑恩方，喬鋼. H序沃爾什變換及其在水聲擴頻通信中的應用[J]. 聲學技術, 2007, 26(3): 477-482.Huang Xiao-ping, Sang En-fang, and Qiao Gang. Fast H-order Walsh transform and its applications in underwater acoustic spread-spectrum communication[J]. Technical Acoustics, 2007, 26(3): 477-482.

[6] 李何明, 張大興. 一種基于Hadamard變換的快速盲水印算法[J]. 杭州電子科技大學學報, 2009, 29(1): 67-70.Li He-ming and Zhang Da-xing. A blind fast image watermarking method based on Hadamard transform[J].Journal of Hangzhou Dianzi University, 2009, 29(1): 67-70.

[7] 張朝宗，郭志平，張朋. 工業CT技術和原理[M]. 北京: 科學出版社, 2009: 第8章.Zhang Chao-zong, Guo Zhi-ping, and Zhang Peng. Industrial CT Technology and Principle[M]. Beijing: Science Press, 2009,Chapter 8.

[8] 王召巴. 基于面陣CCD相機的高能X射線工業CT技術研究[D].[博士論文], 南京理工大學, 2002.Wang Zhao-ba. Study on high energy X-ray industrial CT technology based on area-CCD[D]. [Ph.D. dissertation],Nanjing University of Technology, 2002.

[9] Zarowski C and Yunik M. Spectral filtering using the fast walsh transform[J]. IEEE Transactions on Acoustics, Speech,and Signal Processing, 1985, 33(4): 1246-1252.

[10] Thomas G and Govindan V K. Computationally efficient filtered-backprojection algorithm for tomographic image reconstruction using walsh transform[J]. Journal of Visual Communication & Image Representation, 2006, 17(3):581-588.