Analyse fréquentielle #3: Implémentation et Test de la DFT

Objectifs

Savoir implémenter la DFT en C/Arduino
Test de la DFT sur Arduino Mega/Due
Analyse des performances temporelles du code sur Mega/Due
Analyse de la DFT d’un signal sinusoïdal
Analyse de la DFT d’une entrée réelle
Prendre consciente du problème du sur-échantillonnage
Savoir extraire les amplitudes de la DFT
Comprendre la notion de la fréquence
Savoir définir la fréquence d’échantillonnage
Savoir la formule de la transformée de Fourier discrète (DFT)
Comprendre la notion du spectre
Etc.

Voir le tuto pour plus de détails

Lecture des données

Données synthétiques

La première étape de l’analyse spectrale consiste à récolter les données d’un capteur ou ensemble des capteurs. Le choix de la résolution binaire (Nb) du convertisseur A/N ainsi la fréquence maximale de la bande analogique est cruciale. On choisit le (Nb) en fonction du niveau du signal par rapport au bruit et la fréquence d’échantillonnage (Fs) en fonction de la rapidité du phénomène physique. La fréquence Fs permet de limiter la fenêtre d’observation du signal dans le domaine fréquentiel. En effet, on ne peut observer que les fréquences dans la plage [-Fs,2, Fs/2] environ. En pratique, on se limite aux à la bande 0 et Fs/2. Le pas dans le domaine fréquentiel (ou la résolution fréquentielle), elle est définie par le nombre d’échantillons (N). Elle est égale à Fs/N, autrement dit le nombre d’acquisitions (ou la fenêtre temporelle du signal). Ci-dessous l’exemple de synthèse d’un signal sinusoïdal.

...
  float f0=10.0, t_i=0.0, ts=0.0;
  float A0=1.0, DC=1.0; 
  for(int i=0; i<N; i++)
  {
    t_i=(float)i/(float)N; 
    Data[i]=DC + A0*sin(2.0*PI*f0*t_i); // + A0*sin(2.0*PI*0.5*f0*t_i)+ A0*sin(2.0*PI*(1.0/3.0)*f0*t_i); // DC, f0, f0/2
    ts=(float)1.0/(float)N; 
  }
  Fs=1/ts; 
...

Données réelles

Le code ci-dessous permet de lire N échantillons de l’entrée A0. La période d’échantillonnage est défini par la variable Ts_us (nombre de cycles de 1 µs). Elle est ensuite remesurée afin d’obtenir une meilleure précision sur la mesure des harmoniques (voir le tuto pour plus de détails).

...
  Ts=(float)Ts_us;
  for(int i=0; i<N; i++)
  {
    T1=micros();
    Data[i]=(float)analogRead(A0);
    Data[i]*=5.0/1023.0; 
    delayMicroseconds(Ts_us);
    Ts+=(float)(micros()-T1); 
    Ts=Ts/2.0;
  }
  Fs=1E6/Ts; 
  Serial.println(Fs);return; 
...

Calcul de la DFT

float DFT(float *data, int Ndata, float *real, float *imag, float fs)

La DFT() retourne le spectre d’un signal. Elle prend en entrée trois paramètres :

*data: Tableau contenant les données du signal de taille (N)
Ndata: Nombre des échantillons du tableau (N)
fs: la fréquence d’échardonnage en Hz

Elle retourne trois sorties :

*real: tableau de la partie réelle du spectre de taille N/2
*imag: tableau de la partie réelle du spectre de taille N/2
float: la résolution fréquentielle

Programme complet

#define Sortie  8       // Non utilisée 
#define N       128     // Données N=2^n, TFD Nf=2^(n-1)
#define Ts_us   1000    // Ts=1/Fs: Approchée (voir le programme) 

/*
 * Objectifs:               
 *  1. Savoir implémenter la DFT
 *  2. Test de la DFT avec un Sinus 
 *  3. Test de la DFT avec ADC (Données réelles) 
 *  4. Analyse de la vitesse Mega/Due 
 */

bool State=false; 
unsigned long T1=0, T2=0; 
float Ts=0.0, Fs=0.0; 
float Data[N], Nfft=N>>1; 
float dF=0.0; 
float RealFFT[N>>1], ImagFFT[N>>1]; 
float sMoy=0.0; 
float Am, ReIm[2];
int n=0; 

void setup()
{
  // Affichage 
  Serial.begin(115200); 
}

void loop()
{ 
  // 1. Lecture des données/ Test de la DFT 
  /*
  float f0=10.0, t_i=0.0, ts=0.0;
  float A0=1.0, DC=1.0; 
  for(int i=0; i<N; i++)
  {
    t_i=(float)i/(float)N; 
    Data[i]=DC + A0*sin(2.0*PI*f0*t_i); // + A0*sin(2.0*PI*0.5*f0*t_i)+ A0*sin(2.0*PI*(1.0/3.0)*f0*t_i); // DC, f0, f0/2
    ts=(float)1.0/(float)N; 
  }
  Fs=1/ts; 
  */
  
  // 1. Lecture des données 
  Ts=(float)Ts_us;
  for(int i=0; i<N; i++)
  {
    T1=micros();
    Data[i]=(float)analogRead(A0);
    Data[i]*=5.0/1023.0; 
    delayMicroseconds(Ts_us);
    Ts+=(float)(micros()-T1); 
    Ts=Ts/2.0;
  }
  Fs=1E6/Ts; 
  //Serial.println(Fs);return; 
  

  // 2. Calcul de la FFT
  T1=micros();
  dF=DFT(Data, N, RealFFT, ImagFFT, Fs);
  T2=micros()-T1;
  Serial.print("Temps FFT(s): ");Serial.println((float)T2/1E6); 
  //Serial.print("Résolution Fréq(Hz): ");Serial.println(dF);
  //return; 

  // 3. Afffichage 
  static bool state=false;
  if(state==false)
  {
    for(int i=0; i<Nfft; i++)
    {
      if (i<Nfft-1) 
      {
        float Amp=sqrt((RealFFT[i]*RealFFT[i])+ (ImagFFT[i]*ImagFFT[i]));
        
        Serial.print((float)i*dF);Serial.print(" \t");
        Serial.print((int)floor(1000.0*Amp/((float)N)));Serial.print(" \t"); //||DFT||/N(mV)
        Serial.print(1000.0*Data[2*i+1]);Serial.print(" \t");
        Serial.println(""); 
      }
    }
    state=!state;
  }
  
}


float DFT(float *data, int Ndata, float *real, float *imag, float fs)
{
  int nfft=floor((float)Ndata/2.0); 
  float df=(fs/2.0)/(float)nfft; 
  float somme_r=0.0, somme_i=0.0;  
  float wt=0.0; 
  
  for (int i=0; i<nfft; i++)
  {
      for (int j=0; j<Ndata; j++)
      {
        wt=-2.0*PI*(float)i*(float)j/Ndata;
        somme_r+=data[j]*cos(wt); 
        somme_i+=data[j]*sin(wt); 
      }      
      real[i]= somme_r; somme_r=0.0;
      imag[i]= somme_i; somme_i=0.0;
  }
  return df; 
}