MULTIPLIERLESS LIFTING BASED FFT VIA FAST HARTLEY TRANSFORM
Taizo Suzuki1 , Seisuke Kyochi2 , Yuichi Tanaka3 , Masaaki Ikehara4 , and Hirotomo Aso5
Faculty of EIS, Univ. of Tsukuba, 2 IME Dept., The Univ. of Kitakyushu,
Graduate School of BASE, TUAT, 4 EEE Dept., Keio Univ., 5 EEE Dept., Nihon Univ., Japan
Email: 1 taizo@cs.tsukuba.ac.jp, 2 s-kyochi@kitakyu-u.ac.jp,
3
ytnk@cc.tuat.ac.jp, 4 ikehara@tkhm.elec.keio.ac.jp, 5 aso@ee.ce.nihon-u.ac.jp
13ABSTRACTThe multiplierless fast Fourier transform (FFT) with dyadic-valued(rational) coefficients is important for many signal processing tools.
The proposed lifting based FFT (L-FFT) based on fast Hartley transform (FHT) has a simpler structure than existing ones because fewer
lifting steps need to be approximated. In addition, it has a structureof real-valued calculation followed by complex-valued parts, thereby
it requires fewer memories for the internal implementation than theconventional FFTs.

In addition, since it has a structure of real-valued calculation followed by complex-valued parts, it requires fewer memories for theinternal implementation than the conventional ones. In the simulations, we evaluate the L-HFFT and numerically show its comparable
performance to the L-SFFT in spite of fewer operations.
Notations: I, J, MT , M† , MN and j are an identity matrix, areversal identity matrix, transpose of a matrix M, conjugate transpose of a matrix M, a matrix M with N × N size and an imaginary
number defined as j 2 = −1, respectively.

Index Terms— Dyadic-valued coefficient, fast Fourier transform (FFT), fast Hartley transform (FHT), lifting structure
2.1. DFT and DHT

1. INTRODUCTION

The (m, n)-element of M -channel DFT matrix FM is defined as [2]The fast Fourier transform (FFT) [1], which is a fast algorithm of
discrete Fourier transform (DFT) [2], is one of the most fundamental and common tools in signal processing [3–5]. Generally, theFFT consists of floating-point coefficients and yields floating-point
output signals even if the input signals are integer-valued. In thiscalculation, floating point multipliers waste computational cost andpower consumption which cannot be neglected. To efficiently implement the FFT, floating-point coefficients are often approximatedto dyadic-valued (rational) coefficients. It is the most important issue to reduce its power consumption. Ideally, such approximation
should be carried out without losing perfect reconstruction (PR) or
near PR property.
Lifting structure [6] is a key technology to efficiently approximate floating-point coefficients by dyadic-valued coefficients while
keeping near PR [7]. Oraintara et al. introduced the lifting factorization of a traditional (split-radix based) FFT and realized a multiplierless lifting based FFT (L-FFT) with dyadic-valued coefficients [8].
However, the structure does NOT consider minimizing the numberof lifting steps for fewer adders and bit shifters, and thus it has redundant operations. Moreover, it requires many memories for theinternal implementation because there are many complex-valued operations everywhere.
In this paper, we propose a multiplierless L-FFT based on thefast Hartley transform (FHT) [9] for one and two dimensional (1-Dand 2-D) signals. The FHT is a fast algorithm of discrete Hartleytransform (DHT) [10] which is a DFT like real-valued transform.
The conventional and proposed multiplierless FFTs are named spritradix based L-FFT (L-‘S’FFT) and FHT based L-FFT (L-‘H’FFT),
respectively, to distinguish them throughout this paper. The L-HFFT,
which has a different structure from the traditional radix based FFTs,
has a simpler structure than that of the L-SFFT. Specifically, it can
reduce the number of lifting steps which need to be approximated.

978-1-4799-0356-6/13/$31.00 ©2013 IEEE

2. REVIEW

5603

1[FM ]m,n = √ exp
M(−j2mnπM)

where 0 ≤ m, n ≤ M −1. The inverse DFT (IDFT) is the conjugate†−1transpose of the DFT matrix, i.e., FM= FM. For simplicity, letkM be defined as M = 2 (k ∈ N).
The DHT [10] is a DFT like real-valued transform. The (m, n)element of M -channel DHT matrix HM , which is the symmetric−1Torthogonal matrix as HM= HM= HM , is defined as
1[HM ]m,n = √M((cos2mnπM)(+ sin2mnπM)).

(1)

These DFT and DHT have many fast algorithms, i.e., FFT [1]and FHT [9], respectively, which are usually used in practical implementation.
2.2. Multiplierless Lifting Structure
The lifting structure [6], also known as the ladder structure, is a special type of lattice structure. It is implemented by cascading elementary matrices ― identity matrices with a single nonzero off-diagonal
element.
Fig. 1 shows a basic lifting step. It is expressed byyj (n) = xj (n), yi (n) = xi (n) + round{T xj (n)}zj (n) = yj (n), zi (n) = yi (n) − round{T yj (n)}where round{·} and T are a rounding operation and a lifting coefficient, respectively. Thus the lifting structure with the rounding

ICASSP 2013

Fig. 1. A basic lifting step (white circles: rounding operations)Fig. 2. An approximation from multiplier with dyadic-valued coefficients to bit-shifter and adder in a lifting structure (≫ β means β
bit-shifter).

operation can achieve integer-to-integer transform. Also, the lifting
and its inverse matrices in this case are represented as[][]−1 []
1 T
1 T
1 −T
and=,0 10 101
respectively.
For a high-speed and circuit-size saving implementation, lifting
coefficients are required to approximate floating-point coefficients tosoftware/hardware-friendly dyadic-valued coefficients such as α/2β
(α, β ∈ N) which can be implemented by only adders and bitshifters [7]. The dyadic-valued representation enables to perform
fast implementation in a software and reduce the circuit size. For
example, a coefficient 25/128 can be operated as252 4 + 23 + 20111== 3 + 4 + 7.
128
27
222Hence, the lifting with its coefficient 25/128 and a rounding operation is replaced to the summation of 3, 4 and 7 bit-shifters illustrated in Fig. 2. It is clear that the PR in lifting is always kept even
if floating-point coefficients are approximated to dyadic-valued coefficients. When a floating-point coefficient is allocated β bit, thedyadic-valued coefficient is presented as α/2β .

Fig. 3. The 16-channel FHT Ĥ16where ĤM is an FHT permutated frequency band order of HM , i.e.,T
HM = PTM ĤM = ĤM
PM where PM is a permutation matrix.
Also, RM and WM are expressed asIM1 0RM = 00[1 IM/2
WM = √2 IM/2

FHT in (1) is simplified by rotation matrices as the following recursive algorithm and its 16-channel case is shown in Fig. 3.
[]
ĤM/20ĤM =RM WM
0ĤM/25604]
IM/2,−IM/2

respectively. The (l, l)-elements of CM2 and SM2 are presented by(3. MULTIPLIERLESS L-FFT VIA FHT (L-HFFT)

3.1. FFT based on FHT

0JM2 SM2
0−JM2 CM2 JM200
10for M1 = M/2 + 1 and M2 = M/4 − 1, and[CM2 ]l,l = cos
In this section, we derive a realization of multiplierless FFT for various signal processing tools by the following procedures.
1. An FFT factorization is derived based on FHT.
2. Several obtained rotation matrices are approximated to multiplierless lifting structures, and the scaling factors are moved
to the last part.
(For 2-D signals, e.g., images)
3. It is completed by using 2-D separable block transform.

0CM20SM2 JM22(l + 1)πM)(and [SM2 ]l,l = sin2(l + 1)πM),

respectively, for 0 ≤ l ≤ M/4 − 2.
A representation of M -channel FFT FM is constructed by only
one FHT HM and the matrix GM with real- and complex-valued
coefficients as follows [9]:T
FM = QM GM HM = QM GM ĤM
PM(2)where

1
0
GM = 0001−j
IM/2−1201+j
JM/2−1200100

1+j
JM/2−1 20

1−j
IM/2−12and QM is a permutation matrix. Also, the inverse of GM is pre†sented by the conjugate transpose of GM , i.e., G−1
M = GM .

3.2. Multiplierless Lifting Factorization and Moving Scaling
Factors to The Last Part
First, the rotation matrices in RM and the complex matrices in GM
are factorized into the lifting structures as[] [][][]
1010cos θsin θ1 sin θ= 1−cos θcos θ−1sin θ − cos θ1−1 01sin θsin θand[1−j 1j2] [j1=110−j][

101+j2][11−1]0,1respectively. Next, floating-point coefficients in the lifting rotation
matrices in RM are approximated to dyadic-valued coefficients.
Note that the lifting steps in GM already have only dyadic-valuedcoefficients.
On the other hand, the rotation matrices with the angles π/4 inTFHT part ĤM
except for RM are presented by[][]
1 1cos π4
sin π41.
ππ = √
sin 4 − cos 4
2 1 −1√Since all of the scaling factors 1/ 2 can be moved to the last part ofFFT, the matrices can be implemented by only ±1. Let H̃M be ĤMremoved its scaling factors. Then, (2) is rewritten asT
FM = QM DM GM H̃MPM(3)where DM is the scaling factor matrix. The resulting FFT is shownin Fig. 4.
3.3. Use of 2-D Separable Block Transform
Note that the structure obtained in the above subsection√products
the non-dyadic-valued (floating-point) scaling factor 1/ 2 whenlog2 M is odd. Moreover, if the structure in Fig. 4 is directly applied to image, the complex-valued signals produced after the firstdimensional (i.e., column or row) implementation will be input tothe second dimensional (i.e., row or column) operation. Hence, thestructure requires many memories for the internal implementation.
To solve the problem, we change the order of the 2-D separable
block transform [7] for the proposed FFT. When we apply a blocktransform matrix FM into a 2-D input signal x in column- androw-wise separately, the 2-D output signal y is expressed byT
y = (FM (FM x)T )T = FM xFM(4)TSince FM is factorized as DM GM H̃M
PM as shown in (3), (4) isrepresented byTy = QM DM GM H̃M
PM xPTM H̃M GTM DTM QTM .
TThis equation means that the 2-D block transform by ĤM
is firstlyapplied after the 2-D block transform by PM , the 2-D block transform by GM is secondly applied, the 2-D block transform by DM
is thirdly applied, and the 2-D block transform by QM is finally applied. Therefore,√ when log2 M is odd, the non-dyadic-valued scaling factor 1/ 2 in 1-D transform can be implemented as the dyadicvalued scaling factor 1/2 in 2-D transform. Moreover, it is clear that
the FFT requires fewer memories for the internal implementation
because it can be completely separated to real- and complex-valued
parts.

5605

Table 1. The FFTs for their transform performance comparison (Thesplit-radix based FFT with dyadic-valued coefficients without lifting
factorization is named ‘S’FFT).
Name
Basement
Lifting Dyadic-Value(ideal)
Split-Radix
NANA
SFFT
Split-RadixNA✓
L-SFFT Split-Radix✓✓
L-HFFT
FHT✓✓
Table 2. The numbers of operations in FFTs with 8 bit word length
allocated coefficients (Lift., Add. and Shift. mean the numbers oflifting steps which require finite word length allocation, adders andbit-shifters, respectively).
ChannelL-SFFT
L-HFFT
M
Lift.
Add.
Shift.
Lift.
Add.
Shift.
16
24
178
106
15
151
73
32
78
518
332
51
424
234
64
216
1366
910
147
1098
652
128
486
3356
2298
387
2689
1667
256
1296
8114
5634
963
6395
4093
512
2430 18710 13292 2307 14836
9718
1024
6552 43116 30692 5379 33765 22503*For simplicity, the scaling parts are NOT included in this table.
4. EXPERIMENTAL RESULTS
The FFTs for their transform performance comparison are shownin Table 1. In Sec. 4.2 and 4.3, we used 256-channel FFTs and256 × 256 8-bit grayscale images Camera and Girl as test images.
4.1. Number of Operations in FFTs
Table 2 shows the number of lifting steps which must be approximated, adders and bit-shifters in the L-SFFT and the L-HFFT approximated by 8 bit allocation. It is clear that the L-HFFT has not
only no multiplier but fewer adders and bit-shifters than the L-SFFT.
It is no wonder that the L-HFFT is constructed by fewer lifting steps
than the L-SFFT. Let the numbers of operations in M -channel LSFFT and L-HFFT be defined as SM and HM , respectively. Theyare recursively determined as follows:3(M − 4)
+ S M + 2S M
2423(M − 4)
HM =+ 2H M
24
SM =where S2 = S4 = H2 = H4 = 0.
4.2. Accuracy of FFT
The images Camera and Girl, which are reconstructed by the ideal
inverse FFT (IFFT) after they are transformed by each of SFFT, LSFFT and L-HFFT, are compared with peak-to-noise ratio (PSNR)
as()
2552
PSNR[dB] = 10 log10MSEFig. 4. The 16-channel L-HFFT

Fig. 5. F-D curves of the images reconstructed by ideal IFFT after
they are transformed by each of SFFT, L-SFFT and L-HFFT: (left)
Camera, (right) Girl.

Fig. 6. F-D curves of the images reconstructed by each of the inverse
transforms of SFFT, L-SFFT and L-HFFT after they are transformedby each of SFFT, L-SFFT and L-HFFT: (left) Camera，(right) Girl4.3. Accuracy of Near PR
Table 3. The comparison examples in the images reconstructed byideal IFFT after they are transformed by each of L-SFFT and LHFFT when M = 256.
Numbers
PSNR[dB]
Adders Shifters Camera
GirlL-SFFT (4 bit alloc.)
5,138
3,090
41.69
44.58
L-HFFT (5 bit alloc.)
5,031
2,983
43.64
46.10

Test images Camera and Girl, which are reconstructed by each of
the inverse transforms of SFFT, L-SFFT and L-HFFT after they aretransformed by each of SFFT, L-SFFT and L-HFFT, are comparedwith PSNR. This simulation clarifies an accuracy of reconstructionof each dyadic-valued FFT because the complete PR is lost due toerror generated by bit-shifts in the last scaling part. Fig. 6 shows
F-D curves at the case. Although the SFFT degrades its PSNRs inlower bit allocation, such degrading is NOT shown in the L-SFFT
and the L-HFFT. Hence, the FFTs with lifting structures always preserve near PR.

where MSE is the mean squared error. This simulation clarifies anapproximation accuracy of the dyadic-valued FFTs to the ideal FFT,
under the rounding operation. Fig. 5 shows F-D (finite word length―distortion) curves at the case. Although the L-HFFT seemingly
shows lower PSNRs, the L-HFFT shows comparable performanceto the L-SFFT in spite of fewer operations in practice as shown inTable 3. Also, the L-HFFT requires fewer memories for the internal implementation than the L-SFFT due to real-valued calculation
followed by complex-valued parts.

5. CONCLUSION

5606

This paper has proposed a realization of multiplierless fast Fouriertransform (FFT). It is composed of fast Hartley transform (FHT) andmultiplierless lifting structures, and has only adders and bit-shifters,i.e., no multipliers. Moreover, it has fewer operations and requires
fewer memories for the internal implementation than the conventional FFTs. In spite of such facts, the proposed FFT verified its
comparable performance to the conventional ones in the simulations.

6. REFERENCES[1] J. W. Cooley and J. W. Tukey, An algorithm for the machine
computation of complex Fourier series, Math. Comput., 1965.
[2] N. Ahmed and K. R. Rao, Orthogonal Transforms for Digital
Signal Processing, Berlin, Germany: Springer-Verlag, 1975.
[3] R. Hardie, “A fast image super-resolution algorithm using an
adaptive Wiener filter,” IEEE Trans. Image Process., vol. 16,
no. 12, pp. 2953–2964, Dec. 2007.
[4] M. Ghazal, A. Amer, and A. Ghrayeb, “Structure-oriented
multidirectional Wiener filter for denoising of image and video
signals,” IEEE Trans. Circuits Syst. Video Technol., vol. 18,no. 12, pp. 1797–1802, Dec. 2008.
[5] D. Humphrey and D. Taubman, “A filtering approach to edge
preserving MAP estimation of images,” IEEE Trans. Image
Process., vol. 20, no. 5, pp. 1234–1248, May 2011.
[6] W. Sweldens, “The lifting scheme: A custom-design construction of biorthogonal wavelets,” Appl. Comput. Harmon. Anal.,
vol. 3, no. 2, pp. 186–200, Apr. 1996.
[7] T. Suzuki and M. Ikehara, “Integer discrete cosine transform
via lossless Walsh-Hadamard transform with structural regularity for low-bit-word-length,” IEICE Trans. Fundamentals,
vol. E93-A, no. 4, pp. 734–741, Apr. 2010.
[8] S. Oraintara, Y. J. Chen, and T. Q. Nguyen, “Integer fastFourier transform,” IEEE Trans. Signal Process., vol. 50, no.
3, pp. 607–618, Mar. 2002.
[9] R. N. Bracewel, “The fast Hartley transform,” Proc. of theIEEE, vol. 72, no. 8, pp. 1010–1018, Aug. 1984.
[10] R. N. Bracewel, “Discrete Hartley transform,” J. Opt. Soc.
Am., vol. 73, no. 12, pp. 1832–1835, Dec. 1983.

5607