Transformación rápida de coseno a través de FFT

15

Quiero implementar la transformación rápida de coseno. Leí en Wikipedia , que hay una versión rápida de DCT que se calcula de manera similar a la FFT. Traté de leer el documento citado de Makhoul *, para las implementaciones FTPACK y FFTW que también se usan en Scipy , pero no pude extraer el algoritmo real. Esto es lo que tengo hasta ahora:

Código 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])

Código 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 

Prueba de 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, "Una transformación rápida del coseno en una y dos dimensiones", IEEE Trans. Acust. Speech Sig. Proc. 28 (1), 27-34 (1980).

Framester
fuente
2
¿Pregunta si su código DCT es correcto o algo así?
Jim Clay
Gracias por tus comentarios. Agregué otra oración al principio. Mi objetivo es implementar el FCT sobre la base de la FFT.
Framester
El IDCT se administra en el cómputo del DCT inverso (IDCT) mediante DCT o IFFT .
Royi

Respuestas:

18

He estado leyendo sobre esto y hay varias formas de hacerlo, usando diferentes tamaños N. Mi Matlab está oxidado, así que aquí están en Python ( Nes la longitud de la señal de entrada x, kes arange(N)= ):[0 0,1,2,...,norte-1]

Tipo 2 DCT usando 4N FFT y sin turnos

La señal se [a, b, c, d]convierte

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

Luego toma la FFT para obtener el espectro

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

luego tira todo menos el primero [A, B, C, D], y listo:

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

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

Tipo 2 DCT usando 2N FFT reflejado (Makhoul)

[a, b, c, d]se convierte [a, b, c, d, d, c, b, a]. Tome la FFT de eso para obtener [A, B, C, D, 0, D*, C*, B*], luego deseche todo pero [A, B, C, D]multiplíquelo por mi-jπk2norte ( desplazamiento de media muestra ) para obtener el DCT:

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

Tipo 2 DCT usando 2N FFT acolchado (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]2mi-jπk2norte

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

Y = fft(y)[:N]

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

Tipo 2 DCT usando N FFT (Makhoul)

[a, b, c, d, e, f][a, c, e, f, d, b][A, B, C, D, C*, B*]2mi-jπk2norte

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

En mi máquina, todos estos tienen aproximadamente la misma velocidad, ya que la generación exp(-1j*pi*k/(2*N))lleva más tiempo que la 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
endolito
fuente
2
Gran respuesta, ¡me ayudó mucho con mi implementación! Nota adicional: El último método "Tipo 2 DCT usando N FFT" todavía funciona correctamente si la longitud de la señal es impar; el último elemento se mueve al elemento del medio. He verificado las matemáticas y el código para este hecho.
Nayuki el
1
@Nayuki ¿Estás generando exp(-1j*pi*k/(2*N))o hay un acceso directo a ese paso?
endolith
Estoy generando exp(-1j*pi*k/(2*N))en mi código , porque es necesario un cambio de un cuarto de muestra para que el mapeo de DCT a DFT funcione. ¿Qué te hace preguntar?
Nayuki
Hola, ¿cómo funcionaría esto para el DCT Tipo III, para calcular el inverso del DCT-II?
Jack H
8

X(norte)

dejar

y(norte)={X(norte),norte=0 0,1,...,norte-1X(2norte-1-norte),norte=norte,norte+1,...,2norte-1

El DCT es dado por

C(k)=Rmi{mi-jπk2norteFFT{y(norte)}}

Entonces básicamente creas el 2norte secuencia de longitud y(norte) donde es la primera mitad X(norte) y la segunda mitad es X(norte)invertido Luego, tome la FFT y multiplique ese vector por una rampa de fase. Finalmente, tome solo la parte real y tendrá el DCT.

Aquí está el código en 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));

Editar:

Nota: La fórmula DCT que está usando es:

C(k)=2norte=0 0norte-1X(norte)cos(πk2norte(2norte+1))

Hay varias formas de escalar la suma para que no coincida exactamente con otras implementaciones. Por ejemplo, MATLAB usa:

C(k)=w(k)norte=0 0norte-1X(norte)cos(πk2norte(2norte+1))

dónde w(0 0)=1norte y w(1 ...norte-1)=2norte

Puede dar cuenta de esto escalando adecuadamente la salida.

Jason B
fuente
1
Se supone que y (n) tiene una longitud N, no una longitud 2N. Así es como obtienes la velocidad de cálculo 4x, calculando DCT de longitud N a partir de la señal de longitud N en lugar de 2N FFT a partir de la señal 2N. fourier.eng.hmc.edu/e161/lectures/dct/node2.html www-ee.uta.edu/dip/Courses/EE5355/Discrete%20class%201.pdf
endolith
0

Para la verdadera computación científica, la cantidad de uso de memoria también es importante. Por lo tanto, el punto N FFT es más atractivo para mí. Esto solo es posible debido a la simetría hermitiana de la señal. La referencia de Makhoul se da aquí. Y en realidad tiene el algoritmo para calcular DCT e IDCT o DCT10 y DCT01.
http://ieeexplore.ieee.org/abstract/document/1163351/

Hasbestein
fuente