Analyse spectrale sur Arduino
1. Introduction
Ce document montre comment analyser le signal délivré par un capteur au moyen de la transformée de Fourier discrète (TFD). Lorsque le signal est périodique (ou quasi périodique), l'analyse doit fournir la fréquence du signal et relever ses harmoniques.
Une analyse spectrale fine nécessite de pouvoir mémoriser un grand nombre d'échantillons afin d'effectuer le calcul de la TFD. Si N est le nombre d'échantillons et Te la période d'échantillonnage, la résolution fréquentielle est l'inverse de la durée totale, c'est-à-dire 1/(NTe).
Le seul Arduino à microcontrôleur 8 bits ayant assez de mémoire RAM pour calculer une TFD de taille raisonnable est l'Arduino MEGA, qui possède 8 ko de RAM, ce qui permet de traiter Ne=256 échantillons.
On se fixe pour objectif d'analyser un signal périodique dont la fréquence fondamentale est inférieure à 100 Hz. Afin de maximiser la résolution fréquentielle, on choisit la fréquence d'échantillonnage la plus basse possible, soit fe=200 Hz. Le signal analysé peut cependant contenir des composantes de fréquences supérieures à 100 Hz, qui pourraient se replier dans la bande de fréquence analysée. Pour éviter le repliement, nous allons effectuer la numérisation avec un sur-échantillonnage d'un facteur 5, c'est-à-dire à la fréquence d'échantillonnage f'e=1000 Hz. Un filtrage passe-bas anti-repliement, de fréquence de coupure 100 Hz, sera appliqué au signal échantillonné, avant de ramener sa fréquence d'échantillonnage à 200 Hz pour faire l'analyse spectrale par TFD.
Le calcul de la TFD est effectuée par l'algorithme FFT, détaillé dans Transformée de Fourier rapide.
Le code Arduino présenté ici effectue la numérisation avec le convertisseur analogique/numérique (ADC) du microcontrôleur, l'échantillonnage étant effectué par un timer. Cela permettra de tester le fonctionnement avec un générateur de fonctions appliquant un signal sur l'entrée A0 de la carte. Ce code fonctionnera tel quel avec un capteur analogique. L'algorithme de filtrage et d'analyse pourra sans difficulté être transposé à l'analyse du signal délivré par un capteur numérique.
2. Traitement du signal
Le signal est échantillonné à la fréquence f'e=1000 Hz puis subit un filtrage passe-bas RIF de fréquence de coupure 100 Hz (environ). Le filtre à réponse impulsionnelle finie (RIF) est calculé par la méthode de fenêtrage de la réponse impulsionnelle infinie, au moyen de la fonction scipy.signal.fir. Le script python suivant effectue ce calcul pour un nombre de coefficients M donné et trace la réponse fréquentielle du filtre. Les coefficients sont enregistrés dans un fichier qui sera inclus par le code Arduino. Ce fichier contient l'initialisation d'un tableau float rif[M]. Il comporte aussi l'initialisation d'un tableau float gain_filtre[128] contenant le gain du filtre (le spectre comporte 128 points de 0 à 1000 Hz).
La fréquence de coupure est choisie inférieure à 100 Hz, de manière à pratiquement annuler le gain à 100 Hz et ainsi empêcher le repliement des raies de fréquence supérieure à 100 Hz.
import scipy.signal
import numpy
from matplotlib.pyplot import *
a = 0.07 # rapport de la fréquence de coupure sur la fréquence d'échantillonnage
M=51
b = scipy.signal.firwin(numtaps=M,cutoff=[a],window='hamming',nyq=0.5)
Ne=256
Q=int(Ne/2)
w,h=scipy.signal.freqz(b,worN=Q*5)
G=numpy.absolute(h)
figure()
subplot(211)
plot(w/(2*numpy.pi),G)
ylabel("G")
grid()
subplot(212)
plot(w/(2*numpy.pi),numpy.unwrap(numpy.angle(h)))
xlabel("f/fe")
ylabel("déphasage")
grid()
sa = "%0.3f"%a
sa = sa.replace(".",",")
f = open("RIF_%d_%s_%d.h"%(Ne,sa,M),"w")
f.write("float rif[%d] = {"%M)
for i in range(M-1):
f.write("%f,"%b[i])
f.write("%f"%b[M-1])
f.write("};\n")
f.write("float gain_filtre[%d] = {"%Q)
for i in range(Q-1):
f.write("%f,"%G[i])
f.write("%f"%G[Q-1])
f.write("};\n")
f.close()
filtre_rif.pdf
Après le filtrage par convolution, la fréquence d'échantillonnage est réduite à fe=200 Hz. La réduction consiste à retenir un échantillon sur cinq et à les stocker dans un tableau de taille Ne=256. En pratique, Les échantillons obtenus à la fréquence f'e sont stockés dans un tampon circulaire de taille M. Tous les 5 échantillons, on calcule le produit de convolution et on stocke le résultat dans le tableau de taille Ne. Lorsque ce tableau est plein, on calcule la TFD.
On remarque qu'il n'est pas nécessaire de calculer le produit de convolution pour les échantillons de la sortie qui ne sont pas retenus lors de la réduction. On pourrait aussi faire un filtrage récursif, qui permet d'obtenir la même efficacité avec moins de coefficients. Dans le cas d'un filtrage récursif, il faudrait cependant calculer la sortie filtrée pour tous les échantillons (pas seulement 1 sur 5). L'avantage en efficacité serait sans doute largement compensé par ces calculs supplémentaires. Pour cette application, le filtrage par convolution est donc une meilleure solution.
Le calcul de la TFD fournit un tableau de nombres complexes de taille Ne dont on analyse la première moitié, pour des indices allant de 1 à Ne/2-1.
On suppose que le signal analysé est périodique ou quasi périodique, c'est-à-dire qu'il comporte dans son spectre une raie bien discernable et éventuellement des raies pour différents harmoniques. Une raie prend la forme d'un ou plusieurs éléments dont le plus grand se trouve approximativement à la fréquence de la raie, d'autant plus précisément que le nombre d'échantillons (Ne) est grand. Voici un exemple de spectre (signal créneau) provenant du test montré plus loin :
[f,A] = numpy.loadtxt("creneau-13Hz-spectre.txt",unpack=True)
figure()
plot(f,A,".")
xlabel("f (Hz)")
ylabel("A (V)")
grid()
spectre-1.pdf
Ce spectre comporte 4 raies (plus celle de fréquence nulle, que l'on ignore). Pour chaque raie, on doit repérer le maximum et les points adjacents dont la valeur est supérieure à un seuil. Ce seuil doit être inférieur à l'amplitude des raies que l'on souhaite détecter mais supérieur au bruit et aux raies de repliement (qui sont quasi absentes grace au fitrage). Une fois la raie la plus intense trouvée (ici le fondamental), on pourra ajuster le seuil à une certaine fraction de cette raie.
L'extraction des raies du spectre se fait par l'algorithme suivant :
- (1) Recherche de l'élément de la TFD dont le module est maximal et supérieur à un seuil.
- (2) Si cet élément est trouvé, la fréquence correspondante est celle d'une raie. Sinon, la recherche s'arrête.
- (3) On élimine la raie en annulant les éléments ajdacents dont le module est supérieur au seuil.
- (4) Si la raie qui vient d'être déterminée est la première (donc la plus intense), on ajuste le seuil à une fraction donnée de la hauteur de cette raie.
- (5) On réitère à l'étape (1).
Soit imax l'indice de l'élément de la TFD de plus grand module. La fréquence de la raie est estimée par :
où Te est la période d'échantillonnage après réduction. La précision de cette estimation est de l'ordre de 1/T où T=NTe est la durée du signal analysé.
L'estimation de la fréquence d'une raie peut être améliorée par interpolation quadratique. Soit b le module de l'élement de la TFD d'indice imax. Le module de l'élément précédent, d'indice imax-1, est noté a, celui de l'élément suivant, d'indice imax+1, est noté c. L'interpolation quadratique consiste à chercher la parabole qui passe par ces trois points. Le maximum de cette parabole donne une estimation du maximum de la raie. On obtient la position du maximum sous la forme d'un décalage compris entre -1/2 et 1/2 :
qui permet de calculer l'indice interpolé (nombre réel) du maximum i=imax+p puis la fréquence associée. Une estimation de la hauteur de la raie est donnée par :
3. Programme Arduino
Lien vers le fichier de définition du filtre calculé plus haut : RIF_256_0,070_51.h.
On commence par importer la définition du filtre RIF puis on définit les différents paramètres constants :
#include "RIF_256_0,070_51.h"
//#define DEBUG // à décommenter pour afficher les échantillons et la TFD
#define FILTRAGE // à commenter pour enlever le filtrage
#define M 51 // nombre de coefficients du RIF
#define INTERPOL 1 // analyse des raies par interpolation quadratique
#define SEUIL_REL 0.08 // seuil relatif à l'harmonique la plus haute
#define P 8
#define Ne 256 // 2^P nombre d'échantillons
#define PERIOD 1000 // période d'échantillonnage en microsecondes
#define REDUC 5 // réduction de la fréquence d'échantillonnage après filtrage
#define VREF 4.96 // tension mesurée sur la sortie AREF
const float duree = PERIOD*1e-6*Ne*REDUC; // durée du signal analysé en secondes
Le nombre d'échantillons est de la forme Ne=2P, une condition imposée par l'algorithme FFT. Cet algorithme nécessite deux tableaux de taille Ne contenant des nombres complexes, soit quatre tableaux contenant des nombres réels (type float). Ces tableaux sont déclarés en variables globales :
volatile float A_re[Ne];
float A_im[Ne];
float B_re[Ne];
float B_im[Ne];
Le mot clé volatile est nécessaire pour les variables utilisées dans les interruptions. Voici la déclaration du tampon circulaire, de l'indice indiquant le dernier échantillon dans ce tampon et d'un compteur pour effectuer la réduction de la fréquence d'échantillonnage :
volatile float circulaire[M]; // liste circulaire pour le filtrage
uint16_t dernier; // indice du dernier échantillon dans la liste circulaire
uint8_t cr; // compteur pour la réduction
La variable suivante contient l'indice du tableau en cours de remplissage :
volatile uint16_t indice;
Enfin quelques variables globales dont l'utilité apparaîtra plus loin :
volatile char acquisition=0;
volatile char convol;
char start = 0;
uint8_t flag;
On définit une structure Raie pour transmettre des informations sur une raie du spectre :
typedef struct raie {
int indice; // indice dans le tableau TFD
float pos; // indice réel extrapolé
float amp; // amplitude
} Raie;
La fonction suivante permet de configurer le timer pour l'échantillonnage. Voir la page Conversion analogique-numérique rapide avec échantillonnage pour l'explication du fonctionnement du timer. La période fournie en argument est en microsecondes.
/* configuration de l'horloge d'échantillonnage */
void timer1_init(uint32_t period) {
uint16_t diviseur[6] = {0,1,8,64,256,1024};
TCCR1A = 0;
TCCR1B = 0;
TCCR1B |= (1 << WGM12);
uint32_t top = (F_CPU/1000000*period);
int clock = 1;
while ((top>0xFFFF)&&(clock<5)) {
clock++;
top = (F_CPU/1000000*period/diviseur[clock]);
}
OCR1A = top;
OCR1B = top >> 1;
TIMSK1 = (1 << OCIE1B); // interruption
TCCR1B |= clock;
}
La fonction suivante permet de configurer le multiplexeur et le convertisseur analogique/numérique du microcontrôleur :
/* Initialisation du convertisseur ADC */
void adc_init(uint8_t multiplex, uint8_t prescaler) {
ADMUX = (1 << REFS0);
ADMUX |= multiplex & 0b00011111;
//ADMUX |= (1 << ADLAR); //left adjust
ADCSRA = 0;
ADCSRA |= (1 << ADEN); // enable ADC
ADCSRA |= (1 << ADATE); // Auto Trigger Enable
ADCSRA |= (1 << ADIE); // interrupt enable
ADCSRA |= prescaler ;
ADCSRB = 0;
ADCSRB |= (multiplex & 0b00100000) >> 2;
ADCSRB |= (1 << ADTS2)|(1 << ADTS0);// trigger source : timer 1 comp B
}
La conversion est déclenchée par l'interruption générée par le timer. Lorsque la conversion est terminée, une interruption est générée et la fonction suivante est exécutée :
ISR(ADC_vect) {
circulaire[dernier] = ADCL | (ADCH << 8);
cr += 1;
if (convol&&(cr==REDUC)) {
cr = 0;
#ifdef FILTRAGE
int16_t k = dernier;
int16_t m = 0;
float accum = 0;
while (m<M) {
accum += circulaire[k]*rif[m];
k -= 1;
if (k==-1) k=M-1;
m += 1;
}
A_re[indice] = accum;
#else
A_re[indice] = circulaire[dernier];
#endif
indice += 1;
if (indice==Ne) {
TCCR1B &= ~0x7;
acquisition = 0;
cli();
}
}
if (cr==M) {convol = 1; cr=0;}
dernier += 1;
if (dernier==M) {
dernier=0;
}
}
Cette fonction stocke dans le tampon circulaire l'échantillon qui vient d'être numérisé. Tous les REDUC échantillons, le produit de convolution du tampon circulaire par les coefficients du filtre est calculé puis stocké dans le tableau A_re. Lorsque ce tableau est plein, les interruptions sont annulées afin d'interrompre l'échantillonnage. Si la macro FILTRAGE n'est pas définie, la réduction se fait sans filtrage. La variable convol, initialement nulle, passe à 1 lorsque le tampon circulaire est rempli : cela permet de démarrer la convolution seulement lorsque M échantillons sont disponibles.
La fonction suivante est appelée lors de l'interruption générée par le timer (en plus de la conversion A/N). La présence de cette fonction est indispensable bien qu'elle ne joue aucun rôle dans le traitement du signal. On s'en sert pour faire alterner le niveau de la sortie 18, ce qui permettra de contrôler le bon déroulement des interruptions avec un oscilloscope.
ISR(TIMER1_COMPB_vect) {
if (flag==0) {
flag = 1;
PORTD &= ~(1<<PORTD3); // sortie 18 sur Arduino MEGA, 3 sur UNO
}
else {
flag = 0;
PORTD |= (1<<PORTD3);
}
}
Les deux fonctions suivantes implémentent l'algorithme FFT, tel qu'il est expliqué dans Transformée de Fourier rapide. La première effectue l'inversion des bits d'un nombre entier codé sur P bits :
int inversion(uint16_t i) {
uint8_t bi;
int8_t b;
uint8_t bits[P];
for (b=0; b<P; b++) {
bi = i & 0b1;
i = i >> 1;
bits[b] = bi;
}
uint16_t facteur = 1;
uint16_t j = 0;
for (b=P-1; b>=0; b--) {
j += bits[b]*facteur;
facteur *= 2;
}
return j;
}
La fonction suivante effectue le calcul de la TFD. L'algorithme FFT nécessite deux tableaux (notés A et B). À chaque niveau de TFD, l'un des deux tableaux contient les TFD de niveau inférieur (TFD à N termes) et l'autre tableau reçoit les nouvelles TFD (TFD à 2N termes). Il faut donc alterner le rôle de A et B (source ou destination) à chaque niveau de TFD. En principe, cela se fait en langage C par un échange des pointeurs sur A et B. Cependant, le compilateur de l'IDE Arduino considère comme incompatibles les types float * (pointeur sur un flottant) et float[256] (pointeur sur un tableau de flottants), ce qui empêche de procéder à un échange de pointeurs sur deux tableaux, à moins de créer un troisième tableau de taille égale aux deux premier, ce qui est ici exclu pour des raisons d'économie de la mémoire. La fonction suivante attribue le rôle de source ou de destination aux tableaux A ou B en fonction de la valeur de parite (0 ou 1). Il est convenu que les échantillons du signal sont au départ dans le tableau A (variables globales A_re et A_im) et que la TFD finale se trouve dans le tableau B (variables globales B_re et B_im). Les valeurs initiales du tableau A sont bien sûr perdues.
void fft() {
uint16_t j;
for (int k=0; k<Ne; k++) {
j = inversion(k);
B_re[j] = A_re[k];
B_im[k] = A_im[k];
}
uint16_t taille = 1;
uint16_t taille_precedente;
uint16_t nombre_tfd = Ne;
uint16_t position;
float phi, W_re, W_im;;
float x,y;
uint8_t parite = 0;
for (uint16_t q=0; q<P; q++) {
taille_precedente = taille;
taille *= 2;
nombre_tfd /= 2;
for (uint16_t m=0; m<nombre_tfd; m++) {
phi = -2*PI/taille;
position = m*taille;
uint16_t i,j,k;
for (i=0; i<taille_precedente; i++) {
W_re = cos(phi*i);
W_im = sin(phi*i);
j = position+taille_precedente+i;
k = position +i;
if (parite==0) {
x = B_re[j];
y = B_im[j];
A_re[k] = B_re[k] + x*W_re-y*W_im;
A_im[k] = B_im[k] + x*W_im+y*W_re;
}
else {
x = A_re[j];
y = A_im[j];
B_re[k] = A_re[k] + x*W_re-y*W_im;
B_im[k] = A_im[k] + x*W_im+y*W_re;
}
}
for (i=taille_precedente; i<taille; i++) {
W_re = cos(phi*i);
W_im = sin(phi*i);
j = position+i;
k = position+i-taille_precedente;
if (parite==0) {
x = B_re[j];
y = B_im[j];
A_re[j] = B_re[k] + x*W_re-y*W_im;
A_im[j] = B_im[k] + x*W_im+y*W_re;
}
else {
x = A_re[j];
y = A_im[j];
B_re[j] = A_re[k] + x*W_re-y*W_im;
B_im[j] = A_im[k] + x*W_im+y*W_re;
}
}
}
parite = ~parite;
}
if (parite) {
for (int k=0; k<Ne; k++) {B_re[k] = A_re[k]; B_im[k] = A_im[k];}
}
}
La fonction suivante recherche la raie d'amplitude maximale dans le spectre (le module de la TFD se trouve dans le tableau B_re). Elle recherche l'élément de la TFD dont le module est le plus grand. Si le module maximal est inférieur au seuil, l'indice renvoyé est -1 (raie.indice). Si le module maximal est supérieur au seuil, la fonction renvoie l'indice du maximum et le module (raie.pos et raie.amp). Si la valeur de interpol est 1, une interpolation quadratique est effectuée. Dans ce cas, raie.indice contient l'indice (entier) du maximum, raie.pos contient l'indice extrapolé (nombre réel) et raie.amp contient l'amplitude extrapolée.
void recherche_max(Raie *raie, float seuil, bool interpol) {
uint16_t imax = 1;
float module;
float module_max = 0;
for (uint16_t k=1; k<Ne/2; k++) {
module = B_re[k];
if (module > module_max) {
imax = k;
module_max = module;
}
}
if (module_max < seuil) {raie->indice = -1; return;}
if (interpol) {
float a,b,c;
a = B_re[imax-1];
b = module_max;
c = B_re[imax+1];
if ((a>seuil)&&(c>seuil)) {
raie->indice = imax;
float p = 0.5*(a-c)/(a-2*b+c);
raie->amp = b-0.25*(a-c)*p;
raie->pos = imax+p;
return;
}
}
raie->indice = imax;
raie->pos = imax;
raie->amp = module_max;
}
La fonction suivante efface une raie du spectre. Il s'agit d'annuler les éléments adjacents à l'élément de module maximal (indice imax) et dont les modules sont supérieurs à un seuil.
void effacer_raie(int imax, float seuil) {
int i = imax;
while ((i>0)&&(B_re[i]>seuil)) {
B_re[i] = 0;
i -= 1;
}
i = imax + 1;
while ((i<Ne/2)&&(B_re[i]>seuil)) {
B_re[i] = 0;
i += 1;
}
}
La fonction setup initialise la variable start à 0, ce qui signifie que la numérisation n'a pas démarré, et configure la sortie 18 en sortie.
void setup() {
Serial.begin(115200);
start = 0;
DDRD |= 1 << PORTD3; // PD3 en sortie
}
Voici enfin la fonction loop. Si la valeur start est 0 (c'est le cas initialement), elle démarre une acquisition en déclenchant le timer et en activant les interruptions. Lorsque l'acquisition est en cours, la variable acquisition a la valeur 1 et la fonction loop ne fait rien. Lorsque la valeur de acquisition devient 1, on procède au calcul de la TFD, suivi d'une normalisation (multiplication par 2/Ne). La boucle do {} while (raie.indice>0) effectue la recherche itérative des raies du spectre dont l'amplitude est supérieure à celle de la raie la plus intense multipliée par SEUIL_REL. On affiche sur la console série l'estimation de la fréquence, de l'amplitude et de l'amplitude en décibel (relative à la plus intense). Lorsque l'analyse est terminée, on initialise la variable start à 0 afin de relancer une acquisition à la prochaine exécution de loop. Si la macro DEBUG est définie, il y a aussi un affichage des échantillons du signal et du spectre (en volts).
void loop() {
delay(100);
if (start==0) {
for (uint32_t i=0; i<Ne; i++) {A_re[i] = 0; A_im[i] = 0;}
start = 1;
flag = 0;
acquisition=1;
indice = 0;
dernier = 0;
cr = 0;
convol = 0;
cli();
timer1_init(PERIOD);
adc_init(0,4);
sei();
}
if (acquisition==0) {
#ifdef DEBUG
Serial.println("----------------------");
for (uint16_t k=0; k<Ne; k++) {
Serial.print((float)k*PERIOD*REDUC*1e-6,4);Serial.print(" "); Serial.println(A_re[k]*VREF/1024);
}
#endif
fft();
for (uint16_t k=0; k<Ne/2; k++) {
#ifdef FILTRAGE
B_re[k] = sqrt(B_re[k]*B_re[k]+B_im[k]*B_im[k])*2.0/Ne/gain_filtre[k];
#else
B_re[k] = sqrt(B_re[k]*B_re[k]+B_im[k]*B_im[k])*2.0/Ne;
#endif
}
#ifdef DEBUG
Serial.println("----------------------");
for (uint16_t k=0; k<Ne/2; k++) {
Serial.print(k/duree,4);Serial.print(" "); Serial.println(B_re[k]*5.0/1024,4);
}
#endif
float seuil = 20; // seuil initial
char premier = 1;
Raie raie;
float f;
float amp_max;
Serial.println("-------------------");
do {
recherche_max(&raie,seuil,INTERPOL);
if (premier) {
premier = 0;
seuil = raie.amp * SEUIL_REL; // réglage du seuil par rapport à l'harmonique de plus grande amplitude
amp_max = raie.amp;
}
if (raie.indice>0) {
effacer_raie(raie.indice,seuil);
f = raie.pos*1/duree;
if (raie.amp>0)
{
Serial.print("f = "); Serial.print(f); Serial.print(" +/-"); Serial.print(0.5/duree);
Serial.print(", A = "); Serial.print(raie.amp*VREF/1024);
Serial.print(" ("); Serial.print(20*log10(raie.amp/amp_max)); Serial.println(" dB)");
}
}
}
while (raie.indice>0);
start = 0;
}
}
4. Test de l'analyseur
Le signal est fourni par un générateur de fonctions. La tension appliquée sur l'entrée analogique A0 ne doit pas descendre en dessous de 0 ni monter au dessus de 5 V (un dépassement trop grand détruirait l'entrée du microcontrôleur). Il faut donc appliquer un décalage (offset) de 2,5 V et ne pas dépasser une amplitude de 2,5 V. Il est préférable d'utiliser le circuit suivant, qui permet d'appliquer le décalage à un signal alternatif et protège l'entrée de l'Arduino contre les dépassements de tension. L'amplificateur linéaire intégré est un TLV 2372 alimenté sous une tension de 5 V par l'arduino. Il s'agit d'un ampli rail to rail, dont la tension de sortie peut varier entre 0 et 5 V. Le décalage de 2,5 V est obtenu par réglage du potentiomètre.
Figure pleine pageVoici un exemple avec un signal créneau d'amplitude de crête à crête 4,5 V et de fréquence 13,0 Hz.
[t,u] = numpy.loadtxt("creneau-13Hz-signal.txt",unpack=True)
[f,A] = numpy.loadtxt("creneau-13Hz-spectre.txt",unpack=True)
figure(figsize=(8,8))
subplot(211)
plot(t,u)
xlabel("t (s)")
ylabel("u (V)")
grid()
subplot(212)
plot(f,A)
xlabel("f (Hz)")
ylabel("A (V)")
grid()
signal-1.pdf
f = 13.25 +/-0.39, A = 2.63 (0.00 dB)
f = 39.06 +/-0.39, A = 0.82 (-10.11 dB)
f = 65.62 +/-0.39, A = 0.56 (-13.39 dB)
f = 91.41 +/-0.39, A = 0.31 (-18.47 dB)
Les harmoniques de fréquences inférieures à 100 Hz sont bien détectés (rang 1, 3, 5 et 7). Compte tenu de la résolution fréquentielle, les fréquences sont correctes. Les amplitudes des harmoniques sont théoriquement -9,5 dB, -14,0 dB et -16,9 dB. L'évaluation des amplitudes est donc plutôt mauvaise, ce qui n'est pas étonnant vu le faible nombre d'échantillons (256). Pour obtenir une bonne détermination des hauteurs de raie, il faut calculer la TFD avec un nombre d'échantillons beaucoup plus grand (plusieurs dizaines de millier) ou compléter le signal par des zéros si la durée de l'analyse est trop courte.
Pour comparaison, voici l'analyse du même signal mais sans filtrage :
[t,u] = numpy.loadtxt("creneau-13Hz-signal-brut.txt",unpack=True)
[f,A] = numpy.loadtxt("creneau-13Hz-spectre-brut.txt",unpack=True)
figure(figsize=(8,8))
subplot(211)
plot(t,u)
xlabel("t (s)")
ylabel("u (V)")
grid()
subplot(212)
plot(f,A)
xlabel("f (Hz)")
ylabel("A (V)")
grid()
signal-2.pdf
f = 13.26 +/-0.39, A = 2.68 (0.00 dB)
f = 39.06 +/-0.39, A = 0.75 (-11.10 dB)
f = 65.62 +/-0.39, A = 0.55 (-13.75 dB)
f = 82.03 +/-0.39, A = 0.30 (-19.13 dB)
f = 91.41 +/-0.39, A = 0.27 (-20.07 dB)
f = 29.69 +/-0.39, A = 0.23 (-21.41 dB)
En plus des raies du signal créneau de fréquences inférieures à 100 Hz, l'analyse détecte deux raies qui proviennent du repliement de raies de fréquences supérieures à 100 Hz. Ce résultat est obtenu avec un seuil relatif de 0,08.
Pour analyser des signaux de fréquence supérieure à 100 Hz, il faut augmenter la fréquence d'échantillonnage. Étant donnée l'impossibilité d'augmenter le nombre d'échantillons (à cause de la taille de la RAM), cela se fait au détriment de la précision absolue de la fréquence. Avec le filtre à 51 coefficients, nous avons pu augmenter la fréquence d'échantillonnage jusqu'à f'e=2000 Hz. Au delà, le calcul du produit de convolution n'a plus le temps de se faire pendant la période d'échantillonnage. Avec un filtre à 21 coefficients et a=0,06, nous avons pu augmenter la fréquence d'échantillonnage jusqu'à f'e=4 kHz, ce qui permet en principe d'analyser des signaux jusqu'à 500 Hz. Voici par exemple deux analyses consécutives d'un signal créneau de fréquence 148,3 Hz :
f = 148.12 +/-1.56, A = 2.02 (0.00 dB)
f = 356.25 +/-1.56, A = 0.24 (-18.47 dB)
f = 171.87 +/-1.56, A = 0.16 (-21.76 dB)
-------------------
f = 148.10 +/-1.56, A = 2.02 (0.00 dB)
f = 353.12 +/-1.56, A = 0.22 (-19.30 dB)
La seconde est correcte mais la première fait apparaître une raie de repliement. Le filtre n'est plus assez sélectif pour éliminer complètement le repliement et il faudrait alors augmenter le seuil de détection des raies.
En conclusion, l'analyse spectrale fonctionne bien sur l'Arduino MEGA, à condition de se limiter à 256 échantillons et d'analyser des signaux de fréquences de signaux de l'ordre de quelques dizaines de Hertz. Les fréquences du fondamental et des harmoniques intenses sont bien déterminées mais la détermination de la hauteur des raies est médiocre en raison du faible nombre d'échantillons. Pour l'analyse de signaux dans le domaine audio, l'arduino MEGA n'est pas assez rapide. Il faudra pour cela utiliser un microcontrôleur 32 bits (par exemple l'Arduino Due) qui possèdera en plus une plus grande mémoire RAM et dont le microprocesseur comportera éventuellement une unité de calcul en virgule flottante.
