Catégories
Algorithme Arduino Matlab Projets traitement du signal Traitement du signal sur matlab

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 :

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

Elle retourne trois sorties :

  1. *real: tableau de la partie réelle du spectre de taille N/2
  2. *imag: tableau de la partie réelle du spectre de taille N/2
  3. 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;
}

On verra dans le prochain tuto comment suivre un ou plusieurs harmoniques du spectre.

Laisser un commentaire