Discussion :

[Snark]

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]

// 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)

[Snark]

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]

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]

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]

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 ?

[acx01b]

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)