Fast Cosine Transform via FFT

Saya ingin menerapkan Fast Cosine Transform. Saya membaca di wikipedia , bahwa ada versi cepat dari DCT yang juga dihitung untuk FFT. Saya mencoba membaca makalah yang dikutip Makhoul *, untuk implementasi FTPACK dan FFTW yang juga digunakan dalam Scipy , tetapi saya tidak dapat mengekstrak algoritma yang sebenarnya. Inilah yang saya miliki sejauh ini:

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

Kode 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 

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

* J. Makhoul, "Transformasi kosinus cepat dalam satu dan dua dimensi," IEEE Trans. Acoust. Pidato Sig. Proc 28 (1), 27-34 (1980).

Saya telah membaca tentang ini dan ada beberapa cara untuk melakukannya, menggunakan ukuran yang berbeda N. saya Matlab adalah berkarat, jadi di sini mereka berada di Python ( Nadalah panjang sinyal input x, kadalah arange(N)= ):$[0, 1, 2, ..., N-1]$

Tipe 2 DCT menggunakan 4N FFT dan tanpa shift

Sinyal [a, b, c, d]menjadi

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

Kemudian ambil FFT untuk mendapatkan spektrum

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

lalu buang semuanya kecuali yang pertama [A, B, C, D], dan Anda selesai:

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

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

Tipe 2 DCT menggunakan cermin 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

Tipe 2 DCT menggunakan bantalan 2N FFT (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

Tipe 2 DCT menggunakan 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

Di mesin saya, ini semua kira-kira kecepatan yang sama, karena menghasilkan exp(-1j*pi*k/(2*N))membutuhkan waktu lebih lama daripada FFT. : D

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

membiarkan

$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}$

DCT kemudian diberikan oleh

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

$2N$$y(n)$$x(n)$$x(n)$

Berikut kode di 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));

Edit:

Catatan: Formula DCT yang digunakan adalah:

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

Ada beberapa cara penskalaan penjumlahan sehingga mungkin tidak cocok persis dengan implementasi lainnya. Misalnya, MATLAB menggunakan:

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

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

Anda dapat menjelaskan ini dengan menskalakan output dengan benar.

Untuk komputasi ilmiah sejati, jumlah penggunaan memori juga penting. Oleh karena itu titik N FFT lebih menarik bagi saya. Ini hanya mungkin karena Simetri Hermit dari sinyal. Referensi Makhoul diberikan di sini. Dan sebenarnya memiliki algoritma untuk menghitung DCT dan IDCT atau DCT10 dan DCT01.
http://ieeexplore.ieee.org/abstract/document/1163351/