Infoscience - École polytechnique fédérale de Lausanne

Dragotti, Pier Luigi

Vetterli, Martin

Wavelet Transform Footprints: Catching Singularities For Compression and Denoising

It is known that in the absence of distortion, the minimum average sampling density for a multiband signal is given by its spectral occupancy [1]. Furthermore, there exist nonuniform sampling patterns of the same average sampling density such that reconstruction is feasible even if the actual spectral support of the multiband signal is unknown [2]. This is called spectrum-blind nonuniform sampling. However, if the samples are distorted, an increased sampling density may lead to superior reconstruction.LCAVLIN

Gastpar, Michael

Bresler, Yoram

On the necessary density for spectrum-blind nonuniform sampling

Horbelt, Stefan

Dugelay, Jean-Luc

Active contours for lipreading - combining snakes with templates

Jordan, Fred

Ebrahimi, Touradj

VIEW-DEPENDENT TEXTURE CODING FOR TRANSMISSION OF VIRTUAL ENVIRONMENT
S. Horbelt, F. Jordan, T. Ebrahimi,
Swiss Federal Institute of Technology, Switzerland
Abstract: We propose a novel technique for coding of
texture and three-dimensional data which takes the
viewer position into account. This allows to transmit
only the most visible parts of information that is
needed to render a virtual scene. The rendering
operation is modelled and studied using filtering and
sampling theory. The technique is applied on real data
and results are given for an operational encoding
system.
Keywords: view-dependent scalability, anisotropic
filtering, DCT filtering, very low bitrate, synthetic
texture coding, back-channel, computer graphics.
1 INTRODUCTION
The development of digital communication has led to
an astounding amount of research in the field of
compression and transmission of video information.
This research has developed several standards for the
coding of digital video, which are now widely used. In
parallel, digital networks like the Internet provide new
means for exchange of visual information. This kind
of network offers the ability for the user to interact
directly with a remote database, which was not the
case in TV broadcasting conditions, for instance. In
parallel, we see the emergence of virtual reality
applications. It turns out that the distributed and
interactive access to synthetic worlds could open
incredible breakthroughs for the development of
visual communication. In order to allow such
interactive access to remote data, the development of
techniques for optimal coding of synthetic 2D/3D data
become of utmost importance. This will lead to the
definition of a new generation of standards [1, 6].
We propose a novel way to obtain an optimal coding
of texture and 3D data. It consists of taking into
account the position of the viewer in the 3D virtual
world in order to transmit only the most visible
information. This approach allows therefore to reduce
greatly the transmitted information between a remote
database and a user. At a given time, only the most
important information is sent, depending on object
geometry and viewpoint displacement.
For instance, if we consider a single plane on which a
texture image is mapped, one can notice that the
information displayed after rendering uses only a
small amount of the data contained in the original
texture. This depends on the distortions caused by the
rendering process. They are due to perspective
distortion, which depend upon orientation, distance or
field-of-view as well as other criteria such as
illumination, shadows, fog effect, motion blur, etc.
Little research has already been performed on this
topic, nevertheless, several works in computer
graphics have been made on model simplification in
order to decrease the computing load of graphical
engine [2]. Some of the ideas used in this field can be
exploited for the optimal coding of data. In particular
mip-mapping techniques [3] offer a solution to texture
multi-resolution. Furthermore, Peter Lindstrom [4, 5]
has remarked that texture resolution could be
decreased by taking into account the viewer position.
2 RENDERING OPERATION
2.1 INTRODUCTION
The rendering operation is the mathematical
transformation which projects on a 2D plane, the
image of the 3D scene. Two kinds of projections may
be used: orthogonal and perspective projections.
This section will show that perspective projection of
an image, e.g. texture, can be approximated by a
linear transformation.
2.2 PERSPECTIVE PROJECTION
As was shown by [5], the perspective projection of a
line of length L is a line of length L’ which can be
written as:
(1) ( ) ( )αα 222 sin4
4
cos'
⋅−
⋅⋅⋅=
LD
D
vLL
with D the distance between viewpoint and line
center, α the angle with the line and v the distance to
the projection plane.
Viewer
Texture imageRendered image
Dv α
Figu
re 3
L
F
i
L’
Figu
re 1Figure 1
(2) Lk
D
vLLL D ⋅=⋅⋅=>> αcos'  , case,In 
Suppose a one-dimensional texture f(x) defined along
the x-axis with x∈[0, L[. The perspective projection of
this one-dimensional pattern f(x) is then a simple
change of scale of this pattern:
(3) ).()(' xkfxf =
More generally, [7] has shown that in the case of
orthogonal projection, the projection of a two-
dimensional texture f(x, y) is:
(4) ),(),(' dycxbyaxfyxf ++=  
 
 
In the following, we use (4) in the case of perspective
projection, using the condition given in (2).
3 FILTERING OPERATION
3.1 INTRODUCTION
We will show in the following that the projection of a
texture image can be considered as a filtering
operation followed by a down-sampling.
3.2 DEPTH EFFECT
Let us consider the case of a texture f(x,y) which is
mapped on a plane, normal to the viewing direction.
x
y
z
f(x,y)
f’(x,y)
Projection plane
View point
Figure 2: depth effect, k=2
This projection is a particular case of (4) with b=c=0
and a=d=k. k is proportional to D: k∝D.
(5) ),(),(' kykxfyxf =
So, if f(x,y) is defined at L∗L discrete points (square
texture), the number N of points which are rendered
on the screen is:
(6) 
2
2
k
LN =
If we call F(u ,v) the Fourier transform of f, and FT
the function which associates F to f, we then have:
(7) ),(1)),('(FT
2 k
v
k
uF
k
yxf =
Which means as noticed by [7] that the frequencies of
the projected version is increased by a factor k in both
dimensions. As a consequence, a low-passed version
of f will have a projection similar to f’. This
corresponds to a usual down-sampling operation.
Figure 2 shows an example where the distance has
been chosen such that only a quarter of the frequency
information is used, compared to the original texture.
3.3 TILTING EFFECT IN ORTHOGONAL
PROJECTION
Let us consider the case of plane tilted around y-axis
and projected on the x-y-plane.
x
y
z
f(x,y)
f’(x,y)
Projection plane
Figure 3: tilting effect, cos α= 0.5
The tilting of the square texture f(x,y) mapped on this
plane is another particular case of (4) which can be
written in orthogonal projection as:
(8) ),(),(' ykxfyxf =
 
  
with b=c=0, d=0 and a=k with k∝ 1/cosα.
3.4 ROTATION EFFECT IN ORTHOGONAL
PROJECTION
Let us consider again the case of a plane normal to
viewing direction (α=0). We rotate the plane by the
angle β, then its FT is known by [9]:
(9) )cossin,sincos(
))cossin,sincos((FT
ββββ
ββββ
⋅+⋅−⋅+⋅=
⋅+⋅−⋅+⋅
vuvuF
yxyxf
x
y
z
f(x,y)
α
β
Normal
View Point
Figure 4: Orientation of the plane
The combination of this rotation with tilting and depth
enables to describe the effect in the frequency domain
of any orientation of the plane.
3.5 CONCLUSION
More generally, for a given viewpoint, we have shown
that the projection of a texture is equivalent to an
anisotropic filtering of the texture image.
Consequently it is possible to decrease the number of
bits necessary to encode a texture image by
transmitting only a downscaled version of the texture.
Filter
Viewing conditions
DecoderFFT Coder
Texture Texture
Figure 5: View-dependent coding and decoding
4 USE OF DCT FOR FILTERING
4.1 INTRODUCTION
In the following, we show how DCT can be used to
approximate the above described filtering operation.
Let H(u,v) be the DFT of the transfer function of a
filter and h(x) its impulse response. Let note Fc(u,v)
the DCT of f(x,y) and Wc(u,v) the DCT of the filtered
texture
 w(x,y).
It has been shown by [11] that if:
(10) ),(),(),( vuHvuFcvuWc ⋅=
(where u,v are the coordinates in DCT and Fourier
space respectively for Fc and H)
Let ),(ˆ yxf  be the four-fold symmetric extension of
the function f(x,y). Then
(11) ),(),(ˆ),( yxhyxfyxw ∗=
assuming that H is real and even:
(12) ),(),(),(),( vuHvuHvuHvuH −=−=−−=
Thus, we can replace a filtering operation by using the
DCT transform of a texture image f(x,y) as noticed by
[8]. For coding efficiency reasons, we will use only
binary filters H(u,v) in the following. These kinds of
filters applied in DCT domain have already been
studied by [12, 8, 10]. This filtering operation can be
compared to a zonal coding [14, 15].
4.2 DEPTH EFFECT
As shown by [10], this operation can be considered as
a simple change of scale. The change of scale using
DCT transform has already been studied by [8, 10,
12]. It is possible to use a filter ),( vuHck  applied on
Fc(u,v) as in (10) with the shape described in Figure
6. In practice we simply compute [5]:
(13) ),(),(),( vuHcvuFcvuWc k⋅=
with
(14) 
]}1,0[],1,0{[},{and
elsewhere ,1),(
1
or  
1
  if ,0),(
−−∈
=
−
>
−
>=
LLvu
vuHc
k
L
v
k
L
uvuHc
k
k
H(u,v)
u
u
v
u
u
v
),( vuHk ),( vuHck
(L-1)/k
(L-1)/k
(L-1)/k
(L-1)/k
Figure 6: Filter design for depth effect
4.3 TILTING EFFECT IN ORTHOGONAL
PROJECTION
As shown by (8) this effect can be seen as a simple
change of scale in one direction. The corresponding
kH  and kHc  filters are defined in Figure 7.
(15) 
]}1,0[],1,0{[},{with
else. ,1),(
1
  if ,0),(
−−∈
=
−
>=
LLvu
vuHc
k
L
uvuHc
k
k
u
u
v
u
v
),( vuHk ),( vuHck
(L-1)/k
(L-1)/k
Figure 7: Filter design for tilting effect
4.4 ROTATION EFFECT IN ORTHOGONAL
PROJECTION
The directional filtering in DCT domain has already
been addressed by [13]. Taking into account equation
(9), it follows that the optimal filter in Fourier domain
should be as shown in Figure 8.
u
v
),( vuHck
(L-1)/k
β
),( vuHk
Frequency along
horizontal axis in
rendered image
?
Figure 8: Optimal filter for rotation effect
Since this is a complex filter, there is no way to
implement it by one operation in the DCT domain.
We have studied two real filters as alternatives to the
optimal filter. These filters are symmetrical as shown
by (12). These two filters are sub-optimal for different
reasons:
• H1 is not removing some invisible frequencies.
• H2 is removing some visible frequencies.
u
v
),( vuHck
β
),( vuHk
u
u
v
Figure 9: Real filter 1
u
v
),( vuHck
β
),( vuHk
u
u
v
Figure 10: Real filter 2
Let us design the two filters to make a quantitative
comparison between them (Figure 11). The plane is
rotated by β=45° and titled by α=60°, which means
that only half of the pixels of original texture are
rendered on the screen. Thus, we design the filters
such that their non-zero surface S covers half of the
filter plane.
u
v
v
),(2 vucH k),(1 vucH k
v
v
u
u
2L
S=L2/2
S=L2/2
Figure 11: Shapes of compared filters.
For each of the filters we have made the following
experiment:
1. DCT of the a texture f(x,y) with size L=128
2. Application of filter 1 or 2
3. IDCT of the texture
4. Texture mapping and rendering
Results are shown in Figure 12. We have compared
these two images to the one we get with the original
texture.
Filter H1 Filter H2
Figure 12: Results using filter H1 or H2
Technique PSNR (dB) Removed DCT
Filter H1 32.2811 8010
Filter H2 36.2233 8128
Table 1: PSNR using filter H1 or H2
Table 1 gives for each filter both the PSNR and the
number of DCT coefficients that have been set to zero
by the used filter.
The conclusion is that for this particular orientation
filter H2 gives much better results. Although H2 is
filtering frequencies which appear horizontal on the
rendered image, this effect is less visible than the
filtering performed by H1.
4.5 GENERALIZATION IN PERSPECTIVE
PROJECTION
The assumption is used that perspective projection can
be approximated locally with orthogonal projection.
We use the rotation filter 2 combined with the depth
filter to define the shrinking mask [auv] using the depth
D and the rotation angles α  and β.
For all shrinking mask coefficients auv with
 v∈[0, L-1] and u∈[0, L-1],
We define the normalized coordinates:
V=(v+0.5)/L, V∈[0,1]
U=(u+0.5)/L, U∈[0,1]
Let define the threshold value T as:
T = (k/D) ⋅ cos α;
The shrinking angle βS is defined as:
 βS  = | arcos(cos(90-β)) |
Then for each aij:
if (U>k/D or V>k/D) then auv = 0
if       (sin βS  = 0  and V> T ) then auv = 0
else if (sin βS  = 1  and U> T ) then auv = 0
else if (U >T- (V-T)/ tan βS ) then auv = 0
5 TEXTURE DATA REDUCTION PRINCIPLE
5.1 INTRODUCTION
We have been focusing on the coding of texture
information that is used for generating photo-realistic
landscapes [5]. This kind of data typically represents a
large amount of space (several megabytes). However
at a given time, only a small portion of it is really
needed for visualization. The filtering described in the
last chapter can be used for reducing the transmitted
texture data mapped onto any kind of 3D object.
5.2 GENERALIZATION TO A REGULAR GRID MESH
The filtering discussed in the last chapter for one
plane is now generalized for a 3D regular grid mesh.
We use, at the decoder side, the knowledge of the 3D
regular mesh on which the texture is mapped. The
texture is split into small 8×8 sized blocks, which are
aligned to the cells of the 3D regular grid mesh. The
DCT is computed for each texture block. Then each
DCT block is filtered independently with the
shrinking mask. It uses the local orientation and depth
of the grid cell, which belongs to this DCT block.
Only the DCT coefficient with a corresponding none-
zero shrinking mask coefficient auv are considered for
the transmission and coding.
5.3 STREAMING TEXTURE DATA
Once the texture data needed for the first frame (we
call the data intra data) has been loaded using the
proposed method, it is then straightforward to send
new DCT coefficients (we called them inter data) as
the observer is moving in the scene. Therefore, it is
possible to update the texture information according
to the new viewer position using a very low bit rate
channel. (Figure 13)
Figure 13: View-dependent compression scheme.
5.4 COLOR COMPONENTS
The colors are compressed by applying the same
technique on the 1:4 down-sampled versions of the U
and V texture components.
5.5 ENCODING THE BITSTREAM
The DCT coefficient remaining after the shrinking are
quantized using the quantization tables of MPEG-2,
zig-zag-scanned and variable length coded (VLC)
using the VLC tables of MPEG-4, as detailed in [16].
Due to the shrinking operation and quantization the
inter blocks often contain only zero valued DCT AC
coefficients. So for each block a one bit entry in the
coded pattern value (CPV) table denotes whether any
none-zero DCT AC coefficients will be sent.
The DC coefficients are DPCM and VLC encoded. Its
predictor uses the DC values of the previous blocks
along the two texture axis.
6 RESULTS
Tests have been done with a flight over a landscape
defined by a 128×128 sized elevation grid. The
texture was an RGB image of size 1024×1024 pixels,
which represents a size of 24Mbits. Figure 14 shows
the bitrate obtained to generate each of the frames of
the sequence. This bitrate does not include the
transmission of the elevation data, which represent a
raw data of 32768 bytes with 16 bits per elevation
point (1 % of the texture data). The first frame
requires transmitting 358242 bits, then the numbers of
bits vary between 1000 and 20000 bits for each frame.
For the first frame the intra compression ratio is 1:11,
for the following frames the inter compression ratio is
1:460 and the overall compression ratio is 1:300.
  



 

 

  
    

    
      
Frame number
N
um
be
r 
of
 b
its
Figure 14: Number of bits transmitted to update
the texture image, as the viewpoint is moving.
7 CONCLUSION
We have proposed a technique that uses a back-
channel from the decoder to the coder. This takes into
account the viewing conditions to specify the data
required by the decoder. It enables to transmit and
update texture data as the viewpoint is changing.
Results show that it is possible to reach very low bit
rate coding while keeping high quality of the rendered
images. This technique has been proposed as a
standard for coding of texture images in the
framework of MPEG-4 standardization activity.
Finally, the technique can be easily extended to take
into account other perceptive criteria, like shadows,
speed or fog and therefore allow reducing even more
the size of transmitted data.
This technique covers a wide range of applications:
Streaming of flight simulation data at very low bit-
rates, augmented reality for IFR conditions flying or
exploration of planet surfaces. Streaming of data for
virtual walks through buildings, eg. in architecture or
virtual multi-player worlds.
We would like to thank Swissair for providing us with
photo-realistic landscape textures and elevation data.
8 REFERENCES
1. Peter K. Doenges, Frederic Jordan, Igor
S.Pandzic, 1996, "SNHC Application Objectives
and Requirements", MPEG-4 technical report
M1316, Tempere, Finland.
2. Funkhouser, T.A. and Sequin, Proceedings of
SIGGRAPH 93, pp. 247-254.
3. Jackie Neider, Tom Davis, Mason Woo, OpenGL
programming guide.
4. Peter Lindstrom, David Koller, Larry F.Hodges,
William Ribarsky, Nick Faust, Proc. of
SIGGRAPH 96, pp. 109-118.
5. Peter Lindstrom, David Koller, Larry F.Hodges,
William Ribarsky, Nick Faust, “Level-of-detail
management for real-time rendering of photo-
textured terrain”, Internal report of Georgia Tech,
http://www.cc.gatech.edu/gvu/research/techreport
s95.html
6. Fred Jordan, Touradj Ebrahimi, 1996, “View-
dependent texture for transmission of virtual
environment”, ISO&IEC JTC1, M1501.
7. Don Pearson, 1990, Signal Processing: Image
communication, vol. 2, no. 4, pp. 377-395.
8. King N.Ngan, 1986, Signal Processing, vol. 11,
pp. 249-263.
9. Kenneth R. Castleman, “Digital Image
Processing”, Prentice-Hall.
10. Stephen A. Martucci, 1996, Proceedings of
ICIP’95, vol. 2, pp. 244-247.
11. Chitprasert B., Rao K.R, 1990, Signal Processing,
vol. 19, no. 3, pp. 233-245.
12. Adant J.M., and al., 1987, Signal Processing, vol.
13, pp. 285-397.
13. Toshio Tohne, Hiromi Koda, Hatsukazu Tanaka,
1996, Systems and Computers in Japan, vol. 27,
no. 3, pp. 94-105.
14. Gonzalez Rafael, Woods Richard, “Digital Image
Processing”, Addison –Wesley.
15. Stefan Horbelt, Jon Crowcroft, 1996, “A hybrid
BTC/ ADCT video codex simulation bench", 7th
Int. Workshop on Packet Video, Brisbane,
Australia.
16. Stefan Horbelt, Fred Jordan, Touradj Ebrahimi,
1997, “View Dependent Texture Coding for
Transmission of Virtual Environment”,
ISO/IECJTC1/SC29/WG11, MPEG96/M1723.


View-dependent texture coding for virtual environments

Cvetkovic, Zoran

Overcomplete expansions and robustness

Goodwin, Michael

Time-frequency signal models for music analysis, transformation and synthesis

Khansari, Masoud

McCanne, Steven

The interaction of source and channel coding in wireless and packet video

Kovacevic, Jelena

Design of multidimensional nonseparable regular filter banks and wavelets

Anastassiou, Dimitris

A multiresolution approach for all-digital HDTV

Duhamel, Pierre

CYCLIC CONVOLUTION OF REAL SEQUENCES: HARTLEY VERSUS FOURIER AND NEW SCHEMES
* **
P. DUHAMEL and M. VETTERLI
* CNET/PAB/RPE 38-40, rue du Général Leclerc, F-92131 Issy-les-Moulineaux (France)
** Ecole Polytechnique Fédérale de Lausanne 16, Chemin de Bellerive, CH-1007 Lausanne (Suisse)
ABSTRACT
Recently, new fast transforms (such as the discrete
Hartley transform in particular) have been proposed
which are best suited for the computation of cyclic
convolution of real sequences. Two approaches using
Fourier or Hartley transforms are first compared,
showing that the recently proposed FFT algorithmsfor real data present a lower arithmetic complexity
than the corresponding DHT-based approach. Impro-
vements are made to both types of algorithms, leading
to different trade offs between arithmetic and struc-
tural complexity. We also present a new Hartley
Transform algorithm with lower arithmetic complexity
than any previously published one.
I. INTRODUCTION
The radix-2 Fast Fourier Transform algorithm was
first explained by Cooley and Tukey in 1965, in a
version on complex data [1]. But, as most of the
data to be treated is real, instead of complex, the need
of a version of the FFT on real data removing extra
calculations from the complex case became soon
apparent, and such a program was quickly published
by Bergland [3].
In fact this program does not seem to have been
widely used, and nany people preferred either to
compute a length 2 - real DFT either by the use
of a length 2n-l complex FFT plus some additional
operations, or to coRlpute two real FFT's at a time
by using one length 2 complex FFT.
More recently,the 2n FFT algorithms with the minimum
known number of both multiplications and additions
were also proposed in a real-data version [6,7,10,11].
In the same time, Bracewell [8] proposed the use
of fast Discrete Hartley Transforms (FDHT), which
is real of nature, as a substitute to the FFT for
cyclic convolutions on real data.
This paper is concerned withn the computation of
cyclic convolutions of length -2 real data through the
use of fast transforms.
The well-known approach of computing real convolu-
tions through FFT algorithms on complex data will not
be considered here, and we shall only consider pro-
grams specially designed for real data, since they
allow to obtain the lowest known number of arithmetic
operations, without increase in programlength.
In the first paragraph, we shall first briefly present
the two basic approaches considered here : convolution
using FFT's, and convolution using FDHTs. This
comparison shows that the former has some advantages
when considering arithmetic complexity,while the latter
has a simpler structure, since the DHT is self-inverse.
The second paragraph will present a scheme for
cyclic convolution using two forward DFT's on real
data, and a small number of additional operations.
This solution has the advantage of a simple structure
since it uses only forward DFT's on real data.
The third paragraph will present a hybrid FF/FDHT
algorithm, allowing to obtain any length -2 FDHT
using only 2 additions more than the best FFT algori-
thms known for real data. In this paragraph, we
also show how to derive very quickly improved FDHT
algorithms from their FFT counterparts.
II. THE INITIAL SCHEMES
11.1. Convolution using FFT's
Since the DFT has the convolution property, the
convolution scheme has the structure given in fig. 1.
When the initial sequences x and h are real,
IXhI , Hh , and Vh} have hermitiai symmetry
= Y*k, and the forward DFT has to be performed
on real data, while the inverse FFT has to be perfor-
med on data with hermitian symmetry. Both kind of
FFTs have been published [7,11] have the same
arithmetic complexity, and can be performed in place.(In fact, an FFT on data with hermitian symmetry can
be derived by reversing the flow-graph of an FFT on
real data).
If the FFT's are performed with the lowest arithmetic
complexity, using FFCT [6] or split-radix [7] algorithms,
this leads to the following number of operations
cony. FFT 2 Mfr + 2 + 3 . (2n1 - 1)l) M =n n
= 2n-l(203 +3
(2) A cony. FFT = 2 Afr + 3 (2 - 1)n n
= 2 (6n - 7) + 5
In these operation counts, it has been taken into
account that, due to the symmetry of both1Xkand
lHki, 1''k} can be obtained with 2 real multiplications
for k=O, and k=N/2, and N12-l complex multiplications(3 real mults ± 3 real adds each) for k=1,---N/2-l.
6. 5. 1
ICASSP 86, TOKYO CH2243-4/86/0000-0229 $L00 © 1986 IEEE 229
This arithmetic complexity is rather low, as shown
in table 1, since it is based on powerful FFT algorithm
on real data, but a disadvantage of this scheme is
that it needs both forward transform on real data and
inverse transform on complex data with hermitian
symmetry. This results in an increase in program
length for the cyclic convolution.
11.2. Convolution using FDHT's
The Discrete Hartley Transform is defined as the
sum of the real and imaginary part of the Fourier
Transform
N-Ih -. 27cnk . 2itnk(3) Xk = L x(cos —fl--— + sin
n=0
As can be seen from the definition, the DHT is
its own inverse, and is a real transform. Furthermore,
the convolution in the time domain corresponds to the
following operation in the Hartley domain
(4) Y = X . H(k) + Xhk H(k)
where Heh(k) (resp. H01(k)) is the even (resp. odd) part
of the Hartley transform of h
n
So, the convolution scheme using Hartley transforms
becomes as shown in fig. (2).
When considering eq.(4), the 'multiplications' in
the Hartley domain seems to need 2 multiplications per
point, but it is easy to see that, by grouping the
computation of Yhk and , this number can be
reduced to 3/2 multiplications per point
H(k) H0(k) X(5)
-H0(k) He(k) Xhk
and eq.(5) is easily recognized as a complex multiplica-
tion, needing only 3 real multiplications and 3 additions.
If we consider also that
= H(o)0 0 e
Y/2 = X12 H(N/2)
One can see that the set of "multiplications" in
the Hartley domain needs exactly the same number of
operations as in the Fourier domain.
The remaining part of the scheme given in fig. (2) to
be evaluated consists in the FDHT.
FDHT was first introduced by Bracewell [8] in a
radix 2, decimation in time version, and Sorensen,
3ones, Burrus and Heideman [9] recently proposed
different types of FDHT algorithms.
Inspection of these algorithms allows the following
conclusions
A radix 2 FDI-JT algorithm needs N-2 additions more
than the corresponding FFT algorithm on real data, and
the same number of multiplications.
The split-radix algorithms, which seems already
to be one of the best compromises for real-data FFT's
(lowest known number of operations, in-place computa-
tion, compact program) gives also the lowest known
number of operations for the FDHT, but requires
2( ) additions more than the corresponding
split-radix FFT algorithm on real data.
If the FDHT in fig.(2) is computed through a split-radix
algorithm, this gives the following arithmetic comple-
xity for the cyclic convolution
(7) M°'" DHT = 21(2n 3) + 3
(8) A°" DHT = 2'(6n-7) + 5 + 4 2fl_l..(J)fl_I
One can see that convolution using FDHT's will
have a simpler structure, even if the FDHT programs
have the same code length as FFT programs, due to
the self-inverse property of the Hartley transform.
This simplicity will be obtained at the cost of a small
increase of the number of additions.
In the following, we shall improve both solutions
III. A MORE REGULAR CONVOLUTION SCHEME
USING FFT's
The improvement allowing to use two forward FFT's in
the convolution scheme is based on some of the
remarks made when studying FDHT's
The convolution property in the Hartley domain,
as given in eq.(4) can be rewritten as follows, by
reversing the roles of Xk and Hk
(9) Y = x(k) . H + x(k) Hhk
But, Xeh(k) being the even part of the Hartley trans-
form of x } is also the real part of the DFT of
and x0k is the corresponding imaginary part.
We can then rewrite eq. (5) as follows
(10) Y = H< ReX
Hhk -H ImX
which is also equivalent to a complex product.
Furthermore, since a DHT is the sum of the real and
imaginary part of a DFT, the last FDHT in fig.(2) can
be replaced by an FFT followed by the sum of the real
and imaginary part of each DFT component, thus resul-
ting in the diagram of fig. (3).
This scheme requires exactly the same number of
operations as the usual convolution scheme based
on FFT's, plus N-2 additions due to the +1- box
(11) A°'" fft = 2'1(6n 5) + 3
230 ICASS? 86, TOKYO
= 2(6n-
-) + 5 - . (1)fl_l
(6)
6. 5. 2
This structure is now as regularas the one using
FDE-IT's, to the price of about 2 /3 extra additions
compared to the FDHT-scheme. It has the advantage
of using only usual FF1 algorithms on real data, that
should become widely used, due to their simplicity
and compactness.
IV. IMPROVED FDI-IT ALGORITHM
FDHT algorithms, as found in [8], and [9] can be
improved by making full use of the close relationship
between Fourier and Hartley transforms In fact,
nearly all the extra additions needed to compute
FDHT's instead of FFT's can be nested inside the
twiddle factors of the FF1, thus resulting in hybrid
FDHT/FFT algorithm.
(12) X< = x(cos a,f1! - i sin = xn W
By multiplying the DFT of lxn by (1+j), we obtain
eq. (13)
(13) (1+j) X = (Re 4 -Im 4) +j (Re 4+ Tm 4)
which means that multiplication of by (l+j)gives us the Hartley transform of X
Let us apply this remark to the basic split-radix
decimation-in-time decomposition of the DFT [10]
as given in (14)
N12-l N/4-1(14) 4< = xw< + w1! rO X4fl•1 WN +
N/4-1
+ (l+j) w< : x4÷3
n=O
eq. (15) is now the basic recursion of the improved
FDL-IT algorithm : it transforms the length-N DHT
into a length N12. DHT plus two length N/4 - DFT's,
and some multiplications by twiddle factors of the
k 3kform (l+j) w N or (1+3) WN
As in the case of FFT'skon real data, the butterflies
with twiddle factors WN have to be computed toge-
-kther with the ones with twiddle factors W . This
needs (N/8 - 1) general 6-mult-l 8 adds btterf lies,
plus one 0-mult-6-adds butterfly (corresponding to
k=0) and one 2-mult-4-adds butterfly (corresponding
to k N/8).
Iteration of these addition counts gives exactly the
same number than was obtained for FFT's on real
data, except for the FDHT of length 4 which needs
2 extra addition.
This leads to
(16) A hart = 2n-l (32n - 5) + 6
n
and the whole convolution with this algorithm will need
(17) Aconv hart. = 2n-I (6n +
The number of multiplications is the same as for
the other schemes.
The convolution scheme using this new fast Hartley
transform algorithm has now the same number of
multiplications, and 4 additions more than the best
known scheme, with the advantage of using only
a single DHT program, to be used twice.
It must be noticed that this has to be paid by an
increase in the program length. The new FDHT algo-
rithm has the same structure as the corresponding
FFT algorithm, some of the butterflies in the FFT
program being replaced by the ones corresponding
to eq.(15). Nevertheless, both programs will be the
same size when using autogen techniques.
Another remark is that the same type of improvement,
if applied to other classical decompositions (eg. radix 2,
4, etc...), improves all existing programs.
V. CONCLUSION
In this paper, we have first presented two usual
schemes for cyclic convolution of real signals, via
FFT or FDHT. Operation counts are given for both
of them using the fastest known algorithms, showing
that the FFT-based convolution has a lower arithmetic
complexity, but a more complex structure.
Improvements are then made to both schemes, since
we propose
- a FF1-based scheme using two forward transform
on real data.
- a FDHT-based algorithm with reduced number
of additions.
All operation counts are summarized in Table I
allowing to choose the best tradeoff between low
arithmetic complexity and structural complexity.
REFERENCES
[1] J.W. COOLEY, J.W. TIJKEY "An algorithm for
machine computation of complex Fourier series".
Mat Comput. 1965, vol 19, pp. 297-301.
[2] R.D. PREUSS : "Very fast computation of
the radix 2 Discrete Fourier Transform".
IEEE Trans. on ASSP, 1982, vol. 30, pp. 595-607.
[3] G.D. BERGLAND "A Fast Fourier Transform
Algorithm for real-valued series". Commun.
ACM, vol. 11, pp. 703-710, Oct. 1968.
[4] H. ZIEGLER. "A fast Fourier Transform
rithm for symmetric real-valued series".
trans. on Audio and Electroacoustics,
AU-20, n°5, pp. 353-356, Dec. 1972.
ICASSP 86, TOKYO 231
= X + j Xhk
N/4- 1
4nk
= WN
x4n+3 WN
n=O
N/2-1
(15) (l+j). 4< = x + j Xh1< = (1+j) WN
n=O
1.1/4-1 4nk+ (l+j) Wr x4n+l WN +
n=O
6. 5. 3
Algo-
IEEE
vol.
{5] J.B. MARTENS : "Discrete Fourier Transform
Algorithms for real-valued sequences". IEEE
Trans. on ASSP, vol. ASSP-32, n°2, pp. 390-396,
April 1984.
[6] M. VETTERLI, H.J. NUSSBAUHER : SimpleFFT and DCT algorithms with reduced number
of operations". Signal Processing, vol. 6, n°4,
pp. 267-278, July 1984.
[7] P. DUHAMEL, H. HOLLMANN "Implementation
of 'split-radix" FFT algorithms for complex, real,
and real-symmetric data". ICASSP'85, Tampa.
[8] R.N. BRACEWELL "The Fast Hartley Transform"
Proc IEEE. Vol 22, n°8, Aug. 1984, pp. 1010-1018.
[9] H.V. SORENSEN, D.L. JONES, C.S. BURRUS,
M.T. HEIDEMAN. "On computing the Discrete
Hartley Transform". IEEE Trans. ASSP, vol.
ASSP-33, n°4, pp. 1231-1238, Oct. 1985.
[10] P. DUHAMEL "Un algorithme de transformation
de Fourier rapide a double base". Ann. Tlkom.
Vol. 40, n°9-lO, Sept-Oct. 1985, pp. 481-494.
[11] M. VETTERLI, H.J. NUSSBAUMER "Algorithmes
de transformation de Fourier et en cosinus
mono et bi-dimensionnels". Ann. Télacom.
Vol. 40, n°9-l0, pp. 466-476, Sept—Oct 1985.
{hn }
* Jhkn
{xk}
{Hkh}
Fig. 3 Modified FFT-based convolution scheme
using only forward FFT's on real data.
Number of additions
{xk}
Fig. I : Usual convolution scheme using FFT's.
{Hkh}
Fig. 2 Convolution scheme based on FDHT's.
N
number
of mult.
initial
FFT — based
initial
FDHT — based
FFT — based
with
improved
FD1-IT— based
2 forward FFT's
16
32
64
128
256
512
1024
2048
15
43
115
291
707
1667
3843
8707
19459
49
141
373
933
2245
5253
12037
27141
60421
53
153
393
977
2329
5425
12377
27825
61785
55
1 55
403
995
2371
5507
12547
28163
62467
53
145
377
937
2249
5257
12041
27145
60425
Table 1 : Number of arithmetic operations
for a cyclic convolution on real data.
6. 5. 4
232 ICASSP 86, TOKYO


Cyclic convolution of real sequences: Hartley versus Fourier and new Schemes

École Polytechnique Fédérale de Lausanne

Not a member yet