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).
Réponses:
N
x
k
arange(N)
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é: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]
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 parDCT 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*]
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Éla source
exp(-1j*pi*k/(2*N))
ou existe-t-il un raccourci vers cette étape?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?laisser
Le DCT est alors donné par
Donc, vous créez essentiellement le2 N 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.
Éditer:
Remarque: La formule DCT utilisée est la suivante:
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:
oùw ( 0 ) = 1N--√ et w ( 1 ... N- 1 ) = 2N--√
Vous pouvez en tenir compte en mettant à l'échelle correctement la sortie.
la source
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/
la source