The Discrete Hartley transform (DHT) is an important tool in digital signal processing. We propose a novel approach to perform DHT. We transform DHT into a form expressed in discrete moments via a modular mapping and truncating Taylor series expansion and present a completely new formula for computing DHT. We extend the use of our systolic array for fast computation of moments without any multiplications, to one that computes DHT with only a few multiplications and without any evaluations of triangular functions. The multiplication number used in our method is O(Nlog2N/log2log2N) superior to O(Nlog 2N) in the conventional FDT. The execution time of the systolic array is only O(Nlog2N/log2log2N) for 1-D DHT and O(N k) for k-D DHT (k⩾2). The systolic array consists of very simple processing elements and hence it implies an easy and potential hardware/VLSI implementation. The approach is also applicable to DHT inverses.published_or_final_versio

Chan, FHY

Lam, FK

Li, HF

Liu, JG

HKU Scholars Hub

Title A novel approach to fast discrete Hartley transformAuthor(s) Liu, JG; Chan, FHY; Lam, FK; Li, HFCitationThe 4th International Symposium on Parallel Architectures,Algorithms, and Networks, Perth/Fremantle, WA, Australia, 23-25June 1999, p. 178-183Issued Date 1999URL http://hdl.handle.net/10722/46134Rights Creative Commons: Attribution 3.0 Hong Kong LicenseA Novel Approach to Fast Discrete Hartley TransformJ. G. LIU*, F. H. Y. CHAN”“,  F. K. LAM**, and H. F. Li”“”*Institute for Pattern Recognition & Artificial IntelligenceHuazhong University of Science & TechnologyState Education Minsitry Laboratory forIrngge Processing & Intelligent Control, Wuhan, P. R. ChinaDepartment of Electrical and Electronic EngineeringTh$*$iversity  of Hong Kong, Hong KongDepartment of Computer ScienceUniversity of Concordia, Montreal, CanadaABSTRACTDiscrete Hartley transform (DHT) is animportant tool in digital signal processing. In thispresent paper, we propose a novel approach toperform DHT. We transform DHT into a formexpressed in discrete moments via a modularmapping and truncating Taylor series expansionand present a completely new formula forcomputing DHT. We extend the use of our systolicarray for fast computation of moments without anymultiplications to one that computes DHT with onlya few multiplications and without any evaluationsof triangular functions. The multiplication numberused  in  our  method is  O(NlogJVAog,logJV)superior to O(NlogJV) in the conventional FDT.The execution time of the systolic array is onlyO(NlogflAog,logJV for 1-D DHT and O(Nk) for k-D DHT(k2). The systolic array consists of verysimple processing elements and hence it implies aneasy and potential hardwareNLS1  implementation.The approach is also applicable to DHT inverses.1: IntroductionDiscrete Hartley transform (DHT) is widelyused in signal processing [l-3]. DHT maps a real-valued sequence to a real-valued spectrum whilepreserving some of the useful properties of the DFTand is a powerful tool as a substitute to the FIT forcomputing cyclic convolutions of real data [4-61.There have been many papers reporting thedevelopment with the fast DHT [7-91.On the other hand, computer vision and imageanalysis have propelled the advancement of fastcomputation of discrete moments (DM) [IO-l 31.These two research threads have been developedindependently. Up to now, no fast algorithmsproposed for DHT have some connections withmoments.In this paper, we first construct the bridgebetween DHT and discrete moments (DM) by amodular mapping and making use of Taylorexpansions and hence we can transform DHT intocomputation involving moments. We also discussthe problems of convergence and errors. Based onour approach to the fast calculation of moments[lo], a new systolic array to perform 1-D DHT ispresented, followed by a complexity analysis.Finally, we give our conclusions.2: Using moments to computingThe discrete Hartley transform (DHT) isdefined for a real valued length-N sequence x(n) (0G&N-l),  by the following equations:N-1X(k)=~x(n)cas(2nnklN)n=0N-l= ~x(n)(cos@mk  I N) + sin(2rmk  I N))O<k<N-  1 (2.1)Let us transform Equation (2.1) to look for anew approach to computing it. For every pair of kand i (i, k=O, 1, 2, . . . . . . . . . N-l), and defining S(k,  i)byS(k,  i)=( n I kn=i  (mod N) OlnlN-1 }then defining xk, i (i, k=O, 1, 2, . . . . . . . . . N-l) byO-7695-0231-8/99 $10.00 0 1999 IEEE178,.:6(j) Ski) f 4X1,, = N-I0 orherwise. + gxr.sin(5.,., + (P + l)n)(2m’/  NY /(2p + 2)i,k=O, I,2 ) . . . . . ...) N-l (2.2)For example, for N=4,  For example,Soe=(O,  1,2,3},  soi= s02=  s03=@;s, o={O],  S11={  11, s, 2=(21, S13={3);&O=(O,  21, %I=@),  f322={1,  3], s23=@;s3o=(O),s31=13},s32={2),s33={1}.xoo=x(0)+x(  1)+x(2)+x(3),  xol=O, x0*=0,  xo3=O;x10=  x(O), x11= x(l), x12= x(2), x13=x(3);X20=X(0)+X(2),  X*1=0,  X22=X(  1)+X(3),  X23=0;X30=X(O),  X31=X(3),  X32= X(2), X33=X(  1).Thus, by using the periodic properties of sinefunctions and cosine functions, Equation (2.1) canbe rewritten as follows:X(k)=~Xk,i(~os(~/N)+sin(2d/N))  OlklN-1  (2.3),iDFor OS&N-  1(2.4)This follows immediately by applying the theoremof extended law of the mean [ 141  to cos(2JQ/N)  andsin(27C/N),  where Ri is Taylor remainder term.Substituting them into Equation (2.3),  yieldsN-IX(k) = x‘,O+  ~xr.,~[(-1)‘(27Ti/  N)z’/(2r)!+li, r=0(-1)‘(2ti/  N)*‘+‘/(2r  + I)!]+ Rp= x,,,+~~(-1)‘(2~)“/N”(2r)!~x~,,iz’+~(-1)‘(2x)*‘+‘/N”t1(2~  + @x,,.i”+‘+  R,= XL.~ + za,mt., + R,r=oWhereOlklN- 1 (2.5)(-1)q(2n)‘q/N’p(2q)! r = 2q(-t)‘(2~)‘P+‘/N2q+‘(2q  + t)! r = 2q + 1o I q I p (2.6)N-lm,,, = 5 xk,, i’ (2.7)N-lR, = ;x,,[COSc_,  + (2p+  l)7C/2)(ti/N)‘p”/(2p  + l)ktsin(tS,“,  + (p + l)n)(k/  N)“”  l(2p + 2)!]O<L,.L,<~/N(2.8)If R, is ignored, we haveX(k)=x,,+ za, m,,, OlklN- 1 (2.9)The absolute value of the error introduced byoverlooking R, is bounded byI maxlx(n)l(N  - 1)[(2~)~~‘~/(2p  + 1)!+(2~)*p’z/(2p  + 2)!]n(2.10)From the inequality above, it is clearly shown thatR, converges to zero very rapidly and uniformly.For example, for max 1 x(n) 1 5256,  (O<n<N-1),  N<2048’ ‘=ll’R  1 <5.l3j(lO-7This error c$ satisfy accuracy requirements ofmost applications by computing only thirty-sixterms.Furthermore, we can prove that the least upperbound of p is not more than 0(log2N/Iog210g2N)  asN tends to infinite[  15, 161.Letf(N,p) = maxlx(r)lfi(N-1)(2n)P  /p!= A(N-1)(2~)~/p!(2. 11)Substituting [210gzN/log210gzN]  for p in Equation(2. 1 l)([x]  denotes an integer closest to x), we getf(N, p) = F(N)  = A(N - ,)(,)‘2’w2N”or2’w2N’ /[21og,  N/log, log, N]!Letthen210gzN/logzlogzN =t,ands oandF(N)IAN(271)  /[t]!.N=2 rilG%!b=%lN)f1 = (log,  N)“’tit=-i$log, N)“‘Ir  = l$log,  N)“‘log, log, N/2log,  N= lii(log, log, N)/ 2(log,  N)“’= lJm_(log,  M)/2M”  = 0then;E F(N) I li_iA(2a)“‘(log,  N)“’ /[I]!=b_mA(2~)“‘(log,  N)“‘/t’e~‘&=liimA(2x)(2ae(log,  N)“‘/t)‘/&=OThe final result is obtained by using the wellknown Stirling formula [ 171.  The computation ofX(k) using the approximation of Equation (2. 9)involves the generation of the r-th order moments ofthe transformed data sequence xki and thenperforming a dot product of these moments with aconstant vector (a,) and an addition with xk,O.  Thusif the moments can be computed very quickly andefficiently, so can be the above procedure.3: An algorithm for computing 1-D DHT179According to Section 2, fast DHT includes thefollowing three steps:1. Computing xki (i=O, 1, 2, . . . . . . . . . N-l; k=O, 1, 2,. . . . . . . , N-l) and a,(r=O,  1, 2, . . . . . . . . . 2p+l);2. Computing mk,r  (r=O, 1, 2, . . . . . . . . . 2p+l;  k=O, 1, 2,. . . . . . . . N-l);3. Computing X(k) (k=O, 1,2, . . . . . . N-l).The procedure of computing xk.i and a, is alsocalled the Preprocessing. The algorithm ofcomputing X(k) can be described below.Input initial p, N, x(O), x(l), . . . . . . . x(N-1)perform Preprocessingfor k=O to N- 1compute Momentzpi(xk,N_i,  xk,  N_2,  . . . . . . . . . xk,i)X(k)=0for r=O  to 2p+lX(k)=X(k)+a,mk.,end forX(k)=X(k)+ Xk.0end forMoment2pdXk.N.It  X k ,  N-2, . . . . . . . . ) Xk.1)  is asubroutine that produces the moments up to 2p+lorder [lo].  Now the complexity of the algorithmcan be given. In the Preprocessing, computing a,and index set S(k,  i) can be completed uponcondition that N and p are given. In real-timesystems, these computations can be preprocessedand are not included in execution procedure. Thecomputation of Xki is N’- yNlgcd(k,N) additionsY;o(greatest common divisor). Computing MomentzpiN times needs (2~+2)(2~+3)(N_2)N/2  additions.The computation of X(k) (k=O, 1, 2, . . . . . . . N-l)i n v o l v e s  (2p+2)N  a d d i t i o n s  a n d  (2p+l)Nmultiplications. S o  t h e r e  a r e  a l t o g e t h e rN’- f N Igcd(k,  N) +(2p+2)(2~+3)(N_2)N/2+L-0(2p+2)N additions and (2p+l)N  multiplications inthe algorithms. If N is a prime number,( N’ - $ N I gcd(k, N) ) disappears from the aboveformula.4: Pre-processing array forcomputing Xb,iIn the general case, we propose to use aspecial linear array to implement xk,i,  (i, k=O, 1, . . .N-l ) ,  Each e lement  x(r )  i s  tagged with  arank=kr(mod  N) before it is sent into the array. Forexample, for k=2 and N=4, the samples become(x(O), O),  (x(I), 2), (x(2),  O),  (x(3), 2). These aresent into the linear array one element at a time. Anelement is percolated from left to right until itreaches a cell in the array whose position isidentical to its rank. When this happens, the valueis accumulated in that cell in case it receivesmultiple values (i.e. Su >l). Figure 1 illustratesI Ian example. A cell in a linear array contains anadder only if the corresponding partition Su has asize greater than 1. In 8 clocks, x2,o  will emerge atthe right hand side, as shown in Figure 1. Ingeneral, it takes 2N clocks for the vector to appearon the right. It is obvious that number of additionsrequired in the preprocessing array for each k isequal to N-N/gcd(k,N). So the total number ofaddition required in the preprocessing array is equalN-I N-lto x(N-Nlgcd(k,N)=N’-xNlgcd(k,N).t-0 I=0x(r) to xk,i is noteworthi: If (k, N)=l, (i.e., k and Nare coprime) then (ki i=O, 1, 2, . . . . . . . N-l } is acomplete system of residues modulo N [ 181  and 1 SuI = 1. It means xk,i is a permutation o It remains toconsider the generation of xk,i . One interestingproperty of the mapping off x(r). Particularly, if Nis prime, then Sm={O,  1, . . . N-l} and Saj=B for j#0, but I Ski I =l for all k#O. In other words, xk,i is apermutation of x(r) whenever k#O. In theremaining case, Y,(O) = yx(r) and Q,j=O for j#O. Sor=0in the case where the sample size is a primenumber, a preprocessor can be used to produce thexk,i via permutations of the sampled sequence x(r),after producing Y,(O) = %x(r) in the first pass.r-05: Systolic array for computing DHTThe general network for implementing DHT isdrawn in Figure 2. It consists of the momentgenerator [lo], the processing arrays and a dotproduct array. The moment generator formed of N-2 (2p+l)-networks  with a row of adder-latch couldbe used to generate the 1-D moments. It receivesxk,i, (i, k=O, 1, . . . N-l) that preprocessing arraysproduce. The dot product array performs the dotproduct and produces X(k) at the only output. Thenumbers in square brackets denote the amount oftime delay to keep synchronous pace in Figure 2.The scheduling of the dataflow  is such that mk (0)*aa moves from left to right, accumulates with mk (l)*al to produce the partial sum and this repeats untilX(k) emerges at the far right of the output. As mk(i)is produced earlier one cycle than mk(i+I) (i=O,l,..., 2p+l),  as shown in Figure 2, so the most right180term mk(p)  takes 2 cycles to merge into X(k). It isjust time cost for the dot product.The total execution time of the computingprocedure is:clock cycles. Where N is for preprocessing, and(2p+2)(N-2)  for computing moments, and 2 for thedot product, and N-l for collecting all X(k), (k=l,2, . . . . . . N-l).Since the DHT is self-inverse, the approach isalso applicable to DHT inverse. The method canalso be extended to to multi-dimensional DHT dueto the separability of the DHT kernel. It can alsobe proved by mathematical induction that theexecution time of the systolic array is O(Nk)  for k-DDHT.6: ConclusionsNew fast algorithms and systolic arrays forone-dimensional DHT and their inverse areproposed and the results are extended to k-dimensional DHT (k22)  , and their computationcomplexities are analyzed. The relationshipsbetween these discrete transforms and discretemoments are presented and the fast discretetransforms are converted into fast discrete momenttransforms. It is a completely new approach toDHT.The methods proposed have the advantagesbelow:1. Many multiplications have been converted intoadditions, and triangular functions have beenreplaced by simple polynomial functions. Itdecreases the computational cost and memoryrequirement. Comparing with the direct methodwith computation complexity 0(N2) and theconventional FHT which originates from FFT with0(Nlog2N), the methods proposed are superior tothe direct method and but inferior to FHT) if thecalculations are executed on a sequential machine(for the example given in Section 2: N<_2k,  logzN<11, p=17).  When N is large enough, p can take avalue much less than 1ogzN  and can still satisfy theaccuracy requirement. In other words, in the caseof N large enough, the number of multiplication inour method is less than that in the conventionalFHT though the number of addition in our methodis still more than that in the conventional FHT;2. The systolic arrays consist of only latches, adder-latches and a few multipliers. It results in easyhardware implementation and is very suitable for areal-time system;1813. It can cope with arbitrary data length and notlimited to data length in the form of 2k or any otherforms. It has a very good accuracy andconvergence feature;4. The multiplication number for 1-D DHT in ourmethod is 0(Nlog2N/log210g2N)  superior to O(Nlog2N)  in the conventional FHT;5. The methods are very easily extended to multi-dimensional DHT and their inverses. Thecomputational complexity is O(nk)  for k-dimensional DHT in the systolic arrays.AcknowledgementThis project has been supported by NationalNatural Science Foundation of China under grant69775011 and by University of Hong Kongresearch grants.111PI[31[41151161[71ReferencesR. N. Bracewell, The Hartley Transform, NewYork: Oxford University Press, 1986.S.-C. Pei and I.-I. Yang, “Computing pseudo-Wigner distribution by the fast Hartleytransform,” IEEE Transactions on SignalProcessing, vol. 40, PP. 2346-2349,September, 1992.N.-C. Hu, H.-I Chang, and 0. K. Ersoy,“Generalized discrete Hartley transforms,”IEEE Transactions on Signal Processing, vol.40, pp. 2931-2940, December, 1992.R. N. Bracewell, “The fast Hartley transform,”Proc.  IEEE, vol. 22, pp. 1010-1018, August,1984.P. Duhamel and M. Vetterli, “ImprovedFourier and Hartley transform algorithns:application to cyclic convolution of realdata,” IEEE Transactions on Acoustics,Speech, and Signal Processing, vol. 34, pp.642-644, 1986.H. V. Sorensen, D. L. Jones, C. S. Burrus, andM. T. Heideman, “On computing the discreteHartley transform,” IEEE Transactions onAcoustics, Speech, and Signal Processing,vol. 33, pp. 1231-1238, 1985.H. S. Hou, “The fast Hartley transformalgorithm,” IEEE Trans. Comput., vol. 36, pp.147-156, February, 1987.0. Buneman, “In-situ bit-reversed ordering forHartley transforms,” IEEE Transactions onAcoustics, Speech, and Signal Processing,vol. 37, pp. 577-580, 1989.[9] 0 .  B u n e m a n , “Multidimensional Hartley “Fast algorithm for computation of momenttransforms,” Proc.  IEEE, vol. 75, pp. 267, invariants,” Pattern Recognition, vol. 20, pp.February, 1987. 634-643, 1987.[lo]  Chan, F. H. Y., Lam, F. K., Li, H. F., and Liu,J. G. “An all adder systolic structure for fastcomputation of moments,” Journal of VLSISignal Process. Vol. 54, pp. 159-175,  1996.[ 1 l] Hatamian, M. “A real time two-dimensionalmoment generating algorithm and its singlechip implementation,” IEEE Trans. Acoust.Speech Signal Process. vol. 34, pp. 546-553,1986.[14]  Olmsted, J. M. H. Real Variables, Appleton-Century-Crofts, Inc. 1956.[15]  J. G. Liu, H. F. Li, F. H. Y. Chan and F. K.Lam, “Fast discrete cosine transform viacomputation of moments,” Journal of VLSISignal Processing Vol. 19, No.3,  257-268,1998.[12]  Li, B.-C., and Shen, J. “Pascal triangletransform approach to the calculation of 3-Dmoments,” Computer Vision, Graphics, andImage Processing: Graphical Model andImage Processing, vol. 54, pp. 301-307, 1992.[ 131 Zakaria, M. F., Vroomen, L. J., Zsombar-Murray P. J. A., and Van Kessel, J. M. H. M.[ 161 J. G. Liu, H. F. Li, F. H. Y. Chan and F. K.Lam, “A novel approach to fast discreteFourier transform,” Journal of Parallel andDistributed Computing, vol. 54, pp. 48-58,1998.[ 17 ] Beckenbach, E. F. and Bellman, B.Inequalities. Springer-Verlag, 196 1.[18]  Rose, H. E. A Course in Number Theory,Oxford University Press Inc., 1994.cell # 3X2,i0 x(1)+x(3) 0 x(0)+x(2) >3 2 1 0Figure 1. The linear array to implement x2.i.182X,N.l I........  X,ll x,0,I I\Ru : latch0 : addeflati0 : multiplierFigure 2. The systolic array for 1-D DHT.183

A novel approach to fast discrete Hartley transform

http://hub.hku.hk/bitstream/10722/46134/1/45501.pdf

A novel approach to fast discrete Hartley transform

Abstract

Similar works

Full text

Available Versions

HKU Scholars Hub