An analysis of real-Fourier domain based adaptive algorithms
implemented with the Hartley transform using cosine-sine
symmetries
Citation for published version (APA):
Raghavan, V., Prabhu, K. M. M., & Sommen, P. C. W. (2005). An analysis of real-Fourier domain based
adaptive algorithms implemented with the Hartley transform using cosine-sine symmetries.
IEEE Transactions
on Signal Processing
,
53
(2), 622-629. https://doi.org/10.1109/TSP.2004.838983
DOI:
10.1109/TSP.2004.838983
Document status and date:
Published: 01/01/2005
Document Version:
Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)
Please check the document version of this publication:
• A submitted manuscript is the version of the article upon submission and before peer-review. There can be
important differences between the submitted version and the official published version of record. People
interested in the research are advised to contact the author for the final version of the publication, or visit the
DOI to the publisher's website.
• The final author version and the galley proof are versions of the publication after peer review.
• The final published version features the final layout of the paper including the volume, issue and page
numbers.
Link to publication
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners
and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.
• Users may download and print one copy of any publication from the public portal for the purpose of private study or research.
• You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal.
If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please
follow below link for the End User Agreement:
www.tue.nl/taverne
Take down policy
If you believe that this document breaches copyright please contact us at:
openaccess@tue.nl
providing details and we will investigate your claim.
Download date: 10. Aug. 2022
622 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 2, FEBRUARY 2005
An Analysis of Real-Fourier Domain-Based Adaptive
Algorithms Implemented with the Hartley Transform
Using Cosine-Sine Symmetries
Vasanthan Raghavan, Student Member, IEEE, K. M. M. Prabhu, Senior Member, IEEE, and
Piet C. W. Sommen, Senior Member, IEEE
Abstract—The least mean squared (LMS) algorithm and its
variants have been the most often used algorithms in adaptive
signal processing. However the LMS algorithm suffers from a
high computational complexity, especially with large filter lengths.
The Fourier transform-based block normalized LMS (FBNLMS)
reduces the computation count by using the discrete Fourier trans-
form (DFT) and exploiting the fast algorithms for implementing
the DFT. Even though the savings achieved with the FBNLMS
over the direct-LMS implementation are significant, the compu-
tational requirements of FBNLMS are still very high, rendering
many real-time applications, like audio and video estimation,
infeasible. The Hartley transform-based BNLMS (HBNLMS) is
found to have a computational complexity much less than, and
a memory requirement almost of the same order as, that of the
FBNLMS. This paper is based on the cosine and sine symmetric
implementation of the discrete Hartley transform (DHT), which is
the key in reducing the computational complexity of the FBNLMS
by 33% asymptotically (with respect to multiplications). The
parallel implementation of the discrete cosine transform (DCT) in
turn can lead to more efficient implementations of the HBNLMS.
Index Terms—Adaptive algorithms, DCT, DST, FBNLMS, FFT,
FHT, frequency domain algorithms, LMS algorithms.
I. INTRODUCTION
T
HE ability of an adaptive filter to operate effectively in
an unknown environment and track time variations of the
input statistics makes it a powerful device. It finds wide applica-
tions in diverse fields such as interference and echo cancellation,
channel equalization, linear prediction, and spectral estimation.
The adaptive filter is implemented in the time domain by
algorithms like the least mean squared (LMS) algorithm [1],
the gradient adaptive lattice [2], or the least squares algorithms
[3], [4]. However, as the number of computations per sample
needed to implement the LMS algorithm grows linearly with the
Manuscript received June 16, 2003; revised January 9, 2004. The associate
editor coordinating the review of this paper and approving it for publication was
Prof. Trac D. Tran.
V. Raghavan is with the Department of Electrical and Computer Engi-
neering, University of Wisconsin-Madison, Madison, WI 53706 USA (e-mail:
raghavan@cae.wisc.edu; vasanth_295@yahoo.com).
K. M. M. Prabhu was with Eindhoven University of Technology, 5600 MB
Eindhoven, The Netherlands. He is now with the Department of Electrical Engi-
neering, Indian Institute of Technology Madras, Chennai 600036, India (e-mail:
prabhu_kmm@hotmail.com).
P. C. W. Sommen is with the Department of Electrical Engineering, Eind-
hoven University of Technology, 5600 MB Eindhoven, The Netherlands (e-mail:
p.c.w.sommen@tue.nl).
Digital Object Identifier 10.1109/TSP.2004.838983
filter order, it is less useful in real-time applications like audio
and video estimation. Many techniques have been proposed to
reduce this burgeoning computation count. Transforming the
LMS algorithm to the frequency domain [5], [6] has been the
classical approach, as we can exploit the various forms of fast
Fourier transform (FFT) algorithms to reduce the number of
computations. With the generalization of the transform domain
LMS algorithm [7], we find that we can use orthogonal trans-
forms to yield better computation counts.
Merched
et al. [8]
1
have introduced the implementation of the
FBNLMS in the Hartley domain (the so-called “HBNLMS”).
This method involves extending the data matrices to circulant
symmetric matrices and then using the Hartley transform to di-
agonalize them. The main motive of [8] is to find the com-
monality between adaptive filters using multidelay concepts and
those using filterbanks, whereas in this paper, our main motive
is the efficient implementation of the HBNLMS algorithm. We
exploit the symmetry of the Hartley transform in being decom-
posed into a cosine and a sine transform of reduced order and the
symmetrical structure of the HBNLMS to significantly reduce
the computational complexity of the FBNLMS. This is possible
entirely because the sine transform of a sequence can be imple-
mented using the cosine transform of its alternate sequence [9].
The asymptotic computational savings with the HBNLMS is
20% with respect to multiplications and additions and 33% with
respect to multiplications alone. Since fast cosine transforms
are commonly used in MPEG coding and are efficiently im-
plemented in special VLSI chips, decomposition of the Hartley
transform to a cosine and a sine transform also has other inherent
advantages that are associated with VLSI designs. We also note
that the HBNLMS requires almost the same order of memory as
that of the FBNLMS implementation.
The Hartley transform has been used in the literature for fil-
tering in [10] and [11], but the authors are not aware of an ef-
ficient implementation of the Hartley transform routines based
on sine and cosine symmetric decompositions in the context of
LMS filters. Merched et al. [8] use the regular Hartley trans-
form in their work without exploiting the symmetries of the
sequences involved. Thus, this work is an efficient implemen-
tation of the proposed Hartley transform-based frequency do-
main LMS. The paper is organized as follows. The Fourier do-
main-based block normalized LMS algorithm is described in
1
We would like to thank one of the reviewers for introducing us to the results
of [8].
1053-587X/05$20.00 © 2005 IEEE
Authorized licensed use limited to: Eindhoven University of Technology. Downloaded on September 11, 2009 at 04:32 from IEEE Xplore. Restrictions apply.
RAGHAVAN et al.: ANALYSIS OF REAL-FOURIER DOMAIN-BASED ADAPTIVE ALGORITHMS 623
Section II. The Hartley transform-based BNLMS (HBNLMS) is
introduced in Section III. Sections IV and V deal with the imple-
mentation issues of the HBNLMS and the computation counts
of the discrete Fourier transform (DFT) and discrete Hartley
transform (DHT)-based algorithms, respectively. The memory
requirements of both the algorithms are compared in Section VI
and concluding remarks are provided in Section VII.
A. Notation
We assume real inputs, real “desired” outputs and we start
with a real
. This will ensure that the weight vectors
over all iterations are real. The operations
, and
refer, respectively, to the transpose, Hermitian, and the
conjugate of
, respectively. The operation refers to ele-
mentwise multiplication of two sequences. The superscripts in
and refer, respectively, to “circulant” and “cir-
culant symmetric.”
, , and refer to the Fourier, the
Hartley, the inverse Fourier and the inverse Hartley matrices,
respectively. X(k) and X_H(k) refer, respectively, to the DFT
and the DHT of the sequence x(n). The computation count of
an algorithm is defined as the number of multiplications and
additions, since in most modern-day digital signal processors
(DSPs), the time needed for implementing an addition is almost
the same as that for a multiplication.
The tap inputs of the LMS filter are denoted by
and the tap weights by .
The output, the desired output, and the error in the output are de-
noted by y(n), d(n), and e(n), respectively. Here, M is the order
of the filter used. The filter output y(n) is given by the convolu-
tion sum
(1)
This output y(n) is used to estimate the desired re-
sponse d(n). Let
and
. We do not use an explicit
dependence of the weight vector and the error vector on the
discrete-time variable “n” as we are dealing with a particular
iteration, and for that iteration, these vectors can be assumed to
be constant. In the block algorithm (where B is the filter block
length), we define the block input
(an M B matrix) as
(2)
Then, the block output
is defined as ,
where
(3)
The estimation error
is defined as
(4)
refers to the adaptation constant in the adaptive algorithm.
II. FBNLMS A
LGORITHM
The Fourier transform-based block normalized least mean
square (FBNLMS) algorithm [12], [13] is an efficient way of im-
plementing the block normalized least mean squares (BNLMS)
algorithm [12]. The weight update equation of the BNLMS al-
gorithm is
(5)
Here,
is an estimate of block variance of . The structure
of
allows us to extend it circulantly as follows:
(6)
where
, , and are chosen appropriately. The first row of
is given as
(7)
Here, the size of
(a square matrix of N N), which is
denoted by N, is constrained by
. Usually,
N is chosen to be a power of two so that the FFTs can be used
to reduce computation (described later in this section). Further
manipulation is made easy if we notice that
is a sub-
matrix of
, where . This
submatrix can then be extracted from the product
by
premultiplication with a M
N windowing matrix defined
as
. The following algebraic manipu-
lations then follow trivially:
(8)
It can be seen that
diag , where we exploit the realness of the input vector
and the fact that Fourier matrices diagonalize circulant matrices
[14]. This, therefore, leads to
(9)
where the symbol
denotes the elementwise multiplication of
the two N-point sequences. The weight update equation of the
BNLMS algorithm then becomes
(10)
can be written as follows:
(11)
where
, and .
Equation (11) can then be rewritten similarly as
(12)
Authorized licensed use limited to: Eindhoven University of Technology. Downloaded on September 11, 2009 at 04:32 from IEEE Xplore. Restrictions apply.
624 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 2, FEBRUARY 2005
Since is a sequence that is well-known even before the start
of iterations, we can assume that it is a constant sequence and
therefore,
is a known sequence (or can be computed of-
fline). We therefore need three forward Fourier transforms, two
inverse Fourier transforms, two elementwise multiplications, a
variance estimator and two windows for each iteration. We do
not concern ourselves with the dynamic behavior (convergence
or tracking properties) of the FBNLMS or the Hartley trans-
form-based algorithm that will be introduced in Section III, be-
cause both the algorithms are in principle the same when real
inputs and weights are used. The main motive of this paper is
to exploit the structure that lies hidden in the vectors used in
HBNLMS, in order to reduce the computation count. However,
we would like to comment that, since the FBNLMS algorithm is
an exact transformation of the BNLMS algorithm, the dynamic
behavior of both of them will be essentially the same [12].
III. H
ARTLEY
TRANSFORM-BASED BNLMS (HBNLMS)
The DHT [15], [16] of a real valued N-point sequence x(n) is
defined as follows:
cas
(13)
where, cas
. A circulant matrix can be
diagonalized by a Fourier matrix, viz.
diag ,
where
is the DFT of the first column of [14]. A circulant
symmetric matrix
can be diagonalized by a Hartley matrix.
diag (14)
where
is the first column of C. The FBNLMS algorithm can
then be implemented by extending all data matrices to be cir-
culant and symmetric and using
to diagonalize them [8]. The
order of all circulant symmetric matrices should satisfy the con-
straint
. We will abuse notation and use
(where N stands for the size of the extended circulant
matrix in the FBNLMS) henceforth. The following equations
can then be seen to be true for HBNLMS [8].
where and
.
(15)
where
, ,
and
is the DHT of the first column of . We then
have
(16)
Since
is a sequence that is well known even before the
start of iterations, we can assume that its transform is a known
sequence. Thus, three Hartley transforms and two inverse
Hartley transforms are needed to compute
in the
frequency domain. These Hartley transforms can be computed
using the more efficient DCT and discrete sine transform (DST)
systems. It might initially appear that the HBNLMS has an
inherent disadvantage because of the order of the transforms
needed here. For every Fourier transform (and inverse) of order
N, we need a Hartley transform (and inverse) of order 2N. The
HBNLMS would be computationally disadvantageous when
compared to the FBNLMS (if the symmetries in the sequences
are not used), as is done in [8], where a Hartley transform
of order 2N is implicitly used. The main motive of this work
is to exploit the symmetries inherent in the sequences to use
the computationally more efficient DCT and DST systems to
recursively compute the needed transforms. In this process, we
reduce the needed computations even below that needed by the
FBNLMS (with
real-FFT’s used for a fair comparison). The
symmetry in the sequences helps us in replacing a 2N-point
Hartley transform with cosine transforms of order of N/2 down
until 1 [17], [18].
IV. I
MPLEMENTATION OF THE
HBNLMS
When real data vectors are used, the HBNLMS is just another
mathematically equivalent way of implementing the FBNLMS.
The basic idea of our implementation is that the recursive imple-
mentation of the DCT-I, DST-I, DCT-II and DST-II involve the
same building blocks with signs changed appropriately to give
the needed transforms. The first column of
, sat-
isfies the following properties:
(17)
This then implies the following about
:
cas
(18)
Authorized licensed use limited to: Eindhoven University of Technology. Downloaded on September 11, 2009 at 04:32 from IEEE Xplore. Restrictions apply.
RAGHAVAN et al.: ANALYSIS OF REAL-FOURIER DOMAIN-BASED ADAPTIVE ALGORITHMS 625
Here, is defined as
(19)
This implies that all we need to do to find the 2N-point Hartley
transform is to find the DCT-I [19] of the sequences
and .
The DCT-I of an N-point sequence is computed recursively
using the DCT-I of the even sampled sequence and the DCT-II
of the odd sampled sequence [19], as follows:
(20)
If A and B are the DCT-I of the even samples and DCT-II of
the odd samples of
and and are the last
points of A and B in the same order, respectively,
for to can be written as follows:
(21)
The DCT-I of
can be calculated using A, B,
and as follows:
(22)
This method is employed recursively to obtain A, B, and thus
. Thus, we will need DCT-II of order N/2 down until 1
to compute
. The block diagram of an N-point DCT-I
transformer is shown in Fig. 1. The restructured block diagram
for computing the DCT-I of an alternate sequence is also shown
in Fig. 1.
can be computed using a similar approach. Let
and .
(23)
Since
can be written as the sum of two symmetric vec-
tors (
and ), can be computed using an N-point
transform as follows:
(24)
where .
From (24), it can be seen that
can be computed from
the DCT-I and DST-I [20] of the two sequences,
and
. Let and represent the DCT-I and
the DCT-II of the even and the odd samples, respectively, of
. The suffix R shall denote reordering. [The reordering
(which is an point sequence) of an N-point se-
quence Y(n) is defined as
for to
.] Then, the 2N-point Hartley transform of is
(25)
The DST-I of a sequence can also be obtained from the
DCT-II of that sequence [9]. The DST-I is computed using the
DST-I and DST-II of the even and the odd samples, respec-
tively. Thus, a DST-I of order N requires DST-II’s of order N/2
down until 1. The DST-II of order L can be computed from the
DCT-II of the alternate sequence as shown below:
(26)
Authorized licensed use limited to: Eindhoven University of Technology. Downloaded on September 11, 2009 at 04:32 from IEEE Xplore. Restrictions apply.