一种应用于数字滤波的沃什-新梅森数快速变换方法

摘要

一种应用于数字滤波的沃什-新梅森数快速变换方法，需要沃什哈达玛变换与新梅森数变换结合而成，对于长度为N的沃什-新梅森数变换，输入数据的长度为N，对输入数据x进行重排得到xr，新梅森数变换矩阵通过同样的方法进行重排得到新的新梅森数变换矩阵，最后把同样是N阶的沃什-哈达玛矩阵与重排后的新梅森数矩阵结合在一起，形成阶数为N的T变换矩阵。WHNMNT的正变换首先经过T矩阵,再经过沃什-哈达玛矩阵便能得到输出数据，逆变换步骤与正变换相同。本发明提供一种有效降低计算复杂度、提升快速性的应用于数字滤波的沃什-新梅森数快速变换方法。

说明书

技术领域

本发明涉及数字滤波技术领域，尤其是一种新梅森数快速变换方法，该方法很好的降低新梅森数变换的计算复杂度。

背景技术

在数字信号处理中，数字滤波(卷积)是重要的技术领域，其运算量通常在信号处理算法中占较大比重，因此快速数字滤波对于很多信号处理算法的实际实现具有重要作用。前人关于快速数字滤波的研究一般是基于快速数学变换的，比如快速傅立叶变换(FFT,fast Fouriertransform)。综上，在数字滤波领域，快速变换的研究具有重要意义，其中快速数论变换在运算中不受量化截断的影响，受到了较多关注。

孙子定理是数论中最重要的定理，也是数论变换的基础，它又称为中国剩余定理（CRT，chinese remainder theorem）。中国南北朝的数学著作《孙子算经》中已出现有关孙子定理的命题。宋朝数学家秦九韶于《数书九章》对该命题做出了完整的解答。1972年，Rader根据数论知识提出用梅森数变换（MNT，Mersense number transform）用于计算卷积，此后数论变换开始被国外学者关注。而Reed于1975年对数论变换的变换域进行了扩展，使数论变换能在复数域内进行，这类变换称为复数数论变换（CNT，complex number transform）。

在详述本发明的快速数论变换之前，我们以模数对数论变换进行分类，大致可以分为以下几种类型：梅森数变换、费马数变换（FNT，Fermat number transform）、伪费马数变换、伪梅森数变换、广义数论变换与新梅森数变换（NMNT,new Mersense number transform）。上述以模数划分的每一类数论变换都存在自身对应的快速变换方法。

新梅森数变换是1995年S.Boussakta与A.G.J.Holt提出的一种最新的快速数论变换，该变换依托梅森数的成果，使得梅森数也能作为变换长度为2的幂次方的数论变换模数。新梅森数变换的提出很好的把针对费马数变换的快速算法与梅森变换结合到一起。它同样能应用于数字滤波、卷积和相关等信号处理运算中。现有的新梅森数变换快速算法有基2蝶形结构快速算法、基4蝶形结构快速算法等2的幂次蝶形算法结构。通过对新梅森数变换的变换矩阵进行类似FFT形式的蝶形化处理，可以使运算量下降。但该方法仍旧需要一定量的乘法器去实现，还存在进一步降低计算复杂度的可能。

发明内容：

为了降低数字滤波中现有新梅森数变换的计算复杂度，本发明提供一种有效降低计算复杂度、提升快速性的应用于数字滤波的沃什-新梅森数快速变换方法，变换长度N为2的幂次方。在输入序列的变换长度小于等于64时，本方法的计算复杂度低于已有的新梅森数变换方法。

为了解决上述技术问题采用的技术方案为：

一种应用于数字滤波的沃什-新梅森数快速变换方法，对于长度为N的沃什-新梅森数变换(WHNMNT，Walsh-Hadamard New Mersensenumber transform)，输入数据序列的长度为N，对输入数据x进行重排得到x_r，重排方式是把序列中每个数据的时间序号进行二进制化表示并逐位取反，而后用该取反值作为相应数据在新序列中的位置(时间序号)。例如，假设某个数据的时间序号是1(001₂)，那么重排后它在新序列中的位置应该是4(100₂)，这里下标2表示数的二进制表示；

如果新梅森变换矩阵用NMNT(N)表示，则本发明的快速变换方法需要对新梅森数变换矩阵进行与前述重排方法相同的列重排。不失一般性，称重排后的变换矩阵为NMNTr(N)。由于沃什哈达玛矩阵具有正交性，除了1/N这个因子外，沃什哈达玛正变换和逆变换相同，有如下公式成立：

{x}＝(1/N)WH(N)WH(N){x}；    (1)

式中WH(N)是沃什哈达玛变换矩阵。{x}表示变换的自变量是序列x。

对于输入信号{x}需要进行重排得到{x_r}，通过下标r表示数据重排，重排方式与数论变换矩阵的列重排一样，正变换后得到的序列X重排为X_r，则得到WHNMNT核心公式：

X_r＝NMNTr(N){x_r}＝(1/N)WH(N)WH(N)NMNTr[N]{x_r}    (2)

于是

X_r＝WH(N)T(N){x_r}    (3)

其中

T(N)＝(1/N)WH(N)NMNTr(N)    (4)

此处的T矩阵由新梅森数变换矩阵的位置重排矩阵与沃什哈达玛矩阵结合而成，其对数据的运算称为T变换，因此WHNMNT首先进行T变换，然后再进行沃什哈达玛变换，这和另一种快速数论变换沃什-费马数快速变换不一样。逆变换采用相同的方式，利用新梅森变换逆矩阵与沃什-哈达玛矩阵结合形成逆T矩阵，再经过沃什哈达玛变换得到{x_r}。

进一步，所述T矩阵具有块对角阵结构：

其中n＝log₂N，且T₀具有固定的值$>T_0=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),$>N等于2的幂次方。

T矩阵元素存在于模M_p整数域内，由于其特别的块对角阵结构使其三分之二的元素都为零，这样能更加方便地用于计算。T矩阵结构的证明如下：

参考文献1：Nibouche O.,Boussakta S.,Darnell M.PipelineArchitectures for Radix-2New Mersenne Number Transform[J].IEEETransactions on Circuits and Systems I:Regular Papers.2009,56(8):1668-1680.即Nibouche O.,Boussakta S.,Darnell M.，基2新梅森数变换中的管道架构，IEEE电路与系统汇刊I:长文，2009，第56卷第8期，1668-1680页。

Nibouche等人指出新梅森数变换矩阵通过重排列能得到：

$>\mathrm{NMNTr}\left(N\right)=\left(\begin{array}{cc}A[N/2]& B[N/2]\\ A[N/2]& -B[N/2]\end{array}\right)---\left(6\right)$>

另一方面，我们知道阶数为N=2ⁿ的沃什哈达玛矩阵可写成：

$>\mathrm{WH}\left(N\right)=\left(\begin{array}{cc}\mathrm{WH}[N/2]& \mathrm{WH}[N/2]\\ \mathrm{WH}[N/2]& -\mathrm{WH}[N/2]\end{array}\right)---\left(7\right)$>

于是

$>\left(\begin{array}{c}T\left(N\right)=(1/N)\mathrm{WH}\left[N\right]\mathrm{NMNTr}\left[N\right]\\ =(1/N)\left(\begin{array}{cc}\mathrm{WH}[N/2]& \mathrm{WH}[N/2]\\ \mathrm{WH}[N/2]& -\mathrm{WH}[N/2]\end{array}\right)\left(\begin{array}{cc}A[N/2]& B[N/2]\\ A[N/2]& -B[N/2]\end{array}\right)\\ =(1/N)\left(\begin{array}{cc}2\mathrm{WH}[N/2]A[N/2]& 0\\ 0& 2\mathrm{WH}[N/2]B[N/2]\end{array}\right)\end{array}\right)---\left(8\right)$>

而从文献2：Boussakta S.,Holt A.G.J.Relationship between theFermat number transform and the Walsh-Hadamard transform[J].IEEProceedings G Circuits,Devices and Systems.1989,136(4):191-204.即Boussakta S.,Holt A.G.J.费马数变换和沃什哈达玛变换之间的关系，IEE期刊G辑电路器件和系统，1989，第136卷第4期，191-204页。我们知道T(N)左上对角部分等于T(N/2)，再结合(7)易得

$>\left(\begin{array}{c}T[N/2]=(2/N)\mathrm{WH}[N/2]A[N/2]=(2/N)\mathrm{WH}[N/2]\mathrm{NMNTr}[N/2]\\ =(1/N)\left(\begin{array}{cc}4\mathrm{WH}[N/4]A[N/4]& 0\\ 0& 4\mathrm{WH}[N/4]B[N/4]\end{array}\right)\end{array}\right)---\left(9\right)$>

把公式(9)代入(8)得到：

$>\left(\begin{array}{c}T\left(N\right)=\left(\begin{array}{cc}T(N/2)& 0\\ 0& (2/N)\mathrm{WH}[N/2]B[N/2]\end{array}\right)\\ =\left(\begin{array}{ccc}(4/N)\mathrm{WH}[N/4]A[N/4]& & 0\\ .& & .\\ .& (4/N)\mathrm{WH}[N/4]B[N/4]& .\\ .& & .\\ 0& & (2/N)\mathrm{WH}[N/2]B[N/2]\end{array}\right)\end{array}\right)---\left(10\right)$>

以此类推进行不断迭代，得公式(4)形式的矩阵。

进一步，上述T矩阵存在如下单元子矩阵结构：

$>T_d=\left(\begin{array}{cc}a& b\\ b& a\end{array}\right)---\left(11\right)$>

该矩阵运算可通过如下公式进行变换：

w₁＝x₁a+x₂b＝x₁(a+b)+b(x₂-x₁)

                                       (12)

w₂＝x₂a+x₁b＝x₂(a+b)+b(x₁-x₂)

其中，x₁与x₂表示T变换前的两个输入数据，w₁与w₂表示经过T变换之后的输出数据，公式中的a+b通过事先计算完成。

本发明的有益效果主要表现在：进一步降低现有新梅森数变换方法的计算复杂度，本发明采用一种新的新梅森数变换结构以实现快速变换。在变换长度小于等于64的情况下，本方法比FFT形式的快速新梅森数变换具有更低的计算复杂度。

附图说明

图1为沃什-新梅森数变换快速算法思路，x_r表示对输入数据x的重排，T表示新梅森数变换的重排矩阵与沃什-哈达玛矩阵结合的矩阵，WHT表示沃什-哈达玛矩阵。

图2为长度为8，模为8191的沃什-新梅森数快速变换的实现结构。T₂(4)表示T矩阵中右下对角矩阵块(大小为4×4)。W(.)表示输入数据经过T变换后的值。X(.)表示经过沃什-新梅森数正变换后输出值。

图3为沃什-新梅森变换的正变换和逆变换仿真图，（a）图为输入数据，（b）图为正变换输出数据，（c）图为经过了正变换和逆变换后的输出数据。

图4为通过三个乘法器与两个加法器实现的T矩阵中单元子矩阵运算基本结构图。这里x₁与x₂表示两个T变换前的输入数据，w₁与w₂表示经过T变换之后的输出数据。a与b为单元子矩阵T_d中的元素。

图5为不同方法实现新梅森数变换所需乘法次数的增长趋势示意图。

具体实施方式

下面结合附图对本发明做进一步说明。

参照图1～图5，一种应用于数字滤波的沃什-新梅森数快速变换方法。本发明采用的快速变换方法过程如下：沃什-新梅森数变换快速方法（WHNMNT）需要沃什哈达玛变换与新梅森数变换结合而成，对于长度为N的沃什-新梅森数变换，输入数据的长度为N，对输入数据x进行重排得到x_r，重排方式是把序列中每个数据的时间序号进行二进制化表示并逐位取反，而后用该取反值作为相应数据在新序列中的位置(时间序号)。例如，假设某个数据的时间序号是1(001₂)，那么重排后它在新序列中的位置应该是4(100₂)，这里下标2表示数的二进制表示。进一步考虑N的值，N应等于2的幂次方。

如果新梅森变换矩阵用NMNT(N)表示，则本发明的快速变换方法要对新梅森数变换矩阵进行与前述重排方法相同的列重排。不失一般性，称重排后的变换矩阵为NMNTr(N)。由于沃什哈达玛矩阵具有正交性，除了1/N这个因子外，沃什哈达玛正变换和逆变换相同，有如下公式成立:

{x}＝(1/N)WH(N)WH(N){x}；    (1)

式中WH(N)是沃什哈达玛变换矩阵。{x}表示变换的自变量是序列x。

对于输入信号{x}需要进行重排得到{x_r}，通过下标r表示数据重排，重排方式与数论变换矩阵的列重排一样，，正变换后得到的序列X重排为X_r，则得到WHNMNT核心公式：

X_r＝NMNTr(N){x_r}＝(1/N)WH(N)WH(N)NMNTr[N]{x_r}    (2)

于是

X_r＝WH(N)T(N){x_r}    (3)

其中

T(N)＝(1/N)WH(N)NMNTr(N)    (4)

此处的T矩阵由新梅森数变换矩阵的位置重排矩阵与沃什哈达玛矩阵结合而成，其对数据的运算称为T变换。因此WHNMNT首先进行T变换，然后再进行沃什哈达玛变换。逆变换采用相同的方式，利用新梅森变换逆矩阵与沃什-哈达玛矩阵结合形成逆T矩阵，再经过沃什哈达玛变换得到{x_r}。

本实施例中，新梅森数变换定义

对输入数据x(m)，m为采样时刻的序号，进行长度为N的新梅森数正变换定义如下：

$>X\left(k\right)={⟨\underset{m=0}{\overset{N-1}{\Sigma }}x\left(m\right)\beta \left(\mathrm{mk}\right)⟩}_{M_p},k=\mathrm{0,1,2},...,N-1---\left(15\right)$>

表示模M_p，M_p是梅森数等于2^p-1。

β(m)＝β₁(m)+β₂(m)    (16)

$>\left(\begin{array}{c}{\beta }_1\left(m\right)={⟨\mathrm{Re}{({\alpha }_1+j{\alpha }_2)}^m⟩}_{M_p}\\ {\beta }_2\left(m\right)={⟨\mathrm{Im}{({\alpha }_1+j{\alpha }_2)}^m⟩}_{M_p}\end{array}\right)---\left(17\right)$>

这里q＝2^p-2。Re(.)和Im(.)表示只对括号中的式子截取其实数部分或者虚数部分。α₁和α₂的阶数为N_d=2^p+1，变换长度N=N_d/d，d为2的幂整数。新梅森数的逆变换定义如下：

$>x\left(m\right)={⟨(1/N)\underset{k=0}{\overset{N-1}{\Sigma }}X\left(k\right)\beta \left(\mathrm{mk}\right)⟩}_{M_p},m=\mathrm{0,1,2},...,N-1---\left(18\right)$>

因为M_p是一个奇数，N是2的幂次方，所以对N求逆总是存在。

本实施例中，沃什-新梅森数变换快速算法（WHNMNT）需要沃什哈达玛变换与新梅森数变换结合而成，对于长度为N的沃什-新梅森数变换，输入数据的长度为N，对输入数据x进行重排得到x_r，重排方式是把序列中每个数据的时间序号进行二进制化表示并逐位取反，而后用该取反值作为相应数据在新序列中的位置(时间序号)。例如，假设某个数据的时间序号是1(001₂)，那么重排后它在新序列中的位置应该是4(100₂)，这里下标2表示数的二进制表示。进一步考虑N的值，N应等于2的幂次方。新梅森变换矩阵用NMNT(N)表示，同样对新梅森数变换矩阵进行与上述方法同样的列重排。则重排后的变换矩阵为NMNTr(N)。把沃什-哈达玛矩阵与NMNTr(N)相结合形成T矩阵。沃什-新梅森数变换快速算法思路如图1所示。

T矩阵总是具有块对角阵结构

其中n＝log₂N，且T₀具有固定的值$>T_0=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),$>N等于2的幂次方。

T矩阵元素存在于模M_p整数域内，由于其特别的块对角阵结构使其三分之二的元素都为零，这样能更加方便地用于计算。T矩阵结构的证明如下：

文献1中，Nibouche等人指出新梅森数变换矩阵通过重排列能得到：

$>\mathrm{NMNTr}\left(N\right)=\left(\begin{array}{cc}A[N/2]& B[N/2]\\ A[N/2]& -B[N/2]\end{array}\right)---\left(6\right)$>

另一方面，我们知道阶数为N=2ⁿ的沃什哈达玛矩阵可写成：

$>\mathrm{WH}\left(N\right)=\left(\begin{array}{cc}\mathrm{WH}[N/2]& \mathrm{WH}[N/2]\\ \mathrm{WH}[N/2]& -\mathrm{WH}[N/2]\end{array}\right)---\left(7\right)$>

于是

$>\left(\begin{array}{c}T\left(N\right)=(1/N)\mathrm{WH}\left[N\right]\mathrm{NMNTr}\left[N\right]\\ =(1/N)\left(\begin{array}{cc}\mathrm{WH}[N/2]& \mathrm{WH}[N/2]\\ \mathrm{WH}[N/2]& -\mathrm{WH}[N/2]\end{array}\right)\left(\begin{array}{cc}A[N/2]& B[N/2]\\ A[N/2]& -B[N/2]\end{array}\right)\\ =(1/N)\left(\begin{array}{cc}2\mathrm{WH}[N/2]A[N/2]& 0\\ 0& 2\mathrm{WH}[N/2]B[N/2]\end{array}\right)\end{array}\right)---\left(8\right)$>

而从文献2我们知道T(N)左上对角部分等于T(N/2)，再结合(7)易得

$>\left(\begin{array}{c}T[N/2]=(2/N)\mathrm{WH}[N/2]A[N/2]=(2/N)\mathrm{WH}[N/2]\mathrm{NMNTr}[N/2]\\ =(1/N)\left(\begin{array}{cc}4\mathrm{WH}[N/4]A[N/4]& 0\\ 0& 4\mathrm{WH}[N/4]B[N/4]\end{array}\right)\end{array}\right)---\left(9\right)$>

把公式(9)代入(8)得到：

$>\left(\begin{array}{c}T\left(N\right)=\left(\begin{array}{cc}T(N/2)& 0\\ 0& (2/N)\mathrm{WH}[N/2]B[N/2]\end{array}\right)\\ =\left(\begin{array}{ccc}(4/N)\mathrm{WH}[N/4]A[N/4]& & 0\\ .& & .\\ .& (4/N)\mathrm{WH}[N/4]B[N/4]& .\\ .& & .\\ 0& & (2/N)\mathrm{WH}[N/2]B[N/2]\end{array}\right)\end{array}\right)---\left(10\right)$>

以此类推进行不断迭代，得公式(4)形式的矩阵。

对长度为8，模为8191的沃什-新梅森数快速变换的实现如图2所示。输入信号x=[-23,22,-13,2,3,24,-20,-1,17,-22,-3,6,24,4,-26,-28],对该信号进行WHNMNT的正变换和逆变换，结果如图3所示，显然同时经过正反变换以后，输出数据和输入数据相同，说明了该变换的可逆性和正交性。

为了进一步快速计算WHNMNT，考虑到T矩阵由单元子矩阵T_d构成，而T_d形式上具有如下结构：

$>T_d=\left(\begin{array}{cc}a& b\\ b& a\end{array}\right)---\left(11\right)$>

该矩阵运算可通过如下公式进行变换：

w₁＝x₁a+x₂b＝x₁(a+b)+b(x₂-x₁)

                                   (12)

w₂＝x₂a+x₁b＝x₂(a+b)+b(x₁-x₂)

其中，x₁与x₂表示T变换前的两个输入数据，w₁与w₂表示经过T变换之后的输出数据，公式中的a+b通过事先计算完成，这样在变换的过程中不需要为此做加法运算。公式(12)的乘法次数是3次，比普通直接乘做法少了1次。

由于数学变换的运算复杂度主要由乘法次数决定，本发明给出WHNMNT的乘法复杂度如下：

乘法次数：$>\underset{i=1}{\overset{{\log }_{2}N-2}{\Sigma }}\frac{3}{4}·{\left(\frac{N}{2^i}\right)}^2---\left(13\right)$>

进一步而言，对于WHNMNT的T矩阵，公式(5)中的T₂块矩阵由4个单元矩阵构成，其中只有两个独立单元矩阵，其它两个是它们的复制。而这两个独立的单元矩阵有一个特征，要么元素b为二次幂，要么a+b为二次幂。因此运用(12)计算的时候，我们可以用移位代替二次幂乘法，这样进一步可以降低乘法次数。利用这一特性可以把T₂块矩阵运算的乘法次数控制在6次，而不考虑这一特性的乘法次数为12次。因此当变换长度N=8时，总的乘法次数为6次，而N大于8时，公式（13）变化如下：

乘法次数：$>\underset{i=2}{\overset{{\log }_{2}N-2}{\Sigma }}\frac{3}{4}·{\left(\frac{N}{2^i}\right)}^2+6---\left(14\right)$>

通过下表对不同的新梅森数快速变换方法复杂度进行比较：

表1不同新梅森数变换方法的乘法运算量比较

当变换长度增加时，不同方法实现新梅森数变换所需乘法次数的增长趋势如图5所示。从图中可以发现，在变换长度小于等于64的情况下，本专利的算法比FFT形式的快速新梅森数变换需要的乘法次数要少，从而降低运算量要求，同时避免FFT算法截断误差的存在。由于N=8～64的数论变换应用比较广泛，而且WHNMNT形式上比较简单，因此，WHNMNT在数字滤波领域有实用价值。

图1为沃什-新梅森数快速变换方法思路，T表示将新梅森数变换重排矩阵与沃什-哈达玛矩阵结合后的矩阵，WHT表示沃什-哈达玛矩阵。WHNMNT首先对输入数据进行T变换，然后再进行哈达玛变换，便能得到变换输出。

图2为长度为8，模为8191的沃什-新梅森数快速变换的实现。T₂(4)表示T矩阵中右下对角矩阵块(大小为4×4)。W(.)表示输入数据经过T变换后的值。X(.)表示经过沃什-新梅森数正变换输出值。由于T₀与T₁中的元素为1和0，所以x(1),x(5),x(3),x(7)直接变成W(1),W(2),W(3),W(4)。

图3为沃什-新梅森变换的正变换和逆变换结果图，（a）图为输入数据，（b）图为正变换输出数据，（c）图为正变换和逆变换后输出数据。图中反映了沃什-新梅森变换的正确性。

图4为通过三个乘法器与两个加法器实现的T矩阵中单元子矩阵运算基本结构图。这里x₁与x₂表示两个T变换前的输入数据，w₁与w₂表示经过T变换之后的输出数据。a与b为单元子矩阵T_d中的元素。原本T矩阵基本运算结构通过四个乘法器与两个加法器实现单元子矩阵运算。而经过变换之后，图中的结构能使乘法次数减少一次。但加法次数则会增加一次，用于求解x₁-x₂，而x₂-x₁只需取反即可。公式中的a+b可以通过事先计算完成，这样在变换过程中不需要加法运算参与。

图5为不同算法实现新梅森数变换所需乘法次数的增长趋势。图中我们可以看出直接做新梅森数变换所需的乘法次数最多，当变换长度小于等于64时，沃什-新梅森数变换需要的乘法次数低于FFT形式的新梅森数变换。