Table des matières

Analyse spectrale sur Arduino Due

1. Introduction

Ce document montre comment analyser le signal délivré par un capteur (par exemple un microphone) 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).

L'Arduino Due possède un microcontrôleur 32 bits (SAM3X8E) et 96 ko de RAM, ce qui permet de traiter Ne = 2048 échantillons et permettra donc d'analyser un signal audio avec une bonne précision relative de fréquence.

On se fixe pour objectif d'analyser un signal périodique dont la fréquence fondamentale est inférieure à 1000Hz. Afin de maximiser la résolution fréquentielle, on choisit la fréquence d'échantillonnage la plus basse possible, soit fe = 2000Hz. Le signal analysé peut cependant contenir des composantes de fréquences supérieures à 1000Hz, 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 = 10000Hz. Un filtrage passe-bas anti-repliement, de fréquence de coupure 1000Hz, sera appliqué au signal échantillonné, avant de ramener sa fréquence d'échantillonnage à 2000Hz 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 = 10kHz puis subit un filtrage passe-bas RIF de fréquence de coupure 1000Hz (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[1024] contenant le gain du filtre (le spectre comporte 1024 points de 0 à 1000 Hz).

La fréquence de coupure est choisie inférieure à 1000Hz, de manière à pratiquement annuler le gain à 1000Hz et ainsi empêcher le repliement des raies de fréquence supérieure à 1000Hz.

import scipy.signal
import numpy
from matplotlib.pyplot import *

a = 0.08 # rapport de la fréquence de coupure sur la fréquence d'échantillonnage
M=61
b = scipy.signal.firwin(numtaps=M,cutoff=[a],window='hamming',nyq=0.5)
Ne = 2048
Q=int(Ne/2)
w,h=scipy.signal.freqz(b,worN=Q*5)
figure()
subplot(211)
plot(w/(2*numpy.pi),numpy.absolute(h))
ylabel("G")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.unwrap(numpy.angle(h)))
xlabel("f/fe")
ylabel("déphasage")
grid()

G=numpy.absolute(h)

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[2048] = {")
for i in range(Q-1):
    f.write("%f,"%G[i])
f.write("%f"%G[Q-1])
f.write("};\n")
f.close()


               
filtre_riffiltre_rif.pdf

Après le filtrage par convolution, la fréquence d'échantillonnage est réduite à fe = 2000Hz. La réduction consiste à retenir un échantillon sur cinq et à les stocker dans un tableau de taille Ne = 2048. 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-139Hz-spectre.txt",unpack=True)
figure(figsize=(12,6))
plot(f,A,".")
xlabel("f (Hz)")
ylabel("A (V)")
grid()

              
spectre-1spectre-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.

Le module de la TFD est normalisé avec le gain du filtre, ce qui permet de corriger l'abaissement du gain à l'approche de la fréquence de coupure.

L'extraction des raies du spectre se fait par l'algorithme suivant :

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 :

Te est la période d'échantillonnage après réduction. La précision de cette estimation est de l'ordre de 1/TT=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_2048_0,080_61.h.

On commence par importer la définition du filtre RIF puis on définit les différents paramètres constants :

ArduinoDueAnalyseSpectrale.ino
#include "RIF_2048_0,080_61.h"
//#define DEBUG // à décommenter pour afficher les échantillons et la TFD
#define FILTRAGE // à commenter pour enlever le filtrage
#define M 61 // 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 11
#define Ne 2048 // 2^P nombre d'échantillons
#define PERIOD 100 // période d'échantillonnage en microsecondes
#define REDUC 5 // réduction de la fréquence d'échantillonnage après filtrage
#define VREF 3.3 // tension mesurée sur la sortie AREF

#define ADC_MR_TRIG1 (1<<1)
uint8_t nchannels;
#define MAX_NCHAN 1
uint8_t channels[MAX_NCHAN];

uint8_t gain[MAX_NCHAN];
uint8_t channels_pins[8] ={7,6,5,4,3,2,1,0};

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. La période fournie en argument est en microsecondes.

void config_timer(uint32_t period) {
  uint32_t ticks = 42*period;
  uint32_t channel = 0;
  uint8_t clock = TC_CMR_TCCLKS_TIMER_CLOCK1; // horloge 84MHz/2=42 MHz
  pmc_enable_periph_clk (TC_INTERFACE_ID + 0*3+channel) ;
  TC_Configure(TC0,channel, TC_CMR_WAVE | TC_CMR_WAVSEL_UP_RC | TC_CMR_ACPA_CLEAR | TC_CMR_ACPC_SET |clock); 
  TC_SetRC(TC0,channel,ticks); 
  TC_SetRA(TC0,channel,ticks >> 1);
  TC_Start(TC0, channel);
}
              

La fonction suivante permet de configurer le multiplexeur et le convertisseur analogique/numérique du microcontrôleur :

void config_adc() {
    int k;
    pmc_enable_periph_clk(ID_ADC);
    adc_init(ADC, VARIANT_MCK, ADC_FREQ_MAX, 0);
    adc_configure_timing(ADC, 0, ADC_SETTLING_TIME_3, 1);
    adc_set_resolution(ADC,ADC_12_BITS);
    NVIC_EnableIRQ (ADC_IRQn) ;   
    ADC->ADC_IDR=0xFFFFFFFF ;  
    ADC->ADC_CHDR = 0xFFFF; 
    ADC->ADC_CHSR = 0x0000; 
    ADC->ADC_WPMR &= ~0x1;
    ADC->ADC_CHER = 0x0000;
    ADC->ADC_SEQR1 = 0x00000000;
    ADC->ADC_CGR = 0x00000000;
    ADC->ADC_COR = 0x00000000;
    for (k=0; k<nchannels; k++) {
      ADC->ADC_SEQR1 |= ((channels[k]&0xF)<<(4*k));
      ADC->ADC_CHER |= (1<<k);
      if (gain[k]==2) ADC->ADC_CGR |= 2 << (2*k);
      if (gain[k]==4) ADC->ADC_CGR |= 3 << (2*k);
      ADC->ADC_COR |= 1 << k; 
    }
    ADC->ADC_MR = (ADC->ADC_MR & 0xFFFFFFF0) | ADC_MR_TRIG1 | ADC_MR_TRGEN | ADC_MR_USEQ;
    ADC->ADC_IDR=~(1<<channels[0]); 
    ADC->ADC_IER=1<<channels[0];
}
              

La conversion est déclenchée par l'interruption générée par le convertisseur ADC lorsqu'une converion est terminée :

void ADC_Handler() {
  circulaire[dernier] = *(ADC->ADC_CDR+channels[0]);
    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) {
          TC_Stop(TC0,0);
          pmc_disable_periph_clk (TC_INTERFACE_ID);
          acquisition = 0;
          
      }
  }
  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.

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é.

void setup() {
  Serial.begin(115200);
  start = 0;
  
}                   
                   

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 et division par le gain du filtre. 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;
        nchannels = 1;
        channels[0] = channels_pins[0];
        gain[0] = 0;
        config_adc();
        config_timer(PERIOD);
   }
  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/4096);
     }
     #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]*VREF/4096,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/4096);
                Serial.print("  ("); Serial.print(20*log10(raie.amp/amp_max)); Serial.println(" dB)");
              }
        }
     } 
     while (raie.indice>0);
     start = 0;
  }
}
                   

4. Test de l'analyseur

L'arduino est relié au PC par le port USB Programming Port.

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,5V et ne pas dépasser une amplitude de 2,5V. 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,5V est obtenu par réglage du potentiomètre.

ampliEntree.svgFigure pleine page

Voici un exemple avec un signal créneau d'amplitude de crête à crête 2,0V et de fréquence 139,1Hz.

[t,u] = numpy.loadtxt("creneau-139Hz-signal.txt",unpack=True)
[f,A] = numpy.loadtxt("creneau-139Hz-spectre.txt",unpack=True)
figure(figsize=(8,8))
subplot(211)
plot(t,u)
xlim(0,0.1)
xlabel("t (s)")
ylabel("u (V)")
grid()
subplot(212)
plot(f,A)
xlabel("f (Hz)")
ylabel("A (V)")
grid()
              
signal-1signal-1.pdf
f = 139.07 +/-0.49,  A = 0.90  (0.00 dB)
f = 417.27 +/-0.49,  A = 0.31  (-9.22 dB)
f = 695.31 +/-0.49,  A = 0.19  (-13.39 dB)
f = 973.63 +/-0.49,  A = 0.16  (-15.07 dB)   
              

Les harmoniques de fréquences inférieures à 1000Hz 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,5dB, -14,0dB et -16,9dB. L'évaluation des amplitudes est donc plutôt mauvaise, ce qui n'est pas étonnant vu le faible nombre d'échantillons (2048).

Pour comparaison, voici l'analyse du même signal mais sans filtrage :

[t,u] = numpy.loadtxt("creneau-139Hz-signal-brut.txt",unpack=True)
[f,A] = numpy.loadtxt("creneau-139Hz-spectre-brut.txt",unpack=True)
figure(figsize=(8,8))
subplot(211)
plot(t,u)
xlim(0,0.1)
xlabel("t (s)")
ylabel("u (V)")
grid()
subplot(212) 
plot(f,A)
xlabel("f (Hz)")
ylabel("A (V)")
grid()
              
signal-2signal-2.pdf
f = 139.50 +/-0.49,  A = 1.01  (0.00 dB)
f = 417.97 +/-0.49,  A = 0.41  (-7.78 dB)
f = 696.29 +/-0.49,  A = 0.25  (-12.11 dB)
f = 974.61 +/-0.49,  A = 0.15  (-16.53 dB)
f = 467.77 +/-0.49,  A = 0.11  (-19.12 dB)
f = 746.09 +/-0.49,  A = 0.10  (-19.70 dB)
f = 189.45 +/-0.49,  A = 0.09  (-20.53 dB)
              

En plus des raies du signal créneau de fréquences inférieures à 1000Hz, l'analyse détecte 3 raies provenant du repliement de raies de fréquences supérieures à 1000Hz. Ce résultat est obtenu avec un seuil relatif de 0,08.

En conclusion, nous obtenons une analyse spectrale satisfaisante pour des signaux dans les basses fréquences du domaine audio (inférieures à 1000Hz). Pour cette fréquence d'échantillonnage de 10kHz, nous avons constaté que le programme cesse de fonctionner correctement lorsque M=80, en raison de l'impossibilité de calculer le produit de convolution pendant la période d'échantillonnage. Il est possible d'accélérer le calcul de la convolution en l'effectuant avec des nombres entiers. L'augmentation de l'efficacité du filtrage (par augmentation de M) n'est toutefois pas nécessaire car le spectre est normalisé avec le gain du filtre.

Voyons comment augmenter la plage de fréquences analysées jusqu'à 2000Hz. La fréquence d'échantillonnage finale doit être fe = 4000Hz et la fréquence d'échantillonnage avant réduction doit être f'e = 20kHz. Si l'on souhaite garder la même précision (absolue) de fréquence, il faut doubler le nombre d'échantillons, soit Ne = 4096 (c'est la plus grande taille possible avec la RAM disponible). Le temps disponible pour calculer le produit de convolution étant deux fois moindre, nous réduisons la taille du filtre à M = 31. Le filtre obtenu est beaucoup moins sélectif et présente donc une diminution du gain importante en dessous de 2000Hz. Dans la mesure ou nous effectuons une normalisation de la TFD par ce gain, la réduction de M devrait être acceptable à condition d'abaisser aussi la fréquence de coupure. Pour faire ces changements, on doit poser dans le script python :

a=0.06
Ne=4096
M=31
                  

Les changements à faire dans le programme arduino sont :

#include "RIF_4096_0,06_31.h
#define M 31
#define P 12
#define Ne 4096
#define PERIOD 50
                  

Voici le résultat de l'analyse du même signal créneau que précédemment, qui permet de détecter les harmoniques jusqu'au rang 13 :

f = 138.75 +/-0.49,  A = 1.09  (0.00 dB)
f = 416.99 +/-0.49,  A = 0.42  (-8.33 dB)
f = 972.66 +/-0.49,  A = 0.18  (-15.85 dB)
f = 695.31 +/-0.49,  A = 0.17  (-16.16 dB)
f = 1250.98 +/-0.49,  A = 0.13  (-18.70 dB)
f = 1806.64 +/-0.49,  A = 0.10  (-21.00 dB)
f = 1528.32 +/-0.49,  A = 0.09  (-21.49 dB)    
                  

Voici le résultat de l'analyse d'un signal créneau de fréquence 435,1Hz :

f = 434.94 +/-0.49,  A = 0.90  (0.00 dB)
f = 1304.93 +/-0.49,  A = 0.32  (-9.10 dB)
                  

Dans ce cas, seul l'harmonique de rang 3 a une fréquence inférieure à 2000Hz. Voici ce que donne l'analyse de ce signal sans filtrage :

f = 434.90 +/-0.49,  A = 0.91  (0.00 dB)
f = 1304.86 +/-0.49,  A = 0.33  (-8.81 dB)
f = 1825.20 +/-0.49,  A = 0.21  (-12.59 dB)
f = 955.08 +/-0.49,  A = 0.16  (-14.89 dB)
f = 84.96 +/-0.49,  A = 0.14  (-16.50 dB)
f = 785.16 +/-0.49,  A = 0.11  (-18.08 dB)
f = 1655.27 +/-0.49,  A = 0.10  (-19.36 dB)
f = 1474.61 +/-0.49,  A = 0.09  (-20.59 dB)
f = 604.49 +/-0.49,  A = 0.07  (-21.82 dB)

                  

Il y a dans ce spectre 7 raies de repliement.

En conclusion, l'Arduino Due permet de faire l'analyse spectrale d'un signal du domaine audio (jusqu'à 2000Hz) avec une précision de l'ordre du hertz. Pour les signaux de très basse fréquence (quelques Hz), il faudra bien sûr réduire la fréquence d'échantillonnage afin de garder une bonne précision relative de la fréquence. Les différents paramètres à ajuster en fonction de l'application visée sont :

Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.