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 à 1000 Hz. Afin de maximiser la résolution fréquentielle, on choisit la fréquence d'échantillonnage la plus basse possible, soit fe=2000 Hz. Le signal analysé peut cependant contenir des composantes de fréquences supérieures à 1000 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=10000 Hz. Un filtrage passe-bas anti-repliement, de fréquence de coupure 1000 Hz, sera appliqué au signal échantillonné, avant de ramener sa fréquence d'échantillonnage à 2000 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=10 kHz puis subit un filtrage passe-bas RIF de fréquence de coupure 1000 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[1024] contenant le gain du filtre (le spectre comporte 1024 points de 0 à 1000 Hz).
La fréquence de coupure est choisie inférieure à 1000 Hz, de manière à pratiquement annuler le gain à 1000 Hz et ainsi empêcher le repliement des raies de fréquence supérieure à 1000 Hz.
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_rif.pdf
Après le filtrage par convolution, la fréquence d'échantillonnage est réduite à fe=2000 Hz. 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-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 :
- (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_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 :
#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,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 2,0 V et de fréquence 139,1 Hz.
[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-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 à 1000 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 (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-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 à 1000 Hz, l'analyse détecte 3 raies provenant du repliement de raies de fréquences supérieures à 1000 Hz. 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 à 1000 Hz). Pour cette fréquence d'échantillonnage de 10 kHz, 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'à 2000 Hz. La fréquence d'échantillonnage finale doit être fe=4000 Hz et la fréquence d'échantillonnage avant réduction doit être f'e=20 kHz. 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 2000 Hz. 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,1 Hz :
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 à 2000 Hz. 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'à 2000 Hz) 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 :
- La fréquence d'échantillonnage.
- Le nombre de coefficients du filtre de convolution et sa fréquence de coupure.
- Le seuil de détection des raies.
