用快速哈达玛变换加速滤波反投影算法的滤波过程

乔志伟 韩 焱 魏学业

①(北京交通大学电子信息工程学院 北京 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快速实现线性卷积。可以用图2表示。

图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.