R29. Discrete Cosine Transform in R

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: