Transformation cosinus rapide via FFT

15

Je veux implémenter la transformation cosinus rapide. J'ai lu sur wikipedia , qu'il existe une version rapide du DCT qui est calculée de la même manière que la FFT. J'ai essayé de lire l'article Makhoul * cité , pour les implémentations FTPACK et FFTW qui sont également utilisées dans Scipy , mais je n'ai pas pu extraire l'algorithme en fait. Voici ce que j'ai jusqu'à présent:

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

Code 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 

Essai 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, "Une transformation cosinus rapide en une et deux dimensions", IEEE Trans. Acoust. Speech Sig. Proc. 28 (1), 27-34 (1980).

Framester
la source
2
Demandez-vous si votre code DCT est correct ou quelque chose?
Jim Clay
Merci pour vos commentaires. J'ai ajouté une autre phrase au début. Mon objectif est de mettre en œuvre le FCT sur la base de la FFT.
Framester

Réponses:

18

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

DCT de type 2 utilisant 4N FFT et sans décalage

Le signal [a, b, c, d]devient

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

Prenez ensuite la FFT pour obtenir le spectre

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

puis jetez tout sauf le premier [A, B, C, D], et vous avez terminé:

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

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

DCT de type 2 utilisant un miroir FFT 2N (Makhoul)

[a, b, c, d] devient [a, b, c, d, d, c, b, a][A, B, C, D, 0, D*, C*, B*][A, B, C, D]ejπk2N

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

DCT de type 2 utilisant un rembourrage FFT 2N (Makhoul)

[a, b, c, d]devient [a, b, c, d, 0, 0, 0, 0]. Prenez la FFT de cela pour obtenir [A, B, C, D, E, D*, C*, B*], puis jetez tout mais [A, B, C, D]et multipliez-le par2ejπk2N

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

Y = fft(y)[:N]

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

DCT de type 2 utilisant N FFT (Makhoul)

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

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

Sur ma machine, ce sont tous à peu près la même vitesse, car la génération exp(-1j*pi*k/(2*N))prend plus de temps que la FFT. :RÉ

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
endolith
la source
2
Excellente réponse, beaucoup aidé à ma mise en œuvre! Remarque supplémentaire: La dernière méthode "Type 2 DCT utilisant N FFT" fonctionne toujours correctement si la longueur du signal est impaire; le dernier élément se déplace vers l'élément du milieu. J'ai vérifié les mathématiques et le code pour ce fait.
Nayuki
1
@Nayuki Générez-vous exp(-1j*pi*k/(2*N))ou existe-t-il un raccourci vers cette étape?
endolith
Je génère exp(-1j*pi*k/(2*N))dans mon code , car un décalage d'un quart d'échantillon est nécessaire pour faire fonctionner le mappage DCT vers DFT. Qu'est-ce qui vous fait demander?
Nayuki
Salut, comment cela fonctionnerait-il pour le DCT de type III, pour calculer l'inverse du DCT-II?
Jack H
8

X(n)

laisser

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

Le DCT est alors donné par

C(k)=Re{e-jπk2NFFT{y(n)}}

Donc, vous créez essentiellement le 2N séquence de longueur y(n) où est la première moitié X(n) et la seconde moitié est X(n)renversé. Ensuite, prenez simplement la FFT et multipliez ce vecteur par une rampe de phase. Enfin, ne prenez que la partie réelle et vous avez le DCT.

Voici le code dans 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));

Éditer:

Remarque: La formule DCT utilisée est la suivante:

C(k)=2n=0N-1X(n)cos(πk2N(2n+1))

Il existe plusieurs façons de mettre à l'échelle la sommation afin qu'elle ne corresponde pas exactement aux autres implémentations. Par exemple, MATLAB utilise:

C(k)=w(k)n=0N-1X(n)cos(πk2N(2n+1))

w(0)=1N et w(1...N-1)=2N

Vous pouvez en tenir compte en mettant à l'échelle correctement la sortie.

Jason B
la source
1
y (n) est censé être de longueur N, pas de longueur 2N. C'est ainsi que vous obtenez la vitesse de calcul 4x, en calculant un DCT de longueur N à partir d'un signal de longueur N au lieu de 2N FFT à partir d'un signal 2N. fourier.eng.hmc.edu/e161/lectures/dct/node2.html www-ee.uta.edu/dip/Courses/EE5355/Discrete%20class%201.pdf
endolith
0

Pour un véritable calcul scientifique, la quantité d'utilisation de la mémoire est également importante. Par conséquent, la FFT à point N est plus attrayante pour moi. Cela n'est possible qu'en raison de la symétrie hermitienne du signal. La référence Makhoul est donnée ici. Et a en fait l'algorithme de calcul DCT et IDCT ou DCT10 et DCT01.
http://ieeexplore.ieee.org/abstract/document/1163351/

Hasbestein
la source