Szybka transformacja kosinusowa przez FFT

Chcę wdrożyć Fast Cosine Transform. Czytałem na wikipedii , że istnieje szybka wersja DCT, która jest obliczana podobnie do FFT. Próbowałem przeczytać cytowany artykuł Makhoul * dla implementacji FTPACK i FFTW, które są również używane w Scipy , ale nie byłem w stanie wyodrębnić faktycznego algorytmu. Oto co mam do tej pory:

Kod FFT:

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])

Kod DCT:

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 

Próba FCT:

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])

*JOT. Makhoul, „Szybka transformacja kosinusowa w jednym i dwóch wymiarach”, IEEE Trans. Acoust. Speech Sig. Proc. 28 (1), 27–34 (1980).

Nxkarange(N)$[0, 1, 2, ..., N-1]$

Typ 2 DCT przy użyciu 4N FFT i bez zmian

Sygnał [a, b, c, d]staje się

[0, a, 0, b, 0, c, 0, d, 0, d, 0, c, 0, b, 0, a].

Następnie weź FFT, aby uzyskać widmo

[A, B, C, D, 0, -D, -C, -B, -A, -B, -C, -D, 0, D, C, B]

następnie wyrzuć wszystko oprócz pierwszego [A, B, C, D]i gotowe:

u = zeros(4 * N)
u[1:2*N:2] = x
u[2*N+1::2] = x[::-1]

U = fft(u)[:N]
return U.real

Typ 2 DCT przy użyciu lustrzanego odbicia 2N FFT (Makhoul)

[a, b, c, d][a, b, c, d, d, c, b, a][A, B, C, D, 0, D*, C*, B*][A, B, C, D]$e^{-j \pi k \over 2 N}$

y = empty(2*N)
y[:N] = x
y[N:] = x[::-1]

Y = fft(y)[:N]

Y *= exp(-1j*pi*k/(2*N))
return Y.real

Typ 2 DCT z wykorzystaniem 2N FFT wyściełany (Makhoul)

[a, b, c, d][a, b, c, d, 0, 0, 0, 0][A, B, C, D, E, D*, C*, B*][A, B, C, D]$2e^{-j \pi k \over 2 N}$

y = zeros(2*N)
y[:N] = x

Y = fft(y)[:N]

Y *= 2 * exp(-1j*pi*k/(2*N))
return Y.real

Typ 2 DCT przy użyciu N FFT (Makhoul)

[a, b, c, d, e, f][a, c, e, f, d, b][A, B, C, D, C*, B*]$2 e^{-j \pi k \over 2N}$

v = empty_like(x)
v[:(N-1)//2+1] = x[::2]

if N % 2: # odd length
    v[(N-1)//2+1:] = x[-2::-2]
else: # even length
    v[(N-1)//2+1:] = x[::-2]

V = fft(v)

V *= 2 * exp(-1j*pi*k/(2*N))
return V.real

Na mojej maszynie są one w przybliżeniu tej samej prędkości, ponieważ generowanie exp(-1j*pi*k/(2*N))trwa dłużej niż FFT. :RE

In [99]: timeit dct2_4nfft(a)
10 loops, best of 3: 23.6 ms per loop

In [100]: timeit dct2_2nfft_1(a)
10 loops, best of 3: 20.1 ms per loop

In [101]: timeit dct2_2nfft_2(a)
10 loops, best of 3: 20.8 ms per loop

In [102]: timeit dct2_nfft(a)
100 loops, best of 3: 16.4 ms per loop

In [103]: timeit scipy.fftpack.dct(a, 2)
100 loops, best of 3: 3 ms per loop

$x(n)$

pozwolić

$y(n) = \Bigl\lbrace \begin{array}{ll} x(n), & n=0,1,...,N-1\\ x(2N - 1 - n), & n=N,N+1,...,2N-1 \end{array}$

Następnie podaje DCT

$C(k) = \mathrm{Re}\Bigl\lbrace e^{-j\pi \frac{k}{2N}} FFT\lbrace y(n)\rbrace\Bigr\rbrace$

Więc w zasadzie tworzysz $2N$ sekwencja długości $y(n)$ gdzie jest pierwsza połowa $x(n)$ a druga połowa to $x(n)$wywrócony. Następnie po prostu weź FFT i pomnóż ten wektor przez rampę fazową. Wreszcie weź tylko prawdziwą rolę i masz DCT.

Oto kod w MATLAB.

function C = fdct(x)
    N = length(x);
    y = zeros(1,2*N);
    y(1:N) = x;
    y(N+1:2*N) = fliplr(x);
    Y = fft(y);
    k=0:N-1;
    C = real(exp(-j.* pi.*k./(2*N)).*Y(1:N));

Edytować:

Uwaga: Formuła DCT, której używa, to:

$C(k) = 2\sum_{n=0}^{N-1} x(n) \cos\left(\frac{\pi k}{2N}(2n+1) \right)$

Istnieje kilka sposobów skalowania podsumowania, więc może nie pasować dokładnie do innych implementacji. Na przykład MATLAB używa:

$C(k) = w(k)\sum_{n=0}^{N-1} x(n) \cos\left(\frac{\pi k}{2N}(2n+1) \right)$

where $w(0) = \sqrt\frac{1}{N}$ and $w(1...N-1) = \sqrt{\frac{2}{N}}$

You can account for this by properly scaling the output.

For true scientific computing, the amount of memory usage is important too. Therefore the N point FFT is more attractive to me. This is only possible due to Hermitian Symmetry of the signal. The reference Makhoul is given here. And actually has the algorithm for calculating DCT and IDCT or DCT10 and DCT01.
http://ieeexplore.ieee.org/abstract/document/1163351/