发明内容:
为了降低数字滤波中现有新梅森数变换的计算复杂度,本发明提供一种有效降低计算复杂度、提升快速性的应用于数字滤波的沃什-新梅森数快速变换方法,变换长度N为2的幂次方。在输入序列的变换长度小于等于64时,本方法的计算复杂度低于已有的新梅森数变换方法。
为了解决上述技术问题采用的技术方案为:
一种应用于数字滤波的沃什-新梅森数快速变换方法,对于长度为N的沃什-新梅森数变换(WHNMNT,Walsh-Hadamard New Mersensenumber transform),输入数据序列的长度为N,对输入数据x进行重排得到xr,重排方式是把序列中每个数据的时间序号进行二进制化表示并逐位取反,而后用该取反值作为相应数据在新序列中的位置(时间序号)。例如,假设某个数据的时间序号是1(0012),那么重排后它在新序列中的位置应该是4(1002),这里下标2表示数的二进制表示;
如果新梅森变换矩阵用NMNT(N)表示,则本发明的快速变换方法需要对新梅森数变换矩阵进行与前述重排方法相同的列重排。不失一般性,称重排后的变换矩阵为NMNTr(N)。由于沃什哈达玛矩阵具有正交性,除了1/N这个因子外,沃什哈达玛正变换和逆变换相同,有如下公式成立:
{x}=(1/N)WH(N)WH(N){x}; (1)
式中WH(N)是沃什哈达玛变换矩阵。{x}表示变换的自变量是序列x。
对于输入信号{x}需要进行重排得到{xr},通过下标r表示数据重排,重排方式与数论变换矩阵的列重排一样,正变换后得到的序列X重排为Xr,则得到WHNMNT核心公式:
Xr=NMNTr(N){xr}=(1/N)WH(N)WH(N)NMNTr[N]{xr} (2)
于是
Xr=WH(N)T(N){xr} (3)
其中
T(N)=(1/N)WH(N)NMNTr(N) (4)
此处的T矩阵由新梅森数变换矩阵的位置重排矩阵与沃什哈达玛矩阵结合而成,其对数据的运算称为T变换,因此WHNMNT首先进行T变换,然后再进行沃什哈达玛变换,这和另一种快速数论变换沃什-费马数快速变换不一样。逆变换采用相同的方式,利用新梅森变换逆矩阵与沃什-哈达玛矩阵结合形成逆T矩阵,再经过沃什哈达玛变换得到{xr}。
进一步,所述T矩阵具有块对角阵结构:
其中n=log2N,且T0具有固定的值N等于2的幂次方。
T矩阵元素存在于模Mp整数域内,由于其特别的块对角阵结构使其三分之二的元素都为零,这样能更加方便地用于计算。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等人指出新梅森数变换矩阵通过重排列能得到:
另一方面,我们知道阶数为N=2n的沃什哈达玛矩阵可写成:
于是
而从文献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)易得
把公式(9)代入(8)得到:
以此类推进行不断迭代,得公式(4)形式的矩阵。
进一步,上述T矩阵存在如下单元子矩阵结构:
该矩阵运算可通过如下公式进行变换:
其中,x1与x2表示T变换前的两个输入数据,w1与w2表示经过T变换之后的输出数据,公式中的a+b通过事先计算完成。
本发明的有益效果主要表现在:进一步降低现有新梅森数变换方法的计算复杂度,本发明采用一种新的新梅森数变换结构以实现快速变换。在变换长度小于等于64的情况下,本方法比FFT形式的快速新梅森数变换具有更低的计算复杂度。
具体实施方式
下面结合附图对本发明做进一步说明。
参照图1~图5,一种应用于数字滤波的沃什-新梅森数快速变换方法。本发明采用的快速变换方法过程如下:沃什-新梅森数变换快速方法(WHNMNT)需要沃什哈达玛变换与新梅森数变换结合而成,对于长度为N的沃什-新梅森数变换,输入数据的长度为N,对输入数据x进行重排得到xr,重排方式是把序列中每个数据的时间序号进行二进制化表示并逐位取反,而后用该取反值作为相应数据在新序列中的位置(时间序号)。例如,假设某个数据的时间序号是1(0012),那么重排后它在新序列中的位置应该是4(1002),这里下标2表示数的二进制表示。进一步考虑N的值,N应等于2的幂次方。
如果新梅森变换矩阵用NMNT(N)表示,则本发明的快速变换方法要对新梅森数变换矩阵进行与前述重排方法相同的列重排。不失一般性,称重排后的变换矩阵为NMNTr(N)。由于沃什哈达玛矩阵具有正交性,除了1/N这个因子外,沃什哈达玛正变换和逆变换相同,有如下公式成立:
{x}=(1/N)WH(N)WH(N){x}; (1)
式中WH(N)是沃什哈达玛变换矩阵。{x}表示变换的自变量是序列x。
对于输入信号{x}需要进行重排得到{xr},通过下标r表示数据重排,重排方式与数论变换矩阵的列重排一样,正变换后得到的序列X重排为Xr,则得到WHNMNT核心公式:
Xr=NMNTr(N){xr}=(1/N)WH(N)WH(N)NMNTr[N]{xr} (2)
于是
Xr=WH(N)T(N){xr} (3)
其中
T(N)=(1/N)WH(N)NMNTr(N) (4)
此处的T矩阵由新梅森数变换矩阵的位置重排矩阵与沃什哈达玛矩阵结合而成,其对数据的运算称为T变换。因此WHNMNT首先进行T变换,然后再进行沃什哈达玛变换。逆变换采用相同的方式,利用新梅森变换逆矩阵与沃什-哈达玛矩阵结合形成逆T矩阵,再经过沃什哈达玛变换得到{xr}。
本实施例中,新梅森数变换定义
对输入数据x(m),m为采样时刻的序号,进行长度为N的新梅森数正变换定义如下:
表示模Mp,Mp是梅森数等于2p-1。
β(m)=β1(m)+β2(m) (16)
这里Re(.)和Im(.)表示只对括号中的式子截取其实数部分或者虚数部分。α1和α2的阶数为Nd=2p+1,变换长度N=Nd/d,d为2的幂整数。新梅森数的逆变换定义如下:
因为Mp是一个奇数,N是2的幂次方,所以对N求逆总是存在。
本实施例中,沃什-新梅森数变换快速算法(WHNMNT)需要沃什哈达玛变换与新梅森数变换结合而成,对于长度为N的沃什-新梅森数变换,输入数据的长度为N,对输入数据x进行重排得到xr,重排方式是把序列中每个数据的时间序号进行二进制化表示并逐位取反,而后用该取反值作为相应数据在新序列中的位置(时间序号)。例如,假设某个数据的时间序号是1(0012),那么重排后它在新序列中的位置应该是4(1002),这里下标2表示数的二进制表示。进一步考虑N的值,N应等于2的幂次方。新梅森变换矩阵用NMNT(N)表示,同样对新梅森数变换矩阵进行与上述方法同样的列重排。则重排后的变换矩阵为NMNTr(N)。把沃什-哈达玛矩阵与NMNTr(N)相结合形成T矩阵。沃什-新梅森数变换快速算法思路如图1所示。
T矩阵总是具有块对角阵结构
其中n=log2N,且T0具有固定的值N等于2的幂次方。
T矩阵元素存在于模Mp整数域内,由于其特别的块对角阵结构使其三分之二的元素都为零,这样能更加方便地用于计算。T矩阵结构的证明如下:
文献1中,Nibouche等人指出新梅森数变换矩阵通过重排列能得到:
另一方面,我们知道阶数为N=2n的沃什哈达玛矩阵可写成:
于是
而从文献2我们知道T(N)左上对角部分等于T(N/2),再结合(7)易得
把公式(9)代入(8)得到:
以此类推进行不断迭代,得公式(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矩阵由单元子矩阵Td构成,而Td形式上具有如下结构:
该矩阵运算可通过如下公式进行变换:
其中,x1与x2表示T变换前的两个输入数据,w1与w2表示经过T变换之后的输出数据,公式中的a+b通过事先计算完成,这样在变换的过程中不需要为此做加法运算。公式(12)的乘法次数是3次,比普通直接乘做法少了1次。
由于数学变换的运算复杂度主要由乘法次数决定,本发明给出WHNMNT的乘法复杂度如下:
乘法次数:
进一步而言,对于WHNMNT的T矩阵,公式(5)中的T2块矩阵由4个单元矩阵构成,其中只有两个独立单元矩阵,其它两个是它们的复制。而这两个独立的单元矩阵有一个特征,要么元素b为二次幂,要么a+b为二次幂。因此运用(12)计算的时候,我们可以用移位代替二次幂乘法,这样进一步可以降低乘法次数。利用这一特性可以把T2块矩阵运算的乘法次数控制在6次,而不考虑这一特性的乘法次数为12次。因此当变换长度N=8时,总的乘法次数为6次,而N大于8时,公式(13)变化如下:
乘法次数:
通过下表对不同的新梅森数快速变换方法复杂度进行比较:
表1 不同新梅森数变换方法的乘法运算量比较
当变换长度增加时,不同方法实现新梅森数变换所需乘法次数的增长趋势如图5所示。从图中可以发现,在变换长度小于等于64的情况下,本专利的算法比FFT形式的快速新梅森数变换需要的乘法次数要少,从而降低运算量要求,同时避免FFT算法截断误差的存在。由于N=8~64的数论变换应用比较广泛,而且WHNMNT形式上比较简单,因此,WHNMNT在数字滤波领域有实用价值。
图1为沃什-新梅森数快速变换方法思路,T表示将新梅森数变换重排矩阵与沃什-哈达玛矩阵结合后的矩阵,WHT表示沃什-哈达玛矩阵。WHNMNT首先对输入数据进行T变换,然后再进行哈达玛变换,便能得到变换输出。
图2为长度为8,模为8191的沃什-新梅森数快速变换的实现。T2(4)表示T矩阵中右下对角矩阵块(大小为4×4)。W(.)表示输入数据经过T变换后的值。X(.)表示经过沃什-新梅森数正变换输出值。由于T0与T1中的元素为1和0,所以x(1),x(5),x(3),x(7)直接变成W(1),W(2),W(3),W(4)。
图3为沃什-新梅森变换的正变换和逆变换结果图,(a)图为输入数据,(b)图为正变换输出数据,(c)图为正变换和逆变换后输出数据。图中反映了沃什-新梅森变换的正确性。
图4为通过三个乘法器与两个加法器实现的T矩阵中单元子矩阵运算基本结构图。这里x1与x2表示两个T变换前的输入数据,w1与w2表示经过T变换之后的输出数据。a与b为单元子矩阵Td中的元素。原本T矩阵基本运算结构通过四个乘法器与两个加法器实现单元子矩阵运算。而经过变换之后,图中的结构能使乘法次数减少一次。但加法次数则会增加一次,用于求解x1-x2,而x2-x1只需取反即可。公式中的a+b可以通过事先计算完成,这样在变换过程中不需要加法运算参与。
图5为不同算法实现新梅森数变换所需乘法次数的增长趋势。图中我们可以看出直接做新梅森数变换所需的乘法次数最多,当变换长度小于等于64时,沃什-新梅森数变换需要的乘法次数低于FFT形式的新梅森数变换。