I want to implement the Fast Cosine Transform. I read on wikipedia, that there is a fast version of the DCT which is similarly computed to the FFT. I tried to read the cited Makhoul* paper, for the FTPACK and FFTW implementations that are also used in Scipy, but I were not able to extract the actually algorithm. This is what I have so far:
def fft(x): if x.size ==1: return x N = x.size x0 = my_fft(x[0:N:2]) x1 = my_fft(x[0+1:N:2]) k = numpy.arange(N/2) e = numpy.exp(-2j*numpy.pi*k/N) l = x0 + x1 * e r = x0 - x1 * e return numpy.hstack([l,r])
def dct(x): k = 0 N = x.size xk = numpy.zeros(N) for k in range(N): for n in range(N): xn = x[n] xk[k] += xn*numpy.cos(numpy.pi/N*(n+1/2.0)*k) return xk
def my_fct(x): if x.size ==1: return x N = x.size x0 = my_fct(x[0:N:2]) # have to be set to zero? x1 = my_fct(x[0+1:N:2]) k = numpy.arange(N/2) n = # ??? c = numpy.cos(numpy.pi/N*(n+1/2.0)*k) l = x0 #??? r = x0 #??? return numpy.hstack([l,r])
*J. Makhoul, "A fast cosine transform in one and two dimensions," IEEE Trans. Acoust. Speech Sig. Proc. 28 (1), 27-34 (1980).