#include "complexe.h"
//
 
void Complexe::multi_lambda(float lambda)
{ 
	a *= lambda;
	b *= lambda;
}
 
 
Complexe Complexe::operator+ (Complexe c)
{
	Complexe r;
	r.a=this->a+c.a;
	r.b=this->b+c.b;
	return r;
}
 
 
Complexe Complexe::operator- (Complexe c)
{
    Complexe r;
    r.a=this->a-c.a;
    r.b=this->b-c.b;
    return r;
}
 
 
Complexe Complexe::operator* (Complexe c)
{
	Complexe r;
    r.a=this->a * c.a - this->b * c.b;
    r.b=this->a * c.b + this->b * c.a;
	return r;
}
 
 
 
/* rappel: fonction en pascal
function produit_cplx(c1,c2:complexe):complexe;
 begin
   produit_cplx[1]:=c1[1]*c2[1]-c1[2]*c2[2];
   produit_cplx[2]:=c1[1]*c2[2]+c1[2]*c2[1];
 end;
*/
 
 

 
#ifndef COMPLEXE_H
#define COMPLEXE_H
// Nombres complexes de la forme a+ib
 
//
class Complexe  
{
 
public:
 
    float a; // partie reelle
    float b; // partie imaginaire
    void multi_lambda(float lambda);
    Complexe operator+ (Complexe c);
    Complexe operator- (Complexe c);
    Complexe operator* (Complexe c); //equivaut a une rotation dans le plan (produit des modules, somme des arguments)
 
 
};
 
 
 
 
#endif
 

/*
  Programme écrit par Silicium628
  ce logiciel est libre et open source
 
 
12 juin 2014
 
cette version dérive du projet  "FTT"
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.
afin de l'implémenter sur un uC ATmega8
*/
 
#include "mainwindow.h"
#include "ui_mainwindow.h"
#include "math.h"
 
QString version = "2.0";
 
QColor couleur_ecran = QColor::fromRgb(3, 41, 11, 255);
QColor couleur_ligne = QColor::fromRgb(167, 2, 0, 255);
QColor couleur_trace = QColor::fromRgb(0, 210, 0, 255);
QColor couleur_texte = QColor::fromRgb(255, 255, 0, 255);
QColor couleur_signature = QColor::fromRgb(255, 255, 0, 255);
 
QPen pen_ligne(couleur_ligne, 1, Qt::SolidLine);
QPen pen_trace(couleur_trace, 1, Qt::SolidLine);
QPen pen_reticule(couleur_ligne, 1, Qt::SolidLine);
 
 
typedef float complexe[2];
 
//complexe ech[84]; // nb_ech echantillons  1024    SUPPRIME ! et oui !! j'utilise le même tableau pour tout !
complexe tab_X[84];   // nb_ech valeurs traitées  1024
complexe tab_W[42];    // nb_ech/2              512
 
 
 
float frequence1=400; // Hz
float frequence2=400; // Hz
uint nb_etapes=6; //8
uint nb_ech = pow (2,nb_etapes); // nombre d'échantillons = 2 puissance(nb_etapes)
float frq_echantillonnage=8*frequence1; // Hz
float Te; // période d'échantillonage
 
//int pas_balayage=40;
int hamming = true;
int second_Frq = false;
int bz = true;
uint forme;
 
 
 
 
 
MainWindow::MainWindow(QWidget *parent) :
    QMainWindow(parent)
 
{
    setupUi(this);
    setWindowTitle("Transformee de Fourier Rapide FFT (no objects) - version " + version);
 
    scene = new QGraphicsScene(this);
  //  scene->setBackgroundBrush(couleur_ecran);
    graphicsView1->setScene(scene);
 //   graphicsView1->setGeometry(0,0,530,430);
    graphicsView1->setGeometry(140,30,84,48);
 
    groupe_reticule = new QGraphicsItemGroup();
    scene->addItem(groupe_reticule);
    groupe_trace = new QGraphicsItemGroup();
    scene->addItem(groupe_trace);
 
    effacer_trace();
    RAZ_tableau_echantillons();
    remplir_tableau_echantillons();
 //   tracer_signal();
    calcul_tableau_W();
  //  bit_reverse_tableau_X();
    calcul_fourier();
    tracer_tableau_valeurs();
 
}
 
 
 
 
void MainWindow::effacer_graduation()
{
    foreach( QGraphicsItem *item, scene->items( groupe_reticule->boundingRect() ) )
    {
        if( item->group() == groupe_reticule )  { delete item; }
    }
}
 
 
void MainWindow::effacer_trace()
{
    foreach( QGraphicsItem *item, scene->items( groupe_trace->boundingRect() ) )
    {
        if( item->group() == groupe_trace )  { delete item; }
    }
}
 
 
MainWindow::~MainWindow()
{
 
}
 
 
void MainWindow::tracer_graduations()
{
    float i,x,y;
    float nb_grad_max;
    float intervalle; // séparant les graduations
    float fi;
    QString sti;
    uint decal = 3; // leger decallage vertical pour ne pas masquer la trace
 
    rectangle = new QGraphicsRectItem(0,decal, 511, 400);
    rectangle->setPen(couleur_ligne);
    groupe_trace->addToGroup(rectangle);
 
// lignes verticales
// ici calcul de l'intervalle (en pixels) qui correspond exactement à 100Hz (en fréquence)
 
    intervalle = 2000.0 * nb_ech  * 2.0 / (frq_echantillonnage);
    nb_grad_max = 512 / intervalle;
 
    for (i=0; i<=nb_grad_max; i++)
    {
        x = intervalle * i;
 
        ligne1 = new QGraphicsLineItem(x, decal, x, 400 + decal);
        ligne1->setPen(pen_reticule);
        groupe_reticule->addToGroup(ligne1);
 
        fi = x * frq_echantillonnage  / nb_ech /4;
 
        sti.setNum(fi);
        texte_frq = new QGraphicsTextItem(sti);
        texte_frq->setDefaultTextColor(couleur_texte);
        texte_frq->setPos(x,380);
        if(x<480) // évite que l'écriture ne déborde du cadre a droite
        {
            groupe_reticule->addToGroup(texte_frq);
        }
    }
 
// lignes horizontales
    intervalle = 50;
    nb_grad_max = 400 / intervalle;
 
    for (i=0; i<=nb_grad_max; i++)
    {
        y = 400 + decal - intervalle * i;
 
        ligne1 = new QGraphicsLineItem(0, y, 511, y);
        ligne1->setPen(pen_reticule);
        groupe_reticule->addToGroup(ligne1);
    }
 
    texte_frq = new QGraphicsTextItem("Silicium628");
    texte_frq->setDefaultTextColor(couleur_signature);
    texte_frq->setPos(5,5);
    groupe_reticule->addToGroup(texte_frq);
 
    scene->addItem(groupe_reticule);
}
 
 
void MainWindow::tracer_signal()
{
    uint offset_y = 24;
    float echelle =168/nb_ech;
    uint n;
    float x, y, memo_x, memo_y;
 
    segment_trace = new QGraphicsLineItem(0, offset_y, 512, offset_y);
    segment_trace->setPen(couleur_ligne);
    groupe_trace->addToGroup(segment_trace);
 
    effacer_trace();
 
    x=0;
    y=offset_y;
    for(n=0; n<64; n++)
    {
        memo_x = x;
        memo_y = y;
        x =  n;
 
        y = offset_y - 10 * tab_X[n][1]; // partie réelle
 
        segment_trace = new QGraphicsLineItem(memo_x ,memo_y, x, y);
        segment_trace->setPen(pen_trace);
        groupe_trace->addToGroup(segment_trace);
    }
    scene->addItem(groupe_trace);
}
 
 
 
 
void MainWindow::tracer_tableau_valeurs()
{
    uint offset_y = 46; // 600
    float echelle =168/nb_ech;
    uint n;
    float x, y, memo_x, memo_y;
 
    int z;
 
    if (bz) {z=4;} else {z=1;}
 
 
 
    segment_trace = new QGraphicsLineItem(0, offset_y, 84, offset_y);
    segment_trace->setPen(couleur_ligne);
    groupe_trace->addToGroup(segment_trace);
 
 
 
    x=0;
    y=offset_y;
    for(n=0; n<nb_ech/2; n++)
    {
        memo_x = x;
        memo_y = y;
        x =  echelle * n;
 
     //  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
 
//version sans la racine carrée sqrt (qui pose des problèmes sur ATmega avr-gcc):
       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
 
 
 
        segment_trace = new QGraphicsLineItem(memo_x ,memo_y, x, y);
        segment_trace->setPen(pen_trace);
        groupe_trace->addToGroup(segment_trace);
    }
    scene->addItem(groupe_trace);
}
 
 
 
void MainWindow::RAZ_tableau_echantillons()
{
   uint n;
   for(n=0; n < nb_ech; n++)
   {
       tab_X[n][1] = 0; // partie reelle
       tab_X[n][2] = 0; // partie imaginaire
   }
}
 
 
 
uint bit_reversal(uint num, uint nb_bits)
{
    uint r = 0, i, s;
    if ( num > 1<< nb_bits) { return 0; }
 
    for (i=0; i<nb_bits; i++)
    {
        s = (num & (1 << i));
        if(s) { r |= (1 << ((nb_bits - 1) - i)); }
    }
    return r;
}
 
 
 
void MainWindow::remplir_tableau_echantillons()
{
// création du signal à analyser
    uint n, r, n_max;
    float x;
    Te = 1.0 / frq_echantillonnage;
 
    if (bz==true) {n_max = nb_ech /4; } else {n_max = nb_ech; }
 
    for(n=0; n < n_max; n++) // nb d'échantillons
    {
        x=  cos(2.0 * M_PI * frequence1 * n * Te );
        if (second_Frq) {x+=cos(2 * M_PI * frequence2 * n * Te); }
        if (hamming)  { x = x * (1- cos (2* M_PI * n / n_max ));}
        r=bit_reversal(n,nb_etapes);
        tab_X[r][1] = x; // partie reelle
        tab_X[r][2] = 0; // partie imaginaire
     }
}
 
 
 
 
/*
void MainWindow::bit_reverse_tableau_X()
{
 
// recopie les échantillons en les classant dans l'ordre 'bit reverse'
    uint n,r;
 
    for(n=0; n < nb_ech; n++) // nb d'échantillons
    {
        r=bit_reversal(n,nb_etapes);
        tab_X[n][1] = ech[r][1];
        tab_X[n][2] = ech[r][2];
    }
}
 
*/
 
void MainWindow::calcul_tableau_W()
{
// calcul et memorisation dans un tableau des twiddle factors
    uint n;
    float x;
 
    for(n=0; n<(nb_ech/2-1); n++)
    {
        x=2.0*M_PI * n / nb_ech;
 
        tab_W[n][1] = cos(x);  // partie reelle
        tab_W[n][2] = -sin(x); // partie imaginaire
    }
}
 
 
 
void MainWindow::calcul_fourier()
{
    complexe produit; // voir le type "complexe" déclaré plus haut
 
    uint decal =2;
 
    uint etape, e1, e2, li, w, ci;
    uint li_max=nb_ech;
 
    e2=1;
 
    for (etape=1; etape<=nb_etapes; etape++)
    {
        e1=e2; //(e1 evite d'effectuer des divisions e2/2 plus bas)
        e2=e2+e2;
 
        for (li=0; li<li_max; li+=1)
        {
            ci=li & (e2-1);  // ET bit à bit
            if (ci>(e1-1))
            {
                w=li_max/e2*(li & (e1 -1));  // ET bit à bit  calcul du numéro du facteur de tripatouillages W
 
                produit[1]=tab_W[w][1]*tab_X[li][1]-tab_W[w][2]*tab_X[li][2];
                produit[2]=tab_W[w][1]*tab_X[li][2]+tab_W[w][2]*tab_X[li][1];
 
                tab_X[li][1]=tab_X[li-e1][1]-produit[1];
                tab_X[li][2]=tab_X[li-e1][2]-produit[2];
 
                tab_X[li-e1][1]=tab_X[li-e1][1]+produit[1];
                tab_X[li-e1][2]=tab_X[li-e1][2]+produit[2];
            }
        }
    }
}
 
 
 
 
 
void MainWindow::on_spinBox_frq_valueChanged(int arg1) // SELECTEUR DE FREQUENCE 1
{
   effacer_trace();
   frequence1 = arg1;
   RAZ_tableau_echantillons();
   remplir_tableau_echantillons();
//   tracer_signal();
   calcul_tableau_W();
 //  bit_reverse_tableau_X();
   calcul_fourier();
    tracer_tableau_valeurs();
}
 
 
void MainWindow::on_spinBox_frq_2_valueChanged(int arg1) // SELECTEUR DE FREQUENCE 2
{
    effacer_trace();
    frequence2 = arg1;
    RAZ_tableau_echantillons();
    remplir_tableau_echantillons();
 //   tracer_signal();
    calcul_tableau_W();
  //  bit_reverse_tableau_X();
    calcul_fourier();
    tracer_tableau_valeurs();
}
 
 
 
void MainWindow::on_checkBox1_toggled(bool checked) // HAMMING
{
   effacer_trace();
   if (checked) {hamming = true;} else {hamming = false;}
   RAZ_tableau_echantillons();
   remplir_tableau_echantillons();
//   tracer_signal();
   calcul_tableau_W();
 //  bit_reverse_tableau_X();
   calcul_fourier();
   tracer_tableau_valeurs();
 
}
 
 
void MainWindow::on_checkBox2_toggled(bool checked)  // SECONDE FREQUENCE
{
    effacer_trace();
    if (checked) {second_Frq = true;} else {second_Frq = false;}
    RAZ_tableau_echantillons();
    remplir_tableau_echantillons();
 //   tracer_signal();
    calcul_tableau_W();
  //  bit_reverse_tableau_X();
    calcul_fourier();
    tracer_tableau_valeurs();
}
 
 
void MainWindow::on_checkBox3_toggled(bool checked)
{
    effacer_trace();
    if (checked) {bz = true;} else {bz = false;}
    RAZ_tableau_echantillons();
    remplir_tableau_echantillons();
 //   tracer_signal();
    calcul_tableau_W();
  //  bit_reverse_tableau_X();
    calcul_fourier();
    tracer_tableau_valeurs();
}
 
 
 
 
 
 

/***************************
FFT.cpp
par Silicium628
logiciel Libre Open Source
****************************/
 
#include <avr/io.h>
#include <stdint.h> 
#include <stdlib.h>
#include "timeout.h"   // comprend : #define F_CPU 16000000UL  // 16 MHz
#include "LCD_5110.cpp"
#include <math.h>
 
 
const char *version = "2.5";  // ( signal entrant 10kHz à 90kHz )
 
typedef float complexe[2];
 
complexe tab_X[64];  // nb_ech valeurs traitées
complexe tab_W[32];  // nb_ech/2
 
float frequence1; // Hz
float frq_echantillonnage; // Hz
 
const uint8_t nb_etapes=6; 
uint8_t nb_ech = pow (2,nb_etapes); // nombre d'échantillons = 2 puissance(nb_etapes)
 
char hamming;
char bourrage;
char mode_affi;
 
//const uint8_t hamming = false;
//const uint8_t bourrage = true;
 
 
void RAZ_tableau_echantillons()
{
   uint8_t n;
   for(n=0; n < nb_ech; n++)
   {
       tab_X[n][1] = 0; // partie reelle
       tab_X[n][2] = 0; // partie imaginaire
   }
}
 
 
void init_ports (void)  // ports perso
// 0 = entree, 1=sortie   ;   les 1 sur les pins en entrees activent les R de Pull Up (tirage à VCC)
{
  PORTA = 0b00000111;
  DDRA  = 0b01111000;  //PA7 en entree (ADC7 - entrée signal); PA0, PA1, PA2 jumpers à GND.
 
  PORTB = 0b00000000;
  DDRB |= 0b11111111;
 
  PORTC = 0b00000000;
  DDRC =   0b11111111;
 
  PORTD = 0b00000000; 
  DDRD =  0b11111111;
}
 
/*
The ADC module contains a prescaler, which generates an acceptable ADC clock fre-
quency from any CPU frequency above 100 kHz. The prescaling is set by the ADPS bits
in ADCSRA
*/
 
void InitADC (void)
{
//  SFIOR  = 0b10000000; // ADTS2,1,0 = 100 -> CAN declenché par Timer0 Overflow - p:218
  ADMUX  = 0b01000111;   //Bit 7:6 = ADC Reference Selection = 01 -> AVCC ; //Bits 0:4 - Analog Channel Selection Bits et gain éventuel - p214
// ici Select pin ADC7 using MUX  avec ref tension interne = AREF externe - pas de gain.  
 
/** CHOISIR UNE SEULE LIGNE CI-DESSOUS pour ADCSRA */
 
// parametrage gamme 8-90 kHz
  ADCSRA = 0b10000010;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/4 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 4-56 kHz
//  ADCSRA = 0b10000011;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/8 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 2-31 kHz
//  ADCSRA = 0b10000100;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/16 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 1-15 kHz
//  ADCSRA = 0b10000101;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 300 Hz-8000 Hz
//  ADCSRA = 0b10000110;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 200 Hz- 4300 Hz
//  ADCSRA = 0b10000111;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
}
 
 
 
 
 
uint8_t bit_reversal(uint8_t num, uint8_t nb_bits)
{
    uint8_t r = 0, i, s;
    if ( num > 1<< nb_bits) { return 0; }
 
    for (i=0; i<nb_bits; i++)
    {
        s = (num & (1 << i));
        if(s) { r |= (1 << ((nb_bits - 1) - i)); }
    }
    return r;
}
 
 
 
/***********************************************************************
Remarques :  
- 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,  
(par exemple en 2014 une Nvidia GeForce GTX680 dont le processeur tourne à à 2GHz- 2Go - DDR5 ou une Nvidia Quadro 6000 )
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.
***********************************************************************/  
 
void acquisition()
{  
  uint16_t acqH;
  uint16_t tab_acq[64];
  uint8_t n, r, n_max;
  float Te; // période d'échantillonage
  float x;
  Te = 1.0 / frq_echantillonnage;
  uint16_t valeur_lue;  
 
  if (bourrage) {n_max = nb_ech /4; } else {n_max = nb_ech; }
    //n_max  = 16;
 
    for(n=0; n < n_max-1; n++)
    {    
        ADCSRA |= _BV(ADSC);    //Start conversion - resolution 10bits
        while (ADCSRA & _BV(ADSC) ) {;}    // attend la fin de la converstion
        tab_acq[n] = ADCL;
        acqH = ADCH; // passage à 16 bits pour pouvoir décaller
        tab_acq[n] += acqH << 8;
    }    
    for(n=0; n < n_max-1; n++)
    {
      valeur_lue =tab_acq[n];
      if (hamming)  { x = x * (1- cos (2* M_PI * n / n_max ));} 
      x= (valeur_lue - 512.0) / 100.0;    
      r=bit_reversal(n,nb_etapes);
      tab_X[r][1] = x; // partie reelle      
      tab_X[r][2] = 0; // partie imaginaire
    }
}
 
 
 
void remplir_tableau_echantillons()
{
// création du signal à analyser
    uint8_t n, r, n_max;
    float Te; // période d'échantillonage
    float x;
    Te = 1.0 / frq_echantillonnage;
 
    if (bourrage==true) {n_max = nb_ech /4; } else {n_max = nb_ech; }
 
    for(n=0; n < n_max-1; n++) // nb d'échantillons
    {
        x=  cos(2.0 * M_PI * frequence1 * n * Te );
        if (hamming)  { x = x * (1- cos (2* M_PI * n / n_max ));}
        r=bit_reversal(n,nb_etapes);
        tab_X[r][1] = x; // partie reelle
        tab_X[r][2] = 0; // partie imaginaire
     }
}
 
 
 
void calcul_tableau_W()
{
// calcul et memorisation dans un tableau des twiddle factors
    uint8_t n;
    float x;
 
    for(n=0; n<(nb_ech/2-1); n++)
    {
        x=2.0*M_PI * n / nb_ech;
 
        tab_W[n][1] = cos(x);  // partie reelle
        tab_W[n][2] = -sin(x); // partie imaginaire
    }
}
 
 
 
void calcul_fourier()
{
    complexe produit; // voir le type "complexe" déclaré plus haut
 
    uint8_t decal =2;
 
    uint8_t etape, e1, e2, li, w, ci;
    uint8_t li_max=nb_ech;
 
    e2=1;
 
    for (etape=1; etape<=nb_etapes; etape++)
    {
        e1=e2; //(e1 evite d'effectuer des divisions e2/2 plus bas)
        e2=e2+e2;
 
        for (li=0; li<li_max; li+=1)
        {
            ci=li & (e2-1);  // ET bit à bit
            if (ci>(e1-1))
            {
                w=li_max/e2*(li & (e1 -1));  // ET bit à bit  calcul du numéro du facteur de tripatouillages W
 
                produit[1]=tab_W[w][1]*tab_X[li][1]-tab_W[w][2]*tab_X[li][2];
                produit[2]=tab_W[w][1]*tab_X[li][2]+tab_W[w][2]*tab_X[li][1];
 
                tab_X[li][1]=tab_X[li-e1][1]-produit[1];
                tab_X[li][2]=tab_X[li-e1][2]-produit[2];
 
                tab_X[li-e1][1]=tab_X[li-e1][1]+produit[1];
                tab_X[li-e1][2]=tab_X[li-e1][2]+produit[2];
            }
        }
    }
}
 
 
 
 
 
void tracer_signal()
{
    float a, S;
    uint8_t  r,  memo_x, memo_y, n, x, y;
 
    x=0;
    y=24;
    for(n=0; n<32; n++)
    {
        memo_x = x;
        memo_y = y;
        x = 2*n;
        r=bit_reversal(n, nb_etapes);
        if (hamming ==0) {a = 10.0;} else {a=5.0;}
        S = 20.0 - a * tab_X[r][1]; // partie réelle
        if ( (abs(S) )< 48 )  {y= lrint(S);} else {y=1;}
 
/** choisir une des trois lignes ci-dessous **/    
        //buffer_drawPixel(x, y);    
        //buffer_trace_baton_V(x, y);
        buffer_trace_segment( memo_x, memo_y, x, y);
    }
    affi_buffer(84);
}
 
 
 
void tracer_tableau_valeurs()
{
  uint8_t offset_y = 0;
  uint8_t n ;
  uint8_t memo_x, memo_y, x, y;
  float A;
 
  uint8_t z;
 
    if (bourrage) {z=15;} else {z=1;}
 
    x=0;
    y=offset_y;
    for(n=0; n<32; n++)
    {
        memo_x = x;
        memo_y = y;
 
        x = 2*n;
        //x = n;
 
        A = z *( tab_X[n][1] * tab_X[n][1] + tab_X[n][2] * tab_X[n][2]);
        A = sqrt(A);
        if ( (abs(A) )< 48 )  {y= lrint(A);} else {y=1;}
 
 
/** choisir une des trois lignes ci-dessous **/    
      //  buffer_drawPixel(x, 47-y);    // la plus rapide, mais peu visible...
      //  buffer_trace_baton_V(x, 1 + y);
        buffer_trace_segment( memo_x, memo_y, x, y);
    }
   affi_buffer(84);
}
 
 
 
void test_out(uint8_t etat)
{
  if (etat == 1) { PORTD |= 0b00000010;} else {PORTD &= 0b11111101;}
}
 
 
 
int main(void)
{
 
  init_ports();
  InitADC();
//  InitINTs();
 
  _delay_ms(100);
 
// define the 5 LCD Data pins: SCE, RST, DC, DATA, CLK
    lcd_init(&PORTB, PB1, &PORTB, PB0, &PORTB, PB2, &PORTB, PB3, &PORTB, PB4);
 
  // frequence1 = 100.0; // Hz
   frq_echantillonnage=3200.0; // Hz
 
 
    lcd_goto_LC(0,1);
    lcd_puts2("FTT ATmega32");
 
    _delay_ms(1000);
    lcd_clear();
 
    calcul_tableau_W();
 
    for(;;)
    {
      //frequence1 = 310.0; // Hz
 
      for(uint8_t n=1; n<23; n++)
      {
        mode_affi = (PINA & 0b00000001) == 0;
        hamming   = (PINA & 0b00000010) != 0;
        bourrage  = (PINA & 0b00000100) != 0;
 
        //frequence1 += 30;
 
        RAZ_tableau_echantillons();
        acquisition();
      //  remplir_tableau_echantillons();      
 
        if (mode_affi) 
        {
          buffer_RAZ();
          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
          // _delay_ms(1);
        } 
 
        if (! mode_affi) 
        {  
 
          calcul_fourier();  // 100..120ms
          buffer_RAZ();    // 1.4ms
 
          test_out(1);
            tracer_tableau_valeurs(); // 250 ms
          test_out(0);
        }
        lcd_goto_LC(5,0);
 
 
        if (! mode_affi) {lcd_puts2(" FFT   ");} else {lcd_puts2(" Signal");}  // 7 ms
 
      }
   }
}
 

/***************************
FFT.cpp
par Silicium628
logiciel Libre Open Source
****************************/
 
#include <avr/io.h>
#include <stdint.h> 
#include <stdlib.h>
#include "timeout.h"   // comprend : #define F_CPU 16000000UL  // 16 MHz
#include "LCD_5110-628.cpp"
#include <math.h>
 
 
const char *version = "4.2";  // ( signal entrant 10kHz à 90kHz )
 
//typedef float complexe[2];
typedef int16_t complexe_f[2];  // entier signé pour calculs (rapides) en virgule fixe 
 
//complexe tab_X[64];  // nb_ech valeurs traitées
//complexe tab_W[32];  // nb_ech/2
complexe_f tab_X[64];  // nb_ech valeurs traitées
complexe_f tab_W[32];  // nb_ech/2
 
float frequence1; // Hz
float frq_echantillonnage; // Hz
 
const uint8_t nb_etapes=6; 
uint8_t nb_ech = pow (2,nb_etapes); // nombre d'échantillons = 2 puissance(nb_etapes)
 
char hamming;
char bourrage;
char mode_affi;
 
 
void RAZ_tableau_echantillons()
{
   uint8_t n;
   for(n=0; n < nb_ech; n++)
   {
       tab_X[n][1] = 0; // partie reelle
       tab_X[n][2] = 0; // partie imaginaire
   }
}
 
 
void init_ports (void)  // ports perso
// 0 = entree, 1=sortie   ;   les 1 sur les pins en entrees activent les R de Pull Up (tirage à VCC)
{
  PORTA = 0b00000111;
  DDRA  = 0b01111000;  //PA7 en entree (ADC7 - entrée signal); PA0, PA1, PA2 jumpers à GND.
 
  PORTB = 0b00000000;
  DDRB |= 0b11111111;
 
  PORTC = 0b00000000;
  DDRC =   0b11111111;
 
  PORTD = 0b00000000; 
  DDRD =  0b11111111;
}
 
/*
The ADC module contains a prescaler, which generates an acceptable ADC clock fre-
quency from any CPU frequency above 100 kHz. The prescaling is set by the ADPS bits
in ADCSRA
*/
 
void InitADC (void)
{
//  SFIOR  = 0b10000000; // ADTS2,1,0 = 100 -> CAN declenché par Timer0 Overflow - p:218
  ADMUX  = 0b01000111;   //Bit 7:6 = ADC Reference Selection = 01 -> AVCC ; //Bits 0:4 - Analog Channel Selection Bits et gain éventuel - p214
// ici Select pin ADC7 using MUX  avec ref tension interne = AREF externe - pas de gain.  
 
/** CHOISIR UNE SEULE LIGNE CI-DESSOUS pour ADCSRA */
 
// parametrage gamme 8-90 kHz
//  ADCSRA = 0b10000010;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/4 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 4-56 kHz
//  ADCSRA = 0b10000011;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/8 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 2-31 kHz
//  ADCSRA = 0b10000100;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/16 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 1-15 kHz
//  ADCSRA = 0b10000101;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 300 Hz-8000 Hz
//  ADCSRA = 0b10000110;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
// parametrage gamme 200 Hz- 4300 Hz
  ADCSRA = 0b10000111;  // bit7: Activate ADC; ADPS2..0->Prescaler=f/32 ; bit5 = ADATE (ADC Auto Trigger Enable)->non utilisé ici - p:216 et 218
 
}
 
 
 
 
 
uint8_t bit_reversal(uint8_t num, uint8_t nb_bits)
{
    uint8_t r = 0, i, s;
    if ( num > 1<< nb_bits) { return 0; }
 
    for (i=0; i<nb_bits; i++)
    {
        s = (num & (1 << i));
        if(s) { r |= (1 << ((nb_bits - 1) - i)); }
    }
    return r;
}
 
 
 
/***********************************************************************
Remarques :  
- 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,  
(par exemple en 2014 une Nvidia GeForce GTX680 dont le processeur tourne à à 2GHz- 2Go - DDR5 ou une Nvidia Quadro 6000 )
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.
***********************************************************************/  
 
void acquisition()
{  
  uint16_t acqH;
  uint16_t tab_acq[64];
  uint8_t n, r, n_max;
  float x;
  uint16_t valeur_lue;  
 
  if (bourrage) {n_max = nb_ech /4; } else {n_max = nb_ech; }
    //n_max  = 16;
 
    for(n=0; n < n_max-1; n++)
    {    
        ADCSRA |= _BV(ADSC);    //Start conversion - resolution 10bits
        while (ADCSRA & _BV(ADSC) ) {;}    // attend la fin de la converstion
        tab_acq[n] = ADCL;
        acqH = ADCH; // passage à 16 bits pour pouvoir décaller
        tab_acq[n] += acqH << 8;
    }    
    for(n=0; n < n_max-1; n++)
    {
      valeur_lue =tab_acq[n];
      x= (valeur_lue - 512.0) / 100.0;   
      if (hamming)  { x = x * (1- cos (2* M_PI * n / n_max ));}       
      r=bit_reversal(n,nb_etapes);
      tab_X[r][1] = x * 16; // partie reelle  - * 16 pour obtenir un nombre en virgule fixe
      tab_X[r][2] = 0; // partie imaginaire
    }
}
 
 
 
void remplir_tableau_echantillons()
{
// création du signal à analyser
    uint8_t n, r, n_max;
    float Te; // période d'échantillonage
    float x;
    Te = 1.0 / frq_echantillonnage;
 
    if (bourrage==true) {n_max = nb_ech /2; } else {n_max = nb_ech; }
 
    for(n=0; n < n_max-1; n++) // nb d'échantillons
    {
        x=  cos(2.0 * M_PI * frequence1 * n * Te );
        if (hamming)  { x = x * (1- cos (2* M_PI * n / n_max ));}
        r=bit_reversal(n,nb_etapes);
        tab_X[r][1] = x; // partie reelle
        tab_X[r][2] = 0; // partie imaginaire
     }
}
 
 
 
void calcul_tableau_W()
{
// calcul et memorisation dans un tableau des twiddle factors
    uint8_t n;
    //float x;
    float x, cosx, sinx;
 
    for(n=0; n<(nb_ech/2-1); n++)
    {
        x=2.0*M_PI * n / nb_ech;
				cosx =  16 * cos(x);
				sinx = -16 * sin(x);
 
        tab_W[n][1] =  round(cosx);  // partie reelle
        tab_W[n][2] =  round(sinx);  // partie imaginaire
    }
}
 
 
 
void calcul_fourier()
{
    complexe_f produit; // voir le type "complexe" déclaré plus haut
 
    uint8_t etape, e1, e2, li, w, ci;
    uint8_t li_max=nb_ech;
 
    e2=1;
 
    for (etape=1; etape<=nb_etapes; etape++)
    {
        e1=e2; //(e1 evite d'effectuer des divisions e2/2 plus bas)
        e2=e2+e2;
 
        for (li=0; li<li_max; li+=1)
        {
            ci=li & (e2-1);  // ET bit à bit
            if (ci>(e1-1))
            {
                w=li_max/e2*(li & (e1 -1));  // ET bit à bit  calcul du numéro du facteur de tripatouillages W
 
                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
                produit[2]=((tab_W[w][1]*tab_X[li][2]) + (tab_W[w][2]*tab_X[li][1])) >> 4;
 
                tab_X[li][1]=tab_X[li-e1][1]-produit[1];
                tab_X[li][2]=tab_X[li-e1][2]-produit[2];
 
                tab_X[li-e1][1]=tab_X[li-e1][1]+produit[1];
                tab_X[li-e1][2]=tab_X[li-e1][2]+produit[2];
            }
        }
    }
}
 
 
 
void tracer_signal()
{
    int32_t a, S;
    uint8_t  r,  memo_x, memo_y, n, x, y;
 
    x=0;
    y=24;
    for(n=0; n<32; n++)
    {
        memo_x = x;
        memo_y = y;
        x = 2*n;
        r=bit_reversal(n, nb_etapes);
        if (hamming ==0) {a = 10;} else {a=5;}
        S = (a * tab_X[r][1]); // partie réelle
        S =  S >> 4;
        S = 20 - S;
        if ( (abs(S) )< 48 )  {y= lrint(S);} else {y=1;}
 
/** choisir une des trois lignes ci-dessous **/    
        //buffer_drawPixel(x, y);    
        //buffer_trace_baton_V(x, y);
        buffer_trace_segment( memo_x, memo_y, x, y);
    }
    affi_buffer(84);
}
 
 
 
void tracer_tableau_valeurs()
{
  uint8_t offset_y = 0;
  uint8_t n ;
  uint8_t memo_x, memo_y, x, y;
  uint32_t A;
 
 
  uint8_t z;
 
    if (bourrage) {z=15;} else {z=5;}
 
    x=0;
    y=offset_y;
 
    for(n=0; n<31; n++)  // batons =3.5 ms ; segments = 86 ms
    {
        memo_x = x;
        memo_y = y;
 
        //x = n;        
        x = 2*n; // largeur de la trace = 4 px (1 baton de 3px + 1 séparation de 1px)
 
        A = ((tab_X[n][1] * tab_X[n][1]) + (tab_X[n][2] * tab_X[n][2]));
        A >>=8; // div/256; // 256 = 16 * 16 (compense la multiplication en virgule fixe puis retour à une valeur ordinaire)
        A *= z;
        A = sqrt(A);
 
        if ( (abs(A) )< 48 )  {y= lrint(A);} else {y=47;}
 
/** choisir une des trois lignes ci-dessous; ATTENTION : les rapidités ne sont pas les mêmes !! **/    
      //  buffer_drawPixel(x, 47-y);    // la plus rapide, mais peu visible...
        buffer_trace_baton_V(x, 1 + y);
       // buffer_trace_segment( memo_x, memo_y, x, y);
    }
 
		affi_buffer(84); // 12.6 ms
 
}
 
 
 
 
 
int main(void)
{
 
  init_ports();
  InitADC();
//  InitINTs();
 
  _delay_ms(100);
  lcd_init();
 
  // frequence1 = 100.0; // Hz
  // frq_echantillonnage=3200.0; // Hz
 
 
  lcd_goto_LC(0,1);
  lcd_puts2("FTT ATmega32");
 
  _delay_ms(1000);
  lcd_clear();
 
  calcul_tableau_W();
 
  lcd_goto_LC(5,0);  
  if (! mode_affi) {lcd_puts2("TF Silicium628  ");} else {lcd_puts2(" Signal");}  // 7 ms  
 
    for(;;)
    {
      //frequence1 = 310.0; // Hz
 
      for(uint8_t n=1; n<23; n++)
      {
        mode_affi = (PINA & 0b00000001) == 0;
        hamming   = (PINA & 0b00000010) != 0;
        bourrage  = (PINA & 0b00000100) != 0;
 
        //frequence1 += 30;
 
 
				RAZ_tableau_echantillons();
				acquisition();  // 7ms
 
      //  remplir_tableau_echantillons();      
 
        if (mode_affi) 
        {
          buffer_RAZ();
          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
        } 
 
        if (! mode_affi) 
        {  
          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)
          buffer_RAZ();    // 1.4 ms
 
				  tracer_tableau_valeurs(); // 16ms 
        }
 
      }
   }
}