[Signal] Implémentation filtres Passe-Bas en C

**Snark** · 26/12/2008, 16h20

Bonjour à tous,

Pour des besoins professionnel, je dois filtrer un signal issue de capteurs à l'aide d'un filtre passe-bas. Malheureusement, j'ai un peu de mal à mettre ça en place, malgré toutes les explications que j'ai trouvé sur le net.

Voici en gros, mon programme. Y'a un truc que j'ai pas dû comprendre:

Code C :

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
 
// Creation des buffers
double* hsupp = new double[L];
double* filter = new double[L];
double* result = new double[M];
double* fft_signal = new double[L];
double* fft_filter = new double[L];
double* fft_result = new double[M];
 
// Creation du filtre passe-bas (avec un algo trouvé sur le net)
for (int i=-(L-1.0)/2.0; i<=(L-1.0)/2.0; i++, j++)
	hsupp[j] = i;
for (int i=0; i<L; i++)
	filter[i] = (2.0*fc/Fs)*sinc(2.0*fc*hsupp[i]/Fs);
 
// Calcul des FFT signal et filtre
memcpy(fft_signal, x, M*sizeof(double));
LabCore::rsrfft(fft_signal, OFFSET);
 
memcpy(fft_filter, filter, L*sizeof(double));
LabCore::rsrfft(fft_filter, OFFSET);
 
// Convolution des 2 FFT (à mon avis, c'est ici qu'il y a un problème !)
for (int i=0; i<M; i++)
	fft_result[i] = fft_filter[i] * fft_signal[i];
 
// Calcul de la transformé inverse
memcpy(result, fft_result, M*sizeof(double));
LabCore::irsrfft(result, OFFSET);
 
// Normalisation de FFT (juste pour l'affichage)
fft_signal = normalize(fft_signal, L);
fft_filter = normalize(fft_filter, L);
fft_result = normalize(fft_result, L);
for (int i=0; i<L; i++)
	printf("\"%.09f\";\"%.09f\";\"%.09f\";\"%.09f\";\"%.09f\";\"%.09f\"\n", x[i], hideal[i], result[i], fft_signal[i], fft_filter[i], fft_result[i]);

Au final, j'ai un truc qui n'a rien à voir avec mon signal, et qui, en plus, est devenu symétrique par le centre (ce qui n'est pas du tout le cas de mon signal d'origine). Si on pouvait me dire là où ça cloche

Merci beaucoup.
Snark.

PS: oui, le code est un peu dégueux, mais c'est juste du test...

**pseudocode** · 27/12/2008, 12h51

Envoyé par Snark

// Convolution des 2 FFT (à mon avis, c'est ici qu'il y a un problème !)

Non, ça c'est bon. Le produit de convolution dans le domaine temporel se traduit en une multiplication dans le domaine fréquentiel.

Au final, j'ai un truc qui n'a rien à voir avec mon signal, et qui, en plus, est devenu symétrique par le centre (ce qui n'est pas du tout le cas de mon signal d'origine). Si on pouvait me dire là où ça cloche

Vu ton résultat, je dirais que tu n'as pas fait de zéro-padding (=rajouter des zéros dans tes 2 tableaux pour avoir une taille qui est une puissance de 2)

**acx01b** · 27/12/2008, 13h00

salut

- combien vaut fc dans ton exemple, et Fs ?

pour le zero padding je ne suis pas convaincu, la lib que tu utilises doit le gérer non ?
- combien valent L et M ? il devraient être égaux logiquement

- essaye un (vrai) filtre passe bas de ce type:

Code :

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

1
2
3
for (int i=0; i<L; i++) fft_filter[i] = 0;
for (int i=0; i<L*0.1; i++)
	fft_filter[i] = fft_filter[L-1-i] = 1;

(0.1 <--> fc = Pi/5)

**Snark** · 29/12/2008, 14h45

Envoyé par pseudocode

Non, ça c'est bon. Le produit de convolution dans le domaine temporel se traduit en une multiplication dans le domaine fréquentiel.

Ok, mais j'avais quand même un doute

Vu ton résultat, je dirais que tu n'as pas fait de zéro-padding (=rajouter des zéros dans tes 2 tableaux pour avoir une taille qui est une puissance de 2)

J'utilise directement des tableaux en puissances de 2, pour le moment. Histoire de pas me compliquer la vie.

pour le zero padding je ne suis pas convaincu, la lib que tu utilises doit le gérer non ?

Je n'utilise aucune libs externe. La FFT est directement incluse dans mon prog.

- combien valent L et M ? il devraient être égaux logiquement

Exacte, juste un reste d'essais

- essaye un (vrai) filtre passe bas de ce type:

La FFT de mon filtre et bonne, donc je ne pense pas que ça change grand chose.

Bon, étant donnée que je ne m'en sort pas, je fourni directement mes sources, avec un peu de chances, l'un de vous pourra me débloquer... Ce doit être un truc tout con :
http://blary.jason.free.fr/download/ffttest.tar.bz2

C'est compilable directement sans aucune dépendances externe. Si vous être sous nux, il suffit de taper "make".
Ça renvoie sur la sortie standard un fichier CSV. On peut l'ouvrir avec un tableur pour pouvoir tracer les courbes.

Merci beaucoup, car, là, je sèche vraiment

Et j'en ai besoin.
Snark

**acx01b** · 31/12/2008, 14h38

salut

j'ai trouvé le problème

rajoute cette fonction:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
void fftconvolution(double *a, double *b, double *c, int n)
{
    a[0] = b[0]*c[0];
    a[n/2] = b[n/2]*c[n/2];
    for (int i = 1; i < n/2; i++)
    {
        // multiplication complexe
        a[i] = b[i]*c[i] - b[n-i]*c[n-i];
        a[n-i] = b[i]*c[n-i] + b[n-i]*c[i];
    }
}

et dans ton main:

Code :

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

1
2
// Convolution
fftconvolution(fft_result,fft_signal,fft_filter,M);

au fait ton normalize il n'a pas trop lieu d'être là, en plus il y a un dépassement de tableau dans cette fonction

**Snark** · 16/01/2009, 17h36

Yop,

Je répond un peu tard, désolé. Merci beaucoup pour ton aide.
En effet, c'était une erreur con, j'avais oublié comment on multipliais deux complexe.

Par contre, pour x[0], je suis d'accord pour la multiplication, car ça correspond à l'offset global du signal, mais pour x[n/2], je ne comprend pas vraiment à quoi il correspond.

Bon, avec une convolution correct, ça fonctionne un poil mieux, mais c'est toujours étrange. Après la transformé inverse, j'ai la première moitié du signal à la fin, et la seconde moitié au début.

Il retrouve son sens correct si je corrige la convolution comme cela (c'est pas propre du tout, je sais....) :

Code C :

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
 
void fftconvolution(double *a, double *b, double *c, int n)
{
	a[0] = b[0]*c[0];
	a[n/2] = b[n/2]*c[n/2];
	for (int i=1; i<n/2; i++)
	{
		// multiplication complex
		a[i] = b[i]*c[i] - b[n-i]*c[n-i];
		a[n-i] = b[i]*c[n-i] + b[n-i]*c[i];
	}
 
	// Correction pour avoir le sens correct
	for (int i=1; i<n; i+=2)
		a[i] = -a[i];
}

Mais du coup, j'obtiens des harmoniques supplémentaires sur la fin et le début du signal (voir images). Ça ne m'arrange pas vraiment.

Avant filtrage:

Après filtrage:

Si quelqu'un à une idée de l'origine de mon problème (bon, le soucis d'harmoniques, je pense savoir que ça vient de ma correction... vue que j'ai juste mis un offset sur le signal), je le remercierai beaucoup, voilà 15jrs que je suis là dessus, et j'avance pas.

Je fourni toujours mon programme (version corrigé) si ça peut aider quelqu'un (enfin, surtout moi

)
http://blary.jason.free.fr/download/ffttest.tar.bz2

Snark.

**Snark** · 19/01/2009, 10h29

Hello,

C'est encore moi avec mon problème que je profite pour faire un up.

Bon, j'ai changé un peu la technique de convolution, et ça met un peu plus en évidence mon problème. En faîte, c'est dû à un « effet de bord » à cause du zero padding lors de la convolution:

Y aurait-t-il un moyen pour l'atténuer au maximum ?

**pseudocode** · 19/01/2009, 11h33

C'est moi ou alors tu n'as pas utilisé de fonction de fenêtrage dans ton code ?

**Snark** · 19/01/2009, 14h42

Envoyé par pseudocode

C'est moi ou alors tu n'as pas utilisé de fonction de fenêtrage dans ton code ?

Hein ?! Quoi ?

Je crois qu'il me manque une connaissance là.

En cherchant, j'ai trouvé ça : http://matthieu-brucher.developpez.c...lgo/fft/#LII-B
Par contre, j'ai un peu de mal à être certain de comprendre à quoi ça sert, ce qu'il appelle les « lobes secondaires » c'est les trucs en plus sur les bords ?
Et ça s'utilise juste comme faisant une convolution avec mon signal AVANT de refaire une convolution avec mon filtre ?

J'apprends sur le tas, donc désolé de mes lacunes...

Snark.

**FR119492** · 19/01/2009, 14h45

Salut!
Physiquement, un filtre est un biporte constitué d'un certain nombre d'éléments (résistances, capacités, sources commandées). Le filtre passe-bas le plus simple (du premier ordre) comporte une résistance en série avec une des bornes d'entrée et une capacité en parallèle avec les bornes de sortie. Il te suffira donc d'imaginer le filtre qui aura les caractéristiques désirées, d'écrire ses équations différentielles et d'intégrer ces équations par une méthode pas à pas comme celle de Runge-Kutta.
Jean-Marc Blanc

**pseudocode** · 19/01/2009, 15h17

Envoyé par Snark

Hein ?! Quoi ?

Je crois qu'il me manque une connaissance là.

Ohhhh tu va découvrir les joies du "spectral leakage".

Et ça s'utilise juste comme faisant une convolution avec mon signal AVANT de refaire une convolution avec mon filtre ?

Oui.

Mais grâce la magie de la convolution, tu peux intégrer la fonction de fenêtrage dans ton filtre : tu multiplies terme à terme ton filtre et la fonction de fenêtrage (centrée au milieu du filtre). Ensuite tu utilises ton filtre comme d'habitude (zeropadding, FFT, ...)

Code C :

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

1
2
3
4
5
6
7
8
9
10
 
// Creation du filtre passe-bas 
...
 
// multiplication par la fenetre
for (int i=0; i<L; i++)
	filter[i] *= window(i,L);
 
// Calcul des FFT signal et filtre
...

**acx01b** · 21/01/2009, 15h30

salut
je crois que surtout vu que ce qui t'intéresse c'est une simple convolution temporelle, autant le voir du point de vue temporel: ton filtre est trop grand, il a trop de coefficients

pour t'en convaincre, fais un signal
signal[] = {1,2,-1,-1,0,0,0,........................}
et calcule signal convolué par ton filtre
puis calcule signal convolué par un filtre constitué de très peu de coefficients par exemple filter[] = {0.25,0.5,1,0.5,0.25,0,0,............}

la solution serait alors de calculer les coefficients d'un filtre passe bas d'ordre < 10 (butterworth ou elliptique par exemple) sous matlab .

**j.p.mignot** · 25/01/2009, 17h56

voir le petit soft ci-joint.
Il est basé sur l'appelet

http://www.dsptutor.freeuk.com/remez...terDesign.html

**popof85** · 09/02/2009, 18h16

Bonsoir.
Je suis aussi à la recherche d'une manière de programmer un filtre passe-bas(je dois analyser une courbe de type sinusoîdal). Je débute en langage C et je n'ai aucune idée de la manière de programmer ce filtre...
Est-ce que je peux réutiliser les mêmes commandes (en changeant les fréquences) que celles citées plus haut?
Merci d'avance!

**FR119492** · 10/02/2009, 00h16

Salut!

je dois filtrer un signal

Je reviens sur mon intervention précédente. Il y a une chose que les "traiteurs de signal" ont souvent la fâcheuse habitude d'oublier; un signal est un phénomène physique et le filtrage est un traitement physique qu'on lui applique. Le problème posé doit donc être traité comme tel.
Jean-Marc Blanc

**pseudocode** · 10/02/2009, 09h47

Envoyé par FR119492

Je reviens sur mon intervention précédente. Il y a une chose que les "traiteurs de signal" ont souvent la fâcheuse habitude d'oublier; un signal est un phénomène physique et le filtrage est un traitement physique qu'on lui applique. Le problème posé doit donc être traité comme tel.

Tu veux dire simuler le fonctionnement d'un circuit RC ?

http://www.dspguide.com/ch13/3.htm

**popof85** · 10/02/2009, 11h46

D'un point de vue physique, je comprends bien le fonctionnement et l'utilité d'un filtre.. Mais lorsqu'il s'agit de le traduire en langage C ....

**henderson** · 10/02/2009, 13h42

Salut !
Dans un filtre RC, le condensateur joue le role d'un amortisseur, comme s'opposant à toute variation raide du signal et ce, en fonction de sa constante de temps (et aussi de sa manière à lui de se charger ... ...).
Tu peux en modéliser le mécanisme d'échantillon en échantillon via un traitement par dichotomie, s'en approchant ou s'en éloignant.

soit p la valeur de l'échantillon (passé et calculé)
soit e la valeur de l'échantillon lu
soit r le coefficient de l'amortissement

On commence avec p = 0
La valeur pour chaque échantillon lu est :

p = e + (e-p) * r

Lorsque r = 0.5 le traitement est dichotomique (selon le principe de Zénon d’Élée)

La fréquence de coupure du filtre (si tenté que celà en soit un) resterait à calculer (à partir de r et de la fréquence d'échantillonnage)
Mais, n'étant pas mathématicien dans l'âme ...

A plus !

**FR119492** · 10/02/2009, 13h58

Salut!

Si tu notes u1 la tension d'entrée de ton filtre et u2 sa tension de sortie, tu as:

i = C * du2/dt
u2 = u1 - R * i

Tu élimines i entre ces deux équations:

u2 = u1 - R*C * du2/dt

D'où tu tires:

du2/dt = (u1 - u2) / (R*C)

que tu intègres, par exemple par Runge-Kutta.

Jean-Marc Blanc

**pseudocode** · 10/02/2009, 15h01

Envoyé par popof85

D'un point de vue physique, je comprends bien le fonctionnement et l'utilité d'un filtre.. Mais lorsqu'il s'agit de le traduire en langage C ....

Bon, je vais essayer de faire une implémentation au raz-des-paquerettes en java.

1. Fabriquons un joli signal sinusoidale bruité.

signal = 128*cos(5*w) + 32*bruit_entre_-1_et_1

Code java :

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

1
2
3
4
5
6
7
8
9
int N=800;
int[] signal = new int[N];
for(int i=0;i<N;i++) {
	double t = (double)i/N;
	double w = 2*Math.PI*t;
	int input = (int) (128*Math.cos(5*w));
	int noise = (int) (32*(2*Math.random()-1));
	signal[i] = input+noise;
}

2. on calcule la transformée de fourier complexe du signal, par la methode brutale

Code java :

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

1
2
3
4
5
6
7
8
9
10
11
12
// brutforce Fourier Transform of the 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);
		ImSig[f]-=signal[i]*Math.sin(f*w);
	}
	ReSig[f]/=(N/2);
	ImSig[f]/=(N/2);
}

3. on crée un filtre complexe, en calculant manuellement sa TF (cf. http://www.dspguide.com/ch13/3.htm)

Code java :

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

1
2
3
4
5
6
7
8
9
10
11
12
// Fourier Transform of the filter H(t) = alpha*Math.exp(-alpha*t)
// ==>
// Re(H(w)) = alpha^2/(alpha^2+w^2)
// Im(H(w)) = -alpha*w/(alpha^2+w^2)
double alpha=40.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);
}

4. on effectue le filtrage, c'est à dire une multiplication (complexe) terme a terme entre le signal et le filtre

Code java :

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

1
2
3
4
5
6
7
// output=filter*signal (warning: multiplication of two complex)
int[] ReOutput = new int[N/2];
int[] ImOutput = new int[N/2];
for(int f=0;f<N/2;f++) {
	ReOutput[f]=(int)(ReSig[f]*ReFil[f]-ImSig[f]*ImFil[f]);
	ImOutput[f]=(int)(ReSig[f]*ImFil[f]+ImSig[f]*ReFil[f]);
}

5. On repasse dans le domaine temporel, en calculant la transformée de Fourier INVERSE du signal filtré.

Code java :

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

1
2
3
4
5
6
7
8
// brutforce inverse-Fourier Transform of the output (real part only)
int[] output = new int[N];
for(int f=0;f<(N/2);f++) {
	for(int i=0;i<N;i++) {
		double w = 2*Math.PI*(double)i/N;
		output[i]+=ReOutput[f]*Math.cos(f*w)-ImOutput[f]*Math.sin(f*w);
	}
}

Et voila. C'est fini.

En rouge: le signal (input+noise)
En vert: le signal filtré

[Signal] Implémentation filtres Passe-Bas en C

Traitement du signal

Discussions similaires

Partager

Partager