Mise en oeuvre d'un filtre FIR

Théorie des Filtres FIR : Filtrage Numérique FIR

ALGORITHME

Branchements

Pour caractériser un filtre, il sera nécessaire d’envoyer en entrée un signal comportant toutes les fréquences ( bruit blanc ), et de tracer le spectre du signal en sortie.

Réalisation d’un Filtre FIR

Pour commencer, nous allons mettre en oeuvre un filtre FIR avec des coefficients fixes, déterminés par un script python.

Simulation ( Python )

Q1. A l’aide du script ci dessous, générer les coefficients d’un filtre FIR de taille 64, de fréquence de coupure 4500Hz, et de fréquence d’échantillonnage 44100 HZ.

REMARQUE : les coefficients sont mis en forme dans le fichier FIR.txt.

gene_coeffs_FIR.py

gene_coeffs_FIR.py

#!/usr/bin/python
#jupyter-qtconsole

import numpy as np
from scipy.signal import dlti
from control import *
from scipy import signal
from scipy.signal import freqz, get_window
import matplotlib.pyplot as plt
import sys

#=================================================
#  		CAHIER DES CHARGES
#=================================================
sample_rate=44100.0
Filter_Order = 63 # low pass filter -> ODD | high pass filter -> EVEN
Filter_Type='low' # low / high / bandpass / bandstop
cutoff_1_hz = 0   # 0.0 if low or high pass filter
cutoff_2_hz = 4500.0
numtaps = Filter_Order + 1  # Length of the filter (number of coefficients, i.e. the filter order + 1)
#=================================================
#       CALCUL DES COEFFICIENTS
#=================================================
nyq_rate = sample_rate / 2.
norm_freq1=cutoff_1_hz/nyq_rate
norm_freq2=cutoff_2_hz/nyq_rate

if Filter_Type=='low':
   fir_coeff = signal.firwin(numtaps, norm_freq2, window='hamming')
elif Filter_Type=='high':
   fir_coeff = signal.firwin(numtaps, norm_freq2,pass_zero=False, window='hamming')
elif Filter_Type=='bandpass':
   fir_coeff = signal.firwin(numtaps, [norm_freq1, norm_freq2],pass_zero=False,window='hamming')
elif Filter_Type=='bandstop':   
   fir_coeff = signal.firwin(numtaps, [norm_freq1, norm_freq2],window='hamming')

print(fir_coeff)

#=================================================
#       TRACES
#=================================================
plt.subplot(211)
f, h = signal.freqz(fir_coeff, fs=sample_rate)
h_dB = 20*np.log10(abs(h))
plt.plot(f,h_dB)

plt.subplot(212)
h_Phase = np.unwrap(np.arctan2(np.imag(h),np.real(h)))
plt.subplot(212)
plt.plot(f,h_Phase)

# MISE EN FORME
plt.subplot(211)
plt.ylabel("Module $H(e^{j\omega})\ en\ dB$")
plt.xlabel("Frequence [Hz]")
plt.xlim([0, sample_rate/2])
plt.grid()

plt.subplot(212)
plt.ylabel("Phase $H(e^{j\omega})\ en\ rad$")
plt.xlabel("Frequence [Hz]")
plt.xlim([0, sample_rate/2])
plt.grid()

plt.show()
#=================================================
#      MISE EN FORME DANS UN FICHIER
#=================================================
fd = open("FIR.txt", "a")
fd.truncate(0)
print("float h_%s_%d_%d__f32[%d]=\n{" % (Filter_Type,cutoff_1_hz,cutoff_2_hz,numtaps),file=fd)
for i in range(numtaps):
    print("%f" % fir_coeff[i],end =",",file=fd)
print("\n};\n", file=fd)

print("int16_t h_%s_%d_%d__q15[%d]=\n{" % (Filter_Type,cutoff_1_hz,cutoff_2_hz,numtaps),file=fd)
for i in range(numtaps):
    print("%d" % int(fir_coeff[i]*32768),end =",", file=fd)
print("\n};\n", file=fd)

fd.close()
#============================================================

Mise en Oeuvre

Vieillissement des Echantillons

La réalisation d’un filtrage FIR consiste à effectuer un produit scalaire entre un vecteur de coefficients h et un vecteur d’échantillons x afin de calculer un échantillon de sortie y.
Le vecteur x doit être mis à jour tous les Te afin d’accueillir un nouvel échantillon d’entrée.

Calcul des échantillons de sortie

Copier les coefficients générés par le script python ci-dessus dans le projet STM32 (fichier _FIR_coeff.h), et effectuer la modification suivante :
REMARQUE : les coefficients de filtre utilisés dans cette première application sont au format float.

main.c

//##################################################
#define N_FILTER 64 					// A MODIFIER
#define FILTER_COEFFS h_low_0_4500__f32 // A MODIFIER
//##################################################

float32_t x[N_FILTER];

void BSP_AUDIO_SAI_Interrupt_CallBack()
{
 float32_t yn = 0.0;
 int i = 0;

  BSP_LED_On(LED1);

  x[0] = (float32_t)(rx_sample_L);
 
  for (i=0 ; i<N_FILTER ; i++) yn += FILTER_COEFFS[i]*x[i];   // Calcul nouvel échantillon de sortie
  for (i=(N_FILTER-1) ; i>0 ; i--) x[i] = x[i-1]; // vieillissement échantillons
    
  tx_sample_L = (int16_t)(yn);
  tx_sample_R = tx_sample_L ; // traitement mono

  BSP_LED_Off(LED1);

  return;
}

Test du Filtre et mesure du temps de calcul

Pour caractériser le filtre ( Analyse spectrale ) , mettre en entrée un bruit blanc (bruit contenant toutes les fréquences à amplitude égale) :

Faire lire le fichier suivant à VLC ( vérifier que le paramètre AUDIO –> AUDIO DEVICE correspond bien à la carte son ).

Bruit Blanc

Enregistrer la sortie de la carte avec audacity :

[AUDACITY] Après enregistrement, sélectionner la partie à analyser, et faire ANALYSE –> TRACER LE SPECTRE

Q2. Relever le spectre des échantillons de bruit blanc filtré, et noter le temps pour réaliser le filtrage.

Q3. Implémenter un filtre de type Passe Haut, et mettre en évidence la fréquence de Shannon avec l’analyse spectrale.

Amélioration du Temps de Calcul

Adressage Circulaire

Certains architectures de processeurs proposent des registres afin de gérer des buffers circulaires (pour pointer l’échantillon le plus vieux et le plus jeune).
Cela permet d’éviter une boucle pour faire vieillir le vecteur x.
Le microcontrôleur STM32F746 ne propose pas de gestion de l’adressage circulaire dans son architecture, aussi la manipulation d’un index pour désigner l’échantillon le plus vieux prend au moins autant de temps que de faire vieillir le vecteur x.

L’architecture du STM32 ne propose pas de registres spécifiques permettant de gérer efficacement l’adressage circulaire.
Le temps de calcul ne sera donc pas amélioré avec cette méthode.

Capacité SIMD (Single Instruction Multiple Data)

Dans le cas d’un calcul avec des échantillons et des coefficients entiers 16 bits (q15), il est possible de récupérer en mémoire 32 bits correspondants à deux échantillons, et de traiter alors ces deux échantillons simultanément.
Cela permet de diviser par deux le nombre d’itérations de boucle si l’on traite 2 échantillons simultanément (pour le moment ce n’est pas le cas).

Extrait de arm_fir_fast_q15.c:

/* Read the first two samples from the state buffer:  x[n-N], x[n-N-1] */
x0 = *__SIMD32(px)++;
/* Read the first two coefficients using SIMD:  b[N] and b[N-1] coefficients */
c0 = *__SIMD32(pb)++;
/* acc0 +=  b[N] * x[n-N] + b[N-1] * x[n-N-1] */
acc0 = __SMLAD(x0, c0, acc0);

Utilisation des fonctions CMSIS DSP

Effectuer la modification du code ci-dessous. Attention, FILTER_COEFFS est un tableau de coefficients de filtre entiers (q15).

arm_fir_instance_q15 ARM_FIR_Q15;
int16_t state_q15[N_FILTER+1];
 
void BSP_AUDIO_SAI_Interrupt_CallBack()
{
 int16_t x_L = 0;
 int16_t x_R = 0;
 int16_t y_L = 0;
 int16_t y_R = 0;

  BSP_LED_On(LED1);

  x_L = rx_sample_L;
  x_R = rx_sample_R;

  // SIGNAL PROCESSING ALGORITHM
  arm_fir_fast_q15( &ARM_FIR_Q15, &x_L,&y_L, 1);

  tx_sample_L =(int16_t)y_L ;
  tx_sample_R = tx_sample_L;


  BSP_LED_Off(LED1);

  return;
}

int main(void)
{
...
arm_fir_init_q15(&ARM_FIR_Q15,N_FILTER,FILTER_COEFFS,state_q15,1);
stm32f7_wm8994_init(...);
while(1)
{...}
}

Q4. Vérifier le bon fonctionnement du filtre avec une analyse spectrale et noter le temps de calcul.
Quel peut être l’ordre maximal de filtre que l’on peut atteindre ?

Utilisation du DMA

Rappel : Le pipeline dans un processeur

Problème du Rebouclage

Le rebouclage dans un programme ( lié à un for/while ou test ) est source de perte de temps ; en effet cela casse le pipeline.
Un traitement par blocs permet de calculer plusieurs échantillons de sortie lors d’une itération de boucle :

  acc0 =  b[numTaps-1] * x[n-numTaps-1] + b[numTaps-2] * x[n-numTaps-2] + b[numTaps-3] * x[n-numTaps-3] +...+ b[0] * x[0]
  acc1 =  b[numTaps-1] * x[n-numTaps]   +   b[numTaps-2] * x[n-numTaps-1] + b[numTaps-3] * x[n-numTaps-2] +...+ b[0] * x[1]
  acc2 =  b[numTaps-1] * x[n-numTaps+1] + b[numTaps-2] * x[n-numTaps]   +   b[numTaps-3] * x[n-numTaps-1] +...+ b[0] * x[2]
  acc3 =  b[numTaps-1] * x[n-numTaps+2] + b[numTaps-2] * x[n-numTaps+1] + b[numTaps-3] * x[n-numTaps]   +...+ b[0] * x[3]

( cf arm_fir_fast_q15.c )

La mise en oeuvre de ce traitement nécessite l’utilisation d’un DMA (Direct Memory Access), qui fait le lien entre le périphérique I2S et la mémoire, afin que le processeur puisse procéder à des calculs simultanément.
Cela implique une certaine latence entre les échantillons d’entrée et ceux de sortie.

En régime permanent :

Remplacer le fichier main.c par main_dma.c ( dans STM32CUBEIDE, faire clic droit sur main.c -> Properties : C/C++Build : cocher ‘Exclude ressource from build’, et décocher cette option pour le fichier main_dma.c )

Attention, FILTER_COEFFS est un tableau de coefficients de filtre au format float.

Q5. Mesurer le temps de calcul permettant désormais de calculer PING_PONG_BUFFER_SIZE échantillons de sortie dans la fonction process_buffer.

Q6. Définir un filtre passe bande d’ordre 256 permettant de valider le fonctionnement du DMA.

Calcul des coefficients du Filtre dans la STM32

REMARQUE : Reprendre le fichier main.c et non main_dma.c pour la suite.

Le Keyboard Tracking consiste à adapter la fréquence de coupure du filtre en fonction de la note.

Cela implique de calculer les coefficients du filtre dans le microcontrôleur.

Q7. Compléter la fonction FIR_calc_coeff_f32 ( dans le fichier FIR_filter.c ) afin de calculer les coefficients d’un filtre FIR passe bas en fonction d’une fréquence de coupure souhaitée.

Q8. Réaliser un synthétiseur à synthèse soustractive avec ce filtre FIR, le paramètre k doit être réglable avec un potentiomètre du clavier maître.

Synthese Sonore Mise en oeuvre d'un filtre IIR