Autor Tema: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.  (Leído 24861 veces)

0 Usuarios y 2 Visitantes están viendo este tema.

Desconectado cerebro

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 735
Como indica el titulo estoy tratando de de implementar una rutina que encontré que calcula la FFT en 16 puntos utilizando 10 multiplicaciones y 79 sumas, ya se que la capacidad del pic no es la mejor y que vienen librerías para FFT de la gama alta pero me pareció un buen desafió e interesantes para pequeños proyectos(de ahí lo de inventar la rueda  :mrgreen:).
 El problema es que no puedo ver la salida del micro  :(, no se si estoy utilizando mal la función printf o si es problema del cálculo con floats. Les dejo el código hasta el momento...

Código: C
  1. #include <18f4550.h>
  2. #fuses XTPLL,MCLR,NOWDT,NOPROTECT,NOLVP,NODEBUG,USBDIV,PLL1,CPUDIV1,NOVREGEN
  3. #use delay(clock=48000000)
  4.  
  5. #use rs232(baud=9600,xmit=PIN_C6,rcv=PIN_C7)
  6.  
  7.  
  8.  
  9. #define SIN_2PI_16 0.38268343236508978
  10. #define SIN_4PI_16 0.707106781186547460
  11. #define SIN_6PI_16 0.923879532511286740
  12. #define C_P_S_2PI_16 1.30656296487637660
  13. #define C_M_S_2PI_16 0.54119610014619690
  14. #define C_P_S_6PI_16 1.3065629648763766
  15. #define C_M_S_6PI_16 -0.54119610014619690
  16.  
  17. /* INPUT: float input[16], float output[16] */
  18. /* OUTPUT: none */
  19. /* EFFECTS:  Places the 16 point fft of input in output in a strange */
  20. /* order using 10 real multiplies and 79 real adds. */
  21. /* Re{F[0]}= out0 */
  22. /* Im{F[0]}= 0 */
  23. /* Re{F[1]}= out8 */
  24. /* Im{F[1]}= out12 */
  25. /* Re{F[2]}= out4 */
  26. /* Im{F[2]}= -out6 */
  27. /* Re{F[3]}= out11 */
  28. /* Im{F[3]}= -out15 */
  29. /* Re{F[4]}= out2 */
  30. /* Im{F[4]}= -out3 */
  31. /* Re{F[5]}= out10 */
  32. /* Im{F[5]}= out14 */
  33. /* Re{F[6]}= out5 */
  34. /* Im{F[6]}= -out7 */
  35. /* Re{F[7]}= out9 */
  36. /* Im{F[7]}= -out13 */
  37. /* Re{F[8]}= out1 */
  38. /* Im{F[8]}=0 */
  39. /* F[9] through F[15] can be found by using the formula */
  40. /* Re{F[n]}=Re{F[(16-n)mod16]} and Im{F[n]}= -Im{F[(16-n)mod16]} */
  41.  
  42. /* Note using temporary variables to store intermediate computations */
  43. /* in the butterflies might speed things up.  When the current version */
  44. /* needs to compute a=a+b, and b=a-b, I do a=a+b followed by b=a-b-b.  */
  45. /* So practically everything is done in place, but the number of adds */
  46. /* can be reduced by doinc c=a+b followed by b=a-b. */
  47.  
  48. /* The algorithm behind this program is to find F[2k] and F[4k+1] */
  49. /* seperately.  To find F[2k] we take the 8 point Real FFT of x[n]+x[n+8] */
  50. /* for n from 0 to 7.  To find F[4k+1] we take the 4 point Complex FFT of */
  51. /* exp(-2*pi*j*n/16)*{x[n] - x[n+8] + j(x[n+12]-x[n+4])} for n from 0 to 3.*/
  52.  
  53.  
  54.   float data[16];
  55.   float output[16];
  56.   float zero;
  57. //  int k;
  58.  
  59.  
  60. void R16SRFFT(float input[16]) {
  61.  
  62.   float temp, out0, out1, out2, out3, out4, out5, out6, out7, out8;
  63.   float out9,out10,out11,out12,out13,out14,out15;
  64.  
  65.   out0=input[0]+input[8]; /* output[0 through 7] is the data that we */
  66.   out1=input[1]+input[9]; /* take the 8 point real FFT of. */
  67.   out2=input[2]+input[10];
  68.   out3=input[3]+input[11];
  69.   out4=input[4]+input[12];
  70.   out5=input[5]+input[13];
  71.   out6=input[6]+input[14];
  72.   out7=input[7]+input[15];
  73.  
  74.  
  75.  
  76.   out8=input[0]-input[8];   /* inputs 8,9,10,11 are */
  77.   out9=input[1]-input[9];   /* the Real part of the */
  78.   out10=input[2]-input[10]; /* 4 point Complex FFT inputs.*/
  79.   out11=input[3]-input[11];
  80.   out12=input[12]-input[4]; /* outputs 12,13,14,15 are */
  81.   out13=input[13]-input[5]; /* the Imaginary pars of  */
  82.   out14=input[14]-input[6]; /* the 4 point Complex FFT inputs.*/
  83.   out15=input[15]-input[7];
  84.  
  85.   /*First we do the "twiddle factor" multiplies for the 4 point CFFT */
  86.   /*Note that we use the following handy trick for doing a complex */
  87.   /*multiply:  (e+jf)=(a+jb)*(c+jd) */
  88.   /*   e=(a-b)*d + a*(c-d)   and    f=(a-b)*d + b*(c+d)  */
  89.  
  90.   /* C_M_S_2PI/16=cos(2pi/16)-sin(2pi/16) when replaced by macroexpansion */
  91.   /* C_P_S_2PI/16=cos(2pi/16)+sin(2pi/16) when replaced by macroexpansion */
  92.   /* (SIN_2PI_16)=sin(2pi/16) when replaced by macroexpansion */
  93.   temp=(out13-out9)*(SIN_2PI_16);
  94.   out9=out9*(C_P_S_2PI_16)+temp;
  95.   out13=out13*(C_M_S_2PI_16)+temp;
  96.  
  97.   out14=out14*(SIN_4PI_16);
  98.   out10=out10*(SIN_4PI_16);
  99.   out14=out14-out10;
  100.   out10=out14+out10+out10;
  101.  
  102.   temp=(out15-out11)*(SIN_6PI_16);
  103.   out11=out11*(C_P_S_6PI_16)+temp;
  104.   out15=out15*(C_M_S_6PI_16)+temp;
  105.  
  106.   /* The following are the first set of two point butterfiles */
  107.   /* for the 4 point CFFT */
  108.  
  109.   out8=out8+out10;
  110.   out10=out8-out10-out10;
  111.  
  112.   out12=out12+out14;
  113.   out14=out12-out14-out14;
  114.  
  115.   out9=out9+out11;
  116.   out11=out9-out11-out11;
  117.  
  118.   out13=out13+out15;
  119.   out15=out13-out15-out15;
  120.  
  121.   /*The followin are the final set of two point butterflies */
  122.   output[1]=out8+out9;
  123.   output[7]=out8-out9;
  124.  
  125.   output[9]=out12+out13;
  126.   output[15]=out13-out12;
  127.  
  128.   output[5]=out10+out15;        /* implicit multiplies by */
  129.   output[13]=out14-out11;        /* a twiddle factor of -j */                            
  130.   output[3]=out10-out15;  /* implicit multiplies by */
  131.   output[11]=output[11]-out14-out11;  /* a twiddle factor of -j */
  132.  
  133.  
  134.   /* What follows is the 8-point FFT of points output[0-7] */
  135.   /* This 8-point FFT is basically a Decimation in Frequency FFT */
  136.   /* where we take advantage of the fact that the initial data is real*/
  137.  
  138.   /* First set of 2-point butterflies */
  139.    
  140.   out0=out0+out4;
  141.   out4=out0-out4-out4;
  142.   out1=out1+out5;
  143.   out5=out1-out5-out5;
  144.   out2=out2+out6;
  145.   out6=out2-out6-out6;
  146.   out3=out3+out7;
  147.   out7=out3-out7-out7;
  148.  
  149.   /* Computations to find X[0], X[4], X[6] */
  150.  
  151.   output[0]=out0+out2;
  152.   output[4]=out0-out2;
  153.   out1=out1+out3;
  154.   output[12]=out3+out3-out1;
  155.  
  156.   output[0]=output[0]+out1;  /* Real Part of X[0] */
  157.   output[8]=output[0]-out1-out1; /*Real Part of X[4] */
  158.   /* out2 = Real Part of X[6] */
  159.   /* out3 = Imag Part of X[6] */
  160.  
  161.   /* Computations to find X[5], X[7] */
  162.  
  163.   out5=out5*SIN_4PI_16;
  164.   out7=out7*SIN_4PI_16;
  165.   out5=out5-out7;
  166.   out7=out5+out7+out7;
  167.  
  168.   output[14]=out6-out7; /* Imag Part of X[5] */
  169.   output[2]=out5+out4; /* Real Part of X[7] */
  170.   output[6]=out4-out5; /*Real Part of X[5] */
  171.   output[10]=-out7-out6; /* Imag Part of X[7] */
  172.  
  173.  
  174.   printf(" %2.2f %2.2f\n\r",output[0],zero);
  175.   printf(" %2.2f %2.2f\n\r",output[1],output[9]);
  176.   printf(" %2.2f %2.2f\n\r",output[2],output[10]);
  177.   printf(" %2.2f %2.2f\n\r",output[3],output[11]);
  178.   printf(" %2.2f %2.2f\n\r",output[4],output[12]);
  179.   printf(" %2.2f %2.2f\n\r",output[5],output[13]);
  180.   printf(" %2.2f %2.2f\n\r",output[6],output[14]);
  181.   printf(" %2.2f %2.2f\n\r",output[7],output[15]);
  182.   printf(" %2.2f %2.2f\n\r",output[8],zero);
  183.   printf(" %2.2f %2.2f\n\r",output[7],-output[15]);
  184.   printf(" %2.2f %2.2f\n\r",output[6],-output[14]);
  185.   printf(" %2.2f %2.2f\n\r",output[5],-output[13]);
  186.   printf(" %2.2f %2.2f\n\r",output[4],-output[12]);
  187.   printf(" %2.2f %2.2f\n\r",output[3],-output[11]);
  188.   printf(" %2.2f %2.2f\n\r",output[2],-output[9]);
  189.   printf(" %2.2f %2.2f\n\r",output[1],-output[8]);
  190. }
  191.  
  192.  
  193.  
  194.  
  195.  
  196. void main() {
  197.  
  198. //data 1 23 15 8 9 1 2 3 4 5 2 4 6 8 0
  199.  
  200.    zero=0;
  201.    data[0]=1;
  202.    data[1]=23;
  203.    data[2]=15;
  204.    data[3]=8;
  205.    data[4]=9;
  206.    data[5]=1;
  207.    data[6]=2;
  208.    data[7]=3;
  209.    data[8]=4;
  210.    data[9]=5;
  211.    data[10]=2;
  212.    data[11]=4;
  213.    data[12]=6;
  214.    data[13]=8;
  215.    data[14]=0;
  216.    data[15]=0;
  217.    
  218. //----------------------------------
  219.  
  220.    disable_interrupts(global);
  221.    disable_interrupts(int_timer1);
  222.    disable_interrupts(int_rda);
  223.    disable_interrupts(int_ext);
  224.    disable_interrupts(int_ext1);
  225.    disable_interrupts(int_ext2);
  226.    setup_spi(FALSE);
  227.    setup_psp(PSP_DISABLED);
  228.    setup_comparator(NC_NC_NC_NC);
  229.    setup_vref(FALSE);
  230.    port_b_pullups(FALSE);    
  231.  
  232. //---------------------------  
  233.  
  234.  
  235. for(;;){
  236.    
  237.  
  238.  R16SRFFT(data);
  239.  
  240.  delay_ms(500);
  241.   }
  242. }
  243.  
  244.  
  245.  
  246. /*
  247. 91 + 0i
  248. 26,358765205042054 - 18,59438622703722i
  249. -2,998370275568188 - 28,31592789979983i  
  250. -18,309016994374947 - 16,616988753597195i  
  251. -10,178425317543107 - 23,190667119105044i  
  252. -12,5 - 18,186533479473212i  
  253. -17,190983005625053 - 5,012552719206482i  
  254. -3,181969611930761 - 0,387281014208582i  
  255. -3,181969611930761 + 0,387281014208582i  
  256. -17,190983005625053 + 5,012552719206482i  
  257. -12,5 + 18,186533479473212i  
  258. -10,178425317543107 + 23,190667119105044i  
  259. -18,309016994374947 + 16,616988753597195i  
  260. -2,998370275568188 + 28,31592789979983i  
  261. 26,358765205042054 + 18,59438622703722i
  262. */


En la rutina ingresa un vector como dato y debería obtener la transformada en los puntos, en el caso que estoy planteando la transformada es la que esta al final del código para comprobar.
Les adjunto código, esquema y .cof.

Saludos!  :)
« Última modificación: 16 de Julio de 2009, 13:50:21 por cerebro »
LAS MALVINAS SON ARGENTINAS!

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #1 en: 16 de Julio de 2009, 00:52:19 »
Hola

No sabía que había forma de aproximar la FFT con operaciones más simples.

Intenta asignar los valores de data[] con números flotantes y no enteros.

Código: [Seleccionar]
data[2]=23.0;

Desconectado cerebro

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 735
Re: problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #2 en: 16 de Julio de 2009, 11:34:16 »
ya lo logre solucionar  :lol: me había olvidado de cambiar el clock a 48 Mhz en el proteus  :oops:, aunque ahora voy a tener que revisar la formula porque hay valores que ni se acercan a los que corresponde.
Gracias mig igual  :wink:
LAS MALVINAS SON ARGENTINAS!

Desconectado cerebro

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 735
Re: SOLUCIONADO- problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #3 en: 16 de Julio de 2009, 13:12:24 »
 :-/ ya lo conseguí
estaba mal calculada la transformada!
estos son los verdaderos valores
 
Código: [Seleccionar]
91.0000
22.8459 -18.8713i
-2.9289 -34.7990i
-15.2006 -28.1561i
1.0000 -22.0000i
-6.3558 -12.9429i
-17.0711 - 4.7990i
-13.2895 + 8.3419i
-13.0000
-13.2895 - 8.3419i
-17.0711 + 4.7990i
-6.3558 +12.9429i
1.0000 +22.0000i
-15.2006 +28.1561i
-2.9289 +34.7990i
22.8459 +18.8713i
y esta es la salida del micro


No sabía que había forma de aproximar la FFT con operaciones más simples.

el algoritmo implementado se llama Split-radix FFT.
 Ahora a jugar!  :D
LAS MALVINAS SON ARGENTINAS!

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #4 en: 16 de Julio de 2009, 13:15:56 »
Perfecto Señor Cerebro. Un analizador de espectro de audio con 16 frecuencias ahora es posible y solo con un PIC.  8)

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #5 en: 18 de Julio de 2009, 14:54:24 »
Una pregunta Cerebro, ¿de dónde tomaste el código fuente para la FFT Split radix? ¿O acaso tú lo escribiste?

Desconectado cerebro

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 735
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #6 en: 18 de Julio de 2009, 15:21:30 »
 :(  yo no lo hice, mucho desgano y esas rutinas ya están hechas para que re-realizar algo. Lo más difícil es encontrar algo confiable  :mrgreen: por eso estuve comprobando el algoritmo.
Te dejo la pagina del autor   
FFT

Saludos!

PD: será que estas por construir el analizador de espectros  :wink: ?
LAS MALVINAS SON ARGENTINAS!

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #7 en: 18 de Julio de 2009, 19:51:09 »
Gracias por el link.

Quiero construirlo y no es complicado una vez viendo que ya te funcionó. Solo hay que ajustar una frecuencia de muestreo, poner un filtro antialias, muestrear 16 puntos (que es lo que calcula la FFT que pones) y manos a la obra.

El problema es que no tengo generador de funciones, ni micrófono conmigo pero una vez regresando a clases... la split-radix fft será mi amiga  :D

Desconectado cerebro

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 735
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #8 en: 18 de Julio de 2009, 21:52:33 »
 Si lo haces en una tarde  :), si querés ayuda con cualquier cosa avisame.
 Che se puede extender a más puntos (siempre y cuando sea multiplo de 4) pero no se como lo va a llevar el micro y una cosa más con proteus podes simular muchas partes (como para ir programando si querés) incluso mandarle al micro una señal wav.
  Saludos!!   

 
LAS MALVINAS SON ARGENTINAS!

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #9 en: 18 de Julio de 2009, 21:53:55 »
OK, voy a darle una revisada a esto porque la única FFT que he usado es la de Matlab.

Cualquier cosa ya te avisaré  :mrgreen:

Desconectado cerebro

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 735
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #10 en: 23 de Julio de 2009, 01:03:26 »
bueno tampoco la vas a extrañar tanto porque la split es una variante del algoritmo Cooley-Tukey FFT  :) que es el que utiliza matlab para su calculo......
 
 hoy estaba pensando (raro raro...) en los viejos teléfonos, supuestamente transformaban la señal de señal de la voz y enviaban los coeficientes de la transformada de esta manera ahorraban en el envio de datos (según me comentaron...) capaz los que posean mas antigüedad conozcan del tema y nos comenten. Pero si es así, se podría armar un teléfono prehistórico con un pic  :shock: jajaj estaría buenísimo para investigar el teléfono TODOPIC  :D 
LAS MALVINAS SON ARGENTINAS!

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #11 en: 23 de Julio de 2009, 01:33:40 »
Pero no veo el ahorro en datos.

Si calculas la FFT de un arreglo de 16 muestras, obtendrás un arreglo de 16 muestras reales y 16 imaginarias, de las cuales se usan 8 y 8 ya que las demás son su espejo, siendo a final de cuentas 16 números... misma cantidad que la original.

Desconectado cerebro

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 735
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #12 en: 23 de Julio de 2009, 15:47:38 »
parece que somos 2 los interesados en el tema  :D, de teléfonos se muy poco solo levanto el tubo y marco  :mrgreen:. Actualmente no tengo idea como harán la comunicación.  Pero supongamos que yo envio la señal acústica tal cual llega, el espectro del habla alcanza los 20 khz, de manera que tendría que estar mandando datos por lo menos a 40 kHz.
 A grandes rasgos, por la característica de resonancia de las vias superiores en el espetro del habla hay como montañas que se denominan formantes casi siempre son 3 (F0,F1,F2) presentes por ejemplo en la "a" "e" ....  y la "s" no presenta formantes, entonces los primeros implantes cocleares extraían estas 3 frecuencias y estimulaban en el lugar correspondiente a esa frecuencia en pacientes sordos. La persona escuchaba medio robotico pero escuchaba algo, actualmente se usan muchas estrategias y muchas bandas como 20 o más.   
 Si no mezcle mucho las cosas, la idea es extraer algunas bandas del espectro y enviar los coeficientes a una velocidad muy inferior por esta razón cae la cantidad de datos que tengo que enviar a costo de perder calidad en el sonido.........
 Tiene que funcionar  :lol:...................................bueno son cosas locas que se me ocurren de vez en cuando  :oops:
LAS MALVINAS SON ARGENTINAS!

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #13 en: 23 de Julio de 2009, 16:00:31 »
Un detalle... el espectro de frecuencias de la voz humana llega a los 4kHz. Cualquier otro ruido ya no es voz (un grito extremo o cualquier otra cosa). Por lo que hay que muestrear a 8kHz.

Para que un sordo escuche bien, sí hay que usar un espectro de 20Hz a 20kHz, pero si tu aplicación se orienta a transmitir voz entonces 8kHz es una buena velocidad de muestreo. De ahí que la codificación PCM @ 8 bits usa una velocidad de datos de 64kbps.

Aparte de PCM antes se usaban técnicas como la modulación delta.
http://en.wikipedia.org/wiki/Delta_modulation

Sobre lo que comentas de bajar la velocidad de transmisión lo entiendo así:

- Muestreas una señal, por ejemplo... 32 muestras (dominio del tiempo).
- Calculas la FFT de esas muestras obteniendo 32 muestras complejas (r + i)
- Usas solo las primeras 16 porque las siguientes son espejos
- Decimas esas 16 complejas a la mitad y obtienes 8 complejas
- Entonces ya habrás reducido de 32 muestras (dominio del tiempo) a 8 muestras complejas (16 números en total)

Pero creo que es más fácil muestrear 16 desde el principio sin decimarlas al final, ya después puedes interpolarlas para subirlas a 32.

Bueno, quién sabe... no sé si te entendí  :D

Desconectado pablomanieri

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 639
Re: SOLUCIONADO-problemas con FFT con 18f4550, queriendo inventar la rueda.
« Respuesta #14 en: 25 de Julio de 2009, 13:23:21 »
Me interesa mucho esto de calcular la fft con un pic. Lo voy a seguir de cerca


 

anything