Type-II/III DCT/DST algorithms with reduced number of arithmetic operations

0.22 MB·English

by  Xuancheng Shao

1
Type-II/III DCT/DST algorithms
with reduced number of arithmetic operations
Xuancheng Shao and Steven G. Johnson*
N previous bestDCT-II Newalgorithm
Abstract—Wepresentalgorithmsforthediscretecosinetrans-
16 114 112
form (DCT) and discrete sine transform (DST), of types II
9 and III, that achieve a lower count of real multiplications and 32 290 284
64 706 686
0 additionsthanpreviouslypublishedalgorithms, withoutsacrific-
128 1666 1614
0 ing numerical accuracy. Asymptotically, the operation count is
256 3842 3708
2 reduced from 2Nlog N +O(N) to 17Nlog N +O(N) for a
2 9 2 512 8706 8384
power-of-two transform size N. Furthermore, we show that an
n 1024 19458 18698
a additionalN multiplicationsmaybesavedbyacertainrescaling 2048 43010 41266
J oftheinputsoroutputs,generalizingawell-knowntechniquefor 4096 94210 90264
N = 8 by Arai et al. These results are derived by considering
9 the DCT to be a special case of a DFT of length 4N, with TABLEI
2 certainsymmetries,andthenpruningredundantoperationsfrom FLOPCOUNTS(REALADDS+MULTS)OFPREVIOUSBESTDCT-IIAND

a recent improved fast Fourier transform algorithm (based on a OURNEWALGORITHM
]
A recursive rescaling of the conjugate-pair split radix algorithm).
The improved algorithms for the DCT-III, DST-II, and DST-III
N
follow immediately from the improved count for the DCT-II.
s. Index Terms—discrete cosine transform; fast Fourier trans- of 34Nlog N+O(N) [1],[17].Similarly,the lowest-known
9 2
c form; arithmetic complexity flopcountfortheDCT-IIofsize N =2m >1waspreviously
[
2Nlog N N + 2 for a unitary normalization (with the
2 −
2 additive constant depending on normalization) [6], [7], [12],
I. INTRODUCTION
v
[13], [18]–[26], and could be achieved by starting with the
0 Inthispaper,wedescriberecursivealgorithmsforthetype-
split-radix FFT and discarding redundant operations [6], [7].
5
IIandtype-IIIdiscretecosineandsinetransforms(DCT-IIand
1 (ManyDCTalgorithmswithanunreportedorlargerflopcount
DCT-III,andalsoDST-IIandDST-III),ofpower-of-twosizes,
3 have also been described [27]–[40].) Based on our new FFT
0 that requirefewer totalreal additionsand multiplicationsthan algorithm, the flop counts for the various DCT types were
7 previouslypublishedalgorithms(withanasymptoticreduction reduced using a code generator [10], [41] that automatically
0 of about 5.6%), without sacrificing numerical accuracy. Our
pruned the redundant operations from an FFT with a given
/
s DCT and DST algorithms are based on a recently published symmetry, but neither an explicit algorithm nor a general
c fast Fourier transform (FFT) algorithm, which reduced the
formula for the flop count were presented except for DCT-I
:
v operation count for the discrete Fourier transform (DFT) [1].Inthispaper,weusethesamestartingpointto“manually”
i compared to the previous-best split-radix algorithm [1]. The
X derive a DCT-II algorithm by pruning redundant operations
gains in this new FFT algorithm, and consequently in the
from a real-even FFT, and give the general formula for the
r
a new DCTs and DSTs, stem from a recursive rescaling of new flop count (for N =2m >1):
the internalmultiplicativefactorswithin an algorithmcalled a
“conjugate-pair”split-radixFFT[2]–[5]soastosimplifysome 17 17 1
Nlog N N ( 1)log2Nlog N
ofthemultiplications.InordertoderiveaDCTalgorithmfrom 9 2 − 27 − 9 − 2
this FFT, we simply consider the DCT-II to be a special case 7 3
+ ( 1)log2N + . (1)
of a DFT with real input of a certain even symmetry, and 54 − 2
discard the redundant operations [1], [6]–[10]. The DCT-III,
The first savings over the previous record occur for N =16,
DST-II, and DST-III have identical operation counts to the
andaresummarizedinTableIforseveralN.Wealsoconsider
DCT-II of the same size, since the algorithms are related by
theeffectofnormalizationonthis flopcount:theabovecount
simple transpositions, permutations, and sign flips [11]–[13].
was for a unitary transform, but slightly different counts are
Since 1968, the lowest total count of real additions and
obtained by other choices. Moreover, we show that a further
multiplications, herein called “flops” (floating-point opera-
N multiplicationscanbesavedbyindividuallyrescalingevery
tions), for the DFT of a power-of-two size N = 2m was
output of the DCT-II (or input of the DCT-III). In doing so,
achievedbythesplit-radixalgorithm,with4Nlog N 6N+8
2 − we generalize a result by Arai et al., who showed that eight
flops for N > 1 [6], [8], [14]–[16]. This count was re-
multiplications could be saved by scaling the outputs of a
cently surpassed by new algorithms achieving a flop count
DCT-II of size N = 8 [42], a result commonly applied to
JPEG compression [43].
* Department of Mathematics, Massachusetts Institute of Technology,
CambridgeMA02139. If we merely wished to show that the DCT-II/III could be
2
computed in 17Nlog N + O(N) operations, we could do
9 2
so by using well-known techniques to re-express the DCT in
terms of a real-input DFT of length N plus O(N) pre/post-
processing operations [29], [30], [33], [44]–[46], and then
substituting our new FFT that requires 17Nlog N +O(N)
9 2
operations for real inputs [1]. However, with FFT and DCT
algorithms,thereisgreatinterestinobtainingnotonlythebest
possible asymptotic constant factor, but also the best possible
exact countofarithmeticoperations(which,fortheDCT-IIof
power-of-twosizes,hadnotchangedbyevenoneoperationfor Fig. 1. A DCT-II of size N =8 (open dots, x0,...,x7) is equivalent to
over20years).Ourresult(1)isintendedasanewupperbound asize-4N DFTviainterleaving withzeros (black dots)andextending inan
even(squares)periodic (gray)sequence.
on this (still unknown) minimum exact count, and therefore
we have done our best with the O(N) terms as well as the
atosyamchpiteovtiec(1co)nisstacnlot.seIltyturerlnasteodutto, iwneflla-cktn,otwhantaolugrorailtghomrsithfomr where ωN = e−2Nπi is an Nth primitive root of unity. In
order to relate this to the DCT-II, it is convenientto choose a
expressingtheDCT-IIintermsofareal-inputFFTofthesame
different normalization for the latter transform:
length, but it does not seem obviousa priori that this is what
one obtains by pruning our FFT for symmetric data. N−1 π 1
C =2 x cos n+ k . (4)
Inthefollowingsections,wefirstreviewhowaDCT-IImay k n
(cid:20)N (cid:18) 2(cid:19) (cid:21)
beexpressedasaspecialcaseofaDFT,andhowtheprevious nX=0
optimum flop count can be achieved by pruning redundant Thisnormalizationisnotunitary,butitismoredirectlyrelated
operations from a conjugate-pair split-radix FFT. Then, we totheDFTandthereforemoreconvenientforthedevelopment
briefly review the new FFT algorithm presented in [1], and of algorithms. Of course, any fast algorithm for C trivially
k
derive the new DCT-II algorithm. We follow by considering yields a fast algorithm for CII, although the exact count of
k
the effect of normalization and scaling. Finally, we consider requiredmultiplicationsdependsonthenormalization,anissue
the extension of this algorithm to algorithms for the DCT-III, we discuss in more detail in section V.
DST-II,andDST-III.We close withsome concludingremarks In order to derive C from the DFT formula, one can use
k
about future directions. We emphasize that no proven tight the identity 2cos(πℓ/N)=ω2ℓ +ω4N−2ℓ to write:
4N 4N
lowerboundontheDCT-IIflopcountiscurrentlyknown,and
wemakenoclaimthateq.(1)isthelowestpossible(although N−1 π 1
C =2 x cos n+ k
we have endeavored not to miss any obvious optimizations). k n (cid:20)N (cid:18) 2(cid:19) (cid:21)
nX=0
II. DCT-II IN TERMS OF DFT N−1 N−1
= x ω(2n+1)k+ x ω(4N−2n−1)k (5)
Various forms of discrete cosine transform have been de- n 4N n 4N
nX=0 nX=0
fined, corresponding to different boundary conditions on the
4N−1
transform [47]. Perhaps the most common form is the type- = x˜ ωnk,
n 4N
II DCT, used in image compression [43] and many other
nX=0
applications. The DCT-II is typically defined as a real, or-
wherex˜ isareal-evensequenceoflength4N (i.e.x˜ =
thogonal (unitary), linear transformation by the formula (for n 4N−n
x˜ ), defined as follows for 0 n<N:
k =0,...,N 1): n ≤
−
x˜ =x˜ =0, (6)
2n 4N−2n−2
N−1
2 δ π 1
CII = − k,0 x cos n+ k , (2)
k r N n (cid:20)N (cid:18) 2(cid:19) (cid:21) x˜2n+1 =x˜4N−(2n+1) =xn. (7)
nX=0
for N inputs x and N outputs CII, where δ is the Thus, the DCT-II of size N is preciselythe first N outputsof
n k k,0
Kroneckerdelta(=1fork =0and=0otherwise).However, aDFTofsize4N,ofreal-eveninputs,wheretheeven-indexed
we wish to emphasize in this paper that the DCT-II (and, inputs are zero.
indeed,alltypesofDCT)canbeviewedasspecialcasesofthe This is illustrated by an example, for N = 8, in figure 1.
discrete Fourier transform (DFT) with real inputs of a certain TheeightinputsoftheDCTareshownasopendots,whichare
symmetry, and where only a subset of the outputs need be interleaved with zeros (black dots) and extended to an even
computed. This (well known) viewpoint is fruitful because it (squares) periodic (gray dots) sequence of length 4N = 32
meansthatanyFFT algorithmforthe DFT leadsimmediately corresponding to the DFT. (The type-II DCT is distinguished
to a corresponding fast algorithm for the DCT-II simply by from the other DCT types by the fact that it is even about
discarding the redundant operations [1], [6]–[10]. boththe leftand rightboundariesof the originaldata, andthe
The discrete Fourier transform of size N is defined by symmetry points fall halfway in between pairs of the original
data points.) We will refer, below, to this figure in order to
N−1
X = x ωnk, (3) illustrate what happens when an FFT algorithm is applied to
k n N
this real-symmetric zero-interleaved data.
nX=0
3
Algorithm 1 Standard conjugate-pair split-radix FFT algo- Algorithm 2 Standard conjugate-pair split-radix FFT algo-
rithm of size N. Special-case optimizations for k = 0 and rithm of size N as in Algorithm 1, but specialized for the
k =N/8,aswellasthebasecases,areomittedforsimplicity. case of real inputs. Special-case optimizations for k = 0 and
function X splitfft (x ): k =N/8,aswellasthebasecases,areomittedforsimplicity.
k=0..N−1 ← N n
U splitfft (x ) function X rsplitfft (x ):
k2=0...N/2−1 ← N/2 2n2 k=0..N/2 ← N n
Z splitfft (x ) U rsplitfft (x )
k4=0...N/4−1 ← N/4 4n4+1 k2=0...N/4 ← N/2 2n2
Z′ splitfft (x ) Z rsplitfft (x )
k4=0...N/4−1 ← N/4 4n4−1 k4=0...N/8 ← N/4 4n4+1
for k =0 to N/4 1 do Z′ rsplitfft (x )
X U + ωk−Z +ω−kZ′ fokr4=k0=...N0/8to←N/8 do N/4 4n4−1
endXXXkkkkfo+++←rNN3N///424k←←←UU(cid:0)Ukkk+−N+NN/(cid:0)/4ωk4−Nk+Ziik(cid:0)N(cid:0)ω+ωNkNkωZkZN−k(cid:1)kk−Z−k′ω(cid:1)ωN−N−kkZZk′k′(cid:1)(cid:1) XXXXkNk+/←N4−/Uk4k←←+UUU(cid:0)NNω∗∗//Nk44−−Zkkkω−−+kZiiω(cid:0)(cid:0)N−ωω+kNNkkZωZZk′−(cid:1)kkk−−Zωω′ NN−−∗kkZZkk′′(cid:1)(cid:1)∗
N/2−k ← k − N k N k
end for (cid:0) (cid:1)
III. CONJUGATE-PAIRFFT AND DCT-II
Although the previous minimum flop count for DCT-II in the ordinary split-radix [7], [48] and in the conjugate-pair
algorithmscanbederivedfromtheordinarysplit-radixFFTal- split-radix[1],byeliminatingredundantoperations,toachieve
gorithm[6],[7](anditcanalsobederivedinotherways[12], a flop count of 2Nlog N 4N +6. Specifically, one need
2 −
[13],[18]–[26]),herewewilldothesamethingusingavariant only compute the outputs for 0 k N/2 (where the k =0
≤ ≤
dubbed the “conjugate-pair” split-radix FFT. This algorithm and k =N/2 outputs are purely real), and the corresponding
was originally proposed to reduce the number of flops [2], algorithmisshowninAlgorithm2.Again,forsimplicitywedo
but was later shown to have the same flop count as ordinary not show special-case optimizations for k =0 and k =N/8,
split-radix [3]–[5]. It turns out, however, that the conjugate- soitmayappearthattheX ,X ,andX outputsare
N/8 N/4 3N/8
pairalgorithmexposessymmetriesinthemultiplicativefactors computedtwice.ToderiveAlgorithm2,oneexploitstwofacts.
thatcanbeexploitedtoreducetheflopcountbyanappropriate First, the sub-transforms operate on real data and therefore
rescaling[1],whichwe willemployinthefollowingsections. have conjugate-symmetricoutput. Second, the twiddle factors
in Algorithm1 have some redundancybecause of the identity
A. Conjugate-pair split-radix FFT ωN±(N/4−k) =∓iωN∓k,whichallowsustosharetwiddlefactors
between k and N/4 k in the original loop.
Starting with the DFT of equation (3), the decimation-in- −
time conjugate-pairFFT splits it intothree smaller DFTs:one
of size N/2 of the even-indexedinputs, and two of size N/4: B. Fast DCT-II via split-radix FFT
Before we derive our fast DCT-II with a reduced flop
N/2−1 N/4−1
count, we first derive a DCT-II algorithm with the same flop
X = ωn2kx +ωk ωn4kx
k N/2 2n2 N N/4 4n4+1 countaspreviouslypublishedalgorithms,butstartingwiththe
nX2=0 nX4=0 conjugate-pair split-radix algorithm. This algorithm will then
N/4−1
be modified in section IV-B, below, to reduce the number of
+ω−k ωn4kx , (8)
N N/4 4n4−1 multiplications.
nX4=0 In eq. (5), we showed that a DCT-II of size N, with inputs
where the indices 4n4 ± 1 are computed modulo N. [In xn and outputs Ck, can be expressed as the first N outputs
contrast, the ordinary split-radix FFT uses x4n4+3 for the of a DFT of size N˜ =4N of real-eveninputsx˜n. [Therefore,
third sum (a cyclic shift of x4n4−1), with a corresponding N in eq. (8) above and in the surrounding text are replaced
multiplicative “twiddle” factor of ω3k.] This decomposition below by N˜.] Now, consider what happenswhen we evaluate
N
is repeated recursively (until base cases of size N = 1 or this DFT via the conjugate-pairsplit decompositionin eq. (8)
N = 2, not shown, are reached), as shown by the pseudo- and Algorithm 2. In this algorithm,we computethree smaller
code in Algorithm 1. Here, we denote the results of the three DFTs. First, U is the DFT of size N˜/2 = 2N of the x˜
k 2n
sub-transformsofsizeN/2,N/4,andN/4byUk,Zk,andZk′, inputs,but these are all zero andso Uk is zero. Second,Zk is
respectively. The number of flops required by this algorithm, the DFT of size N˜/4=N of the x˜ inputs(where4n+1
4n+1
after certain simplifications (common subexpression elimina- is evaluated modulo N˜), which by eq. (7) correspond to the
tion and constant folding, and special-case simplifications of original data x by the formula:
n
the constants for k = 0 and k = N/8) and not counting
x 0 n<N/2
d4Natal-oingd2eNpe−nd6eNnt +op8e,raitdieonntsiclailketothoerdicnoamrypustpaltiitornadoifx.ωNk, is y˜n =x˜4n+1 =(cid:26) x2N−21n−2n N/≤2≤n<N , (9)
In the following sections, we will have to exploit further where we have denoted them by y˜ (0 n < N) for
n
≤
simplifications for the case where the inputs x are real. In convenience. That is, Z is the real-input DFT of the even-
n k
this case, X =X∗ (where denotes complex conjugation) indexedelementsofx followedbytheodd-indexedelements
k −k ∗ n
and one can save slightly more than half of the flops, both of x in reverse order. For example, this is shown for N =8
n
4
this algorithm [33] only because they used a suboptimal real-
inputDFT for the subtransform(split-radixbeing then almost
unknown).Using a real-inputsplit-radix DFT to compute Z ,
k
theflopcountforZ is2Nlog N 4N+6fromabove.Toget
k 2 −
the total flop count for Algorithm 3, we need to add N/2 1
−
general complex multiplications by 2ωk (6 flops each) plus
4N
two real multiplications (2 flops), for a total of 2Nlog N
2 −
N +2 flops. This matches the best-known flop count in the
literature (where the +2 can be removed merely by choosing
a different normalization as discussed in section V).
cFoign.ju2g.ate-TpahierFDFCTT-oIfItohfeNN˜ ==382peoxitnentsdexdnda(otapex˜nnd(osqtsu)arisesc)ofmropmutfiedg,urien1a, IV. NEWFFT AND DCT-II
via the DFT Zk of the circled points x˜4n+1, corresponding to the even-
indexed xn followed bytheodd-indexed xn. In this section, we first review the new FFT algorithm
introducedinourpreviousworkbasedonarecursiverescaling
Algorithm 3 Fast DCT-II algorithm, matching previous best of the conjugate-pair FFT [1], and then apply it to derive a
flopcount,derivedfromAlgorithm2 bydiscardingredundant fast DCT-II algorithm as in the previous section.
operations.
function C splitdctII (x ): A. New FFT
k=0..N−1 ← N n
for n=0 to N/2 1 do Basedontheconjugate-pairsplit-radixFFTfromsectionIII,
−
y˜n x2n a new FFT algorithm with a reduced number of flops can
←
y˜N−1−n x2n+1 be derived by scaling the subtransforms [1]. We will not
←
end for reproduce the derivation here, but will simply summarize the
Zk=0...N/2 ←rsplitfftN(y˜n) results. In particular, the original conjugate-pair split-radix
C 2Z Algorithm 1 is split into four mutually recursive subroutines,
0 0
←
for k =1 to N/2 1 do each ofwhich hasthe same split-radixstructure butcomputes
−
C 2 ωk Z a DFT scaled by a different factor. These algorithms are
k ← ℜ 4N k
CN−k ←−(cid:0)2ℑ ω4k(cid:1)NZk shown in Algorithm 4 for the case of complex inputs, and
end for (cid:0) (cid:1) in Algorithm 5 specialized for real inputs, in both of which
C √2Z the scaling factors are combined with the twiddle factors ωk
N/2 ← N/2 N
to reduce the total number of multiplications compared to
Algorithms 1 and 2, respectively. Again, special cases for
k = 0 and k = N/8, as well as the N = 1,2 base cases
in figure 2 with the circled points corresponding to the y˜ ,
n
and obvious simplifications such as common-subexpression
which are clearly the even-indexed x followed by the odd-
n
indexedx inreverse.ThesecondDFT,Z′,ofsizeN˜/4=N elimination, are not shown for simplicity.
n k The key aspect of these algorithmsis the scale factor s ,
isredundant:itistheDFTofx˜ ,butbytheevensymmetry N,k
4n−1
where the subtransformscompute the DFT scaled by 1/s
of x˜ this is equal to x˜ , and therefore Z′ = Z∗ (the ℓN,k
n 4(−n)+1 k k for ℓ = 1,2,4. (The ℓ = 0,1,2 subroutines are condensed
complex conjugate of Z ).
k
into one in Algorithms 4–5, by making ℓ a variable, but in
Therefore, a fast DCT-II is obtained merely by computing
practice they need to be implemented separately in order to
a single real-input DFT of size N to obtain Z = Z∗ ,
k N−k exploit the special cases of the scale factor [1]. The ℓ = 4
and then combining this according to Algorithm 2 to obtain
case is written separately because it is factorized differently.)
C = X for k = 0...N 1. In particular, in the loop of
k k − This scale factor is defined for N = 2m by the following
Algorithm2, all but the X and X terms correspondto
k N˜/4−k recurrence, where k =k mod N:
subscripts N and are not needed. Moreover, we obtain 4 4
≥
X =ωkZ +ω−kZ∗ =2 ωk Z , 1 for N 4
k N˜ k N˜ k ℜ 4N k ≤
XN−k =−i(cid:16)ωNk˜Zk−ωN−˜kZk∗(cid:17)∗(cid:0)=−2ℑ(cid:1)(cid:0)ω4kNZk(cid:1). sN=2m,k = ssNN//44,,kk44csoins((22ππkk44//NN)) foothrekr4w≤iseN/8(1.0)
The resulting algorithm, including the special-case optimiza- 
Thisdefinitionhastheproperties:s =1,s =s ,
tionsfork =0 (whereω0 =1 andZ is real)and k =N/2 N,0 N,k+N/4 N,k
4N 0 and s = s . Also, s > 0 and decays rather
(where ω4NN/2 = (1−i)/√2 and ZN/2 is real), is shown in slowlyN,wNi/t4h−kN: sNN,k,khas an ΩN(,Nk log4cos(π/5)) asymptotic
Algorithm 3.
lower bound[1]. When these scale factors are combinedwith
In fact, Algorithm 3 is equivalent to an algorithm derived the twiddle factors ωk, we obtain terms of the form
N
in a different way by various authors [30], [33] to express a
s
DCT-IIintermsofareal-inputDFTofthesamelength.Here, t =ωk N/4,k, (11)
N,k N s
weseethatthisalgorithmisexactlyequivalenttoaconjugate- N,k
pairsplit-radixFFToflength4N.Previousauthorsobtaineda whichisalwaysacomplexnumberoftheform 1 itan2πk
± ± N
suboptimalflopcount 5Nlog N+O(N)fortheDCT-IIusing or cot2πk i and can therefore be multiplied with two
2 2 ± N ±
5
Algorithm 4 New FFT algorithm [1] of length N (divisible Algorithm5NewFFTalgorithmofsizeN (divisibleby4),as
by4).Thesub-transformsnewfftSℓ (x) forℓ=0are scaled in Algorithm4 butspecializedforthe case of realinputs.The
by s , respectively, while ℓ=0Nis the final6unscaled DFT sub-transformsrnewfftSℓ (x) forℓ=0 arescaled bys ,
ℓN,k N 6 ℓN,k
(s =1).Special-caseoptimizationsfork =0andk =N/8 respectively,whileℓ=0isthefinalunscaledDFT(s =1).
0,k 0,k
are not shown. Special-case optimizations for k = 0 and k = N/8 are not
function X newfftSℓ (x ): shown.
computeks=D0.F.NT−/1s← for ℓ=N0,1n,2 function X rnewfftSℓ (x ):
U{ neℓNw,kfftS2ℓ (x ) } computeks=D0.F.NT/2/ s← for ℓ=N0,1,n2
ZZfokkk′r244===k000=.........NNN0///244to−−−111N←←←/4nn−eeww1ffffdottSSN1N1N///244((xx244nnn244+−11)) ZZU{kkk′442===000.........NNN///884 ←←←rrrnnneeewwwℓNffffff,ktttSSS1N1N2Nℓ///442(((xxx442nnn442+−)11))}
X U + t Z +t∗ Z′ (s /s ) for k =0 to N/8 do
k ← k (cid:16) N,k k N,k k(cid:17)· N,k ℓN,k X U + t Z +t∗ Z′ (s /s )
XXkk++NN//42 ←UUkk+−N/(cid:16)4tN,kZk+t∗N,kZk′(cid:17)·(sN,k/sℓN,k) Xkk+←N/4k←U(cid:16)N∗N/4,−kkk N,k k(cid:17)· N,k ℓN,k
Xk++−3Nii(cid:16)(cid:16)/4tt←NN←,,kkZZUkkk+−−Ntt/∗N∗N4,,kkZZkk′′(cid:17)(cid:17)··((ssNN,,kk//ssℓℓNN,,kk++NN//44)) XXN−−/4−ii(cid:16)(cid:16)ktt←NN,,kkUUZZN∗kk/−−4−ttk∗N∗N,,kkZZkk′′(cid:17)(cid:17)∗··(s(sNN,k,k//ssℓNℓN,k,+k+NN/4/)4)
end for N/2−k ← k ∗
function Xk=0..N−1 ←newfftS4N(xn): −(cid:16)tN,kZk+t∗N,kZk′(cid:17) ·(sN,k/sℓN,k)
computes DFT / s end for
4N,k
U{Zk2=0...N/2−1 ←nneewwffff}ttSS2N1/2((xx2n2) ) funcctoimonpuXteks=D0.F.NT/2/ s←rnewfftS4N(xn):
Zfokk′r44==k00=......NN0//44to−−11N←←/4new1ffdotSN1N//44(x44nn44+−11) ZU{k2=0...N/4 ←rrnneeww4Nffff,ktt}SS12N/2((xx2n2) )
− k4=0...N/8 ← N/4 4n4+1
Xk ←hUk+(cid:16)tN,kZk+t∗N,kZk′(cid:17)i·(sN,k/s4N,k) fZok′r4=k0=...N0/8to←Nr/n8edwofftS1N/4(x4n4−1)
X U t Z +t∗ Z′
k+N/2 ←h (ksN−,k(cid:16)/sN4,Nk,k+kN/2N),k k(cid:17)i Xk ←hUk+(cid:16)tN,kZk+t∗N,kZk′(cid:17)i·(sN,k/s4N,k)
Xk+N/4 ←h·Uk+N/4−i(cid:16)tN,kZk−t∗N,kZk′(cid:17)i Xk+N/4 ←hU(N∗s/4−/ks−i(cid:16)tN,kZ)k−t∗N,kZk′(cid:17)i
(s /s ) N,k 4N,k+N/4
· N,k 4N,k+N/4 · ∗
Xk+3N/4 ←h·U(skN+,Nk//s44+N,ik(cid:16)+t3NN,/k4Z)k−t∗N,kZk′(cid:17)i XN/4−k ←h·U(NsN/4,−k/ks−4Ni,N(cid:16)t/N4−,kkZ)k−t∗N∗,kZk′(cid:17) i
end for X U∗ t Z +t∗ Z′
N/2−k ←h k −(cid:16) N,k k N,k k(cid:17) i
(s /s )
N,k 4N,N/2−k
·
end for
fewerrealmultiplicationsthanarerequiredtomultiplybyωk.
N
Because of the symmetry s = s , it is possible to
N,N/4−k N,k
specialize Algorithm 4 for real data in the same way as in
Algorithm 2, because the scale factor preserves the conjugate ThenumberM(N)ofrealmultiplicationssavedoverordinary
symmetry Xk = XN∗−k of the outputs of all subtransforms split radix is given by:
[1].1
The resulting flop count, for arbitrary complex data x , is
n
then reduced from 4Nlog2N −6N +8 to M(N)= 2Nlog N 38N +2log N
9 2 − 27 2
34 124
Nlog N N 2log N 2 16
9 2 − 27 − 2 + 9(−1)log2Nlog2N − 27(−1)log2N. (13)
2 16
( 1)log2Nlog N + ( 1)log2N +8. (12)
− 9 − 2 27 −
Inparticular,assumingthatcomplexmultiplicationsareimple- If the data are purely real, it was shown that M(N)/2
mentedasfourrealmultiplicationsandtworealadditions,then multiplications are saved over the corresponding real-input
the savings are purely in the number of real multiplications. split-radixalgorithm[1].Theseflopcountsaretocomputethe
original, unscaled definition of the DFT. If one is allowed to
1In Algorithm 5, we have also utilized various symmetries of the scale scaletheoutputsbyanyfactordesired,thenscalingby1/s
factors, such as s2N,N/4−k = s2N,k+N/4, as described in our previous [corresponding to newfftS1 (x ) in Algorithm 4], savesNa,nk
work [1], to make explicit the shared multiplicative factors between the N n
different terms. additional MS(N) M(N) 0 multiplications for complex
− ≥
6
Algorithm6NewDCT-IIalgorithmofsizeN =2m,basedon as for Algorithm 3 because the products ωk s can be
4N N,k
discarding redundant operations from the new FFT algorithm precomputed. Therefore, we save a total of M (N)/2 flops
S
of size 4N, achieving new record flop count. comparedtothepreviousbestflopcountof2Nlog N N+2,
2 −
function C newdctII (x ): resultingintheflopcountofeq.(1).Thisformulamatchesthe
k=0..N−1 N n
←
for n=0 to N/2 1 do flop count that was achieved by an automatic code generator
−
y˜ x in Ref. [1].
n 2n
←
y˜ x Because the new DCT algorithm is mathematically merely
N−1−n 2n+1
←
end for a special set of inputsfor the underlyingFFT, it will have the
Z rnewfftS1 (y˜ ) same favorable logarithmic error bounds as discussed in the
k=0..N/2 ← N n
C 2Z previous section.
0 0
←
for k =1 to N/2 1 do
−
Ck ←2ℜ ω4kNsN,kZk V. NORMALIZATIONS
CN−k ←−(cid:0)2ℑ ω4kNsN(cid:1),kZk In the above definition of DCT-II, we use a scale factor
end for (cid:0) (cid:1) “2” in front of the summation in order to make the DCT-
C √2Z
N/2 N/2 II directly equivalent to the corresponding DFT. But in some
←
circumstances,itisusefultomultiplybyotherfactors,anddif-
ferentnormalizationsleadtoslightlydifferentflopcounts.For
inputs, where: example,onecouldusetheunitarynormalizationfromeq.(2),
whichreplaces2by 2/N or 1/N (fork=0)andrequires
2 20
M (N)= Nlog N N thesamenumberofflpops.Ifonpeusestheunitarynormalization
S 9 2 − 27 multiplied by √N, then one saves two multiplications in the
2 7
+ ( 1)log2Nlog N ( 1)log2N +1. (14) k = 0 and k = N/2 terms (in this normalization, C0 = Z0
9 − 2 − 27 − and C = Z ) for both Algorithm 3 and Algorithm 6.
N/2 N/2
For real inputs, one similarly saves M (N)/2 flops for (This leads to a commonly quoted 2Nlog N N formula
rnewfftS1 (x ) operating on real x in ASlgorithm 5. for the previous flop count.) 2 −
N n n
Althoughthedivisionbya cosineorsineinthescalefactor It is also common to compute a DCT-II with scaled out-
may, at first glance, seem as if it may exacerbate floating- puts, e.g. for the JPEG image-compression standard where
point errors, this is not the case. Cooley–Tukey based FFT an arbitrary scaling can be absorbed into a subsequent quan-
algorithms have O(√logN) root-mean-square error growth, tization step [43], and in this case the scaling can save 8
and O(logN) error bounds, in finite-precision floating-point multiplications [42] over the 42 flops required for an unitary
arithmetic (assuming that the trigonometric constants are 8-point DCT-II. Since our newfftS1 (x ) attempts to be the
N n
precomputed accurately) [49]–[51], and the new FFT is no optimally scaled FFT, we should be able to derive this scaled
different [1]. The reason for this is straightforward: it never DCT-II by using rnewfftS1 (x )/2 for our length-4N DFT
N n
addsorsubtractsscaledandunscaledvalues.Instead,wherever instead of rnewfftS0 (x), and again pruning the redundant
N
theoriginalFFTwouldcomputea+b,thenewFFTcomputes operations. Doing so, we obtain an algorithm identical to
s (a+b) for some fixed scale factor s. Algorithm6exceptthat2ωk s isreplacedbyt ,which
· 4N N,k 4N,k
requires fewer multiplications, and also we now obtain C =
0
Z and C = Z . This algorithm computes C /2s
B. Fast DCT-II from new FFT 0 N/2 N/2 k 4N,k
insteadofC .Inparticular,wesaveexactlyN multiplications
k
To derive the new DCT-II algorithmbased on the new FFT
over Algorithm 6, which matches the result by Ref. [42]
oftheprevioussection,wefollowexactlythesameprocessas
for N = 8 but generalizes it to all N. This savings of N
in Sec. III-B. We express the DCT-II of length N in terms of
multiplicationsforascaledDCT-IIalsomatchestheflopcount
a DFT of length 4N, apply the FFT algorithm 5, and discard
that was achievedby an automatic code generatorin Ref. [1].
the redundant operations. As before, the U subtransform
k
is identically zero, Z is the transform of the even-indexed
k VI. FAST DCT-III,DST-II, AND DST-III
elements of x followed by the odd-indexed elements in
n Given any algorithm for the DCT-II, one can immediately
reverse order, and Z′ =Z∗ is redundant.
k N−k obtain algorithms for the DCT-III, DST-II, and DST-III with
TheresultingalgorithmisshowninAlgorithm6,anddiffers
exactlythesamenumberofarithmeticoperations.Inthisway,
fromtheoriginalfastDCT-IIofAlgorithm3inonlytwoways.
any improvement in the DCT-II, such as the one described
First, instead of calling the ordinary split-radix (or conjugate-
in this paper, immediately leads to improved algorithms for
pair) FFT for the subtransform, it calls newfftS (x ).
N n thosetransformsandviceversa.Inthissection,wereviewthe
Second, because this subtransform is scaled by 1/s , the
N,k equivalence between those transforms.
twiddle factor ωk in Algorithm 6 is multiplied by s .
4N N,k
To derive the flop count for Algorithm 6, we merely need
to add the flop count for the subtransform [which saves A. DCT-III
M (N)/2 = 1Nlog N + O(N) flops, from eq. (14), ToobtainaDCT-IIIalgorithmfromaDCT-IIalgorithm,one
S 9 2
compared to ordinary real-input split radix] with the number can exploit two facts. First, the DCT-III, viewed as a matrix,
of flops in the loop, where the latter is exactly the same is simply the transpose of the DCT-II matrix. Second, any
7
size-2 DCT-II size-2 DCT-III Algorithm 7 New DCT-III algorithm of size N =2m, based
x 2 C = 2(x + x) network C 2 x = 2C + 2 C on network transposition of Algorithm 6, with the same flop
0 0 0 1 transposition 0 0 0 1
2 2 2 2 count.
function D newdctIII (x ):
k=0..N−1 N n
←
x C = 2(x x) C x = 2C 2 C Z 2x
1 2 1 0 1 1 2 1 0 1 0 ← 0
for n=1 to N/2 1 do
Z 2ω−ns −(x ix )
Fig. 3. Example ofnetwork transposition forasize-2 DCT-II(left). When n ← 4N N,n n− N−n
thelinearnetworkistransposed(alledgesarereversed),theresultingnetwork Z Z∗
computes thetransposed-matrix operation: asize-2DCT-III(right). endNfo−rn ← n
Z √2x
N/2 N/2
←
y˜ rnewfftS(T)1(Z )
linear algorithm can be viewed as a linear network, and the k ← N n
for k =0 to N/2 1 do
transpose operation is computed by the network transposition
−
CT y˜
of this algorithm [52]. To review, a linear network represents 2k ← k
CT y˜
the algorithm by a directed acyclic graph, where the edges 2k+1 ← N−1−k
end for
represent multiplications by constants and the vertices repre-
sent additions of the incoming edges. Network transposition
simply consists of reversing the direction of every edge. We
DFT of the same size [30], [33]. In order to minimize the
prove below that the transposed network requires the same
flop count once the scale factors are included, however, it is
numberofadditionsandmultiplicationsfornetworkswith the
important that this real-output DFT operate on scaled inputs
samenumberofinputsandoutputs,andthereforetheDCT-III
ratherthanscaledoutputs,soitisintrinsicallydifferent(trans-
can be computed in the same number of flops as the DCT-
posed) from Algorithm 4. The real-output (scaled, conjugate-
II by the transposed algorithm.The DCT-III correspondingto
symmetric input) version of rnewfftS1 (x ) can be derived
the transpose of eq. (4) is N n
bynetworktranspositionof the real-inputcase (since network
N−1 π 1 transposition interchanges inputs and outputs). Equivalently,
CTk =2 xncos(cid:20)Nn(cid:18)k+ 2(cid:19)(cid:21). (15) whereasthernewfftS1N(xn)wasascaled-outputdecimation-
nX=0 in-time algorithm specialized for real inputs, the transposed
Since the DCT-III is the transpose of the DCT-II, a fast DCT- algorithmrnewfftS(T)1(x )isascaled-inputdecimation-in-
N n
III algorithm is obtained by network transposition of a fast frequencyalgorithmspecializedforrealoutputs.Wedonotlist
DCT-II algorithm. For example, network transposition of a rnewfftS(T)1(x ) explicitly here, however; our main point
N n
size-2 DCT-III is shown in figure. 3. is simply to establish that an algorithm for the DCT-III with
A proof that the transposed network requires the same exactly the same flop count (1) follows immediately from the
numberofflopsisasfollows.Clearly,thenumberofmultipli- new DCT-II.
cations,thenumberofedgeswithweight= 1,isunchanged There are also other ways to derive a DCT-III algorithm
6 ±
bytransposition.Thenumberofadditionsisgivenbythesum with the same flop countwithoutusingnetworktransposition,
of indegree 1 for all the vertices, except for the N input of course. For example, one can again consider the DCT-III
−
vertices which have indegree zero. That is, for a set V of to be a real-input DFT of size 4N with certain symmetries
vertices and a set E of edges: (differentfromtheDCT-II’ssymmetries),andpruneredundant
operationsfromAlgorithm5,resultinginadecimation-in-time
# adds=N + [indegree(v) 1]
− algorithm. Such a derivation is useful because it allows one
vX∈V
to derive a scaled-output DCT-III as well, which turns out to
=N + E V , (16)
| |−| | be a subroutine of a new DCT-IV algorithm that we describe
using the fact thatthe sum of the indegree overall verticesis in another manuscript, currently in press [53]. However, this
simplythenumber E ofedges.Becausetheaboveexpression derivation is more complicated than the one above, resulting
| |
is obviously invariant under network transposition as long as in four mutually recursive DCT-III subroutines, and does not
N is unchanged(i.e. same numberof inputs and outputs),the lower the flop count, so we omit it here.
number of additions is invariant.
More explicitly, a fast DCT-III algorithm derived from the B. DST-II and DST-III
transpose of our new Algorithm 6 is outlined in Algorithm 7. The DST-II is defined, using a unitary normalization, as:
WhereasfortheDCT-IIwecomputedthereal-inputDFT(with
N−1
2 δ π 1
cthoennjupgoastet--psyromcmesesterdictohuetpcuotm)polfextheourtepaurtrsanZgked=inpZuN∗ts−ky˜,nnaonwd SkII =r −Nk,N nX=0xnsin(cid:20)N (cid:18)n+ 2(cid:19)k(cid:21) (17)
we do the reverse. We first preprocess the inputs to form
for k = 1,...,N (not 0,...,N 1). As in section II, it is
a complex array Z = Z∗ , then perform a real-output, −
k N−k more convenient to develop algorithms for an unnormalized
scaled-inputDFTtoobtainy˜ ,andfinallyassigntheevenand
k variation:
oddelementsoftheresultCT fromthefirstandsecondhalves
k N−1
π 1
of y˜k. Without the scale factors, this is equivalent to a well- S =2 x sin n+ k (18)
k n
knownalgorithmtoexpressaDCT-IIIintermsofarealoutput (cid:20)N (cid:18) 2(cid:19) (cid:21)
nX=0
8
similar to our normalization of C and CT above. Although by arithmetic [10]. Nevertheless, minimal-arithmetic hard-
k k
we could derive fast algorithms for S directly by treating it coded DCTs of small sizes are often used in audio and
k
as a DFT of length 4N with odd symmetry, interleaved with image compression[43],[47],and the availability of any new
zeros,anddiscardingredundantoperationssimilartoabove,it algorithm with regular structure amenable to implementation
turns out there is a simpler technique. The DST-II is exactly leads to rich new opportunitiesfor performanceoptimization.
equivalent to a DCT-II in which the outputs are reversed and Similar arithmetic improvementsalso occur for other types
every other input is multiplied by 1 [11]–[13]: ofDCTandDST[1],aswellasforthemodifiedDCT(MDCT,
−
closely related to a DCT-IV), and the explicit description of
N−1
π 1
S =2 ( 1)nx cos n+ k (19) thosealgorithmsisthesubjectofanothermanuscriptcurrently
N−k n
nX=0 − (cid:20)N (cid:18) 2(cid:19) (cid:21) inpress[53].Wehavealreadyshownthatthenewexactcount
fortheDCT-Iis2Nlog N 3N+2log N+5 M(2N)/4=
fork =0,...,N 1.ItthereforefollowsthataDST-IIcanbe 2 − 2 −
− 17Nlog N +O(N), where M(N) is given by eq. (13) [1].
computed with the same number of flops as a DCT-II of the 9 2
However, no new exact count for arbitrary N = 2m has yet
same size, assuming that multiplications by 1 are free—the
− beenpublishedfortheDCT-IVandrelatedtransforms.Again,
reason for this is that sign flips can be absorbed at no cost
it immediately follows that the asymptotic cost of a DCT-IV
by converting additions into subtractions or vice versa in the
(and hence an MDCT) is reduced to 17Nlog N + O(N)
subsequent algorithmic steps. Therefore, our new flop count 9 2
simply by applying known algorithms to express a DCT-
(1) immediately applies to the DST-II.
IV in terms of a DCT-II of size N [12] (although this
Similarly, a DST-III is given by the transpose of the DST-
algorithm has large numerical errors [10]) or in terms of
II (which is the inverse, for the unitary normalization). In
two DCT-II/III transforms of size N/2 [19]. However, again
unnormalized form, this is
we wish to establish a new lowest-known upper bound for
N π 1 the exact count of operations of the DCT-IV, and not just
ST =2 x sin n k+ (20)
k n (cid:20)N (cid:18) 2(cid:19)(cid:21) the asymptotic constant factor. (It turns out that the DCT-IV
nX=1 algorithm obtained by pruning redundant computations from
for k = 0,...,N 1. This must have the same number of our new FFT, this time of symmetric data with size 8N [10],
−
flops as the DST-II by the network transposition argument is closely related to the algorithm by Wang et al. [19], just
above.Alternatively,it can also be obtainedfrom the DCT-III as as the algorithm for the DCT-II in this paper was closely
by reversing the inputs and multiplying every other outputby related to a known algorithm derived by other means [30],
1 [12]: [33].)
−
This work was supported in part by a grant from the
N−1
π 1
ST =( 1)k x cos n k+ , (21) MIT Undergraduate Research Opportunities Program. The
k − nX=0 N−n (cid:20)N (cid:18) 2(cid:19)(cid:21) authorsarealsogratefultoM.Frigo,co-authorofFFTWwith
S. G. Johnson [10], for many helpful discussions.
which can be obtained from CT with the same number of
k
flops.
REFERENCES
VII. CONCLUSION [1] S. G. Johnson and M. Frigo, “A modified split-radix FFT with fewer
arithmetic operations,” IEEE Trans. Signal Processing, vol. 55, no. 1,
We have shown that a improved count of real additions
pp.111–119, 2007.
and multiplications (flops), eq. (1), can be obtained for the [2] I. Kamar and Y. Elcherif, “Conjugate pair fast Fourier transform,”
DCT/DST types II/III by pruning redundant operations from Electron. Lett.,vol.25,no.5,pp.324–325,1989.
a recent improved FFT algorithm. We have also shown how [3] R. A. Gopinath, “Comment: Conjugate pair fast Fourier transform,”
Electron. Lett.,vol.25,no.16,p.1084,1989.
to save N multiplications by rescaling the outputs of a DCT- [4] H.-S. Qian and Z.-J. Zhao, “Comment: Conjugate pair fast Fourier
II (or the inputs of a DCT-III), generalizing a well-known transform,”Electron.Lett.,vol.26,no.8,pp.541–542,1990.
technique for N = 8 [42]. However, we do not claim that [5] A.M.KrotandH.B.Minervina,“Comment:ConjugatepairfastFourier
transform,”Electron.Lett.,vol.28,no.12,pp.1143–1144,1992.
ournew countis optimal—itmay be thatfurthergainscan be
[6] M. Vetterli and H. J. Nussbaumer, “Simple FFT and DCT algorithms
obtained by, for example, extending the rescaling technique with reduced number of operations,” Signal Processing, vol. 6, no. 4,
of [1] to greater generality. In particular, as has been pointed pp.267–278, 1984.
[7] P.Duhamel,“Implementation of“split-radix” FFTalgorithms forcom-
out by other authors [8], any improvement in the arithmetic
plex,real,andreal-symmetricdata,”IEEETrans.Acoust.,Speech,Signal
complexity or flop counts for the DFT immediately leads to Processing,vol.34,no.2,pp.285–295, 1986.
corresponding improvements in DCTs/DSTs, and vice versa. [8] P.DuhamelandM.Vetterli, “FastFouriertransforms:atutorial review
and a state of the art,” Signal Processing, vol. 19, pp. 259–299, April
Our most important result, we believe, is the fact that there
1990.
suddenlyappearstobenewroomforimprovementinproblems [9] R. Vuduc and J. Demmel, “Code generators for automatic tuning of
that had seen no reductions in flop counts for many years. numericalkernels:experiences withFFTW,”inProc.Semantics,Appli-
cation, andImplementation ofCodeGenerators Workshop, (Montreal),
The question of the minimal operation counts for basic
Sept.2000.
lineartransformationssuchasDCTsandDSTsisoflongstand- [10] M. Frigo and S. G. Johnson, “The design and implementation of
ing theoretical interest. The practical impact of a few percent FFTW3,”Proc.IEEE,vol.93,no.2,pp.216–231,2005.
[11] Z.Wang,“Afastalgorithmforthediscretesinetransformimplemented
improvement in flop counts, admittedly, is less clear because
by the fast cosine transform,” IEEE Trans. Acoust., Speech, Signal
computationtime on modernhardwareisoftennotdominated Processing,vol.30,no.5,pp.814–815, 1982.
9
[12] S. C. Chan and K. L. Ho, “Direct methods for computing discrete [40] G. Plonka and M. Tasche, “Fast and numerically stable algorithms
sinusoidal transforms,” IEE Proceedings F, vol. 137, no. 6, pp. 433– for discrete cosine transforms,” Lin. Algebra and its Appl., vol. 394,
442,1990. pp.309–345, 2005.
[13] P. Lee and F.-Y. Huang, “Restructured recursive DCT and DST algo- [41] M. Frigo, “A fast Fourier transform compiler,” in Proc. ACM SIG-
rithms,”IEEETrans.SignalProcessing,vol.42,no.7,pp.1600–1609, PLAN’99 Conference on Programming Language Design and Imple-
1994. mentation(PLDI),vol.34,(Atlanta,Georgia),pp.169–180,ACM,May
[14] R. Yavne, “An economical method for calculating the discrete Fourier 1999.
transform,”inProc.AFIPSFallJointComputerConf.,vol.33,pp.115– [42] Y.Arai,T.Agui,andM.Nakajima,“AfastDCT-SQschemeforimages,”
125,1968. Trans.IEICE,vol.71,no.11,pp.1095–1097, 1988.
[15] P. Duhamel and H. Hollmann, “Split-radix FFT algorithm,” Electron. [43] W.B.PennebakerandJ.L.Mitchell,JPEGStillImageDataCompres-
Lett.,vol.20,no.1,pp.14–16,1984. sionStandard. NewYork:VanNostrandReinhold, 1993.
[16] J.B.Martens,“Recursivecyclotomicfactorization—anewalgorithmfor [44] P. N. Swarztrauber, “Vectorizing the FFTs,” in Parallel Computations
calculatingthediscreteFouriertransform,”IEEETrans.Acoust.,Speech, (G.Rodrigue, ed.),pp.52–83,NewYork:Academic Press,1982.
SignalProcessing,vol.32,no.4,pp.750–761, 1984. [45] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling,
[17] T.LundyandJ.VanBuskirk,“AnewmatrixapproachtorealFFTsand Numerical Recipes in C: The Art of Scientific Computing. New York,
convolutionsoflength2k,”Computing,vol.80,no.1,pp.23–45,2007. NY:Cambridge Univ.Press,2nded.,1992.
[18] B. G. Lee, “A new algorithm to compute the discrete cosine trans- [46] C.vanLoan,ComputationalFrameworksfortheFastFourierTransform.
form,”IEEETrans.Acoust., Speech, Signal Processing, vol. 32, no.6, Philadelphia: SIAM,1992.
pp.1243–1245, 1984. [47] K.R.Rao andP.Yip, Discrete CosineTransform: Algorithms, Advan-
[19] Z. Wang, “On computing the discrete Fourier and cosine transforms,” tages, Applications. Boston,MA:Academic Press,1990.
IEEE Trans. Acoust., Speech, Signal Processing, vol. 33, no. 4, [48] H.V.Sorensen,D.L.Jones,M.T.Heideman,andC.S.Burrus,“Real-
pp.1341–1344, 1985. valuedfastFouriertransformalgorithms,”IEEETrans.Acoust.,Speech,
[20] N. Suehiro and M. Hatori, “Fast algorithms for the DFT and other SignalProcessing,vol.35,no.6,pp.849–863, 1987.
sinusoidaltransforms,”IEEETrans.Acoust.,Speech,SignalProcessing, [49] W.M.GentlemanandG.Sande,“FastFouriertransforms—forfunand
vol.34,no.3,pp.642–644, 1986. profit,”Proc.AFIPS,vol.29,pp.563–578, 1966.
[21] H. S. Hou, “A fast algorithm for computing the discrete cosine trans- [50] J. C. Schatzman, “Accuracy of the discrete Fourier transform and the
form,”IEEETrans.Acoust.,Speech,SignalProcessing,vol.35,no.10, fast Fourier transform,” SIAM J. Scientific Computing, vol. 17, no. 5,
pp.1455–1461, 1987. pp.1150–1166, 1996.
[51] M.TascheandH.Zeuner, HandbookofAnalytic-Computational Meth-
[22] F.ArguelloandE.L.Zapata,“Fastcosinetransformbasedonthesuc-
odsinAppliedMathematics,ch.8,pp.357–406. BocaRaton,FL:CRC
cessivedoublingmethod,”Electronicletters,vol.26,no.19,pp.1616–
Press,2000.
1618,1990.
[52] R. E. Crochiere and A. V. Oppenheim, “Analysis of linear digital
[23] C. W. Kok, “Fast algorithm for computing discrete cosine transform,”
networks,”Proc.IEEE,vol.63,no.4,pp.581–595, 1975.
IEEETrans.SignalProcessing, vol.45,no.3,pp.757–760,1997.
[53] X.ShaoandS.G.Johnson,“Type-IVDCT,DST,andMDCTalgorithms
[24] J. Takala, D. Akopian, J. Astola, and J. Saarinen, “Constant geometry
with reduced numbers of arithmetic operations,” Signal Processing,
algorithmfordiscretecosinetransform,”IEEETrans.SignalProcessing,
2008. Inpress.
vol.48,no.6,pp.1840–1843, 2000.
[25] M. Pu¨schel and J. M. F. Moura, “The algebraic approach to the
discrete cosine and sine transforms and their fast algorithms,” SIAM
J.Computing, vol.32,no.5,pp.1280–1316, 2003.
[26] M. Pu¨schel, “Cooley–Tukey FFT like algorithms for the DCT,” in
Proc.IEEEInt’lConf.Acoustics,Speech,andSignalProcessing,vol.2,
(HongKong),pp.501–504,April2003.
[27] N. Ahmed, T. Natarajan, and K. R. Rao, “Discrete cosine transform,”
IEEETrans.Comput.,vol.23,no.1,pp.90–93,1974.
[28] R.M.Haralick,“Astorageefficientwaytoimplementthediscretecosine
transform,”IEEETrans.Comput.,vol.25,pp.764–765,1976.
[29] W.-H. Chen, C. H. Smith, and S. C. Fralick, “A fast computational
algorithm for the discrete cosine transform,” IEEE Trans. Commun.,
vol.25,no.9,pp.1004–1009, 1977.
[30] M.J.NarasimhaandA.M.Peterson,“Onthecomputationofthediscrete
cosine transform,”IEEETrans.Commun.,vol.26,no.6,pp.934–936,
1978.
[31] B. D. Tseng and W. C. Miller, “On computing the discrete cosine
transform,”IEEETrans.Comput.,vol.27,no.10,pp.966–968,1978.
[32] L. P. Yaroslavskii, “Shifted discrete Fourier transforms,” Problems of
Information Transmission,vol.15,no.4,pp.324–327, 1979.
[33] J.Makhoul,“Afastcosinetransforminoneandtwodimensions,”IEEE
Trans. Acoust., Speech, Signal Processing, vol. 28, no. 1, pp. 27–34,
1980.
[34] H. S. Malvar, “Fast computation of the discrete cosine transform and
the discrete Hartley transform,” IEEE Trans. Acoust., Speech, Signal
Processing,vol.35,no.10,pp.1484–1485,1987.
[35] W. Li, “A new algorithm to compute the DCT and its inverse,” IEEE
Trans.SignalProcessing,vol.39,no.6,pp.1305–1313,1991.
[36] G. Steidl and M. Tasche, “A polynomial approach to fast algorithms
for discrete Fourier-cosine and Fourier-sine transforms,” Math. Comp.,
vol.56,no.193,pp.281–296,1991.
[37] E.FeigandS.Winograd,“Onthemultiplicative complexityofdiscrete
cosinetransforms,”IEEETrans.Info.Theory,vol.38,no.4,pp.1387–
1391,1992.
[38] J.AstolaandD.Akopian,“Architecture-oriented regularalgorithmsfor
discrete sine and cosine transforms,” IEEE Trans. Signal Processing,
vol.47,no.4,pp.1109–1124, 1999.
[39] Z.Guo, B.Shi, and N.Wang, “Twonew algorithms based onproduct
system for discrete cosine transform,” Signal Processing, vol. 81,
pp.1899–1908, 2001.