A Butterfly Structured Design of The Hybrid Transform
Coding Scheme
Jingning Han, Yaowu Xu, and Debargha Mukherjee
Google Inc., 1950 Charleston Road, Mountain View, CA 94043
Emails: {jingning,yaowu,debargha}@google.com
Abstract—The hybrid transform coding scheme that alternates
amongst the asymmetric discrete sine transform (ADST) and the discrete
cosine transform (DCT) depending on the boundary prediction condi-
tions, is an efficient tool for video and image compression. It optimally
exploits the statistical characteristics of prediction residual, thereby
achieving significant coding performance gains over the conventional
DCT-based approach. A practical concern lies in the intrinsic conflict
between transform kernels of ADST and DCT, which prevents a butterfly
structured implementation for parallel computing. Hence the hybrid
transform coding scheme has to rely on matrix multiplication, which
presents a speed-up barrier due to under-utilization of the hardware,
especially for larger block sizes. In this work, we devise a novel ADST-
like transform whose kernel is consistent with that of DCT, thereby
enabling butterfly structured computation flow, while largely retaining
the performance advantages of hybrid transform coding scheme in terms
of compression efficiency. A prototype implementation of the proposed
butterfly structured hybrid transform coding scheme is available in the
VP9 codec repository.
I. INTRODUCTION
Transform coding is a central component in video and image
compression. Many research efforts have been devoted to optimize the
transform kernel to fully exploit signal correlation for compression
gains. A recent approach that jointly optimized spatial prediction
and the choice of the subsequent transform for video and image
compression was developed in [1] and [2], where it was shown that
the optimal Karhunen-Loeve transform (KLT) given available, partial
boundary information is well approximated by a close relative of the
discrete sine transform (DST), with basis vectors that tend to vanish at
the known boundary and maximize energy at the unknown boundary.
The overall intra coding scheme thus switches between this variant
of DST named asymmetric DST (ADST), and the conventional DCT,
depending on the prediction direction and boundary information. This
adaptive prediction-transform approach, namely hybrid transform
coding scheme, was experimentally shown to significantly outperform
the DCT-based intra-frame prediction-transform coding.
On the hardware design side, the transform module typically
contributes a large portion of codec computational complexity, and
hence a butterfly structured implementation that allows parallel com-
puting via single instruction multiple data (SIMD) operations [3] is
highly desirable. The design of butterfly structured discrete cosine
transform can trace back to 1970’s (e.g., [4].) A recent development
on fast transform using integer transform was proposed in [5], where
it approximates the DCT transform element-wisely using a matrix
whose entries are all small integers. The computation hence only
involves additions and shifts. Similar principle was also applied to the
ADST in [2], both targetting at 4×4 block size. In fact, methods along
this line are typically limited to smaller transform dimensions. When
it comes to transform dimension of 8 × 8 or above, it is in general
difficult to find an orthogonal matrix whose elements are small
integers and approximates the DCT closely. Hence, its advantages in
simple computation will gradually disappear as transform size grows.
It is noteworthy that larger block size transforms provides higher
transform coding gains for stationary signal and are experimentally
proved to contribute compression efficiency in various video codecs.
Challenges arise in the design of fast ADST, and hence hybrid
transform coding scheme, of any block sizes. The original ADST
kernel was derived as sin(
n (2k−1)π
2N+1
) in [1], where N is the
block dimension, n and k denote the time and frequency indexes,
respectively, both ranging from 1 to N. The DCT kernel, on the
other hand, is of form cos(
(2n−1)(k−1)π
2N
). The butterfly structured
implementations of sinusoidal transforms exist if and only if the
denominator of the kernel argument, i.e., (2N +1) for the ADST and
2N for DCT, is a composite number (and ideally can be decomposed
into product of small integers). For this reason, most block-based
video (and image) codecs are designed to make the block size power
of two, e.g., N = 4, 8, 16, etc, for efficient computation of DCT
transformation. It, however, makes the original ADST not capable of
fast implementation. For example, when N = 8, (2N + 1) turns out
to be 17, which is a prime number that precludes the possibility of
butterfly structure.
This work resolves this intrinsic conflict between DCT and the
original ADST by designing a new variant of ADST whose kernel
is of the form sin(
(2n−1)(2k−1)π
4N
). Clearly the denominator of the
kernel argument, 4N , is consistent with that of DCT, in that if 2N
is a power of two, so is 4N. Therefore, it can be implemented
in a butterfly structure, and is henceforth referred to as btf-ADST.
A prototype butterfly structured implementation of btf-ADST was
initially provided in [6]. The btf-ADST is a basis-wise approximation
to the original ADST, and well preserves the its superior coding
performance given the boundary condition. We hence use this btf-
ADST to replace the original ADST in the hybrid transform coding
scheme. The overall scheme thus selects the appropriate 1-D trans-
forms amongst the btf-ADST and DCT depending on the prediction
direction to form the 2-D transformation. Further note that while there
are many variants of butterfly design for the btf-ADST and DCT,
most of them targeted only on reducing the number of multiplica-
tions without consideration on the rounding effects of intermediate
steps
1
, which potentially affects the accuracy of the computation
hence incurring round-trip error. In our implementation, we use the
structure that has more rotation (which involves more multiplication)
steps in the first few stages, so as to reduce the accuracy impact
due to rounding error. It is experimentally demonstrated that our
proposed approach, in conjunction with SIMD, significantly reduces
the runtime of hybrid transform coding scheme in terms of CPU
cycles, while largely retaining its compression gains.
II. SPATIAL PREDICTION AND TRANSFORM CODING
We revisit the mathematical theory that derived the original ADST,
in the context 1-D first-order Gauss-Markov model, given partial
1
In practice, all the computations are performed in the integer format for
speed reasons.
prediction boundary [1], which leads to our btf-ADST proposed in
this work.
Consider a zero-mean, unit variance, first-order Gauss-Markov
sequence
x
k
= ρx
k−1
+ e
k
, (1)
where ρ is the correlation coefficient, and e
k
is a white Gaussian noise
process with variance 1 − ρ
2
. Let x = [x
1
, x
2
, ··· , x
N
]
T
denote
the random vector to be encoded given x
0
as the available (one-
sided) boundary. The superscript T denotes matrix transposition. The
recursion (1) translates into the following set of equations:
x
1
= ρx
0
+ e
1
x
2
− ρx
1
= e
2
.
.
.
x
N
− ρx
(N−1)
= e
N
, (2)
or in compact notation:
Qx = b + e (3)
where
Q =
1 0 0 0 . . .
−ρ 1 0 0 . . .
0 −ρ 1 0 . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 . . . 0 −ρ 1
, (4)
and b
= [ρx
0
, 0, ··· , 0]
T
and e = [e
1
, e
2
, ··· , e
N
]
T
capture the
boundary information and innovation process, respectively. It can be
shown that Q is invertible, and thus:
x = Q
−1
b + Q
−1
e , (5)
where the superscript −1 indicates matrix inversion. As expected, the
“boundary response” or prediction, Q
−1
b, in (5) satisfies
Q
−1
b = [ρx
0
, ρ
2
x
0
, ··· , ρ
N
x
0
]
T
. (6)
The prediction residual
y = Q
−1
e (7)
is to be compressed and transmitted, which motivates the derivation
of its KLT. The autocorrelation matrix of y is given by:
R
yy
= E{yy
T
} = Q
−1
E{ee
T
}(Q
T
)
−1
= (1 −ρ
2
)Q
−1
(Q
T
)
−1
.
(8)
Thus, the KLT for y is a unitary matrix that diagonalizes
Q
−1
(Q
T
)
−1
, and hence also the more convenient:
P
1
= Q
T
Q =
1 + ρ
2
−ρ 0 0 . . .
−ρ 1 + ρ
2
−ρ 0 . . .
0 −ρ 1 + ρ
2
−ρ . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 . . . −ρ 1 + ρ
2
−ρ
0 . . . 0 −ρ 1
.
(9)
Although P
1
is Toeplitz, note that the element at the bottom right
corner is different from all the other elements on the principal
diagonal, i.e., it is not 1+ρ
2
. This irregularity complicates an analytic
derivation of the eigenvalues and eigenvectors of P
1
. As a subterfuge,
we approximate P
1
with
ˆ
P
1
=
1 + ρ
2
−ρ 0 0 . . .
−ρ 1 + ρ
2
−ρ 0 . . .
0 −ρ 1 + ρ
2
−ρ . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 . . . −ρ 1 + ρ
2
−ρ
0 . . . 0 −ρ 1 + ρ
2
− ρ
(10)
which is obtained by replacing the bottom-right corner element with
1 + ρ
2
− ρ. The approximation clearly holds for ρ → 1, which is
indeed a common approximation that describes the spatial correlation
of video/image signals. The unitary matrix T
S
that diagonalizes
ˆ
P
1
,
and hence an approximation for the required KLT of y, can be shown
as the following relative of the common DST:
[T
S
]
j,i
=
2
√
2N + 1
sin
(2j − 1)iπ
2N + 1
(11)
where j, i ∈ {1, 2, ··· , N} are the frequency and time indexes of
the transform kernel respectively. Needless to say, the constant matrix
T
S
is independent of the statistics of the innovation e
k
, and can be
used as an approximation for KLT when boundary information x
0
is
available.
III. BUTTERFLY STRUCTURED VARIANT OF ADST
A key observation of the above derived ADST is that the rows
of T
S
(i.e., basis functions of the transform) possess smaller values
in the beginning (closer to the known boundary), and larger values
towards the other end. For instance, consider the row with j = 1
(i.e., the basis function with the lowest frequency). In the case where
N 1, the first sample (i = 1) is
2
√
2N+1
sin
π
2N+1
≈ 0, while the
last sample (i = N ) takes the maximum value
2
√
2N+1
sin
Nπ
2N+1
≈
2
√
2N+1
. This effectively exploits the fact that pixels closer to the
known boundary are better predicted and hence have statistically
smaller variance than those at far end. It inspires our search for a
unitary sinusoidal transform that resembles the compression perfor-
mance of the ADST, to overcome the intricacy of butterfly design
of ADST and hence hybrid transform coding for parallel computing.
We hence devise a new variant of DST as an alternative:
[T
btf
]
j,i
=
r
2
N
sin
(2j − 1)(2i − 1)π
4N
!
, (12)
where j, i ∈ {1, 2, ··· , N } denote the frequency and time indexes
respectively. Clearly, it also possesses the property of asymmetric
basis function, but has the denominator of kernel argument, 4N,
consistent with that of DCT, thereby allowing the butterfly structured
implementation. We refer to it as btf-ADST henceforth.
A butterfly structure of btf-ADST was initially provided in [6],
where it assumed no precision loss in the intermediate steps that
involve multiplications by irrational numbers. In practice, all these
computations are performed in the integer format, which inevitably
incurs rounding effects accumulated through every stage. To minimize
the round-trip error, we modify the structure to make the initial stages
consist more multiplications, so that the rounding errors are less
magnified. In keeping the conventions used in [6], let I and D
N
denote the re-ordering operations:
I =
1
···
···
1
D
N
=
1
−1
1
···
−1
. (13)
Let P
J
be the permutation matrix that move the first half of the vector
entries to the even-numbered position, and the second half entries to
the odd-numbered position but in a reversed order:
P
J
=
1 0 ···
0 ··· 0 1
0 1 ···
0 ··· 1 0
··· ··· ···
, (14)
where J is the height of the matrix. It formulates a second permuta-
tion:
H
N
= P
N
P
N/2
P
N/2
···
P
4
I
4
P
4
I
4
···
P
4
I
4
P
4
I
4
. (15)
Similarly the permutation operator Q
J
moves the odd-numbered
entries to be in reversed order:
Q
J
=
1 0 ···
0 ··· 0 1
0 0 1 ···
0 ··· 1 0
··· ··· ···
0 1 0 ··· 0
. (16)
Let J = log
2
N, we define the following building blocks that
formulate the butterfly structure.
Type 1 Translational Operators: matrices U
N
(j), j =
1, 2, ··· , J − 1 are defined as
U
N
(j) =
B(j)
B(j)
···
B(j)
, (17)
where
B(j) =
I
2
j
I
2
j
I
2
j
−I
2
j
. (18)
Type 2 Rotational Operators: V
N
(j), j = 1, 2, ··· , J − 1, are
block diagonal matrices:
V
N
(j) =
I
2
j
E(j)
I
2
j
···
E(j)
, (19)
where E(j) = diag{T
1/2
j+1
, T
5/2
j+1
, ··· , T
(2
j+1
−3)/2
j+1
} and
T
r
=
cos rπ sin rπ
sin rπ −cos rπ
. (20)
As a special case,
V
N
(J) =
T
1/4N
T
5/4N
···
T
(2N−3)/4N
. (21)
Given the above established building blocks, the btf-ADST can be
decomposed as:
T
btf
= D
N
· H
T
N
· V
N
(1) · U
N
(1) · V
N
(2)
···U
N
(J − 1) · V
N
(J) · Q
N
· I
N
, (22)
which directly translates into a butterfly graph. The theoretical coding
performance of the two ADST variants and DCT compared against
the KLT will be provided next, followed by the experimental codec
evaluation.
IV. QUANTITATIVE ANALYSIS
We quantitatively evaluate the performance of the btf-ADST,
original ADST, and DCT, against the KLT (of y in Sec. II) in terms
of coding gains [7] under the assumed signal model, at different
correlation coefficient values.
Let the prediction residual, y, be transformed to
z = Ay = [z
1
, z
2
, ..., z
N
]
T
(23)
with an N ×N unitary matrix A. The objective of the encoder is to
distribute a fixed number of bits to the different elements of z such
that the average distortion is minimized. This bit-allocation problem
is addressed by water filling algorithm of [7]. Under assumptions such
as a Gaussian source, high-quantizer resolution, negligible quantizer
overload, and with non-integer bit-allocation allowed, it can be shown
that the minimum distortion (mean squared error) obtainable is
proportional to the geometric mean of the transform domain sample
variances σ
2
z
i
, i.e.,
D
A
∝ (
N
Y
i=1
σ
2
z
i
)
1/N
, (24)
where for Gaussian source the proportionality coefficient is indepen-
dent of the transform A. These variances can be obtained as the
diagonal elements of the autocorrelation matrix of z:
R
zz
= E[zz
T
] = AE[yy
T
]A
T
= (1 −ρ
2
)AP
−1
1
A
T
(25)
where we have used (8) and (9). The coding gain in dB of any
transform A is now defined as:
G
A
= 10log
10
(D
I
/D
A
) . (26)
Here I is the N ×N identity matrix, and hence D
I
is the distortion
resulting from direct quantization of the untransformed vector y.
The coding gain, G
A
thus provides a comparison of the average
distortion incurred with and without the transformation A. Note that
for any given A (including the btf-ADST, ADST, DCT, and KLT
of y), computing R
zz
, and hence σ
2
z
i
, does not require making any
approximations for P
1
. When A is the KLT of y, R
zz
is a diagonal
matrix (with diagonal elements equal to the eigen values of P
−1
1
),
the transform coefficients z
i
are uncorrelated, and the coding gain
reaches its maximum.
Fig. 1 compares btf-ADST, ADST and DCT in terms of their
coding gains, relative to KLT, specifically it depicts G
T
btf
− G
KLT
,
G
T
s
− G
KLT
, and G
T
c
− G
KLT
, versus the correlation coefficient
ρ. Clearly the original ADST well approximates KLT at various
values of the correlation coefficient ρ. The maximum gap between
ADST and KLT, or the maximum loss of optimality, is less than 0.05
dB. The proposed btf-ADST closely resembles the performance of
ADST and KLT, at a maximum loss of 0.15 dB. In comparison, DCT
performs poorly (by about 0.55 dB loss) for the practically relevant
case of high correlation (ρ ≈ 0.95). At low correlation (ρ → 0),
the autocorrelation matrix of the prediction residual, R
yy
≈ I, and
hence any unitary matrix, including ADST and DCT, will function
as a KLT. The block length used in obtaining the results of Fig. 1 is
N = 8.
0 0.2 0.4 0.6 0.8 1
−0.6
−0.5
−0.4
−0.3
−0.2
−0.1
0
0.1
inter pixel correlation coefficient
coding gains (dB)
KLT
btf−ADST
original ADST
DCT
Fig. 1. Theoretical coding gains of btf-ADST, ADST, and DCT, relative
to KLT, all applied to the prediction residuals, plotted versus the inter pixel
correlation coefficient; the block dimension is 8 × 8.
V. EXPERIMENTAL RESULTS
The proposed btf-ADST was employed to replace the original
ADST in the hybrid transform coding scheme. This btf-ADST/DCT
hybrid transform coding scheme was implemented in the VP9 codec
[8]. We first evaluated its compression efficiency as compared to
the original ADST/DCT hybrid transform and the conventional 2D-
DCT coding scheme. Fig. 2 demonstrates the rate-distortion per-
formance comparison for sequence harbour at CIF resolution.
Clearly, btf-ADST/DCT closely tracks the original ADST/DCT, while
both significantly outperform the conventional 2D-DCT approach.
Similar results were observed over a wide varieties of sequences and
resolutions.
We next consider the computational efficiency of the proposed
btf-ADST/DCT hybrid transform. The btf-ADST enables the but-
terfly structured implementation, which in conjunction with SIMD
operations provides significant speed-up as compared to the original
ADST operating via matrix multiplication. We compare the runtime
of the btf-ADST/DCT and the original ADST/DCT hybrid transform
schemes, in terms of the average CPU cycles, as shown in Fig. 3.
The implementation was using streaming SIMD extension 2 (SSE2)
and the experiments were running on a 64-bit platform. Evidently the
proposed btf-ADST allows efficient hardware utilization and thereby
substantial codec speed-up, while closely resembling the compression
gains of the original ADST.
VI. CONCLUSIONS
This work devised a novel variant of ADST transform whose kernel
approximates the original ADST basis-wisely and is consistent with
the DCT kernel, thereby enabling the butterfly structured implemen-
tation. The proposed scheme allows efficient hardware utilization for
significant codec speed-up, while largely retaining the advantageous
compression performance of hybrid transform coding scheme.
450 500 550 600 650 700 750 800 850 900
36
37
38
39
40
41
42
43
bit−rate (kbit/s)
PSNR (dB)
2D−DCT
ADST/DCT
btf−ADST/DCT
Fig. 2. Coding performance comparison of 2D-DCT, ADST/DCT hybrid
transform, and btf-ADST/DCT hybrid transform for sequence harbour at
CI F resolution.
4 6 8 10 12 14 16
0
2000
4000
6000
8000
10000
12000
14000
transform dimension
CPU cycles
ADST/DCT
btf−ADST/DCT
Fig. 3. Computational complexity in terms of CPU cycles. The btf-ADST
was implemented using streaming SIMD extension 2 and the experiments
were running on a 64-bit platform.
REFERENCES
[1] J. Han, A. Saxena, and K. Rose, “Towards jointly optimal spatial
prediction and adaptive transform in video/image coding,” IEEE Proc.
ICASSP, pp. 726–729, Mar. 2010.
[2] J. Han, A. Saxena, V. Melkote, and K. Rose, “Jointly optimized spatial
prediction and block transform for video and image coding,” IEEE Trans.
on Image Processing, vol. 21, pp. 1874–1884, April 2012.
[3] J. L. Hennessy and D. A. Patterson, Computer Architecture-a quantitative
approach, 4th Ed., Mogan Kaufmann, 2007.
[4] W.-H. Chen, C. Smith, and S. Fraclick, “A fast computational algorithm
for the discrete cosine transform,” IEEE Trans. on Communications, vol.
25, pp. 1004–1009, Sep. 1977.
[5] H. S. Malvar, A. Hallapuro, M. Karczewicz, and L. Kerofsky, “Low-
complexity transform and quantization in H.264/AVC,” IEEE Trans. on
Circuits and Systems for Video Technology, vol. 13, pp. 598–603, July
2003.
[6] Z. Wang, “Fast algorithms for the discrete W transform and for the
discrete Fourier transform,” IEEE Trans. on Acoustics, Speech, and Signal
Proc., vol. 32, no. 4, pp. 803–816, Aug. 1984.
[7] N. S. Jaynat and P. Noll, “Digital coding of waveforms,” Engle-wood
Cliffs, NJ: Prentice-Hall, 1984.
[8] open source, http://www.webmproject.org/code/.