Transformée de Fourier sur ATmega32 et Arduino Mega2560

La rapidité d'exécution du calcul est essentielle lorsqu'il s'agit de traiter un signal en temps réel. Elle déterminera la fréquence maximale qui puisse être analysée.

1 La partie théorique

Je vous invite dans un premier temps, si ce n'est déjà fait, à consulter l'étude analytique de la Transformée de Fourier Rapide sur mes pages pleines de jolies équations mathématiques commentées :


2 Le programme en C++ et Qt4

La page principale :

CODE SOURCE C++ Qt4
  1. /*
  2.   Programme écrit par Silicium628
  3.   ce logiciel est libre et open source
  4.  
  5.   */
  6. #include "mainwindow.h"
  7. #include "ui_mainwindow.h"
  8. #include "math.h"
  9. #include "complexe.cpp"
  10.  
  11. QString version = "2.4";
  12.  
  13. QColor couleur_ecran = QColor::fromRgb(3, 41, 11, 255);
  14. QColor couleur_ligne = QColor::fromRgb(167, 2, 0, 255);
  15. QColor couleur_trace = QColor::fromRgb(0, 210, 0, 255);
  16. QColor couleur_texte = QColor::fromRgb(255, 255, 0, 255);
  17. QColor couleur_signature = QColor::fromRgb(255, 255, 0, 255);
  18.  
  19. QPen pen_ligne(couleur_ligne, 1, Qt::SolidLine);
  20. QPen pen_trace(couleur_trace, 1, Qt::SolidLine);
  21. QPen pen_reticule(couleur_ligne, 1, Qt::SolidLine);
  22.  
  23. Complexe ech[1024]; // nb_ech echantillons
  24. Complexe tab_X[1024]; // nb_ech valeurs traitées
  25. Complexe tab_W[512]; // nb_ech/2
  26.  
  27.  
  28. float frequence1=400; // Hz
  29. float frequence2=400; // Hz
  30. float frequence3=100; // Hz
  31. uint nb_etapes=8;
  32. uint nb_ech = pow (2,nb_etapes); // nombre d'échantillons = 2 puissance(nb_etapes)
  33. float frq_echantillonnage=16*frequence1; // Hz
  34. float Te; // période d'échantillonage
  35.  
  36. //int pas_balayage=40;
  37. int hamming = false;
  38. int second_Frq = false;
  39. int bz = false;
  40. int modul = false;
  41.  
  42. MainWindow::MainWindow(QWidget *parent) :
  43. QMainWindow(parent)
  44.  
  45. {
  46. setupUi(this);
  47. setWindowTitle("Transformee de Fourier Rapide FFT - version " + version);
  48.  
  49. scene = new QGraphicsScene(this);
  50. scene->setBackgroundBrush(couleur_ecran);
  51. graphicsView1->setScene(scene);
  52. graphicsView1->setGeometry(0,0,530,430);
  53.  
  54. groupe_reticule = new QGraphicsItemGroup();
  55. scene->addItem(groupe_reticule);
  56. groupe_trace = new QGraphicsItemGroup();
  57. scene->addItem(groupe_trace);
  58.  
  59. calcul_tableau_W();
  60. tracer_graduations();
  61. }
  62.  
  63.  
  64.  
  65.  
  66. void MainWindow::effacer_graduation()
  67. {
  68. foreach( QGraphicsItem *item, scene->items( groupe_reticule->boundingRect() ) )
  69. {
  70. if( item->group() == groupe_reticule ) { delete item; }
  71. }
  72. }
  73.  
  74.  
  75. void MainWindow::effacer_trace()
  76. {
  77. foreach( QGraphicsItem *item, scene->items( groupe_trace->boundingRect() ) )
  78. {
  79. if( item->group() == groupe_trace ) { delete item; }
  80. }
  81. }
  82.  
  83.  
  84. MainWindow::~MainWindow()
  85. {
  86.  
  87. }
  88.  
  89.  
  90. void MainWindow::tracer_graduations()
  91. {
  92. float i,x,y;
  93. float nb_grad_max;
  94. float intervalle; // séparant les graduations
  95. float fi;
  96. QString sti;
  97. uint decal = 3; // leger decallage vertical pour ne pas masquer la trace
  98.  
  99. rectangle = new QGraphicsRectItem(0,decal, 511, 400);
  100. rectangle->setPen(couleur_ligne);
  101. groupe_trace->addToGroup(rectangle);
  102.  
  103. // lignes verticales
  104. // ici calcul de l'intervalle (en pixels) qui correspond exactement à 100Hz (en fréquence)
  105.  
  106. intervalle = 2000.0 * nb_ech * 2.0 / (frq_echantillonnage);
  107. nb_grad_max = 512 / intervalle;
  108.  
  109. for (i=0; i<=nb_grad_max; i++)
  110. {
  111. x = intervalle * i;
  112.  
  113. ligne1 = new QGraphicsLineItem(x, decal, x, 400 + decal);
  114. ligne1->setPen(pen_reticule);
  115. groupe_reticule->addToGroup(ligne1);
  116.  
  117. fi = x * frq_echantillonnage / nb_ech /4;
  118.  
  119. sti.setNum(fi);
  120. texte_frq = new QGraphicsTextItem(sti);
  121. texte_frq->setDefaultTextColor(couleur_texte);
  122. texte_frq->setPos(x,380);
  123. if(x<480) // évite que l'écriture ne déborde du cadre a droite
  124. {
  125. groupe_reticule->addToGroup(texte_frq);
  126. }
  127. }
  128.  
  129. // lignes horizontales
  130. intervalle = 50;
  131. nb_grad_max = 400 / intervalle;
  132.  
  133. for (i=0; i<=nb_grad_max; i++)
  134. {
  135. y = 400 + decal - intervalle * i;
  136.  
  137. ligne1 = new QGraphicsLineItem(0, y, 511, y);
  138. ligne1->setPen(pen_reticule);
  139. groupe_reticule->addToGroup(ligne1);
  140. }
  141.  
  142. texte_frq = new QGraphicsTextItem("Silicium628");
  143. texte_frq->setDefaultTextColor(couleur_signature);
  144. texte_frq->setPos(5,5);
  145. groupe_reticule->addToGroup(texte_frq);
  146.  
  147. scene->addItem(groupe_reticule);
  148. }
  149.  
  150.  
  151. void MainWindow::tracer_signal()
  152. {
  153. uint offset_y = 100;
  154. float echelle =1024/nb_ech;
  155. uint n;
  156. float x, y, memo_x, memo_y;
  157.  
  158. segment_trace = new QGraphicsLineItem(0, offset_y, 512, offset_y);
  159. segment_trace->setPen(couleur_ligne);
  160. groupe_trace->addToGroup(segment_trace);
  161.  
  162. effacer_trace();
  163.  
  164. x=0;
  165. y=offset_y;
  166. for(n=0; n<nb_ech; n++)
  167. {
  168. memo_x = x;
  169. memo_y = y;
  170. x = echelle * n /2;
  171.  
  172. y = offset_y - 20 * ech[n].a;
  173.  
  174. segment_trace = new QGraphicsLineItem(memo_x ,memo_y, x, y);
  175. segment_trace->setPen(pen_trace);
  176. groupe_trace->addToGroup(segment_trace);
  177. }
  178. scene->addItem(groupe_trace);
  179. }
  180.  
  181.  
  182.  
  183. void MainWindow::tracer_tableau_valeurs()
  184. {
  185. uint offset_y = 350; // 600
  186. float echelle =1024/nb_ech;
  187. uint n;
  188. float x, y, memo_x, memo_y;
  189.  
  190. int z;
  191.  
  192. if (bz) {z=4;} else {z=1;}
  193.  
  194. segment_trace = new QGraphicsLineItem(0, offset_y, 512, offset_y);
  195. segment_trace->setPen(couleur_ligne);
  196. groupe_trace->addToGroup(segment_trace);
  197.  
  198.  
  199.  
  200. x=0;
  201. y=offset_y;
  202. for(n=0; n<nb_ech/2; n++)
  203. {
  204. memo_x = x;
  205. memo_y = y;
  206. x = echelle * n;
  207.  
  208. y = offset_y - z * sqrt( tab_X[n].a * tab_X[n].a + tab_X[n].b * tab_X[n].b );
  209.  
  210. segment_trace = new QGraphicsLineItem(memo_x ,memo_y, x, y);
  211. segment_trace->setPen(pen_trace);
  212. groupe_trace->addToGroup(segment_trace);
  213. }
  214. scene->addItem(groupe_trace);
  215. }
  216.  
  217.  
  218.  
  219. void MainWindow::RAZ_tableau_echantillons()
  220. {
  221. uint n;
  222. for(n=0; n < nb_ech; n++)
  223. {
  224. ech[n].a = 0; // partie reelle
  225. ech[n].b = 0; // partie imaginaire
  226. }
  227. }
  228.  
  229.  
  230.  
  231.  
  232. void MainWindow::remplir_tableau_echantillons()
  233. {
  234. // création du signal à analyser
  235. uint n, n_max;
  236. float x;
  237.  
  238. Te = 1.0 / frq_echantillonnage;
  239.  
  240. if (bz==true) {n_max = nb_ech /4; } else {n_max = nb_ech; }
  241.  
  242. for(n=0; n < n_max; n++) // nb d'échantillons
  243. {
  244. x= cos(2.0 * M_PI * frequence1 * n * Te );
  245.  
  246. if (modul) {x *= 1 + 0.8 * cos(2 * M_PI * frequence3 * n * Te);} // modulation d'amplitude
  247.  
  248. if (second_Frq) {x+=cos(2 * M_PI * frequence2 * n * Te); }
  249. if (hamming) { x = x * (1- cos (2* M_PI * n / n_max ));}
  250. ech[n].a = x; // partie reelle
  251. ech[n].b = 0; // partie imaginaire
  252. }
  253. }
  254.  
  255.  
  256. //void MainWindow::remplir_tableau_echantillons()
  257. //{
  258. //// création du signal à analyser
  259. // uint n;
  260. // float x;
  261.  
  262. // Te = 1.0 / frq_echantillonnage;
  263.  
  264. // for(n=0; n < nb_ech; n++) // nb d'échantillons
  265. // {
  266. // if ((n >= 50) and (n<=55)) {ech[n].a = 30;} else {ech[n].a = 0;}
  267.  
  268. // }
  269. //}
  270.  
  271.  
  272. uint bit_reversal(uint num, uint nb_bits)
  273. {
  274. uint r = 0, i, s;
  275. if ( num > (1<< nb_bits)) { return 0; }
  276.  
  277. for (i=0; i<nb_bits; i++)
  278. {
  279. s = (num & (1 << i));
  280. if(s) { r |= (1 << ((nb_bits - 1) - i)); }
  281. }
  282. return r;
  283. }
  284.  
  285.  
  286.  
  287. void MainWindow::bit_reverse_tableau_X()
  288. {
  289.  
  290. // recopie les échantillons en les classant dans l'ordre 'bit reverse'
  291. uint n,r;
  292.  
  293. for(n=0; n < nb_ech; n++) // nb d'échantillons
  294. {
  295. r=bit_reversal(n,nb_etapes);
  296. tab_X[n] = ech[r];
  297. }
  298. }
  299.  
  300.  
  301.  
  302. void MainWindow::calcul_tableau_W()
  303. {
  304. // calcul et memorisation dans un tableau des twiddle factors
  305. uint n;
  306. float x;
  307.  
  308. for(n=0; n<(nb_ech/2-1); n++)
  309. {
  310. x=2.0*M_PI * n / nb_ech;
  311.  
  312. tab_W[n].a = cos(x); // partie reelle
  313. tab_W[n].b = -sin(x); // partie imaginaire
  314. }
  315. }
  316.  
  317.  
  318.  
  319. void MainWindow::calcul_fourier()
  320. {
  321. Complexe produit; // voir la classe "Complexe" : complexe.h et complexe.ccp
  322.  
  323. uint etape, e1, e2, li, w, ci;
  324. uint li_max=nb_ech;
  325.  
  326. e2=1;
  327.  
  328. for (etape=1; etape<=nb_etapes; etape++)
  329. {
  330. e1=e2; //(e1 evite d'effectuer des divisions e2/2 plus bas)
  331. e2=e2+e2;
  332.  
  333. for (li=0; li<li_max; li+=1)
  334. {
  335. ci=li & (e2-1); // ET bit à bit
  336. if (ci>(e1-1))
  337. {
  338. w=li_max/e2*(li & (e1 -1)); // ET bit à bit calcul du numéro du facteur de tripatouillages W
  339.  
  340. produit = tab_W[w] * tab_X[li]; // le twiddle factor est lu en memoire; le produit est une mutiplication de nb complexes. Voir "complexe.cpp"
  341.  
  342. tab_X[li]=tab_X[li-e1]-produit; // concerne la ligne basse du croisillon; soustraction complexe; Voir "complexe.cpp"
  343. tab_X[li-e1]=tab_X[li-e1] + produit; // concerne la ligne haute du croisillon; addition complexe; Voir "complexe.cpp"
  344. }
  345. }
  346. }
  347. }
  348.  
  349.  
  350.  
  351. void MainWindow::on_btn2_clicked() // BOUTON "FOURIER"
  352. {
  353. effacer_trace();
  354. calcul_tableau_W();
  355. RAZ_tableau_echantillons();
  356. remplir_tableau_echantillons();
  357. tracer_signal();
  358. bit_reverse_tableau_X();
  359. calcul_fourier();
  360. tracer_graduations();
  361. tracer_tableau_valeurs();
  362. }
  363.  
  364.  
  365.  
  366. void MainWindow::on_spinBox_frq_valueChanged(int arg1) // SELECTEUR DE FREQUENCE 1
  367. {
  368. effacer_trace();
  369. frequence1 = arg1;
  370. RAZ_tableau_echantillons();
  371. remplir_tableau_echantillons();
  372. tracer_signal();
  373. calcul_tableau_W();
  374. bit_reverse_tableau_X();
  375. calcul_fourier();
  376. // tracer_graduations();
  377. tracer_tableau_valeurs();
  378. }
  379.  
  380.  
  381. void MainWindow::on_spinBox_frq_2_valueChanged(int arg1) // SELECTEUR DE FREQUENCE 2
  382. {
  383. effacer_trace();
  384. frequence2 = arg1;
  385. RAZ_tableau_echantillons();
  386. remplir_tableau_echantillons();
  387. tracer_signal();
  388. calcul_tableau_W();
  389. bit_reverse_tableau_X();
  390. calcul_fourier();
  391. // tracer_graduations();
  392. tracer_tableau_valeurs();
  393. }
  394.  
  395.  
  396. void MainWindow::on_spinBox_frq_3_valueChanged(int arg1)
  397. {
  398. effacer_trace();
  399. frequence3 = arg1;
  400. RAZ_tableau_echantillons();
  401. remplir_tableau_echantillons();
  402. tracer_signal();
  403. calcul_tableau_W();
  404. bit_reverse_tableau_X();
  405. calcul_fourier();
  406. // tracer_graduations();
  407. tracer_tableau_valeurs();
  408. }
  409.  
  410.  
  411. void MainWindow::on_checkBox1_toggled(bool checked) // HAMMING
  412. {
  413. effacer_trace();
  414. if (checked) {hamming = true;} else {hamming = false;}
  415. RAZ_tableau_echantillons();
  416. remplir_tableau_echantillons();
  417. tracer_signal();
  418. calcul_tableau_W();
  419. bit_reverse_tableau_X();
  420. calcul_fourier();
  421. // tracer_graduations();
  422. tracer_tableau_valeurs();
  423.  
  424. }
  425.  
  426.  
  427. void MainWindow::on_checkBox2_toggled(bool checked) // SECONDE FREQUENCE
  428. {
  429. effacer_trace();
  430. if (checked) {second_Frq = true;} else {second_Frq = false;}
  431. RAZ_tableau_echantillons();
  432. remplir_tableau_echantillons();
  433. tracer_signal();
  434. calcul_tableau_W();
  435. bit_reverse_tableau_X();
  436. calcul_fourier();
  437. // tracer_graduations();
  438. tracer_tableau_valeurs();
  439. }
  440.  
  441.  
  442. void MainWindow::on_checkBox3_toggled(bool checked)
  443. {
  444. effacer_trace();
  445. if (checked) {bz = true;} else {bz = false;}
  446. RAZ_tableau_echantillons();
  447. remplir_tableau_echantillons();
  448. tracer_signal();
  449. calcul_tableau_W();
  450. bit_reverse_tableau_X();
  451. calcul_fourier();
  452. // tracer_graduations();
  453. tracer_tableau_valeurs();
  454. }
  455.  
  456. void MainWindow::on_checkBox4_clicked(bool checked)
  457. {
  458. effacer_trace();
  459. if (checked) {modul = true;} else {modul = false;}
  460. RAZ_tableau_echantillons();
  461. remplir_tableau_echantillons();
  462. tracer_signal();
  463. calcul_tableau_W();
  464. bit_reverse_tableau_X();
  465. calcul_fourier();
  466. // tracer_graduations();
  467. tracer_tableau_valeurs();
  468. }
  469.  
  470.  
  471.  

3 -

Le code de la classe "Complexe" : complexe.ccp

CODE SOURCE C++ Qt4
  1.  
  2. #include "complexe.h"
  3. //
  4.  
  5. void Complexe::multi_lambda(float lambda)
  6. {
  7. a *= lambda;
  8. b *= lambda;
  9. }
  10.  
  11.  
  12. Complexe Complexe::operator+ (Complexe c)
  13. {
  14. Complexe r;
  15. r.a=this->a+c.a;
  16. r.b=this->b+c.b;
  17. return r;
  18. }
  19.  
  20.  
  21. Complexe Complexe::operator- (Complexe c)
  22. {
  23. Complexe r;
  24. r.a=this->a-c.a;
  25. r.b=this->b-c.b;
  26. return r;
  27. }
  28.  
  29.  
  30. Complexe Complexe::operator* (Complexe c)
  31. {
  32. Complexe r;
  33. r.a=this->a * c.a - this->b * c.b;
  34. r.b=this->a * c.b + this->b * c.a;
  35. return r;
  36. }
  37.  
  38.  
  39.  
  40. /* rappel: fonction en pascal
  41. function produit_cplx(c1,c2:complexe):complexe;
  42.  begin
  43.   produit_cplx[1]:=c1[1]*c2[1]-c1[2]*c2[2];
  44.   produit_cplx[2]:=c1[1]*c2[2]+c1[2]*c2[1];
  45.  end;
  46. */
  47.  
  48.  

4 -

Header file : complexe.h

CODE SOURCE C++ Qt4
  1.  
  2. #ifndef COMPLEXE_H
  3. #define COMPLEXE_H
  4. // Nombres complexes de la forme a+ib
  5.  
  6. //
  7. class Complexe
  8. {
  9.  
  10. public:
  11.  
  12. float a; // partie reelle
  13. float b; // partie imaginaire
  14. void multi_lambda(float lambda);
  15. Complexe operator+ (Complexe c);
  16. Complexe operator- (Complexe c);
  17. Complexe operator* (Complexe c); //equivaut a une rotation dans le plan (produit des modules, somme des arguments)
  18.  
  19.  
  20. };
  21.  
  22.  
  23.  
  24.  
  25. #endif
  26.  

5 Saisie d'écran de l'exécutable :

Le tracé du haut représente le signal temporel tel qu'il est échantillonné. Celui du bas représente sa transformée de Fourier (spectre fréquentiel). Bien que le nombre de points soit identique dans les deux cas (512 points), les unités ne sont pas les mêmes, des millisecondes en haut et des Hertz en bas.

Les courbes sont mises à jour automatiquement ("instantanément" dirons nous, vu la puissance de nos ordinateurs) lorsqu'on modifie les fréquences à la souris.

Le programme peut calculer la transformée de la somme de deux signaux sinusoïdaux de fréquences différentes.
La technique de "bourrage de zéros" (dont je n'ai pas parlé) est prise en compte. Utilisée conjointement avec la fenêtre de Hamming pour un signal constitué de deux sinusoïdes, elle permet d'obtenir des effets "remarquables" lorsque les deux fréquences se rapprochent puis se croisent... Je vous laisse découvrir cela.

Des modifications simples de ce programme peuvent permettre d’analyser d'autres formes d'ondes.

6 Documents

7 Version en C (n'utilisant pas la class 'complexes')

Il m'a été demandé s'il était possible de programmer la FFT en langage C plutôt qu'en C++, en particulier sans faire appel à la class (objet) 'Complexe'. ça tombe bien, cela faisait partie de mes projets afin de pouvoir à terme implémenter la FFT sur un microcontrôleur ATmega.

Voici donc le code dans lequel j'ai remplacé les manipulations de l'objet appelé 'complexe' par des opérations arithmétiques classiques sur les parties réelles et imaginaires et les tableaux de nombres complexes par des tableaux de deux nombres réels (partie réelle et imaginaire).

Le code obtenu est moins élégant, il perd en lisibilité, mais il constitue le premier pas vers la portabilité sur ATmega que je publierai prochainement. Ce qui reste du code en C++ et Qt4 ne concerne que l'affichage sur PC, sur l'ATmega ce sera remplacé par quelques opérations de lecture-écriture sur fichiers ou acquisition directe en entrée et sortie sur affichage (Module LCD 84*48 Nokia 5110 par exemple)

8 -

Le programme en C (sans la class 'Complexe') :

CODE SOURCE C++ Qt4
  1. /*
  2.   Programme écrit par Silicium628
  3.   ce logiciel est libre et open source
  4.  
  5.  
  6. 12 juin 2014
  7.  
  8. cette version dérive du projet "FTT"
  9. Il a pour but de tester la faisabilité de programmer la transformée en C et non en C++, c'est à dire sans faire appel à la programmation orientée objet.
  10. afin de l'implémenter sur un uC ATmega8
  11. */
  12.  
  13. #include "mainwindow.h"
  14. #include "ui_mainwindow.h"
  15. #include "math.h"
  16.  
  17. QString version = "2.0";
  18.  
  19. QColor couleur_ecran = QColor::fromRgb(3, 41, 11, 255);
  20. QColor couleur_ligne = QColor::fromRgb(167, 2, 0, 255);
  21. QColor couleur_trace = QColor::fromRgb(0, 210, 0, 255);
  22. QColor couleur_texte = QColor::fromRgb(255, 255, 0, 255);
  23. QColor couleur_signature = QColor::fromRgb(255, 255, 0, 255);
  24.  
  25. QPen pen_ligne(couleur_ligne, 1, Qt::SolidLine);
  26. QPen pen_trace(couleur_trace, 1, Qt::SolidLine);
  27. QPen pen_reticule(couleur_ligne, 1, Qt::SolidLine);
  28.  
  29.  
  30. typedef float complexe[2];
  31.  
  32. //complexe ech[84]; // nb_ech echantillons 1024 SUPPRIME ! et oui !! j'utilise le même tableau pour tout !
  33. complexe tab_X[84]; // nb_ech valeurs traitées 1024
  34. complexe tab_W[42]; // nb_ech/2 512
  35.  
  36.  
  37.  
  38. float frequence1=400; // Hz
  39. float frequence2=400; // Hz
  40. uint nb_etapes=6; //8
  41. uint nb_ech = pow (2,nb_etapes); // nombre d'échantillons = 2 puissance(nb_etapes)
  42. float frq_echantillonnage=8*frequence1; // Hz
  43. float Te; // période d'échantillonage
  44.  
  45. //int pas_balayage=40;
  46. int hamming = true;
  47. int second_Frq = false;
  48. int bz = true;
  49. uint forme;
  50.  
  51.  
  52.  
  53.  
  54.  
  55. MainWindow::MainWindow(QWidget *parent) :
  56. QMainWindow(parent)
  57.  
  58. {
  59. setupUi(this);
  60. setWindowTitle("Transformee de Fourier Rapide FFT (no objects) - version " + version);
  61.  
  62. scene = new QGraphicsScene(this);
  63. // scene->setBackgroundBrush(couleur_ecran);
  64. graphicsView1->setScene(scene);
  65. // graphicsView1->setGeometry(0,0,530,430);
  66. graphicsView1->setGeometry(140,30,84,48);
  67.  
  68. groupe_reticule = new QGraphicsItemGroup();
  69. scene->addItem(groupe_reticule);
  70. groupe_trace = new QGraphicsItemGroup();
  71. scene->addItem(groupe_trace);
  72.  
  73. effacer_trace();
  74. RAZ_tableau_echantillons();
  75. remplir_tableau_echantillons();
  76. // tracer_signal();
  77. calcul_tableau_W();
  78. // bit_reverse_tableau_X();
  79. calcul_fourier();
  80. tracer_tableau_valeurs();
  81.  
  82. }
  83.  
  84.  
  85.  
  86.  
  87. void MainWindow::effacer_graduation()
  88. {
  89. foreach( QGraphicsItem *item, scene->items( groupe_reticule->boundingRect() ) )
  90. {
  91. if( item->group() == groupe_reticule ) { delete item; }
  92. }
  93. }
  94.  
  95.  
  96. void MainWindow::effacer_trace()
  97. {
  98. foreach( QGraphicsItem *item, scene->items( groupe_trace->boundingRect() ) )
  99. {
  100. if( item->group() == groupe_trace ) { delete item; }
  101. }
  102. }
  103.  
  104.  
  105. MainWindow::~MainWindow()
  106. {
  107.  
  108. }
  109.  
  110.  
  111. void MainWindow::tracer_graduations()
  112. {
  113. float i,x,y;
  114. float nb_grad_max;
  115. float intervalle; // séparant les graduations
  116. float fi;
  117. QString sti;
  118. uint decal = 3; // leger decallage vertical pour ne pas masquer la trace
  119.  
  120. rectangle = new QGraphicsRectItem(0,decal, 511, 400);
  121. rectangle->setPen(couleur_ligne);
  122. groupe_trace->addToGroup(rectangle);
  123.  
  124. // lignes verticales
  125. // ici calcul de l'intervalle (en pixels) qui correspond exactement à 100Hz (en fréquence)
  126.  
  127. intervalle = 2000.0 * nb_ech * 2.0 / (frq_echantillonnage);
  128. nb_grad_max = 512 / intervalle;
  129.  
  130. for (i=0; i<=nb_grad_max; i++)
  131. {
  132. x = intervalle * i;
  133.  
  134. ligne1 = new QGraphicsLineItem(x, decal, x, 400 + decal);
  135. ligne1->setPen(pen_reticule);
  136. groupe_reticule->addToGroup(ligne1);
  137.  
  138. fi = x * frq_echantillonnage / nb_ech /4;
  139.  
  140. sti.setNum(fi);
  141. texte_frq = new QGraphicsTextItem(sti);
  142. texte_frq->setDefaultTextColor(couleur_texte);
  143. texte_frq->setPos(x,380);
  144. if(x<480) // évite que l'écriture ne déborde du cadre a droite
  145. {
  146. groupe_reticule->addToGroup(texte_frq);
  147. }
  148. }
  149.  
  150. // lignes horizontales
  151. intervalle = 50;
  152. nb_grad_max = 400 / intervalle;
  153.  
  154. for (i=0; i<=nb_grad_max; i++)
  155. {
  156. y = 400 + decal - intervalle * i;
  157.  
  158. ligne1 = new QGraphicsLineItem(0, y, 511, y);
  159. ligne1->setPen(pen_reticule);
  160. groupe_reticule->addToGroup(ligne1);
  161. }
  162.  
  163. texte_frq = new QGraphicsTextItem("Silicium628");
  164. texte_frq->setDefaultTextColor(couleur_signature);
  165. texte_frq->setPos(5,5);
  166. groupe_reticule->addToGroup(texte_frq);
  167.  
  168. scene->addItem(groupe_reticule);
  169. }
  170.  
  171.  
  172. void MainWindow::tracer_signal()
  173. {
  174. uint offset_y = 24;
  175. float echelle =168/nb_ech;
  176. uint n;
  177. float x, y, memo_x, memo_y;
  178.  
  179. segment_trace = new QGraphicsLineItem(0, offset_y, 512, offset_y);
  180. segment_trace->setPen(couleur_ligne);
  181. groupe_trace->addToGroup(segment_trace);
  182.  
  183. effacer_trace();
  184.  
  185. x=0;
  186. y=offset_y;
  187. for(n=0; n<64; n++)
  188. {
  189. memo_x = x;
  190. memo_y = y;
  191. x = n;
  192.  
  193. y = offset_y - 10 * tab_X[n][1]; // partie réelle
  194.  
  195. segment_trace = new QGraphicsLineItem(memo_x ,memo_y, x, y);
  196. segment_trace->setPen(pen_trace);
  197. groupe_trace->addToGroup(segment_trace);
  198. }
  199. scene->addItem(groupe_trace);
  200. }
  201.  
  202.  
  203.  
  204.  
  205. void MainWindow::tracer_tableau_valeurs()
  206. {
  207. uint offset_y = 46; // 600
  208. float echelle =168/nb_ech;
  209. uint n;
  210. float x, y, memo_x, memo_y;
  211.  
  212. int z;
  213.  
  214. if (bz) {z=4;} else {z=1;}
  215.  
  216.  
  217.  
  218. segment_trace = new QGraphicsLineItem(0, offset_y, 84, offset_y);
  219. segment_trace->setPen(couleur_ligne);
  220. groupe_trace->addToGroup(segment_trace);
  221.  
  222.  
  223.  
  224. x=0;
  225. y=offset_y;
  226. for(n=0; n<nb_ech/2; n++)
  227. {
  228. memo_x = x;
  229. memo_y = y;
  230. x = echelle * n;
  231.  
  232. // y = offset_y - z * sqrt( tab_X[n][1] * tab_X[n][1] + tab_X[n][2] * tab_X[n][2] ); // a=parties réelles, b=parties imaginaires
  233.  
  234. //version sans la racine carrée sqrt (qui pose des problèmes sur ATmega avr-gcc):
  235. y = offset_y - z * ( tab_X[n][1] * tab_X[n][1] + tab_X[n][2] * tab_X[n][2] ) / 8; // a=parties réelles, b=parties imaginaires
  236.  
  237.  
  238.  
  239. segment_trace = new QGraphicsLineItem(memo_x ,memo_y, x, y);
  240. segment_trace->setPen(pen_trace);
  241. groupe_trace->addToGroup(segment_trace);
  242. }
  243. scene->addItem(groupe_trace);
  244. }
  245.  
  246.  
  247.  
  248. void MainWindow::RAZ_tableau_echantillons()
  249. {
  250. uint n;
  251. for(n=0; n < nb_ech; n++)
  252. {
  253. tab_X[n][1] = 0; // partie reelle
  254. tab_X[n][2] = 0; // partie imaginaire
  255. }
  256. }
  257.  
  258.  
  259.  
  260. uint bit_reversal(uint num, uint nb_bits)
  261. {
  262. uint r = 0, i, s;
  263. if ( num > 1<< nb_bits) { return 0; }
  264.  
  265. for (i=0; i<nb_bits; i++)
  266. {
  267. s = (num & (1 << i));
  268. if(s) { r |= (1 << ((nb_bits - 1) - i)); }
  269. }
  270. return r;
  271. }
  272.  
  273.  
  274.  
  275. void MainWindow::remplir_tableau_echantillons()
  276. {
  277. // création du signal à analyser
  278. uint n, r, n_max;
  279. float x;
  280. Te = 1.0 / frq_echantillonnage;
  281.  
  282. if (bz==true) {n_max = nb_ech /4; } else {n_max = nb_ech; }
  283.  
  284. for(n=0; n < n_max; n++) // nb d'échantillons
  285. {
  286. x= cos(2.0 * M_PI * frequence1 * n * Te );
  287. if (second_Frq) {x+=cos(2 * M_PI * frequence2 * n * Te); }
  288. if (hamming) { x = x * (1- cos (2* M_PI * n / n_max ));}
  289. r=bit_reversal(n,nb_etapes);
  290. tab_X[r][1] = x; // partie reelle
  291. tab_X[r][2] = 0; // partie imaginaire
  292. }
  293. }
  294.  
  295.  
  296.  
  297.  
  298. /*
  299. void MainWindow::bit_reverse_tableau_X()
  300. {
  301.  
  302. // recopie les échantillons en les classant dans l'ordre 'bit reverse'
  303.   uint n,r;
  304.  
  305.   for(n=0; n < nb_ech; n++) // nb d'échantillons
  306.   {
  307.   r=bit_reversal(n,nb_etapes);
  308.   tab_X[n][1] = ech[r][1];
  309.   tab_X[n][2] = ech[r][2];
  310.   }
  311. }
  312.  
  313. */
  314.  
  315. void MainWindow::calcul_tableau_W()
  316. {
  317. // calcul et memorisation dans un tableau des twiddle factors
  318. uint n;
  319. float x;
  320.  
  321. for(n=0; n<(nb_ech/2-1); n++)
  322. {
  323. x=2.0*M_PI * n / nb_ech;
  324.  
  325. tab_W[n][1] = cos(x); // partie reelle
  326. tab_W[n][2] = -sin(x); // partie imaginaire
  327. }
  328. }
  329.  
  330.  
  331.  
  332. void MainWindow::calcul_fourier()
  333. {
  334. complexe produit; // voir le type "complexe" déclaré plus haut
  335.  
  336. uint decal =2;
  337.  
  338. uint etape, e1, e2, li, w, ci;
  339. uint li_max=nb_ech;
  340.  
  341. e2=1;
  342.  
  343. for (etape=1; etape<=nb_etapes; etape++)
  344. {
  345. e1=e2; //(e1 evite d'effectuer des divisions e2/2 plus bas)
  346. e2=e2+e2;
  347.  
  348. for (li=0; li<li_max; li+=1)
  349. {
  350. ci=li & (e2-1); // ET bit à bit
  351. if (ci>(e1-1))
  352. {
  353. w=li_max/e2*(li & (e1 -1)); // ET bit à bit calcul du numéro du facteur de tripatouillages W
  354.  
  355. produit[1]=tab_W[w][1]*tab_X[li][1]-tab_W[w][2]*tab_X[li][2];
  356. produit[2]=tab_W[w][1]*tab_X[li][2]+tab_W[w][2]*tab_X[li][1];
  357.  
  358. tab_X[li][1]=tab_X[li-e1][1]-produit[1];
  359. tab_X[li][2]=tab_X[li-e1][2]-produit[2];
  360.  
  361. tab_X[li-e1][1]=tab_X[li-e1][1]+produit[1];
  362. tab_X[li-e1][2]=tab_X[li-e1][2]+produit[2];
  363. }
  364. }
  365. }
  366. }
  367.  
  368.  
  369.  
  370.  
  371.  
  372. void MainWindow::on_spinBox_frq_valueChanged(int arg1) // SELECTEUR DE FREQUENCE 1
  373. {
  374. effacer_trace();
  375. frequence1 = arg1;
  376. RAZ_tableau_echantillons();
  377. remplir_tableau_echantillons();
  378. // tracer_signal();
  379. calcul_tableau_W();
  380. // bit_reverse_tableau_X();
  381. calcul_fourier();
  382. tracer_tableau_valeurs();
  383. }
  384.  
  385.  
  386. void MainWindow::on_spinBox_frq_2_valueChanged(int arg1) // SELECTEUR DE FREQUENCE 2
  387. {
  388. effacer_trace();
  389. frequence2 = arg1;
  390. RAZ_tableau_echantillons();
  391. remplir_tableau_echantillons();
  392. // tracer_signal();
  393. calcul_tableau_W();
  394. // bit_reverse_tableau_X();
  395. calcul_fourier();
  396. tracer_tableau_valeurs();
  397. }
  398.  
  399.  
  400.  
  401. void MainWindow::on_checkBox1_toggled(bool checked) // HAMMING
  402. {
  403. effacer_trace();
  404. if (checked) {hamming = true;} else {hamming = false;}
  405. RAZ_tableau_echantillons();
  406. remplir_tableau_echantillons();
  407. // tracer_signal();
  408. calcul_tableau_W();
  409. // bit_reverse_tableau_X();
  410. calcul_fourier();
  411. tracer_tableau_valeurs();
  412.  
  413. }
  414.  
  415.  
  416. void MainWindow::on_checkBox2_toggled(bool checked) // SECONDE FREQUENCE
  417. {
  418. effacer_trace();
  419. if (checked) {second_Frq = true;} else {second_Frq = false;}
  420. RAZ_tableau_echantillons();
  421. remplir_tableau_echantillons();
  422. // tracer_signal();
  423. calcul_tableau_W();
  424. // bit_reverse_tableau_X();
  425. calcul_fourier();
  426. tracer_tableau_valeurs();
  427. }
  428.  
  429.  
  430. void MainWindow::on_checkBox3_toggled(bool checked)
  431. {
  432. effacer_trace();
  433. if (checked) {bz = true;} else {bz = false;}
  434. RAZ_tableau_echantillons();
  435. remplir_tableau_echantillons();
  436. // tracer_signal();
  437. calcul_tableau_W();
  438. // bit_reverse_tableau_X();
  439. calcul_fourier();
  440. tracer_tableau_valeurs();
  441. }
  442.  
  443.  
  444.  
  445.  
  446.  
  447.  

9 -

21 juin 2014:
Je viens de supprimer le tableau des valeurs d'entrées, j'utilise directement le tableau de travail. (84x2x4 = 672 octets de gagnés en RAM une fois le programme implémenté sur un ATmega. C'est l'ATmega32 (2k SRAM) qui va être content !!! Sans ça j'aurais dû utiliser un ATmega128 (8k RAM) parce que j'ai besoin de pas mal de RAM également pour gérer efficacement l'afficheur LCD Nokia 5110. A ce propos je vous signale que j'ai fini (hier) cette implémentation sur ATmega et que ça fonctionne parfaitement. Il me reste plus qu'à entrer les valeurs avec le convertisseur A/N, à toiletter un peu le source obtenu et à la publier ici, ce qui ne devrait pas tarder, promis !

J'ai supprimé une vidéo ici parce que Dailymotion m'y incrustait une pub au début. Ce site n'est pas un panneau publicitaire !

10 Implémentation de la Transformée de Fourier Rapide sur ATmega32

11 -

Le premier jump permet de choisir entre visualiser le signal d'entrée ou bien sa transformée de Fourier.
Le deux suivants permettent d'appliquer (ou non) le hamming et le bourrage de zéros.
Le signal d'entrée doit avoir une amplitude de 1V environ et une fréquence comprise entre 100 Hz et 700 Hz.
Je vais pouvoir augmenter cette fréquence par simple changement des paramètres du CONV A/N de l'ATmega. A suivre donc...

12 -

22 juin 2014:

Nouvel affichage avec une courbe fine en place des bâtons verticaux;
La courbe est obtenue pour le calcul avec une fenêtre d'acquisition rectangulaire (au lieu de la fenêtre de hamming utilisée précédemment pour la vidéo). On distingue nettement les lobes secondaires de part et d'autre du pic central, malgré la (très) faible résolution de l'afficheur Nokia 5110 (84x48 pixels).

C'est une expérience qui repose sur une démonstration mathématique relativement corsée mais qui reste réalisable par tout un chacun pour une somme modique, de l'ordre de 5 euro en commandant l'ATmega et le LCD sur eBay. (Je me permets de donner ce mauvais conseil parce qu'il est difficile de trouver des ATmega32 chez l'épicier du coin).

13 Le programme (firmware) en C pour l'ATmega

CODE SOURCE pour ATmega32
  1. /***************************
  2. FFT.cpp
  3. par Silicium628
  4. logiciel Libre Open Source
  5. ****************************/
  6.  
  7. #include <avr/io.h>
  8. #include <stdint.h>
  9. #include <stdlib.h>
  10. #include "timeout.h" // comprend : #define F_CPU 16000000UL // 16 MHz
  11. #include "LCD_5110.cpp"
  12. #include <math.h>
  13.  
  14.  
  15. const char *version = "2.5"; // ( signal entrant 10kHz à 90kHz )
  16.  
  17. typedef float complexe[2];
  18.  
  19. complexe tab_X[64]; // nb_ech valeurs traitées
  20. complexe tab_W[32]; // nb_ech/2
  21.  
  22. float frequence1; // Hz
  23. float frq_echantillonnage; // Hz
  24.  
  25. const uint8_t nb_etapes=6;
  26. uint8_t nb_ech = pow (2,nb_etapes); // nombre d'échantillons = 2 puissance(nb_etapes)
  27.  
  28. char hamming;
  29. char bourrage;
  30. char mode_affi;
  31.  
  32. //const uint8_t hamming = false;
  33. //const uint8_t bourrage = true;
  34.  
  35.  
  36. void RAZ_tableau_echantillons()
  37. {
  38. uint8_t n;
  39. for(n=0; n < nb_ech; n++)
  40. {
  41. tab_X[n][1] = 0; // partie reelle
  42. tab_X[n][2] = 0; // partie imaginaire
  43. }
  44. }
  45.  
  46.  
  47. void init_ports (void) // ports perso
  48. // 0 = entree, 1=sortie ; les 1 sur les pins en entrees activent les R de Pull Up (tirage à VCC)
  49. {
  50. PORTA = 0b00000111;
  51. DDRA = 0b01111000; //PA7 en entree (ADC7 - entrée signal); PA0, PA1, PA2 jumpers à GND.
  52.  
  53. PORTB = 0b00000000;
  54. DDRB |= 0b11111111;
  55.  
  56. PORTC = 0b00000000;
  57. DDRC = 0b11111111;
  58.  
  59. PORTD = 0b00000000;
  60. DDRD = 0b11111111;
  61. }
  62.  
  63. /*
  64. The ADC module contains a prescaler, which generates an acceptable ADC clock fre-
  65. quency from any CPU frequency above 100 kHz. The prescaling is set by the ADPS bits
  66. in ADCSRA
  67. */
  68.  
  69. void InitADC (void)
  70. {
  71. // SFIOR = 0b10000000; // ADTS2,1,0 = 100 -> CAN declenché par Timer0 Overflow - p:218
  72. ADMUX = 0b01000111; //Bit 7:6 = ADC Reference Selection = 01 -> AVCC ; //Bits 0:4 - Analog Channel Selection Bits et gain éventuel - p214
  73. // ici Select pin ADC7 using MUX avec ref tension interne = AREF externe - pas de gain.
  74.  
  75. /** CHOISIR UNE SEULE LIGNE CI-DESSOUS pour ADCSRA */
  76.  
  77. // parametrage gamme 8-90 kHz
  78. ADCSRA = 0b10000010; // bit7: Activate ADC; ADPS2..0->Prescaler=f/4 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  79.  
  80. // parametrage gamme 4-56 kHz
  81. // ADCSRA = 0b10000011; // bit7: Activate ADC; ADPS2..0->Prescaler=f/8 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  82.  
  83. // parametrage gamme 2-31 kHz
  84. // ADCSRA = 0b10000100; // bit7: Activate ADC; ADPS2..0->Prescaler=f/16 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  85.  
  86. // parametrage gamme 1-15 kHz
  87. // ADCSRA = 0b10000101; // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  88.  
  89. // parametrage gamme 300 Hz-8000 Hz
  90. // ADCSRA = 0b10000110; // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  91.  
  92. // parametrage gamme 200 Hz- 4300 Hz
  93. // ADCSRA = 0b10000111; // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  94.  
  95. }
  96.  
  97.  
  98.  
  99.  
  100.  
  101. uint8_t bit_reversal(uint8_t num, uint8_t nb_bits)
  102. {
  103. uint8_t r = 0, i, s;
  104. if ( num > 1<< nb_bits) { return 0; }
  105.  
  106. for (i=0; i<nb_bits; i++)
  107. {
  108. s = (num & (1 << i));
  109. if(s) { r |= (1 << ((nb_bits - 1) - i)); }
  110. }
  111. return r;
  112. }
  113.  
  114.  
  115.  
  116. /***********************************************************************
  117. Remarques :
  118. - si l'on veut augmenter la vitesse d'acquisition et maîtriser précisément
  119. sa durée d'exécution il y a lieu de remplacer la fonction "acquisition()" par
  120. une routine écrite en assembleur et faire appel à une interruption pour le timing.
  121. - si l'on veut de plus augmenter la vitesse du traitement de la FFT, il faudra
  122. remplacer la fonction "calcul_fourier()" par une routine en assembleur, et en fait
  123. dans ce cas on envisagera de réécrire tout le programme en assembleur.
  124. - si on désire une plus grande rapidité encore, on pourra se tourner vers
  125. un microcontrôleur ARM cadencé à... 700 MHz (voire plus) et on pense alors bien sûr au Raspberry Pi !
  126. - si l'on veut calculer encore plus vite on peut utiliser le processeur d'une carte graphique,
  127. (par exemple en 2014 une Nvidia GeForce GTX680 dont le processeur tourne à à 2GHz- 2Go - DDR5 ou une Nvidia Quadro 6000 )
  128. ou au moins le FPU intégré dans les microprocesseurs de nos ordinateurs.
  129. - et si... on exige encore plus de vitesse on se tournera vers un Field Programmable Gate Arrays (FPGA)
  130. ou processeur de signal numérique (DSP) programmé en assembleur.
  131. voir :
  132. - http://fr.wikipedia.org/wiki/Processeur_de_signal_num%C3%A9rique
  133. - http://www.orfony.fr/bdtech/DSP.html
  134. - ADSP-TS001 (Analog Device)
  135. - Stratix FPGA (Altera)
  136.  
  137. Toutefois cette implémentation sur un simple ATmega à 4 euros n'a pas la prétention
  138. d’afficher des performances époustouflantes, mais elle permet de COMPRENDRE et de maîtriser
  139. la transformée de Fourier et de la voir fonctionner.
  140. ***********************************************************************/
  141.  
  142. void acquisition()
  143. {
  144. uint16_t acqH;
  145. uint16_t tab_acq[64];
  146. uint8_t n, r, n_max;
  147. float Te; // période d'échantillonage
  148. float x;
  149. Te = 1.0 / frq_echantillonnage;
  150. uint16_t valeur_lue;
  151.  
  152. if (bourrage) {n_max = nb_ech /4; } else {n_max = nb_ech; }
  153. //n_max = 16;
  154.  
  155. for(n=0; n < n_max-1; n++)
  156. {
  157. ADCSRA |= _BV(ADSC); //Start conversion - resolution 10bits
  158. while (ADCSRA & _BV(ADSC) ) {;} // attend la fin de la converstion
  159. tab_acq[n] = ADCL;
  160. acqH = ADCH; // passage à 16 bits pour pouvoir décaller
  161. tab_acq[n] += acqH << 8;
  162. }
  163. for(n=0; n < n_max-1; n++)
  164. {
  165. valeur_lue =tab_acq[n];
  166. if (hamming) { x = x * (1- cos (2* M_PI * n / n_max ));}
  167. x= (valeur_lue - 512.0) / 100.0;
  168. r=bit_reversal(n,nb_etapes);
  169. tab_X[r][1] = x; // partie reelle
  170. tab_X[r][2] = 0; // partie imaginaire
  171. }
  172. }
  173.  
  174.  
  175.  
  176. void remplir_tableau_echantillons()
  177. {
  178. // création du signal à analyser
  179. uint8_t n, r, n_max;
  180. float Te; // période d'échantillonage
  181. float x;
  182. Te = 1.0 / frq_echantillonnage;
  183.  
  184. if (bourrage==true) {n_max = nb_ech /4; } else {n_max = nb_ech; }
  185.  
  186. for(n=0; n < n_max-1; n++) // nb d'échantillons
  187. {
  188. x= cos(2.0 * M_PI * frequence1 * n * Te );
  189. if (hamming) { x = x * (1- cos (2* M_PI * n / n_max ));}
  190. r=bit_reversal(n,nb_etapes);
  191. tab_X[r][1] = x; // partie reelle
  192. tab_X[r][2] = 0; // partie imaginaire
  193. }
  194. }
  195.  
  196.  
  197.  
  198. void calcul_tableau_W()
  199. {
  200. // calcul et memorisation dans un tableau des twiddle factors
  201. uint8_t n;
  202. float x;
  203.  
  204. for(n=0; n<(nb_ech/2-1); n++)
  205. {
  206. x=2.0*M_PI * n / nb_ech;
  207.  
  208. tab_W[n][1] = cos(x); // partie reelle
  209. tab_W[n][2] = -sin(x); // partie imaginaire
  210. }
  211. }
  212.  
  213.  
  214.  
  215. void calcul_fourier()
  216. {
  217. complexe produit; // voir le type "complexe" déclaré plus haut
  218.  
  219. uint8_t decal =2;
  220.  
  221. uint8_t etape, e1, e2, li, w, ci;
  222. uint8_t li_max=nb_ech;
  223.  
  224. e2=1;
  225.  
  226. for (etape=1; etape<=nb_etapes; etape++)
  227. {
  228. e1=e2; //(e1 evite d'effectuer des divisions e2/2 plus bas)
  229. e2=e2+e2;
  230.  
  231. for (li=0; li<li_max; li+=1)
  232. {
  233. ci=li & (e2-1); // ET bit à bit
  234. if (ci>(e1-1))
  235. {
  236. w=li_max/e2*(li & (e1 -1)); // ET bit à bit calcul du numéro du facteur de tripatouillages W
  237.  
  238. produit[1]=tab_W[w][1]*tab_X[li][1]-tab_W[w][2]*tab_X[li][2];
  239. produit[2]=tab_W[w][1]*tab_X[li][2]+tab_W[w][2]*tab_X[li][1];
  240.  
  241. tab_X[li][1]=tab_X[li-e1][1]-produit[1];
  242. tab_X[li][2]=tab_X[li-e1][2]-produit[2];
  243.  
  244. tab_X[li-e1][1]=tab_X[li-e1][1]+produit[1];
  245. tab_X[li-e1][2]=tab_X[li-e1][2]+produit[2];
  246. }
  247. }
  248. }
  249. }
  250.  
  251.  
  252.  
  253.  
  254.  
  255. void tracer_signal()
  256. {
  257. float a, S;
  258. uint8_t r, memo_x, memo_y, n, x, y;
  259.  
  260. x=0;
  261. y=24;
  262. for(n=0; n<32; n++)
  263. {
  264. memo_x = x;
  265. memo_y = y;
  266. x = 2*n;
  267. r=bit_reversal(n, nb_etapes);
  268. if (hamming ==0) {a = 10.0;} else {a=5.0;}
  269. S = 20.0 - a * tab_X[r][1]; // partie réelle
  270. if ( (abs(S) )< 48 ) {y= lrint(S);} else {y=1;}
  271.  
  272. /** choisir une des trois lignes ci-dessous **/
  273. //buffer_drawPixel(x, y);
  274. //buffer_trace_baton_V(x, y);
  275. buffer_trace_segment( memo_x, memo_y, x, y);
  276. }
  277. affi_buffer(84);
  278. }
  279.  
  280.  
  281.  
  282. void tracer_tableau_valeurs()
  283. {
  284. uint8_t offset_y = 0;
  285. uint8_t n ;
  286. uint8_t memo_x, memo_y, x, y;
  287. float A;
  288.  
  289. uint8_t z;
  290.  
  291. if (bourrage) {z=15;} else {z=1;}
  292.  
  293. x=0;
  294. y=offset_y;
  295. for(n=0; n<32; n++)
  296. {
  297. memo_x = x;
  298. memo_y = y;
  299.  
  300. x = 2*n;
  301. //x = n;
  302.  
  303. A = z *( tab_X[n][1] * tab_X[n][1] + tab_X[n][2] * tab_X[n][2]);
  304. A = sqrt(A);
  305. if ( (abs(A) )< 48 ) {y= lrint(A);} else {y=1;}
  306.  
  307.  
  308. /** choisir une des trois lignes ci-dessous **/
  309. // buffer_drawPixel(x, 47-y); // la plus rapide, mais peu visible...
  310. // buffer_trace_baton_V(x, 1 + y);
  311. buffer_trace_segment( memo_x, memo_y, x, y);
  312. }
  313. affi_buffer(84);
  314. }
  315.  
  316.  
  317.  
  318. void test_out(uint8_t etat)
  319. {
  320. if (etat == 1) { PORTD |= 0b00000010;} else {PORTD &= 0b11111101;}
  321. }
  322.  
  323.  
  324.  
  325. int main(void)
  326. {
  327.  
  328. init_ports();
  329. InitADC();
  330. // InitINTs();
  331.  
  332. _delay_ms(100);
  333.  
  334. // define the 5 LCD Data pins: SCE, RST, DC, DATA, CLK
  335. lcd_init(&PORTB, PB1, &PORTB, PB0, &PORTB, PB2, &PORTB, PB3, &PORTB, PB4);
  336.  
  337. // frequence1 = 100.0; // Hz
  338. frq_echantillonnage=3200.0; // Hz
  339.  
  340.  
  341. lcd_goto_LC(0,1);
  342. lcd_puts2("FTT ATmega32");
  343.  
  344. _delay_ms(1000);
  345. lcd_clear();
  346.  
  347. calcul_tableau_W();
  348.  
  349. for(;;)
  350. {
  351. //frequence1 = 310.0; // Hz
  352.  
  353. for(uint8_t n=1; n<23; n++)
  354. {
  355. mode_affi = (PINA & 0b00000001) == 0;
  356. hamming = (PINA & 0b00000010) != 0;
  357. bourrage = (PINA & 0b00000100) != 0;
  358.  
  359. //frequence1 += 30;
  360.  
  361. RAZ_tableau_echantillons();
  362. acquisition();
  363. // remplir_tableau_echantillons();
  364.  
  365. if (mode_affi)
  366. {
  367. buffer_RAZ();
  368. tracer_signal(); // entrée ATTENTION : il faut faire cet affichage AVANT de calculer la TF parce que les résultats écrasent le tableau d'entrée
  369. // _delay_ms(1);
  370. }
  371.  
  372. if (! mode_affi)
  373. {
  374.  
  375. calcul_fourier(); // 100..120ms
  376. buffer_RAZ(); // 1.4ms
  377.  
  378. test_out(1);
  379. tracer_tableau_valeurs(); // 250 ms
  380. test_out(0);
  381. }
  382. lcd_goto_LC(5,0);
  383.  
  384.  
  385. if (! mode_affi) {lcd_puts2(" FFT ");} else {lcd_puts2(" Signal");} // 7 ms
  386.  
  387. }
  388. }
  389. }
  390.  

14 -

26 juin 2014:
L'acquisition fonctionne : je mets en ligne cette vidéo montrant le traitement du signal en temps réel pour une plage de fréquences de 300Hz à 2kHz environ.

28 juin 2014:
J'ai un peu optimisé la fonction d'acquisition ce qui permet d'augmenter la gamme de fréquences à analyser : elle s'étend maintenant de 10 kHz à 90 kHz (elle était de 300HZ à 2kHz dans la version précédente !). Ce qui est moins ridicule n'est-ce pas ? On est maintenant bien au delà du simple domaine audio.

J'ai supprimé une vidéo ici parce que Dailymotion m'y incrustait une pub au début. Ce site n'est pas un panneau publicitaire !

15 Remarques à propos de la vitesse :

La fréquence max est actuellement égale à 90 kHz.
- si l'on veut augmenter la vitesse d'acquisition et maîtriser précisément sa durée d'exécution il y a lieu de remplacer la fonction "acquisition()" par une routine écrite en assembleur et faire appel à une interruption pour le timing.
- si l'on veut de plus augmenter la vitesse du traitement de la FFT, il faudra remplacer la fonction "calcul_fourier()" par une routine en assembleur, et en fait dans ce cas on envisagera de réécrire tout le programme en assembleur.
- si on désire une plus grande rapidité encore, on pourra se tourner vers un microcontrôleur ARM cadencé à... 700 MHz (voire plus) et on pense alors bien sûr au Raspberry Pi !
- si l'on veut calculer encore plus vite on peut utiliser le processeur d'une carte graphique, ou au moins le FPU intégré dans les microprocesseurs de nos ordinateurs.
- et si... on exige encore plus de vitesse on se tournera vers un Field Programmable Gate Arrays (FPGA) ou processeur de signal numérique (DSP) programmé en assembleur. voir :
- http://fr.wikipedia.org/wiki/Processeur_de_signal_num%C3%A9rique
- http://www.orfony.fr/bdtech/DSP.html
- ADSP-TS001 (Analog Device)
- Stratix FPGA (Altera)

Toutefois cette implémentation sur un simple ATmega à 4 euros n'a pas la prétention d’afficher des performances époustouflantes, mais elle permet de COMPRENDRE et de maîtriser la transformée de Fourier et de la voir fonctionner.

16 Version avec calculs en virgule fixe

J'ai modifié le code en effectuant les calculs de la FFT sur des nombres en virgule fixe (le processeur traite des entiers 16 bits ce qui lui va fort bien, au lieu de nombres à virgule flottante "float" ce qui est très lent).

Résultat : La FFT est calculée en 3,5ms au lieu de 120 ms !!!! Soit 34 fois plus vite ! Voici ce nouveau code. Mais maintenant ce qui reste très lent c'est l'affichage. Mais je vais résoudre ça Mille Millions de Mille Sabords !

CODE SOURCE pour ATmega32
  1. /***************************
  2. FFT.cpp
  3. par Silicium628
  4. logiciel Libre Open Source
  5. ****************************/
  6.  
  7. #include <avr/io.h>
  8. #include <stdint.h>
  9. #include <stdlib.h>
  10. #include "timeout.h" // comprend : #define F_CPU 16000000UL // 16 MHz
  11. #include "LCD_5110-628.cpp"
  12. #include <math.h>
  13.  
  14.  
  15. const char *version = "4.2"; // ( signal entrant 10kHz à 90kHz )
  16.  
  17. //typedef float complexe[2];
  18. typedef int16_t complexe_f[2]; // entier signé pour calculs (rapides) en virgule fixe
  19.  
  20. //complexe tab_X[64]; // nb_ech valeurs traitées
  21. //complexe tab_W[32]; // nb_ech/2
  22. complexe_f tab_X[64]; // nb_ech valeurs traitées
  23. complexe_f tab_W[32]; // nb_ech/2
  24.  
  25. float frequence1; // Hz
  26. float frq_echantillonnage; // Hz
  27.  
  28. const uint8_t nb_etapes=6;
  29. uint8_t nb_ech = pow (2,nb_etapes); // nombre d'échantillons = 2 puissance(nb_etapes)
  30.  
  31. char hamming;
  32. char bourrage;
  33. char mode_affi;
  34.  
  35.  
  36. void RAZ_tableau_echantillons()
  37. {
  38. uint8_t n;
  39. for(n=0; n < nb_ech; n++)
  40. {
  41. tab_X[n][1] = 0; // partie reelle
  42. tab_X[n][2] = 0; // partie imaginaire
  43. }
  44. }
  45.  
  46.  
  47. void init_ports (void) // ports perso
  48. // 0 = entree, 1=sortie ; les 1 sur les pins en entrees activent les R de Pull Up (tirage à VCC)
  49. {
  50. PORTA = 0b00000111;
  51. DDRA = 0b01111000; //PA7 en entree (ADC7 - entrée signal); PA0, PA1, PA2 jumpers à GND.
  52.  
  53. PORTB = 0b00000000;
  54. DDRB |= 0b11111111;
  55.  
  56. PORTC = 0b00000000;
  57. DDRC = 0b11111111;
  58.  
  59. PORTD = 0b00000000;
  60. DDRD = 0b11111111;
  61. }
  62.  
  63. /*
  64. The ADC module contains a prescaler, which generates an acceptable ADC clock fre-
  65. quency from any CPU frequency above 100 kHz. The prescaling is set by the ADPS bits
  66. in ADCSRA
  67. */
  68.  
  69. void InitADC (void)
  70. {
  71. // SFIOR = 0b10000000; // ADTS2,1,0 = 100 -> CAN declenché par Timer0 Overflow - p:218
  72. ADMUX = 0b01000111; //Bit 7:6 = ADC Reference Selection = 01 -> AVCC ; //Bits 0:4 - Analog Channel Selection Bits et gain éventuel - p214
  73. // ici Select pin ADC7 using MUX avec ref tension interne = AREF externe - pas de gain.
  74.  
  75. /** CHOISIR UNE SEULE LIGNE CI-DESSOUS pour ADCSRA */
  76.  
  77. // parametrage gamme 8-90 kHz
  78. // ADCSRA = 0b10000010; // bit7: Activate ADC; ADPS2..0->Prescaler=f/4 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  79.  
  80. // parametrage gamme 4-56 kHz
  81. // ADCSRA = 0b10000011; // bit7: Activate ADC; ADPS2..0->Prescaler=f/8 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  82.  
  83. // parametrage gamme 2-31 kHz
  84. // ADCSRA = 0b10000100; // bit7: Activate ADC; ADPS2..0->Prescaler=f/16 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  85.  
  86. // parametrage gamme 1-15 kHz
  87. // ADCSRA = 0b10000101; // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  88.  
  89. // parametrage gamme 300 Hz-8000 Hz
  90. // ADCSRA = 0b10000110; // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  91.  
  92. // parametrage gamme 200 Hz- 4300 Hz
  93. ADCSRA = 0b10000111; // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  94.  
  95. }
  96.  
  97.  
  98.  
  99.  
  100.  
  101. uint8_t bit_reversal(uint8_t num, uint8_t nb_bits)
  102. {
  103. uint8_t r = 0, i, s;
  104. if ( num > 1<< nb_bits) { return 0; }
  105.  
  106. for (i=0; i<nb_bits; i++)
  107. {
  108. s = (num & (1 << i));
  109. if(s) { r |= (1 << ((nb_bits - 1) - i)); }
  110. }
  111. return r;
  112. }
  113.  
  114.  
  115.  
  116. /***********************************************************************
  117. Remarques :
  118. - si l'on veut augmenter la vitesse d'acquisition et maîtriser précisément
  119. sa durée d'exécution il y a lieu de remplacer la fonction "acquisition()" par
  120. une routine écrite en assembleur et faire appel à une interruption pour le timing.
  121. - si l'on veut de plus augmenter la vitesse du traitement de la FFT, il faudra
  122. remplacer la fonction "calcul_fourier()" par une routine en assembleur, et en fait
  123. dans ce cas on envisagera de réécrire tout le programme en assembleur.
  124. - si on désire une plus grande rapidité encore, on pourra se tourner vers
  125. un microcontrôleur ARM cadencé à... 700 MHz (voire plus) et on pense alors bien sûr au Raspberry Pi !
  126. - si l'on veut calculer encore plus vite on peut utiliser le processeur d'une carte graphique,
  127. (par exemple en 2014 une Nvidia GeForce GTX680 dont le processeur tourne à à 2GHz- 2Go - DDR5 ou une Nvidia Quadro 6000 )
  128. ou au moins le FPU intégré dans les microprocesseurs de nos ordinateurs.
  129. - et si... on exige encore plus de vitesse on se tournera vers un Field Programmable Gate Arrays (FPGA)
  130. ou processeur de signal numérique (DSP) programmé en assembleur.
  131. voir :
  132. - http://fr.wikipedia.org/wiki/Processeur_de_signal_num%C3%A9rique
  133. - http://www.orfony.fr/bdtech/DSP.html
  134. - ADSP-TS001 (Analog Device)
  135. - Stratix FPGA (Altera)
  136.  
  137. Toutefois cette implémentation sur un simple ATmega à 4 euros n'a pas la prétention
  138. d’afficher des performances époustouflantes, mais elle permet de COMPRENDRE et de maîtriser
  139. la transformée de Fourier et de la voir fonctionner.
  140. ***********************************************************************/
  141.  
  142. void acquisition()
  143. {
  144. uint16_t acqH;
  145. uint16_t tab_acq[64];
  146. uint8_t n, r, n_max;
  147. float x;
  148. uint16_t valeur_lue;
  149.  
  150. if (bourrage) {n_max = nb_ech /4; } else {n_max = nb_ech; }
  151. //n_max = 16;
  152.  
  153. for(n=0; n < n_max-1; n++)
  154. {
  155. ADCSRA |= _BV(ADSC); //Start conversion - resolution 10bits
  156. while (ADCSRA & _BV(ADSC) ) {;} // attend la fin de la converstion
  157. tab_acq[n] = ADCL;
  158. acqH = ADCH; // passage à 16 bits pour pouvoir décaller
  159. tab_acq[n] += acqH << 8;
  160. }
  161. for(n=0; n < n_max-1; n++)
  162. {
  163. valeur_lue =tab_acq[n];
  164. x= (valeur_lue - 512.0) / 100.0;
  165. if (hamming) { x = x * (1- cos (2* M_PI * n / n_max ));}
  166. r=bit_reversal(n,nb_etapes);
  167. tab_X[r][1] = x * 16; // partie reelle - * 16 pour obtenir un nombre en virgule fixe
  168. tab_X[r][2] = 0; // partie imaginaire
  169. }
  170. }
  171.  
  172.  
  173.  
  174. void remplir_tableau_echantillons()
  175. {
  176. // création du signal à analyser
  177. uint8_t n, r, n_max;
  178. float Te; // période d'échantillonage
  179. float x;
  180. Te = 1.0 / frq_echantillonnage;
  181.  
  182. if (bourrage==true) {n_max = nb_ech /2; } else {n_max = nb_ech; }
  183.  
  184. for(n=0; n < n_max-1; n++) // nb d'échantillons
  185. {
  186. x= cos(2.0 * M_PI * frequence1 * n * Te );
  187. if (hamming) { x = x * (1- cos (2* M_PI * n / n_max ));}
  188. r=bit_reversal(n,nb_etapes);
  189. tab_X[r][1] = x; // partie reelle
  190. tab_X[r][2] = 0; // partie imaginaire
  191. }
  192. }
  193.  
  194.  
  195.  
  196. void calcul_tableau_W()
  197. {
  198. // calcul et memorisation dans un tableau des twiddle factors
  199. uint8_t n;
  200. //float x;
  201. float x, cosx, sinx;
  202.  
  203. for(n=0; n<(nb_ech/2-1); n++)
  204. {
  205. x=2.0*M_PI * n / nb_ech;
  206. cosx = 16 * cos(x);
  207. sinx = -16 * sin(x);
  208.  
  209. tab_W[n][1] = round(cosx); // partie reelle
  210. tab_W[n][2] = round(sinx); // partie imaginaire
  211. }
  212. }
  213.  
  214.  
  215.  
  216. void calcul_fourier()
  217. {
  218. complexe_f produit; // voir le type "complexe" déclaré plus haut
  219.  
  220. uint8_t etape, e1, e2, li, w, ci;
  221. uint8_t li_max=nb_ech;
  222.  
  223. e2=1;
  224.  
  225. for (etape=1; etape<=nb_etapes; etape++)
  226. {
  227. e1=e2; //(e1 evite d'effectuer des divisions e2/2 plus bas)
  228. e2=e2+e2;
  229.  
  230. for (li=0; li<li_max; li+=1)
  231. {
  232. ci=li & (e2-1); // ET bit à bit
  233. if (ci>(e1-1))
  234. {
  235. w=li_max/e2*(li & (e1 -1)); // ET bit à bit calcul du numéro du facteur de tripatouillages W
  236.  
  237. produit[1]=((tab_W[w][1]*tab_X[li][1]) - (tab_W[w][2]*tab_X[li][2])) >> 4; // le décallage ( /16) because multiplication en virgule fixe
  238. produit[2]=((tab_W[w][1]*tab_X[li][2]) + (tab_W[w][2]*tab_X[li][1])) >> 4;
  239.  
  240. tab_X[li][1]=tab_X[li-e1][1]-produit[1];
  241. tab_X[li][2]=tab_X[li-e1][2]-produit[2];
  242.  
  243. tab_X[li-e1][1]=tab_X[li-e1][1]+produit[1];
  244. tab_X[li-e1][2]=tab_X[li-e1][2]+produit[2];
  245. }
  246. }
  247. }
  248. }
  249.  
  250.  
  251.  
  252. void tracer_signal()
  253. {
  254. int32_t a, S;
  255. uint8_t r, memo_x, memo_y, n, x, y;
  256.  
  257. x=0;
  258. y=24;
  259. for(n=0; n<32; n++)
  260. {
  261. memo_x = x;
  262. memo_y = y;
  263. x = 2*n;
  264. r=bit_reversal(n, nb_etapes);
  265. if (hamming ==0) {a = 10;} else {a=5;}
  266. S = (a * tab_X[r][1]); // partie réelle
  267. S = S >> 4;
  268. S = 20 - S;
  269. if ( (abs(S) )< 48 ) {y= lrint(S);} else {y=1;}
  270.  
  271. /** choisir une des trois lignes ci-dessous **/
  272. //buffer_drawPixel(x, y);
  273. //buffer_trace_baton_V(x, y);
  274. buffer_trace_segment( memo_x, memo_y, x, y);
  275. }
  276. affi_buffer(84);
  277. }
  278.  
  279.  
  280.  
  281. void tracer_tableau_valeurs()
  282. {
  283. uint8_t offset_y = 0;
  284. uint8_t n ;
  285. uint8_t memo_x, memo_y, x, y;
  286. uint32_t A;
  287.  
  288.  
  289. uint8_t z;
  290.  
  291. if (bourrage) {z=15;} else {z=5;}
  292.  
  293. x=0;
  294. y=offset_y;
  295.  
  296. for(n=0; n<31; n++) // batons =3.5 ms ; segments = 86 ms
  297. {
  298. memo_x = x;
  299. memo_y = y;
  300.  
  301. //x = n;
  302. x = 2*n; // largeur de la trace = 4 px (1 baton de 3px + 1 séparation de 1px)
  303.  
  304. A = ((tab_X[n][1] * tab_X[n][1]) + (tab_X[n][2] * tab_X[n][2]));
  305. A >>=8; // div/256; // 256 = 16 * 16 (compense la multiplication en virgule fixe puis retour à une valeur ordinaire)
  306. A *= z;
  307. A = sqrt(A);
  308.  
  309. if ( (abs(A) )< 48 ) {y= lrint(A);} else {y=47;}
  310.  
  311. /** choisir une des trois lignes ci-dessous; ATTENTION : les rapidités ne sont pas les mêmes !! **/
  312. // buffer_drawPixel(x, 47-y); // la plus rapide, mais peu visible...
  313. buffer_trace_baton_V(x, 1 + y);
  314. // buffer_trace_segment( memo_x, memo_y, x, y);
  315. }
  316.  
  317. affi_buffer(84); // 12.6 ms
  318.  
  319. }
  320.  
  321.  
  322.  
  323.  
  324.  
  325. int main(void)
  326. {
  327.  
  328. init_ports();
  329. InitADC();
  330. // InitINTs();
  331.  
  332. _delay_ms(100);
  333. lcd_init();
  334.  
  335. // frequence1 = 100.0; // Hz
  336. // frq_echantillonnage=3200.0; // Hz
  337.  
  338.  
  339. lcd_goto_LC(0,1);
  340. lcd_puts2("FTT ATmega32");
  341.  
  342. _delay_ms(1000);
  343. lcd_clear();
  344.  
  345. calcul_tableau_W();
  346.  
  347. lcd_goto_LC(5,0);
  348. if (! mode_affi) {lcd_puts2("TF Silicium628 ");} else {lcd_puts2(" Signal");} // 7 ms
  349.  
  350. for(;;)
  351. {
  352. //frequence1 = 310.0; // Hz
  353.  
  354. for(uint8_t n=1; n<23; n++)
  355. {
  356. mode_affi = (PINA & 0b00000001) == 0;
  357. hamming = (PINA & 0b00000010) != 0;
  358. bourrage = (PINA & 0b00000100) != 0;
  359.  
  360. //frequence1 += 30;
  361.  
  362.  
  363. RAZ_tableau_echantillons();
  364. acquisition(); // 7ms
  365.  
  366. // remplir_tableau_echantillons();
  367.  
  368. if (mode_affi)
  369. {
  370. buffer_RAZ();
  371. tracer_signal(); // entrée ATTENTION : il faut faire cet affichage AVANT de calculer la TF parce que les résultats écrasent le tableau d'entrée
  372. }
  373.  
  374. if (! mode_affi)
  375. {
  376. calcul_fourier(); // 7 ms avec calculs en virgule fixe (int 16 bits) au lieu de 100..120ms avec calculs en virgule flotante (float 32 bits)
  377. buffer_RAZ(); // 1.4 ms
  378.  
  379. tracer_tableau_valeurs(); // 16ms
  380. }
  381.  
  382. }
  383. }
  384. }
  385.  

17 Version turbo !

01 juillet 2014:

Après avoir diminué drastiquement la durée de calcul de la TF en utilisant un format de nombres en virgule fixe, je me suis attaqué à la partie affichage.

J'ai revu les fonctions de commande de l'afficheur LCD 5110, tout en gardant l'écriture en langage C.

Résultat : La boucle complète comprenant l'acquisition, le calcul FFT et l'affichage s’exécute en 28 ms à peine ! (Au départ on tournait à 400 ms environ).
Et ça change tout ! 28 ms ça permet de faire une analyse visuelle de l'audio en temps réel. Un "Analyser de spectre audio", vous savez avec des petites barres verticales qui dansent en mesure avec la musique. Et ça marche. Je vous posterai ici la vidéo prochainement.

On en voit partout de ces trucs là, dans les logiciels, sur Androïd etc... Oui mais celui-ci on l'a fait entièrement nous même, en comprenant toutes les étapes du raisonnement qui conduit au résultat. Souvenez-vous, on est parti de ça ==>

18 -


19 Documents :


20 Le coin de l'ARDUINO

Moi qui apprécie bien les microcontrôleurs AVR je surveille du coin de l'oeil ce qui se passe au royaume de l'ARDUINO. On y rencontre visiblement pas mal de fous furieux dans mon genre. Mais l'aspect commercial de la chose m'a toujours rebuté. Toutefois avec l'apparition de nombreux clones les choses évoluent vers plus d'ouverture. Ce qui m'a vraiment intéressé ces temps derniers c'est la symbiose réussie entre des cartes puissantes (basées sur des ATmega 2560) et des petits afficheurs TFT couleur au format des téléphones portables (3.5") tactiles qui plus est.

J'ai donc décidé de tester ces bestioles et la transformée de Fourier était le projet idéal pour ça. Je me suis offert un clone SainSmart 2560 + Ecran TFT 320x240 px + carte d'adaptation servant à relier l'un à l'autre.

21 -

22 -

Après une demi journée pour apprivoiser l'animal (et réussir à le programmer en C++ avec Geany comme éditeur en remplacement de l'EDI Arduino qui ne me convenait pas) j'ai pu facilement adapter notre programme de Transformée de Fourier Rapide pour la carte Mega 2560 et l'afficheur TFT. Voici le résultat :

Les deux petits carrés verts à droite ce sont des boutons. si ! si ! lorsqu'on clique dessus avec un stylet, il s'allument et s’éteignent en modifiant au passage les paramètres du programme : celui du haut permet d'appliquer ou non la fenêtre de Hamming sur le signal d'entrée, et le second applique ou non le bourrage de zéros. Je vous montre cela en vidéo ci-dessous.

23 -

Je tiens à préciser que la musique d'illustration de la vidéo ne constitue nullement le signal à analyser (Il est créé en interne par la fonction sinus...) Mais bien entendu je vais activer la partie acquisition du signal d'entrée par le convertisseur ADC de l'ATmega. Comme d'habitude je vous tiens au courant de l'évolution des choses.

24 Le programme (firmware) en C++ pour la carte Mega2560 + TFT

CODE SOURCE pour ATmega2560
  1.  
  2. #include <UTFT.h>
  3. #include <ITDB02_Touch.h>
  4.  
  5. #include <avr/io.h>
  6. #include <stdint.h>
  7. #include <stdlib.h>
  8. //#include "timeout.h" // comprend : #define F_CPU 16000000UL // 16 MHz
  9. #include <math.h>
  10.  
  11. // Declare which fonts we will be using
  12. extern uint8_t SmallFont[];
  13.  
  14. // Arduino Mega:
  15. // -------------------
  16. // Standard Arduino Mega/Due shield : <display model>,38,39,40,41
  17. // CTE TFT LCD/SD Shield for Arduino Mega : <display model>,38,39,40,41
  18.  
  19.  
  20. const char *version = "1.1"; // ( signal entrant 10kHz à 90kHz )
  21.  
  22. UTFT TFT320(ITDB32S,38,39,40,41);
  23. ITDB02_Touch TACTILE(6,5,4,3,2);
  24.  
  25. typedef float complexe[2];
  26.  
  27. complexe tab_X[128]; // nb_ech valeurs traitées
  28. complexe tab_W[64]; // nb_ech/2
  29.  
  30.  
  31. float frequence1; // Hz
  32. float frq_echantillonnage; // Hz
  33.  
  34. const uint8_t nb_etapes=7;
  35. uint8_t nb_ech = pow (2,nb_etapes); // nombre d'échantillons = 2 puissance(nb_etapes)
  36.  
  37. uint8_t hamming;
  38. uint8_t bourrage;
  39. uint8_t mode_affi;
  40.  
  41. uint8_t etat_boutons; // 8 bits definissant chacun l'état d'un bouton, soit 8 boutons en tout
  42.  
  43. void init_ports (void) // ports perso
  44. // 0 = entree, 1=sortie ; les 1 sur les pins en entrees activent les R de Pull Up (tirage à VCC)
  45. {
  46. PORTA = 0b00000111;
  47. DDRA = 0b01111000; //PA7 en entree (ADC7 - entrée signal); PA0, PA1, PA2 jumpers à GND.
  48.  
  49. PORTB = 0b00000000;
  50. DDRB |= 0b11111111;
  51.  
  52. PORTC = 0b00000000;
  53. DDRC = 0b11111111;
  54.  
  55. PORTD = 0b00000000;
  56. DDRD = 0b11111111;
  57. }
  58.  
  59.  
  60.  
  61. /*
  62. The ADC module contains a prescaler, which generates an acceptable ADC clock fre-
  63. quency from any CPU frequency above 100 kHz. The prescaling is set by the ADPS bits
  64. in ADCSRA
  65. */
  66.  
  67. void InitADC (void)
  68. {
  69. // SFIOR = 0b10000000; // ADTS2,1,0 = 100 -> CAN declenché par Timer0 Overflow - p:218
  70. ADMUX = 0b01000111; //Bit 7:6 = ADC Reference Selection = 01 -> AVCC ; //Bits 0:4 - Analog Channel Selection Bits et gain éventuel - p214
  71. // ici Select pin ADC7 using MUX avec ref tension interne = AREF externe - pas de gain.
  72.  
  73. /** CHOISIR UNE SEULE LIGNE CI-DESSOUS pour ADCSRA */
  74.  
  75. // parametrage gamme 8-90 kHz
  76. // ADCSRA = 0b10000010; // bit7: Activate ADC; ADPS2..0->Prescaler=f/4 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  77.  
  78. // parametrage gamme 4-56 kHz
  79. // ADCSRA = 0b10000011; // bit7: Activate ADC; ADPS2..0->Prescaler=f/8 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  80.  
  81. // parametrage gamme 2-31 kHz
  82. // ADCSRA = 0b10000100; // bit7: Activate ADC; ADPS2..0->Prescaler=f/16 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  83.  
  84. // parametrage gamme 1-15 kHz
  85. // ADCSRA = 0b10000101; // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  86.  
  87. // parametrage gamme 300 Hz-8000 Hz
  88. // ADCSRA = 0b10000110; // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  89.  
  90. // parametrage gamme 200 Hz- 4300 Hz
  91. ADCSRA = 0b10000111; // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
  92.  
  93. }
  94.  
  95. void RAZ_tableau_echantillons()
  96. {
  97. uint8_t n;
  98. for(n=0; n < nb_ech; n++)
  99. {
  100. tab_X[n][1] = 0; // partie reelle
  101. tab_X[n][2] = 0; // partie imaginaire
  102. }
  103. }
  104.  
  105. uint8_t bit_reversal(uint8_t num, uint8_t nb_bits)
  106. {
  107. uint8_t r = 0, i, s;
  108. if ( num > 1<< nb_bits) { return 0; }
  109.  
  110. for (i=0; i<nb_bits; i++)
  111. {
  112. s = (num & (1 << i));
  113. if(s) { r |= (1 << ((nb_bits - 1) - i)); }
  114. }
  115. return r;
  116. }
  117.  
  118.  
  119.  
  120. void acquisition()
  121. {
  122. uint16_t acqH;
  123. uint16_t tab_acq[64];
  124. uint8_t n, r, n_max;
  125. float x;
  126. uint16_t valeur_lue;
  127.  
  128. if (bourrage != 0) {n_max = nb_ech; } else {n_max = nb_ech /2; }
  129. //n_max = 16;
  130.  
  131. for(n=0; n < n_max-1; n++)
  132. {
  133. ADCSRA |= _BV(ADSC); //Start conversion - resolution 10bits (0..1024)
  134. while (ADCSRA & _BV(ADSC) ) {;} // attend la fin de la converstion
  135. tab_acq[n] = ADCL;
  136. acqH = ADCH; // passage à 16 bits pour pouvoir décaller
  137. tab_acq[n] += acqH << 8; // ici tab_acq[n] = 0..1024
  138. }
  139. for(n=0; n < n_max-1; n++)
  140. {
  141. valeur_lue =tab_acq[n];
  142. x= (valeur_lue - 512.0) / 100.0; // 512 = 1024/2
  143. if (hamming) { x = x * (1- cos (2* M_PI * n / n_max ));}
  144. r=bit_reversal(n,nb_etapes);
  145. tab_X[r][1] = x * 16; // partie reelle - * 16 pour obtenir un nombre en virgule fixe
  146. tab_X[r][2] = 0; // partie imaginaire
  147. }
  148. }
  149.  
  150.  
  151.  
  152. void remplir_tableau_echantillons()
  153. {
  154. // création du signal à analyser
  155. uint8_t n, r, n_max;
  156. float Te; // période d'échantillonage
  157. float y;
  158. float a;
  159. Te = 1.0 / frq_echantillonnage;
  160.  
  161. if (bourrage != 0) { n_max = nb_ech /4; } else { n_max = nb_ech; }
  162.  
  163. for(n=0; n < n_max-1; n++) // nb d'échantillons
  164. {
  165. y= cos( 2.0 * M_PI * frequence1 * n * Te );
  166. if (hamming) { y = y * (1- cos (2* M_PI * n / n_max ));}
  167.  
  168.  
  169. r=bit_reversal(n,nb_etapes);
  170. y *= 2;
  171.  
  172. tab_X[r][1] = y; // partie reelle
  173. tab_X[r][2] = 0; // partie imaginaire
  174. }
  175. }
  176.  
  177.  
  178.  
  179. void calcul_tableau_W()
  180. {
  181. // calcul et memorisation dans un tableau des twiddle factors
  182. uint8_t n;
  183. //float x;
  184. float x, cosx, sinx;
  185.  
  186. for(n=0; n<(nb_ech/2-1); n++)
  187. {
  188. x=2.0*M_PI * n / nb_ech;
  189.  
  190. tab_W[n][1] = cos(x); // partie reelle
  191. tab_W[n][2] = -sin(x); // partie imaginaire
  192. }
  193. }
  194.  
  195.  
  196.  
  197.  
  198. void calcul_fourier()
  199. {
  200. complexe produit; // voir le type "complexe" déclaré plus haut
  201.  
  202. uint8_t etape, e1, e2, li, w, ci;
  203. uint8_t li_max=nb_ech;
  204.  
  205. e2=1;
  206.  
  207. for (etape=1; etape<=nb_etapes; etape++)
  208. {
  209. e1=e2; //(e1 evite d'effectuer des divisions e2/2 plus bas)
  210. e2=e2+e2;
  211.  
  212. for (li=0; li<li_max; li+=1)
  213. {
  214. ci=li & (e2-1); // ET bit à bit
  215. if (ci>(e1-1))
  216. {
  217. w=li_max/e2*(li & (e1 -1)); // ET bit à bit calcul du numéro du facteur de tripatouillages W
  218.  
  219. produit[1]=(tab_W[w][1]*tab_X[li][1]) - (tab_W[w][2]*tab_X[li][2]);
  220. produit[2]=(tab_W[w][1]*tab_X[li][2]) + (tab_W[w][2]*tab_X[li][1]);
  221.  
  222. tab_X[li][1]=tab_X[li-e1][1]-produit[1];
  223. tab_X[li][2]=tab_X[li-e1][2]-produit[2];
  224.  
  225. tab_X[li-e1][1]=tab_X[li-e1][1]+produit[1];
  226. tab_X[li-e1][2]=tab_X[li-e1][2]+produit[2];
  227. }
  228. }
  229. }
  230. }
  231.  
  232.  
  233.  
  234.  
  235. void tracer_signal()
  236. {
  237.  
  238. uint8_t offset_y = 60;
  239. int32_t a, S;
  240. uint8_t r, memo_x, memo_y, n, x, y;
  241.  
  242. x=0;
  243. y= offset_y;
  244.  
  245. TFT320.setColor(200, 0, 0);
  246. TFT320.drawLine(1,offset_y,319,offset_y);
  247.  
  248. TFT320.setColor(255, 255, 0);
  249. TFT320.setBackColor(100, 100, 100);
  250.  
  251.  
  252.  
  253. for(n=0; n<128; n++)
  254. {
  255. memo_x = x;
  256. memo_y = y;
  257. x = 2*n;
  258. r=bit_reversal(n, nb_etapes);
  259. if (hamming !=0) {a = 10.0;} else {a=20.0;}
  260. S = a * tab_X[r][1]; // partie réelle
  261. if ( (abs(S) )< 240 ) {y= lrint(S);} else {y=1;}
  262.  
  263. /** choisir une des lignes ci-dessous **/
  264.  
  265. //buffer_trace_baton_V(x, y);
  266.  
  267. y= offset_y - y;
  268. TFT320.drawLine(memo_x, memo_y, x, y);
  269.  
  270. }
  271. // affi_buffer(84);
  272. }
  273.  
  274.  
  275.  
  276.  
  277. void tracer_tableau_valeurs()
  278. {
  279. uint8_t offset_x = 100;
  280. uint8_t offset_y = 200;
  281. uint8_t n ;
  282. uint8_t memo_x, memo_y, x, y;
  283. float A;
  284.  
  285. uint8_t z;
  286.  
  287. if (bourrage != 0 ) {z=4;} else {z=1;}
  288.  
  289. x=0;
  290. y= offset_y;
  291.  
  292. TFT320.setColor(200, 0, 0);
  293. TFT320.drawLine(1,offset_y,319,offset_y);
  294.  
  295. TFT320.setColor(0, 180, 200);
  296. TFT320.setBackColor(0, 0, 0);
  297.  
  298.  
  299. for(n=0; n<64; n++)
  300. {
  301. memo_x = x;
  302. memo_y = y;
  303.  
  304. //x = n;
  305. x = offset_x + 2*n;
  306.  
  307. A = ((tab_X[n][1] * tab_X[n][1]) + (tab_X[n][2] * tab_X[n][2]));
  308. A = z * sqrt(A);
  309.  
  310. if ( (abs(A) )< 240 ) {y= lrint(A);} else {y=240;}
  311.  
  312. y= offset_y - y;
  313.  
  314. /** choisir une des trois lignes ci-dessous; ATTENTION : les rapidités ne sont pas les mêmes !! **/
  315. // buffer_drawPixel(x, 47-y); // la plus rapide, mais peu visible...
  316. // buffer_trace_baton_V(x, 1 + y);
  317. // buffer_trace_segment( memo_x, memo_y, x, y);
  318.  
  319. //TFT320.drawPixel(x, 200-y);
  320. TFT320.drawLine(memo_x, memo_y, x, y);
  321. }
  322.  
  323. // affi_buffer(84); // 12.6 ms
  324.  
  325. }
  326.  
  327.  
  328. void trace_entete()
  329. {
  330. TFT320.setColor(64, 64, 64);
  331. TFT320.fillRect(0, 0, 319, 13); // bandeau en haut
  332.  
  333. TFT320.setColor(64, 64, 64);
  334. TFT320.fillRect(0, 226, 319, 239); // bandeau en bas
  335.  
  336. TFT320.setColor(255, 255, 255);
  337. TFT320.setBackColor(64, 64, 64);
  338. TFT320.print("Transformee de Fourier", CENTER, 1);
  339.  
  340. TFT320.setColor(255,255,0);
  341. TFT320.print("ATmega2560 - Silicium628", CENTER, 227);
  342.  
  343. }
  344.  
  345.  
  346. void trace_reticule()
  347. {
  348. // cadre
  349. TFT320.setColor(0, 150, 0);
  350. TFT320.setBackColor(0, 0, 0);
  351. TFT320.drawRect(0, 14, 319, 225);
  352.  
  353. //reticule
  354. TFT320.drawLine(159, 15, 159, 224);
  355. TFT320.drawLine(1, 119, 318, 119);
  356.  
  357. for (int i=9; i<310; i+=10) { TFT320.drawLine(i, 117, i, 121); }
  358. for (int i=19; i<220; i+=10) { TFT320.drawLine(157, i, 161, i); }
  359. }
  360.  
  361.  
  362. void trace_bouton(uint8_t n, uint8_t etat)
  363. {
  364. uint16_t position_x = 290;
  365. uint8_t dy = 25 * n;
  366. if (etat != 0) { TFT320.setColor(0, 255, 0); } else { TFT320.setColor(0, 80, 0); }
  367. TFT320.fillRoundRect(position_x, 20+dy, position_x+20, 40+dy);
  368. TFT320.setColor(255, 255, 255);
  369. TFT320.drawRect(position_x, 20+dy, position_x+20, 40+dy);
  370.  
  371. }
  372.  
  373.  
  374.  
  375. void trace_ecran()
  376. {
  377. TFT320.setColor(0, 0, 0);
  378. TFT320.setBackColor(0, 0, 0);
  379. TFT320.fillRect(0, 14, 319, 225);
  380.  
  381. trace_entete();
  382. trace_reticule();
  383.  
  384. trace_bouton(0, hamming);
  385. trace_bouton(1, bourrage);
  386.  
  387. TFT320.setColor(255, 255, 255);
  388. TFT320.print("H", 265, 25);
  389. TFT320.print("B", 265, 50);
  390.  
  391. }
  392.  
  393.  
  394.  
  395. void setup()
  396. {
  397. randomSeed(analogRead(0));
  398.  
  399. InitADC();
  400.  
  401. TFT320.InitLCD();
  402. TFT320.setFont(SmallFont);
  403. TFT320.clrScr();
  404.  
  405. TACTILE.InitTouch(LANDSCAPE);
  406. TACTILE.setPrecision(PREC_MEDIUM);
  407.  
  408. trace_ecran();
  409. }
  410.  
  411.  
  412. void polling_boutons()
  413. {
  414. uint8_t xt0 = 8; // position horizontale des boutons (alignés verticalement)
  415. uint8_t yt0;
  416. uint8_t ds = 8; // espacement vertical des boutons
  417. uint8_t xt, yt;
  418.  
  419. while (TACTILE. dataAvailable() == true)
  420. {
  421. TACTILE.read();
  422. xt = TACTILE.getX();
  423. yt = TACTILE.getY();
  424.  
  425. TFT320.printNumI(xt, 280, 215, 3, ' ');
  426. TFT320.printNumI(yt, 280, 225, 3, ' ');
  427.  
  428. yt0 = 28; // position verticale (centre) du premier bouton
  429. if ((xt > xt0-ds) && (xt < xt0+ds) && (yt > yt0-ds) && (yt < yt0+ds))
  430. {
  431. etat_boutons ^=0b00000001;
  432. hamming = etat_boutons & 0b00000001;
  433. trace_ecran();
  434. delay(300);
  435. }
  436.  
  437. yt0 = 54; // position verticale (centre) du second bouton
  438. if ((xt > xt0-ds) && (xt < xt0+ds) && (yt > yt0-ds) && (yt < yt0+ds))
  439. {
  440. etat_boutons ^=0b00000010;
  441. bourrage = etat_boutons & 0b00000010;
  442. trace_ecran();
  443. delay(300);
  444. }
  445. }
  446. delay(100);
  447.  
  448. }
  449.  
  450.  
  451. void loop()
  452. {
  453.  
  454. randomSeed(analogRead(0));
  455.  
  456. frequence1 = 500.0; // Hz
  457. frq_echantillonnage=3200.0; // Hz
  458.  
  459. calcul_tableau_W();
  460.  
  461. //frequence1 = 310.0; // Hz
  462.  
  463. // for(uint8_t n=1; n<23; n++)
  464. // {
  465.  
  466. //frequence1 += 30;
  467.  
  468.  
  469. RAZ_tableau_echantillons();
  470. // acquisition(); // 7ms
  471. remplir_tableau_echantillons();
  472. tracer_signal(); // entrée ATTENTION : il faut faire cet affichage AVANT de calculer la TF parce que les résultats écrasent le tableau d'entrée
  473. calcul_fourier(); // 7 ms avec calculs en virgule fixe (int 16 bits) au lieu de 100..120ms avec calculs en virgule flotante (float 32 bits)
  474. tracer_tableau_valeurs(); // 16ms
  475.  
  476. polling_boutons();
  477.  
  478. }
  479.  
  480.  

25 Documents pour l'ARDUINO (ou clone)

Cette archive comprend le fichier principal (.ino) et les bibliothèques incluses (UTFT à ne pas confondre avec FFT, TFT c'est pour l'afficheur TFT) et ITDB02_Touch.cpp pour la partie tactile de l'afficheur)

26 Une application sur ESP32 :

27 -

Liens...

34168