He estado haciendo pruebas calculado el error cuadrático medio entre x y IFFT(FFT(x)) y se obtienen valores inferiores a 1E-30 para tamaños de entrada lo suficientemente grandes (65536 o más). El algoritmo no es muy eficiente en términos de memoria ya que hace varios
new en cada llamada recursiva, pero eso facilitará la comprensión del mismo.El siguiente paso será utilizar esta implementación de la FFT para implementar el algoritmo que dejé aparcado a un lado hace algunas semanas: el de separación de fuentes de sonido mediante discriminación por acimut. El objetivo final es extraer la voz de una grabación comercial estéreo. A ver si sale. Ya iré contando mis logros y mis sinsabores :-).
Messenger holen , ver?ndern Archetypal , Landschaft Mail oder ver?ndern Sie Ihr Konto in der Umgebung Mitglied Center oder Anzeige aus . Die Band unter , dass die Anzahl der Pr?sentationen
einen offenen Geist , das Leben leben , zu lernen jordan 11 bred und zu wachsen, und folgen Sie unserer Firma Mantra, das stolz über unserer Haustür wird angezeigt: " . Keine Louis Vuitton Outlet perfekte Menschen erlaubt"
auferlegten Fristen und innerhalb von Teams , die von gew?hlten Politikern geführt werden gruppiert und Louis Vuitton Outlet Store nach einer bestimmten Aufgabe . "Wir bekommen eine erstaunliche Menge an Arbeit in
Firma umfangreiche Bibliothek ausleihen der pers?nlichen Louis Vuitton Outlet Store Entwicklung Unterrichtsmaterialien und wird ermutigt, ihre Erfahrungen mit anderen zu teilen , die F?rderung der pers?nlichen
überlegenheit keine Mühe jordan retro 4 bred F?rdern gealterte Computer-Kommunikation auf einen neuen Computer, ist keine Synchronisierung zwischen Computer und wollte nicht schwer Eintrag von überall air jordan 11 space jam in der
REWW nimmt diese Konzepte zu Herzen , wenn es um erfolgreiche Immobilien-Investitionen geht. Bedenken Sie, dass wenn space jams jordans Sie Louis Vuitton Outlet ein System der Immobilien-Investment haben , unabh?ngig davon, ob für
deville waches We own a local watches shop in the marketing and sincerely wish to cooperate with every replica watch resellers from all over the world
Si, la verdad que un muy buen resumen.
Por otra parte, he descubierto este blog hace poco y me ha encantado, sobre todo los artículos que tienen que ver con programación y audio.
Muy bueno!
Tío: Eres el puto amo :-). Muchas gracias por el algoritmo, al final con todo el rollo del problema con el servidor no pude pillar el que me dejaste aquí y cogí y adapté otro que encontré en esta web.
De todas formas voy a tomar nota del tuyo. Está muy bien comentado :-). Gracias, crack.
Pues sí se pierde.
Poner [ code ] no ayuda.
¿Y así se pierde la identación
también?
Marchando un pseudocódigo. La variante de FFT no es la típica de Cooley y Tukey, sino la de Pease, además implementa el truquito para trabajar con reales. El problema es que está en Octave, pero cambiando las multiplicaciones de las exponenciales complejas por senos y cosenos de reales y quitando los +1 de los accesos a arrays (bárbara costumbre de Octave siguiendo las pautas de MatLab de no usar el índice 0) te debería servir. El porqué estoy ahora mismo empollado en FFTs y encima uso la variante de Pease es una historia larga de contar.
function Y=peaseFftReal(X)
N = length(X);
r = log2(N);
if (rem(r, 1) != 0)
disp "El tamaño de la entrada ha de ser potencia de 2"
return;
end
% Transformada real de tamaño N con una FFT compleja de tamaño N/2
% Basicamente se usa la parte imaginaria, que sería 0 si se
% trabajase con señales reales, para albergar la mitad de la señal
N = N/2;
r=r-1;
Y = zeros(1,N);
for k=[0:N-1]
kr = bitReverse(k, r);
Y(kr+1) = X(2*k +1) + j*X(2*k+1 +1);
end;
% Fijate que a r, numero de etapas -> log2(N), ya se le resto 1
% porque Y tiene tamaño N/2
for l=[1:r]
% La variante de Pease necesita un array auxiliar
% Es out-of-place, lo calculado no sobreescribe en los mismos
% índices que lo leido.
Yaux = Y;
% Variante de Pease es bastante particular, fíjate que las mariposas
% son idénticas en todas las etapas (la variable l, número de etapa,
% no se considera al formar las direcciones a leer-escribir).
% El cálculo distingue entre las dos mitades del array de datos.
% En la primera mitad, para el dato i, se suman los datos de
% las posiciones 2i y 2i+1
for i = 0:N/2-1
Y(i+1) = Yaux(2*i +1) + Yaux(2*i+1 +1);
end
% La segunda mitad del array se calcula cogiendo los dos
% mismos números, restándolos y multiplicándolos por la
% exponencial compleja discreta que corresponda.
Nl = 2^(r-l+1);
for i = 0:N/2-1
algo = 2^(r-l)-1;
algo=bitand(i, algo);
algo=bitReverse(algo, r-l);
twiddle = e^(-2*pi*j*(algo)/Nl);
Y(i+N/2 +1) = (Yaux(2*i +1) - Yaux(2*i+1 +1))*twiddle;
end
end
% Se deja "como ejercicio para el alumno con inquitudes"
% que una vez accedidos los datos en 2i y 2i+1,
% calcule Ynuevo en i y en i+N/2. Tal como están los búcles
% evidencian la forma tan curiosa de la variante de Pease,
% pero esto fuerza a que se tengan que leer dos veces los datos.
% Y eso es todo.
% A continuación se "desenvuelve" la transformada compleja para
% restaurar la transformada real de tamaño doble
Yaux = Y;
Y(N +1) = real(Y(0 +1)) - imag(Y(0 +1));
Y(0 +1) = ((1-j*e^(-j*0))*Yaux(0 +1) + (1+j*e^(-j*0))*Yaux(0 +1)')/2;
for k=[1:N-1]
Y(k +1) = ((1-j*e^(-j*2*pi*k/(2*N)))*Yaux(k +1) +
(1+j*e^(-j*2*pi*k/(2*N)))*Yaux(N-k +1)')/2;
Y(2*N-k +1) = real(Y(k +1)) -j*imag(Y(k +1));
end
N = 2*N;
r=r+1;
return;
Saludos, de un crack de mentirijilla (que los viernes trabaja poco y bloguea mucho), a un crack de verdad,
Jose.
Joder, tio. Ojalá todos los comentarios fuesen así :-). Lo cierto es que buscaba más bien una implementación más pedagógica que otra cosa. Quiero implementar la separación de fuentes de audio por acimut y el algoritmo que tenía era el que tuvimos que hacer en Señales de 2º (¿recuerdas? :-) ), que no me iba muy bien. Partiendo de esta implementación en Java puedo optimizarlo de forma más sencilla: lo que no me gustaba de muchas implementaciones que vi por ahí era que estaban muy optimizadas, iban rapidísimo, pero eran un poco difíciles de coger o de adaptar (sólo puedes limitarte a copiar y pegar). Gracias por la referencia a la transformada Hartley, miraré a ver si me sirve. En cuanto a lo que comentas sobre la FFTW, está claro que es la mejor opción, pero bueno, me apetecía implementarla yo mismo... gajes del masoquismo... ;-).
Gracias Jose. Eres un crack.
Por comentar algo, un par de apuntes perfectamente ignorables.
1a) Todas las implementaciones sencillas son ineficientes.
1b) Las implementaciones eficientes son ininteligibles.
1c) Corolario a los epígrafes a y b: Si quieres ganar en eficiencia, emplea código ajeno, claro que llamar a una API no es demasiado instructivo.
2) Si lo que vas a implementar corre sobre un ordenador (da igual que sea PC, Mac, Sparc, ...): usa la librería FFTW del M.I.T. (www.fftw.org; o incluso 'apt-get install fftw'). Sobre estas plataformas la ganancia se logra en base a subdividir el problema en subproblemas que presenten los tamaños más eficientes para las jerarquías de cache de la máquina en cuestión: ej: FFT(n=8) -> FFT(2)yFFT(4) ó FFT(4)yFFT(2) ó FFT(8). Y eso es lo que hace ésta librería, para un cierto tamaño del problema corre las posibles descomposiciones y determina cúal es la más eficiente, lo que queda registrado para posteriores ejecuciones de FFTs de ese tamaño (es lo que denomina 'wisdom'). La manera de descomponer las matrices de la DFT de forma automática para los distintos tamaños entra directamente en el terreno de lo indescifrable (álgebra tensorial, productos de Kronecker, ...).
3) Si vas a desarrollar sobre DSPs o FPGAs, hay un proyecto denominado SPIRAL: www.spiral.net, que directamente te da el código que has de compilar para tu arquitectura y tamaño de problema concreto.
4) Si vas a desarrollar para la GBA: n.p.i., pero trata de buscar FFTs de código abierto para el tipo de operandos adecuado, me imagino que coma fija de 16 bits, o bizarreces similares.
5) La precisión que obtienes ahora mismo es seguramente innecesaria y te vaya más rápido si te pasas a simple precisión (debes de estar trabajando con doubles).
6) Si vas a funcionar con datos reales, que es lo más probable, busca cómo realizar FFTs de reales (transformada Hartley: FHT, o los típicos trucos: cómo calcular dos FFTs de 2N puntos con una FFT compleja de tamaño N, o una FFT real de tamaño N con una FFT compleja de tamaño N/2, creo que tienes un capítulo al respecto en el Proakis, Manolakis).
7) Estuve echando un ojo al enlace que dejas de Princeton y es directamente criminal, usa la variante Cooley y Tukey que si tiene alguna ventaja es precisamente que se puede realizar 'in-place', esto es sin consumir absolutamente ninguna memoria extra y la implementan recursiva y duplicando la memoria en cada pasada para evitar hacer un bucle for: Princeton pa'septiembre.
8) Hay disponile pseudo-código, o incluso código en prácticamente cualquier lenguaje: casi la única recomendación es precisamente que no sea recursivo (los matemáticos por no trabajar usan esta forma de expresar el algoritmo, pero tú que has hecho un compilador y sabes lo que es gestionar la pila ...).
Saludos,
Jose.
Lo sentimos. No se permiten nuevos comentarios después de 90 días.
Calendario
