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).
Respuestas:
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 ([ 0 , 1 , 2 , . . . , N- 1 ]
N
es la longitud de la señal de entradax
,k
esarange(N)
= ):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: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 porTipo 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]
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*]
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. :REfuente
exp(-1j*pi*k/(2*N))
o hay un acceso directo a ese paso?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?dejar
El DCT es dado por
Entonces básicamente creas el2 N secuencia de longitud y( n ) donde es la primera mitad x ( n ) y la segunda mitad es x ( n ) 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.
Editar:
Nota: La fórmula DCT que está usando es:
Hay varias formas de escalar la suma para que no coincida exactamente con otras implementaciones. Por ejemplo, MATLAB usa:
dóndew ( 0 ) = 1norte--√ y w ( 1 ... N- 1 ) = 2norte--√
Puede dar cuenta de esto escalando adecuadamente la salida.
fuente
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/
fuente