Filtrons avec Octave

Chipotons un peu avec GNU Octave (MATLAB)... Je m'étais intéressé, il y a quelques mois à la détection des signaux faisant basculer les compteurs électriques bihoraires au moyen d'un simple fil sur le connecteur microphone d'un PC. Je vais utiliser un bout de fichier pour faire du 'traitement de signal' (filtrage numérique).

$ echo 'date "+%g%m%d.%H%M%S"' > now; chmod +x now
$ arecord -f S16_LE -t raw -r 25000 > 25000.`./now`
$ dd if=25000.131204.064758 of=sample.raw bs=1000 count=5 skip=1
$ sox -r 25000 -e signed -b 16 -c 1 sample.raw sample.wav

Et dans Octave,

x = wavread('sample.wav');
y = x - mean(x);   # remove the DC component
plot(y)
print -dpng 'bzzz_50Hz.png'

donne

et

plotspec(y, 1/25000)
print -dpng '50Hz_spectrum.png'

donne

On voit (ou pourrait voir en zoomant avec le bouton de gauche de la souris) que la composante principale du spectre se trouve à +/- 50 Hz mais que l'on trouve des harmoniques avec une magnitude atteignant 10% de la composante principale. Ce n'est probablement pas la bonne manière de faire (?) mais le but est ici de jouer avec les fonctions de filtrage de MATLAB/GNU Octave. Par exemple, éliminons les harmoniques avec un FIR :

b = fir1(120, 0.0035);
z=filter(b, 1, y);
clf;
plotspec(z, 1/25000);
print -dpng 'spectrum_50Hz_filtered_FIR_120_0.0035.png'

donne

La réponse fréquentielle du filtre nous est donnée par :

clf;
freqz(b, 1, 512, 'half', 25000);
print -dpng 'filter_FIR_120_0.0035.png'

Essayons d'isoler l'harmonique à 250 Hz en concevant un FIR avec remez() (firpm() dans MATLAB. On spécifie la réponse fréquentielle approximative que l'on désire au moyen de deux vecteurs (fréquences/réponses) :

f = [0 0.03 0.04 1];
a = [0 1    0    0];
w = [1 1];
b = remez (120, f, a, w);
clf;
freqz(b, 1,512,'half',25000);
print -dpng 'remez.120_0-1-0-0_0-.03-.04-1_freq.png'
z=filter(b, 1, y);
clf;
plot(z);
hold on;
plot(y, 'r');
print -dpng 'remez.120_0-1-0-0_0-.03-.04-1_result.png'

donne

et

Le filtre n'est pas très sélectif mais le 250 Hz ressort bien (en bleu).

---

...En cours de rédaction / y-a un truc qui cloche...
Le 0.0035 passé à 'fir1()' est à comparer avec 50/(25000/2) ('fréquence normalisée'), on coupe un peu dans le 50Hz, cela aurait dû être supérieur à 0.004 (50/12500) (?). La pente du filtre n'est pas extraordinaire, on est à -10 dB à 250 Hz...