Tanggal :August 5, 2020

R29. Discrete Cosine Transform in R

Spread the love

The Discrete Cosine Transform (DCT) leads to a real transform, whereas the fft leads to a complex transform. The relevant equations are given, for example, at Wikipedia

The program below finds the DCT and inverse DCT using fft and ifft. For both fft and ifft, we have to solve a 2N problem. For example, for the DCT, x is joined with its reverse. Besides the fft and ifft, we have to add phase terms and proper scaling constants.

The example x is of length 16. We plot x and its DCT. In R, the negative index indicates the element that we do not want.


# ex29.R

dct <- function(X) {
x <- c(X,rev(X))
N <- 2*length(X)
f0 <- 1/sqrt(2*N)
f1 <- 1/sqrt(N)
Y <- fft(x)
kin <- 0:(N-1)
Mul <- cos(pi/N*kin)-1i*sin(pi/N*kin)
Y <- Y * Mul
Y <- Y * f1
Y[1] <- Y[1]*f0/f1
Re(Y[1:(N/2)])
}

idct <- function(X) {
N <- length(X)
f0 <- 1/sqrt(4*N)
f1 <- 1/sqrt(2*N)
kin <- 0:(N-1)
FTR <- X*cos(.5*pi/N*kin)
FTI <- X*sin(.5*pi/N*kin)
FTR <- FTR/f1
FTI <- FTI/f1
FTR[1] <- FTR[1] * f1/f0
FTI[1] <- FTI[1] * f1/f0
FTR1 <- c(FTR,0,rev(FTR[-1]))
FTI1 <- c(FTI,0,rev(-1*FTI[-1]))
y <- fft(FTR1+1i*FTI1, TRUE)/length(FTR1)
Re(y[1:N])
}

x <- c(0:3,2:4,10:2)
y <- dct(x)
z <- idct(y)
err <- sd(z-x)
cat('err =',err,'n')
plot(y, col = 'blue', type = 'b',
xlab = 'n', ylab = 'x, y',
main = 'x (red), y = dct(x) (blue) ')
points(x, col = 'red', type = 'b')
# err = 1.110223e-16

Output:

Share

Leave a Reply

Your email address will not be published. Required fields are marked *