How to perform a *Fast* DCT (Discrete Cosine Trans

问题:

Using Rprof revealed the dct in the dtt package to be the main offender in a piece of R code that was running quite slowly. Swapping it out for fft in the stats package (which is not the same transformation, but should take the same time to compute) my run time improved dramatically. In fact 2/3 of my Rprof lines were dct calls previously and a mere 3 lines of about 600 were fft calls after making the switch.

Is the dct implementation in the dtt package not done using the fast discrete fourier transform? If so, is there a package that does have one? (I know one could double ones data, then extract coefficients for the dct from those fft coefficients, but a straight fast dct would certainly be nicer, and there really ought to be one somewhere).

回答1:

There appears to be no fast dct, but there is an fft (fast fourier transform) in the stats package, so here is how you could go about getting the fast dct using fft.

Use this at your own risk. I haven't done any serious checking of it. I checked it on a couple of vectors of different sizes and it gave the same results as function dct in package dtt on those. If anyone wants to double check me by comparing it to the output from dct then feel free to do so and post your results.

Take your vector and extend it to a vector twice as long as follows: If your vector is v=(1,2,3) then double the entries to w=(1,2,3,3,2,1). Note the ordering. If your vector is v=(1,2,4,9) then double the entries to w=(1,2,4,9,9,4,2,1)

Let N be the length of your ORIGINAL vector (before you doubled its length).

Then the first N coefficients of .5 * fft(w)/exp(complex(imaginary=pi / 2 / N)*(seq(2*N)-1)) should agree with computing dct(v) except it should be dramatically faster in almost all cases.

Speed considerations. If you prime factor N then the time it takes to compute the fast dct is like the time it takes to do a slow dct for each of those prime factors. So if N is 2^K it is like doing a K different slow dct transforms on a vector of length two, so its really fast. If N is prime (worst case scenario) then there is no speed up at all. The greatest speedup is on vectors that are a power of two in length.

Note: The R code above looks incredibly unfriendly, so let me say what is going on. After doubling the length in the right way, the first N coefficients of the fft you get are almost the right thing. However the coefficients need to be tweaked a bit. Let P stand for e^(pi * i / 2 / N). Leave the first coefficient alone. Divide the second coefficient by P, divide the third by P^2, divide the fourth by P^3, etc... Then divide all the coefficients by 2 (including the first one) to agree with the normalization R uses for the dct.

This should give the same thing as using the dct function in package dtt but be dramatically faster in almost all cases.

回答2:

Disclaimer: I have never used the dtt package and cannot compare my results to its results. My answer is generic regarding DCT/FFT.

John's idea is correct, but he is two-off regarding repeating the vector and has to compensate it by tweaking the coefficients (the e^(pi * i / 2 / N) factor). With the right extension of the original vector, the FFT produces directly correct results(*). To quote from Wikipedia:

"The DCT-I is exactly equivalent (up to an overall scale factor of 2), to a DFT of 2N-2 real numbers with even symmetry. For example, a DCT-I of N=5 real numbers abcde is exactly equivalent to a DFT of eight real numbers abcdedcb (even symmetry), divided by two."

I.e., if we have a vector v = [1, 2, 3, 4, 5] on which we wish to perform the DCT, we should construct a new vector w = [1, 2, 3, 4, 5, 4, 3, 2] and perform the FFT on it. Notice that the first and the last component of v appear only once in w, in the original order!

This works because Fourier transform of an even function (function symmetric around zero) consists purely of real (cosine) coefficients. If we'd construct the vector w by including the whole reversed v, as John suggested, it would be symmetric around -0.5. Due to this tiny shift the Fourier transform would also produce imaginary (sine) coefficients.

(*) The method here produces DCT-I. John's method seems to aim at DCT-II.

回答3:

Both John and Igor are correct in their answers, but I thought they lacked some example code.

library(dtt)

par(mfrow=c(3, 1), mar=c(2, 3, 1, 0.1), mgp=c(2, 0.8, 0))

set.seed(1)
N <- 60
v <- sin(seq(0, pi*2*4, l=N))/2 + 
     sin(seq(0, pi*2*7, l=N))/3 + 
     sin(seq(0, pi*2*15, l=N))/2 + 
     runif(N, -1, 1)/40 +
     runif(N, -1, 1)/40

plot(v, type="o")

DCT-I:

v.dct1 <- dct(v, variant=1)

w <- c(v, v[(N - 1):2])
w.dct1 <- 0.5 * Re(fft(w)[1:N])


plot(v.dct1, type="l", col="#00000088")
abline(h=0, col="#00000044")
lines(w.dct1, col=2, lty=3, lwd=2)
legend("topright", c("dtt::dct", "fft"), bty="n", col=1:2, lwd=1)
legend("topleft", "DCT-I", bty="n")

DCT-II:

v.dct2 <- dct(v, variant=2)

P <- exp(complex(imaginary=pi / 2 / N)*(seq(2*N)-1)) 

w <- c(v, v[N:1])
w.dct2 <- 0.5 * Re(fft(w)[1:N]/P)

plot(v.dct2, type="l", col="#00000088")
abline(h=0, col="#00000044")
lines(w.dct2, col=2, lty=3, lwd=2)
legend("topright", c("dtt::dct", "fft"), bty="n", col=1:2, lwd=1)
legend("topleft", "DCT-II", bty="n")