快速卷积中嵌套算法的设计与实现

i快速卷积中嵌套算法的设计与实现摘 要离散富里叶变换(DFT)和卷积计算在图象、数字信号处理中起着重要的作用，因此对快速算法的研究早就引起人们足够的重视。针对卷积算法的计算进行深入研究，发现在离散卷积计算过程中的计算量会随着输入信号序列的长度而急速增加，传统的卷积计算算法已不能满足要求，本文研究了如何将一维卷积变换成二维卷积或多维卷积，而多维卷积中由包含简单的一维卷积，从而进行嵌套计算。研究了利用短卷积嵌套计算长卷积的算法，最终实现 16 点循环卷积嵌套算法，大幅度减少了卷积的计算量。关键词：卷积,快速,嵌套iiDesign and implementation of the nested algorithm of fast convolutionAbstractDiscrete Fourier Transform (DFT) and convolution calculation plays an important role in the image, a digital signal processing, and therefore the study of fast algorithms have attracted enough attention. Through the In-depth study of the convolution algorithm, we find the calculation of the convolution calculation will increases with the length of the rapid of the input signal sequence. Conventional convolution calculation algorithm cant meet the requirements. This paper studies how the one-dimensional convolution transform into a multi-dimensional convolution or two-dimensional, and multi-dimensional convolution include simple one-dimensional convolution. Using short convolution calculate long convolution, and finally achieving 16 points nested circular convolution algorithm, significantly reducing the computation of convolution.Key Words: convolution;fast,;nested目 录摘 要 .iAbstract.ii第一章 绪论 .11.1 课题研究背景 .11.2 快速卷积算法的发展历史 .31.3 课题研究内容 .4第二章 快速卷积算法运算中的问题 .52.1 数字信号处理中的计算问题 .52.1.1 滤波和相关 .52.1.2 离散傅里叶变换 .82.2 算法序列 .11第三章利用短卷积嵌套计算长卷积算法原理简介 .133.1 二维卷积与多维卷积 .133.2 Agarwal-Cooley 卷积算法 .193.3 分裂算法 .28第四章 快速卷积嵌套算法实现 .344.1 16 点循环卷积算法实现 .344.2 算法性能分析 .35第五章 总结 .37参考文献 .381第一章 绪论1.1课题研究背景随着大规模集成电路、计算机和数字信号处理器（DSP）的快速发展，数字信号处理技术已得到广泛应用。然而在许多具体的应用场合需要处理的数据量和计算量都是巨大的，特别是在某些实时性要求较高的应用领域，对信号的处理速度是苛刻的。当然可以设计研制速度更快的处理器来满足要求，但研制新的处理器需要较大的投入和较长的研制时间，况且还有赖于技术的进步。另一种方案就是设计高效的算法使计算量降到可以接受的水平，实际上也正是出现了著名的离散傅里叶变换（DFT）的快速算法（FFT）后，数字信号处理技术才迅速发展起来并走向实用的。通常我们在表达一个算法的时候是以一种人们容易理解的方式表达出来的，然而这种表达方式计算起来并不高效，需要用较多的乘法和加法次数；另外一种表示方法就是从计算效率的角度来表示一个算法，但这种表达方法往往又会使得表达式变得让人难于理解。例如，我们要计算下面一个表达式的值 A(1-1)acdb如果按照这个表达式计算需要 4 次乘法和 3 次加法，不难看出上面的表达式也可写成(1-2)()Aabcd这样只需要 1 次乘法和 2 次加法就可以了。另一个典型的例子是计算两个复数的乘积(1-3)()()ejfajbcjd其中 直接计算需要 4 次乘法和 2 次加法，如eacbd果把它们写成2(1-4)()()ebacd()()fbacd则需要 3 次乘法和 5 次加法就可以了，通常 DSP 计算乘法要比计算加法来得复杂。如果 c 和 d 都是常数，那么 和 可以预先计算好，这样只需要3 次乘法和 3 次加法。这个复数乘法也可用矩阵表示(1-5)ecdafb按照第二个算法还可以写成(1-6)0110()0()ce adf bc 或写成简洁的形式(1-7)eaCMAfb矩阵 A 与向量 相乘只是一些加法运算，因此称 A 为前加矩阵；矩Tab阵 M 为一对角阵，它代表了算法中的所有乘法运算；C 则称后加矩阵。如果将 作为一个系统的输入数据， 作为系统的输出数据。那么，, ,ef这个算法可描述为：首先将输入数据进行一系列的加法运算，然后是一些乘法运算，最后又是一系列的加法运算。今后我们介绍的快速算法大多可以用这种矩阵结构表示，即它的中间是一个对角矩阵表示乘法运算，两边的矩阵元素是1，0，-1 等简单的数，表示加法运算。我们再来看两个矩阵相乘的问题设 ， ，其中 F 为某一数域，则 C 的元素CAB,mllnmnFC， (1-8)1lijikjcab1,;,ijn 若直接计算，每计算一个 C 中的元素就需要 次乘法和 次加法，因此l(1)l总共需要 次乘法和 次加法。现在我们来改进这个算法，使得乘法lmn()lmn次数几乎减少一半，而加法次数略有增加。利用下面的恒等式3(1-9)1212121()ababab则 (1-10)2 2,12,2,12,21,2,1,221,;,l lijikjikjikjikjl likikjcabababimjn 这里假定 是偶数，否则可在 A 中加入一列 0，在 B 中加入一行 0。上面的l表达式中第一项需要的乘法次数为 ，每一个元素都需要计算，因此共 次2l 2lmn乘法；第二项与 j 无关，每一行只要计算一次，因此共 次乘法；而第三项2l与 i 无关，每一列只要计算一次，因此共 次乘法；这样计算 C 总共需要2ln次乘法。1()2lmnl对每一个元素需要计算的加法次数为 次，对每一行或每一列要计算312l的加法次数为 次，因此加法次数总共为 次。12l 12lmnlmn对于一个大的矩阵乘法次数大约是直接计算的一半。但值得注意的是这种算法对计算误差会比较敏感，特别是采用定点数计算的时候，这是工程技术人员值得注意的问题。但本课程主要讨论快速算法的问题，对计算误差的问题不做进一步的讨论。1.2快速卷积算法的发展历史有限长序列卷积(1-11)1,01,.0yMlmlNxg广泛地应用于相关函数和功率谱的计算以及有限长脉冲响应数字滤波器的4设计和实现等众多的数字信号处理领域。然而,直接计算(l-11)式所需乘法及加法运算量，当卷积结果长度 N 较大时，将随着 量级增长很大。以致在快2(N)O速算法提出之前，(l-11)式的应用受到了很大的限制。1960 年和 1963 年 Good 和 Thomas 分别发表了他们的快速算法，但并没有引起人们的重视。1965 年，Coley 和 Tukcy 提出了一种称之为快速傅里叶变换(FFT)的计算机实现算法。应用该算法计算(l-11)式，其乘法及加法运算量均可降低到。1966 年 Stockham 将它用于卷积的计算，从而也推动的数字信2(Nlog)O号处理技术的快速发展。Cooley-Tukey 算法也被人们广泛第学习和研究，这种算法也就是一般的数字信号教课书中都有介绍的 FFT 算法；1978 年，wiongrad 推出 WFT 算法，使得(l-11)式的运算量进一步得到降低。然而 FFT 及 WFT 都存在两点不足，即：需要计算三角函数以及需要施行复数运算，即使输入值为实序列也是如此，而且当卷积结果长度 N 增大时，所需乘法及加法量仍然很大，难于实现实时相关或实时卷积计算。早在 1972 年，C.M.Rader 在前人关于 Galois 域上的离散傅里叶变换研究的基础上提出了应用默森变换方法，在取默森 M 为模的整数环上计算序列卷积的默森变换(MN)T 算法闭。接着 Agarwal 和 Burrsu 进行了更深入的研究，奠定了数论变换(NTT)的基础。NTT 的应用使得卷积的乘法运算量降低到 ，并且其算法结构便于(N)O硬件实现。然而遗憾的 NTT 受到字长及变换长度的严格限制，因此在应用上就有很大的局限性。1978 年，Nussbaumer 提出了快速多项式变换(FPT)算法,虽然解除了字长及变换长度受限的约束，但其乘法及加法运算仍然保持在的量级上，因而更适合于软件实现。2(Nlog)O1.3课题研究内容(1)研究快速卷积算法的基本原理和实现方法;(2)设计和实现快速卷积算法;(3)设计和推导快速卷积的嵌套算法并对算法性能进行评估;5第 2章 快速卷积算法运算中的问题2.1数字信号处理中的计算问题2.1.1滤波和相关数字信号处理中的主要任务是对信号进行滤波，有两种形式的数字滤波器：有限长冲激响应滤波器（FIR）和无限长冲激响应滤波器（IIR）无限长冲激响应滤波器也称递归滤波器。它们的结构可用 图 2.1.1 和 图 2.1.2 表示。D D D D. . . . .g0g1g2g3gL - 1 d2, d1, d0 s2, s1, s0gL - 2图 2.1.1FIR 滤波器对于 FIR 滤波器其输入输出关系为线性卷积。设 d 为 N 点的序列，滤波器有 L 个系数用 g 表示，即 ， ，那么其011,Nd 011,Lgg输出序列 s 就是， (2-1)10Liiksg,2iL数字信号处理中的另一个问题是计算两个序列的相关 r， (2-2)10Liikrdg0,12iN显然它也可以用卷积来计算，只要将一个序列颠倒就可以了。与线性卷积密切相关的是圆周卷积，设 d 和 g 都是 n 点的序列，圆周卷积定义为， (2-3)1()0nniiks0,1i表示对 n 取模，在含义清楚的情况下下标 n 也可不写。因 时，()n ki6； 时， 于是()ikiki()iknik， (2-4)10iiknikisdgg0,1in注意到当 时，第一项的 ， 是第二项的 ，所以kiiki0nikd， (2-5)1100nniikikinsdgs 0,1i此式表明线性卷积中下标大于等于 n 的部分折叠到下标 0 开始的地方，并将它们叠加起来。卷积还可以用多项式乘法来表示。设， (2-6)10()Niidx10()Liigx以及(2-7)20()LNiisxs则(2-8)()()sxdgx不难验证 的系数就是 d 和 g 的线性卷积。显然 还可以写成()sx ()sx，因此线性卷积是符合交换律的，还可以写成()sxgd(2-9)10Niiksgd圆周卷积也可以用多项式表示，设， (2-10)10()niidx10()niigx以及(2-11)10()niisxs则7(2-12)()mod1()mod1()n nsxxgxg显然圆周卷积也是符合交换律的，对 取模只要简单地将多项式中的n代以 1，而 代以 。 图 2.1.1 是用 FIR 滤波器来得到圆周卷积，nx,0nixix输入序列重复 2 次，在输出的 3n-1 点中，中间的 n 点即为圆周卷积。D D D. . . . .g0g1g2gn - 1dn - 1 d1, d0, dn - 1 , d1, d0gn - 2 10,nxsx中间 n 点为圆周卷积图 2.1.1用 FIR滤波器得到圆周卷积圆周卷积是一种快速计算线性卷积的重要手段。为计算一个长序列的线性卷积，通常将输入序列分段计算，然后将输出序列拼接。对输入数据的分段一般有两种方法：一是数据分段不重叠，这样卷积的输出需要叠加；另一种数据分段重叠，这需要将卷积输出中的混叠部分舍去。假定滤波器系数的长度为 n1 点，数据分段不重叠的方法称重叠相加法（overlap-add） ，将数据分成长度为 n2 点的互不重叠的数据段，记， ， (2-13)2,miind20,1 ,012m 然后，计算 与 ， 的 点圆周卷积，显然计算,iig 2n结果就是两个序列的线性卷积，记为 ,mis， (2-14)11,(),00nnmiikikkksdgd0,1in由于 ，所以2,iinmd(2-15)2 221 1,() ,(),0 0n ni ikmnkmiknmink ksgdgs 8因为 是（ ）点的序列，因此上式中相邻两段之间有 点的,mis12n 1n重叠，需要将它们叠加起来。如果 ，则可能出现两个以上的段重叠，21n所以一般要取 。21分段有重叠的方法称重叠保留法（overlap-save） 。记， ， (2-16)2,miind120,n 1,02m ，当 (2-17)11,(),00niikmikkksgdg112,in该式表明，在每次的计算结果中后面的 n2 个点等于线性卷积，所以我们只要保留后面的 n2 个点，而丢弃前面不正确的 n1-1 个点。最后将每次计算保留下来的 n2 个点的正确值拼接起来就可以了。D D D D. . . . .- 1h1h2hL - 1hL a2, a1, a0 q3, q2, q1, q0h3图 2.1.2 递归滤波器对于 IIR 滤波器其输入输出的关系为 ，这也可用多1LiLikiqahq项式来表示。设， (2-18)10()niiax0()Liihx其中 ，则01h(2-19)()()axQhxr其中 为商多项式， 为剩余多项式，当 n 个数据全部进入移位寄()Qxr存器后，其输出序列就是 的系数，而 系数的值保留在个寄存器中。()x()rx2.1.2离散傅里叶变换9离散傅里叶变换离（DFT）是数字信号处理中的另一个重要的计算问题，它与卷积是密切相关的，设 是复数域或实数域中的 n 维向,0,1ivn量，其一维离散傅里叶变换是复数域中的另一个 n 维向量 V，它定义为， ，其中 10nikkiVv, 2,1jne也可用矩阵的形式表示：(2-20)0 0211 142()2 212(1)(1)1 11nnnnn nvVv 由 V 也可求出原序列 v，即离散傅里叶反变换（IDFT）(2-21)10nikiV证明如下： 因 (2-22)111()000nnnikikl klillVv而(2-23)()1()00nlinkli il于是(2-24)1100nnikliiVv其中(2-25)10ilil事实上由离散傅里叶变换的矩阵形式可以看出10(2-26)10 02142(1)2 212(1)()1 11nnnnn nv Vv 显然(2-27)12 12(1)42(1) 24212(1)() (1)(1)(1)1n nnnnnnn 离散傅里叶变换与圆周卷积有密切的关系，即卷积定理。设序列 e 是序列f 和 g 的圆周卷积， 。 (2-28)1()0niillefg0,1in则， 。 (2-29)kkEFG,其中 分别是 e，f，g 的 DFT。,kkF证明如下：(2-30)1111()() )0000nnnnkl kikiliil ill lkief fGF从而 kkE还可以定义二维和多维离散傅里叶变换，由于多维离散变换的可分离性因此可以转化为一维离散傅里叶变换来计算。而一维离散傅里叶变换也可转化为多维离散傅里叶变换来计算。设 v 是复数域或实数域中的 二维数组，则n其二维离散傅里叶变换为复数域中的 二维数组，定义为n11(2-31)1, ,00,1nikk ii nVv 其中 ， ，也可以用矩阵表示2jne2jne(2-32)ikikv(2-33)21242()1(1)(1)1niknnn (2-34)21242()1(1)(1)niknnn 其反变换为(2-35)1, ,0 0,1nikii kk ivVn 或用矩阵形式表示为(2-35)1ikikvn2.2算法序列通常我们都是用一个表达式或者它的等价形式来表示一个算法的，然而它并没有具体地规定计算过程，它只表达了输入量和输出量之间的关系。如果我们将这种关系用具体的计算步骤表达出来就是算法序列。例如按照 Error! Reference source not found.式来计算两个复数的乘法，我们可以写出下面的计算步骤12(2-36)12344tacbdettcft其中 t 是在计算过程中引入的中间变量，通过这个算法序列首先可以清楚地看出这个计算过程需要 4 次乘法和 2 次加法，其次我们可以方便地将这个算法序列变成计算机程序。如果按照 Error! Reference source not found.式来计算两个复数的乘法，又可写出另外一个算法序列(2-37)123415263465tabdcttabetf这个算法序列中用到了 3 次乘法和 5 次加法。如果 d 和 c 是已知的，则和 可预先计算好，这样加法次数就降为 3 次。2t313第三章利用短卷积嵌套计算长卷积算法原理简介3.1二维卷积与多维卷积二维卷积是对一个二维数组进行的一种运算。可设想为一个二维的滤波过程，其中的一个二维数组是二维滤波器的系数，二维的数据通过它后形成一个二维的输出。图像信号是典型的二维信号，因此二维卷积对图像信号的处理是很有用的。现在给出它的定义，设是 mn 的数组，,01,;0,1ijdmj 是 pq 的数组，ijgp 它们的二维卷积定义为：，1()0mnijikjlklsd0,12;0,12imjn 是一个( m+p-1)(n+q-1)的数组。,2;,ij pjnq 根据这个定义容易推广到更高维的情况，但一般仅限于讨论二维的情况。二维卷积也可用多项式表示，将数组的每一行用多项式表示，这样 d 和 g 就成为具有 m 和 p 个以多项式为元素的向量：，10(),1njiijdxm 10(),1qjiijgxp写成矩阵的形式： 0010(1)011 1()0(1)(1)()() nnmmmddx xd 1400010(1)11 1()0(1)(1)()() qqppppgxgx计算二维卷积的问题就变成了计算 和 两个多项式向量的一维卷积dxg问题： 10()()miiksxx则是 m+p-1 个多项式的向量：()sx20(),1,2nqjiijsxsp计算两个多项式向量的卷积在结构上与一维卷积完全相同，因此前面的算法在这里也适用。所不同的只是原来的乘法变成多项式乘法，而加法变成了多项式加法。其中的每一个多项式乘法又是一个卷积。这个计算过程可用图3.1.1 表示。 ()DxG计算开始()(),0,11kkSxxMmp计算多项式乘积 ()sx重构结束由多项式向量()dx和g和lDlG计算 ,01()llSMnq计算 s重构和g和()kx()k和由 的系数构成的向量 d对每个 k注 ：(),Mmpnq表示计算 m p 和计算 n q 点线性卷积所需的乘法次数 。15图 3.1.1 mn 与 pq 的二维卷积二维卷积还可以写成二元多项式的形式，这时：，10(,)nijijdxydxy10(,)pqijijggxy20(,)(,),mpnqijijs s类似地可用多项式的形式写出两个 mn 维数组的两维圆周卷积：，10(),1njiijdx 10(),1njiijgxm()0()()odmniikksx则是 m 个多项式的向量：()sx10(),1njiijsx或写成二元多项式的形式： (,)(,),mododnsxygdxyxy其计算过程与计算线性卷积类似，如 图 3.1.2，这个计算过程可解释为计算长度为 m 的多项式圆周卷积，而多项式圆周卷积中的乘法是以 为模的多1nx项式乘法，因而是一个长度为 n 的圆周卷积。16()DxG计算 开始 ()()mod10,1nkkSxxM计算多项式乘积 ()sx重构结束由多项式向量()dx和g和lDlG计算 ,01()llSMn计算 s重构和g和()kx()k和由 的系数构成的向量 d对每个 k注 ：(),Mmn表示计算 m 和计算 n 点圆周卷积所需的乘法次数 。图 3.1.2 mn 二维圆周卷积由 图 3.1.2 可知计算一个 mn 的二维圆周卷积需要的乘法次数为：()()M加法次数为： ()()(AmnAnm一个二维数组可以看成是 mn 维的或者是 nm 维的其乘法次数是相同的而加法次数是不同的。例如计算 79 的二维圆周卷积，因为； ，它们的乘法次数为 304，而加法次(7)16,()70M(9)1,()74M数按 79 计算时为 1814，按 97 计算时为 1848。而直接计算需要次乘法和 次加法。29386302这个计算过程还可以从另一个角度解释。计算长度为 m 的多项式圆周卷积，首先要计算 和 ， 。它们是用两个矩阵乘以多项式()kGx()kD,1()M向量 和 得到的，这两个矩阵分别用 和 表示，它们都是gdmBA17的矩阵，下标 m 表示它们是 m 点圆周算法中的 B 矩阵和 A 矩阵。而()Mm和 也可用矩阵表示。gxd记 则：11TnXx，()mijDxAd()mijGBgX然后计算， 。()()o1nkkSGx0,()1kM这意味着要计算矩阵 和 的每一行元素的 n 点圆周卷积。因mijAdijBg此我们将它们的行变成列再乘以 和 于是得到 ，和nTuvijmDAd，将它们的对应元素相乘 ，TuvnijmGBg vvSG； ，这需要 次乘法。于是矩0,1()M 0,1()1vM ()nM阵 的每一列就是多项式向量 中每一个元素的系数。最后将矩阵nuvCSx的列变成行再左乘 就得到了二维圆周卷积的结果 。把这v mCTmuvnCS个计算过程画成流程图，如 图 3.1.3。18计算开始计算结束TuvnijmDAdijGBg计算uvuvS0,1()Mn m TmuvnCS计算二维D F T 开始计算结束计算I D F TuvuvSG0,1m nuivjuvijdijijg1iujvSn图 3.1.3 用 Winograd 算法计算 图 3.1.4 用二维 DFT 计算若用 D 和 G 表示 d 和 g 的二维 DFT，；由于二维 DFT 也符合卷积定理，因,01,;0,1uvuvSmvn 此 S 就是 s 的二维 DFT。也可用矩阵表示：，uivjuvijduivjuvijGg1ijvijsSmn其中：，11()muim 2jme，11()nvjn 2jne其计算流程图如 图 3.1.4。可以看到它们的计算过程是相似的，所不同的19是用二维 DFT 计算时中间需要 mn 次乘法，并且在计算 DFT 和 IDFT 时还需要大量的乘法，而且是复数乘法；而快速卷积计算时中间需要 次乘()Mnm法，而且在第一步和最后一步的计算中是不需要乘法的。现在我们要回到用短卷积计算长卷积的问题上来，基本思路是将一个长的序列变换成二维或多维的数组，从而转化成一个二维或多维的卷积问题，而二维和多维的卷积是可以通过一些短的一维卷积来计算的。3.2 Agarwal-Cooley 卷积算法设 d 和 g 是 n 点的序列，则它们的圆周卷积， 。1()0nniiksgd0,1in如果 ，且 互素，则根据孙子定理可将一维的下标 i 和 k 用二维01n1,的下标 和 来表示。它们的关系为(,)i)k和 0011modin0011modkn0011,;,iNiinn 10011od,;0,kk 其中 满足方程： ，这样 n 点圆周卷积就可改写为01,NNn011011001101()()nnnii ikNikNknksgd令：， ，0110,iNniis01101,iNniid01101,iNniig。,; 于是：2001010101, (),(),nnni ikikksgd这是一个二维圆周卷积，如果按照上一节的方法计算它其乘法次数和加法次数为，01()(Mnn1001)()()AnMAn或 0根据这些公式将会看到短卷积算法中的乘法次数对整个算法的计算量起着更重要的作用。例如我们要计算 504 点的圆周卷积，因为 504978，(7)16()70MA(9)1()MA而 8 点的算法有两种情况和(8)4()(8)2()7按第一种情况计算得到 (50)19645M(504)97887016457628A按第二种情况计算得到 (504)1962348(504)987710428263A可以看到按照第二种 8 点的算法计算时其乘法和加法次数都要比用第一种算法少，虽然第二种算法的加法次数比第一种多了 26 次。因此在选择短卷积算法时一般取乘法次数少的为好。Agarwal-Cooley 卷积算法的性能列于表 3.2.1 中。21表 3.2.1 Agarwal-Cooley 卷积算法的性能乘法次数 加法次数3 85 0卷积长度5 61 82 02 43 03 64 86 07 28 41 2 01 8 02 1 02 4 03 6 04 2 05 0 48 09 51 3 22 0 02 6 63 2 05 6 09 5 01 2 8 01 0 5 62 2 8 03 2 0 03 6 4 81 8 42 3 02 7 24 1 85 0 59 0 01 1 2 01 4 5 02 1 0 03 0 9 65 4 7 07 9 5 81 0 1 7 61 4 7 4 82 0 4 2 02 6 3 0 4每点乘法次数 每点加法次数2 . 1 12 . 5 02 . 3 32 . 6 72 . 6 42 . 7 53 . 3 33 . 6 93 . 8 14 . 6 75 . 2 86 . 1 04 . 4 06 . 3 37 . 6 27 . 2 41 0 . 2 21 1 . 5 01 1 . 3 31 3 . 9 31 4 . 0 31 8 . 7 51 8 . 6 72 0 . 1 42 5 . 0 02 5 . 8 03 0 . 3 93 7 . 9 04 2 . 4 04 0 . 9 74 8 . 6 85 2 . 1 9nMnAn/Mn/An8 4 01 0 0 81 2 6 02 5 2 07 6 8 01 0 0 3 21 2 1 6 02 9 1 8 45 2 7 8 87 1 2 6 59 5 7 4 42 4 1 6 8 09 . 1 49 . 9 59 . 6 51 1 . 5 86 2 . 8 47 0 . 7 07 5 . 9 99 5 . 9 0下面我们用 Agarwal-Cooley 算法把两个短卷积算法表示成一个算法。上一节中在计算二维卷积时首先要计算 ，和 ，TuvnijmDAdTuvnijmGBg我们把这些矩阵写成分块的形式，并重新排列。22记 的列向量为 ， 那么ijd01(),jjmjd,1n01011T TuvnijnmmnTmmDAdAA将矩阵 用它的元素具体写出来就是n0010,11(),0()1,()1,nnnnMnnMnaaA 所以 010()1,00110,1, ()1,TnnMuvmmnnaaDAA 记 ，101,1uuunnmnmnaa ,()u这是一个 的向量，它应该就是矩阵 转置的列向量，现在用()MuvD它构成一个 的列向量，仍然记为 ，则：1n0010,10 011 1(),0()1, ()1,() nnmnmmMnnMnn nmmmaAaAD 记00,1011(),(), (),nnnnmMMMmmnmaaaAAA 23这称为矩阵 和矩阵 的 Kronecker 积或称直积张量积等。它就是用nAm的每一个元素去乘以 从而构成一个 的矩阵。这样就可以nA()Mnm将 写成简洁的形式，我们仍然用 d 表示数据但它已不是原来的二维数组，而D是用原来的二维数组的列向量构成的具有 mn 个元素的一维数组。nmDA对滤波器系数做同样的处理，g 也是用原来的二维数组的列向量构成的一个一维数组，则 nmGBg然后计算 ，或写成矩阵的形式：,01,()1kkSDMM其中 M 是由 G 的元素构成的 的对角矩阵：()()nn，在形式上将它写成：01()1Mmn。这里的 是 的 个列向量构成的一维数组。记 的列向量也STuv() TuvS就是 的行向量为 ， 则： ，并且：uv01,()uuMmS,()1n 01()Mn0101,0101()10,(),()11,()011 nnnTmuvnmmMnnMnnMn ccCSC 24其中 ， 。0101,()1()jjjMnjnmmncCcC ,1jn用它构成一个 的一维数组，我们仍然记它为 s，但它是由原来的二维数组的列向量构成的一维数组。于是： 0010,()1001 11,01, 1,()()MnnmnmnnnnnmmcCcCsS 最后可以将它写成统一的形式： nmnnmsCMAd如果我们用原来二维数组的行向量构成一维数组，则可的到另外一种对称的形式： mnnmns d注意这里两个表达式中的 和 的排列次序是不同的。作为例子我们用 2d点和 3 点的短卷积算法来构造一个 6 点的算法。2 点的 Winograd 短卷积算法为, 000111sGd 00112Gg3 点的 Winograd 短卷积算法为 00 011 122 231120s dG 00112233Gg首先我们将 点的一维数组用孙子定理变换成二维的数组，然后6n25再用这个二维数组的列向量构成一个一维向量。变换关系为 ，0132mod6ii处理的过程示于 图 3.2.1。 0i1012345dd4235d01mo6ii034125dd图 3.2.1 数组整序32 11001112A01 02 33 44 15 26 571101DddD 可以看到输入数据序列的次序是乱的。为了将它调整为自然的次序，我们定义初等矩阵 ：ijP26101101ij ijP 第 行第 行左乘某个矩阵将交换这个矩阵的第 i 和 j 行，右乘则交换第 i 和 j 列，ijP并且 。于是：1ijijp01 02 33 44 15 26 571110011DddD 0123450102110112 dd02345dd由于这个初等矩阵右乘于 A，因此相当于交换 A 矩阵的列。对 G 做同样的处理得到：2701 02 33 44 15 26 57 11103306212GggG 034125gg由于 G 是可以预先计算的，因此可以不交换 B 的列，当然也可进行交换。这样就变成： 01 02 13 24 35 46 571300121gGg 于是： 011 022 133 24 4 35 5 46 6 57 7110121SG dSGd 00 13 24 31 42 55 6710101 2210Sss S 01234567SS28输出数据的次序也是乱的，需要将它调整为自然的次序。只要在 C 矩阵的左边乘一个初等矩阵，即交换 C 矩阵的行。最后的算法写成： 00 11 22 33 44 55 67111020101Gss G023452dd其中：01 02 13 24 35 46 571300121GggG 3.3 分裂算法在 Agarwal-Cooley 算法中，我们将一个长的一维数组转化成一个二维或多维的数组，从而就转化成计算二维或多维卷积的问题。而多维卷积可以通过一维卷积的嵌套来进行计算。例如前面讨论的二维卷积，我们首先将它看成是一维的多项式卷积，而多项式卷积中的每一个乘法又都是一维的卷积，这样就把一个二维的卷积转化成一维的卷积了。对于更多维的情况也是一样，例如要计算三维的卷积，我们可以把它转化成许多二维的卷积，而二维的卷积又可转化为许多一维的卷积。最后只要将多维的数组再还原为一维数组就可以了。所以这个问题归结为计算多维卷积的问题。我们还是以二维卷积来加以讨论。对于二维卷积可以用二元的多项式乘积表示： (,)(,),mod()()sxygdxypxqy29其中，10(,)mnijTijdxydxyXY10(,)pqijTjijgxygxyXY，1mX 1n，010,11,0,nmmndd 00,11,0,1,qppqgg如果 ， ，计算结果是线性卷积。deg()2pxdeg()2qx如果 ，且 ，则计算结果为圆周卷积。,nq()1,np假设 可分解为 u 个互素多项式， 可分解为 v 个互素多项式：()x()x； 。011()()upx 011()()qqx这样就可用计算 ， 的剩余多项式来,ijq,;,ijv 计算这个二元多项式的乘积，即计算： ,()(,),mod()()ijijij ijsxygxypxqy。01;0,1uv 然后用孙子定理重构 。用 和 分别表示计算以 和(,)sxy()iM()j ()ipx为模的多项式乘积所需要的乘法次数，则总的乘法次数为()jqx10()uvijijpq这样就将一个次数较高的二元多项式的乘积分裂为次数较小的二元多项式乘积。而每一个二元多项式的乘积都是一个二维卷积的问题，可以用前面介绍的嵌套的方法计算。分裂的方法是多样的。可以只对 或 分裂，当然也可以都分裂。下()pxq面以一个 20 点的圆周卷积为例说明它的计算过程。首先将一个 20 点的一维数组变换为 45 的二维数组：30100628415739119014mod20idddiii对 g 也做同样的处理，然后计算二元多项式乘积： 45(,)(,),1modsxygxyxy将 分解为 ，而对 不做分解，这51y432141x样就把它分裂成： 400(,)(,),od1sxygdxyxy 3211m1下一步是用孙子定理重构 ，因为(,)sxy432324)15yyy所以 432320 1451(,)1(,) (,)55modsxyysxyyysxx这里用了一次重构，如果将 也分解，它将分裂为更多的多项式乘积，4x那么要做两次重构。因为： 0162845713923 231 4(,)1ddydxyx 所以3101628457139230 1(,)1dddxyx 016284517392 31 1dddxx041641248459979139231 21 31230416412484599791(,)dddydxyxdyyydd 392321046142148459979319xdyyy 对 做同样的处理：(,)gxy0016284517139162841ggggxx2310416412484359979392104614214184235997939(,)gxygygygyxggygygy 因此计算 时是计算 4 点的一维400(,)(,),modsxydxx圆周卷积，它需要 5 次乘法和 15 次加法。而计算 则是计算 4443211(,)(,),11sgyyy点的多项式卷积，因此需要 5 次多项式乘法和 15 次多项式加法。对于模为的乘法需要 9 次乘法和 16 次加法。因此 59=45 次乘法，432yy516+4