Fabrication d'un filtre passe-bas

**Flo Flo** · 29/06/2010, 18h15

Bonjour à tous

Malheureusement pour moi, j'ai un boulot dans lequel j'ai pas mal de choses à faire avec la programmation et c'est un domaine dans lequel je ne suis vraiment pas performant dans la mesure où personne ne m'a jamais donné de cours là dessus. Par conséquent, ma seule chance de progresser et d'arriver ensuite à me débrouiller par moi même est de me tourner vers cette communauté, qui m'a déjà sorti des épines du pied à plusieurs reprises (et je vous en remercie).

Mon problème actuel consiste à créer un filtre passe bas en langage C. En gros, la situation se présente de la manière suivante:

J'ai un fichier d'entrée avec deux colonnes, la première étant du temps et la seconde des mesures (en l'occurrence magnétiques). Il se trouve que mon signal magnétique a une furieuse tendance a danser dans tous les sens et je voudrais le calmer un peu, autrement dit filtrer les fréquences supérieures à environ 50Hz pour ne laisser apparaître que les grandes tendances. Malheureusement, comme il s'agit de données expérimentales, elles ne répondent pas à une équation particulière. Je ne peux donc pas dire que mon signal est du type Y = A * sin (w.t + Phi) + B ou quoi que ce soit d'autre.

Étant complètement débutant en C, je ne sais pas du tout par où aborder le problème. Quelqu'un aurait il des idées sur la marche à suivre?

Je vous remercie tous par avance.

**haraelendil** · 29/06/2010, 18h27

Il s'agit donc juste de faire un filtre passe bas numérique "tout simple", après le type et l'ordre dépend de ton application, des restrictions quand au temps de calcul ou à la puissance de la machine?

Sinon après le problème n'est pas l'implémentation en C mais surtout la conception du filtre en elle même, as-tu des notions/connaissances la dessus?

**Flo Flo** · 30/06/2010, 10h50

Bonjour,

Je te remercie pour ta réponse. Mes connaissances dans ce domaine sont très limitées. On m'en avait vaguement parlé quand j'étais encore étudiant, mais ça remonte à pas mal d'années et quand on ne pratique pas, on oublie

En fait, oui, il s'agit d'un filtre simple, mais par manque de connaissances, aussi bien d'un point de vue théorique que de programmation, je ne sais pas par où aborder le problème.

**haraelendil** · 30/06/2010, 11h01

Pour le début, le plus simple est d'utiliser un simple filtre FIR (filtre à réponse impulsionnelle finie) qui est un filtre du type y(n) = a.s(n) + b.s(n-1) + ...
et ainsi de suite selon l'ordre que tu désire (plus grand sera l'ordre du filtre, plus abrupte sera la coupure, mais ça demande plus de mémoire et de calcul).
Après, quoi mettre dans a et b et ... , je te laisse un peu chercher, sinon c'est pas drôle et ça n'a aucun intérêt.

**jowo** · 30/06/2010, 13h18

Bonjour,

Un logiciel comme MatLab ne serait-il pas préférable à un programme écrit en c pour résoudre ce type de problème?

**jabbounet** · 30/06/2010, 13h35

Salut,

un livre comme celui ci peux t'intéresser.
http://www.nrbook.com/a/bookcpdf.php

par contre il date un peu et est en anglais.

**Flo Flo** · 30/06/2010, 13h41

C'est vraiment cool de recevoir toute cette aide sur ce site.

Jowo, malheureusement, il se trouve que je ne dispose pas de licence Matlab sur mon ordinateur et ça ne devrait pas changer à court terme

. Je ne peux qu'avoir recours au langage C.

Jabbounet, merci pour le lien. J'ai commencé à le lire.

Haraelendil, j'apprécie vraiment ta manière de me venir en aide. C'est parfait. Ca m'oblige à cogiter, ce qui est la meilleure solution.

Il me semble donc que les coefficients dépendent de l'ordre du filtre que je veux obtenir et qu'ils sont définis en fonction du rapport entre la fréquence de coupure (probablement aux alentours de 50 Hz en ce qui me concerne) et la fréquence de Niqvist (Je ne suis pas sûr de l'orthographe), qui elle même vaut 1/T avec T le pas d'échantillonnage. Je suis complètement à l'ouest ou ça va à peu près?

**haraelendil** · 30/06/2010, 15h09

Oui c'est à peu près ça, le nombre de coefficient que tu aura dépendra de l'ordre du filtre, un filtre d'ordre un s'écrit y(n)=a.x(n)+b.x(n-1), un filtre d'ordre deux s'écrit y(n)=a.x(n)+b.x(n-1)+c.x(n-2) et ainsi de suite.

Ensuite la valeur de tes coefficients dépend principalement de trois choses: le type de filtre que tu veux, la fréquence de coupure et la fréquence d'échantillonnage de ton système.

**Flo Flo** · 30/06/2010, 15h41

OK. Voici déjà un premier point d'éclairci.

Pour ce qui concerne le filtre, je pense qu'un ordre 2 serait raisonnable. Quant aux fréquences d'échantillonnage et de coupure, je les connais. Elles sont respectivement de 0.20Hz (c'est à dire un échantillon toutes les 5 secondes) et de 50Hz.

Comme je ne cherche pas à faire un filtre particulièrement perfectionné, je ne pense pas avoir besoin de quoi que ce soit d'autre. J'aurais tendance à commencer le code comme ça:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define BOOL int
#define TRUE 1
#define FALSE 0
main()
{
	char Div_lu[500];
        int iDiv = 100;
	FILE *ID, *OF;
        float time_Div;
        double X, Y, Z, Tot;
        BOOL FIRST=TRUE;
	ID = fopen("Dive.txt", "r");
	OF = fopen("Filtre.txt", "w");
/*Question: Comment arriver à faire comprendre au programme que le signal sur lequel il doit travailler n'est pas défini par une équation, mais qu'il est composé de toutes les valeurs lues dans le tableau de départ?
 
/*Calcul de la transformée de Fourier du signal:
 
int[] ReSig = new int[N/2];
int[] ImSig = new int[N/2];
for(int f=0;f<(N/2);f++) {
	for(int i=0;i<N;i++) {
		double w = 2*Math.PI*(double)i/N;
		ReSig[f]+=signal[i]*Math.cos(f*w);  
/* Ici, il faut définir une variable signal, qui correspond aux valeurs lues dans le fichier d'entrée.
		ImSig[f]-=signal[i]*Math.sin(f*w);
	}
	ReSig[f]/=(N/2);
	ImSig[f]/=(N/2);
}
 
/* Transformée de Fourier du filtre:
 
double alpha=50.0;
double[] ReFil = new double[N/2];
double[] ImFil = new double[N/2];
for(int f=0;f<(N/2);f++) {
	double w = 2*Math.PI*f;
	ReFil[f]=alpha*alpha/(alpha*alpha+w*w);
	ImFil[f]=-w*alpha/(alpha*alpha+w*w);
}

Après, je verrais bien un filtrage en multipliant les parties réelles et imaginaires entre elles (respectivement).

Puis une transformée de Fourier inverse.

J'imagine qu'il y a des tonnes de trucs qui ne vont pas dans ce que j'ai posté ici. C'est un mix que j'ai fait entre le peu que je connais en programmation et des bouts de code récupérés sur le Net. Et il y a des toutes façons des choses que je ne sais pas faire, comme par exemple le coup du signal, qui est issu d'un tableau extérieur.

**haraelendil** · 30/06/2010, 15h53

Déjà un gros problème d'ordre théorique, il t'es impossible de filtrer à 50Hz un signal échantillonné à 0.2 Hz...

Par définition, un signal numérique échantillonné à 0.2Hz ne contient des informations que sur la bande de fréquence de 0 à 0.1Hz (Te/2 => théorème de Shannon-Nyquist), en partant du principe que le signal a été conditionné à l'échantillonnage pour éviter tout repliement de spectre.

Donc, si ton signal est tel quel et que tu ne peux ré-effectuer d'acquisition, tu ne peux pas (et ça n'a pas de sens) le filtrer à 50Hz...
Si tu veux nettoyer ton signal d'une bonne partie des variations que tu vois, tu peux le filtrer à 0.05Hz par exemple (ne sachant pas le domaine dans lequel tu travaille et ce que tu mesure, je ne sais pas si c'est approprié ou pas)

**Flo Flo** · 30/06/2010, 16h10

Oui, maintenant que tu le dis, ça paraît évident. J'ai un peu fumé la moquette, sur ce coup là.

Après ré examen du look du signal, 0.05Hz paraît convenable. Par contre, je ne peux pas faire de ré-acquisiton, parce que, pour tout te dire, ce sont des données magnétiques qui ont été acquises lors d'une campagne en mer il y a plusieurs années. Je fais un nouvel essai:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define BOOL int
#define TRUE 1
#define FALSE 0
main()
{
	char Div_lu[500];
        int iDiv = 100;
	FILE *ID, *OF;
        float time_Div;
        double X, Y, Z, Tot;
        BOOL FIRST=TRUE;
	ID = fopen("Dive.txt", "r");
	OF = fopen("Filtre.txt", "w");
/*Question: Comment arriver à faire comprendre au programme que le signal sur lequel il doit travailler n'est pas défini par une équation, mais qu'il est composé de toutes les valeurs lues dans le tableau de départ?*/
 
/*Calcul de la transformée de Fourier du signal:*/
 
int[] ReSig = new int[N/2];
int[] ImSig = new int[N/2];
for(int f=0;f<(N/2);f++) {
	for(int i=0;i<N;i++) {
		double w = 2*Math.PI*(double)i/N;
		ReSig[f]+=signal[i]*Math.cos(f*w);  
/* Ici, il faut définir une variable signal, qui correspond aux valeurs lues dans le fichier d'entrée.*/
		ImSig[f]-=signal[i]*Math.sin(f*w);
	}
	ReSig[f]/=(N/2);
	ImSig[f]/=(N/2);
}
 
/* Transformée de Fourier du filtre:*/
 
double alpha=0.05;
double[] ReFil = new double[N/2];
double[] ImFil = new double[N/2];
for(int f=0;f<(N/2);f++) {
	double w = 2*Math.PI*f;
	ReFil[f]=alpha*alpha/(alpha*alpha+w*w);
	ImFil[f]=-w*alpha/(alpha*alpha+w*w);
}

**haraelendil** · 30/06/2010, 16h29

Ba maintenant, il faut juste que tu détermine tes trois déficients à appliquer (si on part sur de l'ordre deux), après, pour réellement filtrer le signal, c'est très simple:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
double e[N]; //ton signal
double s[N]; //ton signal filtré
s[0] = 0; s[1] = 0;
for(unsigned int i=2; i<N; i++)
s[i]=a*e[i]+b*e[i-1]+c*e[i-2];

Avec N le nombre d'échantillons que tu as dans ton signal.

Après, comme je l'ai dit, le plus dur reste de trouver les valeurs des coefficients...

Invité · 30/06/2010, 16h46

Le fait que ce signal aille dans tous les sens veut peut-être dire que le vrai signal samplé à la base contenait des fréquence supérieures à 0.1 hz, qui se seraient "repliées" et ont bousillé ton signal. Est-ce qu'on peut voir la gueule de tes données ?

Déjà, dans quelle bande de fréquence pense-tu que le signal qui t'intéresse soit ? La question c'est : y'a t'il vraiment quelque chose à récupérer ?

Ensuite, tu te mets à utiliser une FFT ce qui est overkill pour ton cas et très sujet à erreurs. Je te conseillerais de ne pas utiliser la FFT.

Il existe deux filtres tout simple que tu peux faire :

1-pole lowpass IIR
output[i] = output[i - 1] * A + input[i] * (1 - A)
ou A va déterminer ta fréquence de coupure.
Si A = 0.5 par exemple dans ton cas tu coupe a 0.05hz (cette coupure n'est jamais nette).
Par contre tu va avoir une réponse en phase pas linéaire ce qui concrètement veut dire que les basses fréquences vont prendre du "retard" (est-ce acceptable ?).

1-pole lowpass FIR
Par exemple
output[i] = input[i] * 0.5 + input[i - 1] * 0.5
Ca ne distord pas la phase, il coupe les fréquence avec la même pente que l'autre mais moins proprement.

Invité · 30/06/2010, 16h49

/*Question: Comment arriver à faire comprendre au programme que le signal sur lequel il doit travailler n'est pas défini par une équation, mais qu'il est composé de toutes les valeurs lues dans le tableau de départ?*/

Le programme n'a rien à comprendre, il fait exactement ce que tu lui dit et n'accorde de signification à rien.

**Flo Flo** · 30/06/2010, 17h18

OK,

Merci beaucoup.

Je vais essayer de bidouiller avec ça. Je reviendrai sûrement dans pas très longtemps avec ce que j'aurai essayer de pondre.

A plus

**Flo Flo** · 30/06/2010, 17h26

Excuse moi ponce, je répondais à haraelendil

Bien sûr tu peux voir mes données.

Je vais les poster sous la forme d'un fichier texte.

Il y a pas mal de colonnes dans ce fichier. La première est le temps et les composantes X, Y et Z du champ sont respectivement les colonnes 4, 5 et 6. La dernière correspond au champ total, qui est égal à la moyenne quadratique des 3 composantes

**diogene** · 30/06/2010, 17h29

Visiblement, le problème est mal défini. Cette histoire de 50Hz et de fréquence d'échantillonnage de 0.2Hz le montre.
La première question avant de filtrer le signal parasite est de connaitre ses caractéristiques, au minimum sa fréquence et son étendue spectrale.
Il faudrait commencer par examiner la transformée de Fourier du signal total pour caractériser le signal parasite.

@ponce

Ensuite, tu te mets à utiliser une FFT ce qui est overkill pour ton cas et très sujet à erreurs. Je te conseillerais de ne pas utiliser la FFT.

Il n'utilise pas une FFT mais la formule brute de la transformée de Fourier discrète (DFT) ce qui est très pénalisant effectivement en temps de calcul. Mais il est facile de trouver des FFT sur le Net.
Je ne comprend pas ta remarque sur les erreurs dues à l'utilisation de la DFT.

**Flo Flo** · 30/06/2010, 17h30

Tu vas voir, si tu traces les composantes en fonction du temps, il y a des grandes tendances qui se dégagent (autrement dit à basse fréquence / grande longueur d'onde) et du bruit à fréquence beaucoup plus élevée. C'est ce bruit dont j'aimerais m'affranchir pour avoir un signal net.

PS, JE te remercie pour les exemples de filtres que tu as indiqués

**Flo Flo** · 30/06/2010, 17h35

Ouh là, j'ai tellement d'aide que je n'arrive même plus à poster des messages suffisamment rapidement pour tenir compte de tout ce qui arrive.

Mais je ne vais certainement pas m'en plaindre

Diogene, effectivement, les 50Hz étaient foireux. C'est moi qui ai dit une bêtise, qu'Haerendil a vue et m'a faite corriger. Le bruit semble être autour de 0.1Hz environ

Invité · 01/07/2010, 15h56

Envoyé par diogene

Il n'utilise pas une FFT mais la formule brute de la transformée de Fourier discrète (DFT) ce qui est très pénalisant effectivement en temps de calcul. Mais il est facile de trouver des FFT sur le Net.
Je ne comprend pas ta remarque sur les erreurs dues à l'utilisation de la DFT.

Effectivement c'est une DFT.

J'avais l'impression que Flo Flo veut appliquer son filtre dans le domaine de Fourier d'où ma remarque.

Fabrication d'un filtre passe-bas

C

Discussions similaires

Partager

Partager