Discussion :

[Cyborg]

Bonjour à tous,

Bon, j'ai un soucis. J'ai essayé de programmer un filtre numérique, basé sur un filtre de Butterworth, en appliquant directement le carré Gain du filtre sur mon signal, en utilisant comme librairie de FFT : FFTW.

Mon soucis ... C'est que ça ne fait pas ce que je souhaite. J'ai un signal cardiaque de fréquence environ 60 bpm (au début en tout cas), soit un signal périodique de fréquence 1Hz. Hors, lorsque j'applique mon filtre avec une frequence de coupure à 8Hz, il ne me reste déjà plus rien de mon signal cardiaque !!

Voici mon code, si vous avec une idée du type d'erreur que j'ai pu faire :

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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
static const unsigned int ordre = 3 ;
 
/******************************************************************************
 ***********************    passe_bande_frequentiel    ************************
 ******************************************************************************
 * But :
 *  - filtre un signal par un filtre passe-bande de Butterworth
 * Entrees :
 *  - p_d_donnees : tableau contenant les donnees a filtrer
 *  - ul_nb_donnees : nombre d'echantillons du tableau
 *  - d_frequence_echantillonnage : frequence d'echantillonnage du signal
 *  - d_frequence_coupure_ph : frequence de coupure du passe-haut
 *  - d_frequence_coupure_pb : frequence de coupure du passe-bas
 * Sorties :
 *  - p_d_donnees_filtrees : tableau contenant les donnees apres filtrage
 *  - code_erreur : code d'erreur 
 *      0 : pas d'erreur
 *     -1 : mauvais parametres
 *     -2 : erreur allocation mémoire
 ******************************************************************************/
void FiltragePasseBandeFFTW ( const double * p_d_donnees ,
                              const unsigned long ul_nb_donnees ,
                              const double d_frequence_echantillonnage ,
                              const double d_frequence_coupure_ph ,
                              const double d_frequence_coupure_pb ,
                              double * p_d_donnees_filtrees ,
                              int * code_erreur )
{
  unsigned long i ;
  unsigned long ul_dim_fft = ( unsigned long ) ( ul_nb_donnees / 2 ) + 1 ;
  double d_gain_ph , d_gain_pb ;
  fftw_complex * fc_sortie_fft ;
  fftw_plan fp_plan_fft ;
 
  if ( code_erreur )
    * code_erreur = 0 ;
 
  if ( ( ! p_d_donnees ) || ( ! p_d_donnees_filtrees ) ||
       ( ! d_frequence_echantillonnage ) || ( ! ul_nb_donnees ) ||
       ( ( d_frequence_coupure_ph <= 0 ) && ( d_frequence_coupure_pb <= 0 ) ) ||
       ( 2 * d_frequence_coupure_ph > d_frequence_echantillonnage ) ||
       ( 2 * d_frequence_coupure_pb > d_frequence_echantillonnage ) )
    { // Controle des frequences de coupures ( il faut que Fc<Fe/2 )
      if ( code_erreur )
        * code_erreur = -1 ;
    }
  else
    {
      fc_sortie_fft = ( fftw_complex * ) fftw_malloc ( ul_nb_donnees * sizeof ( fftw_complex ) ) ;
      if ( ! fc_sortie_fft )
        {
          if ( code_erreur )
            * code_erreur = -2 ;
        }
      else
        {
          // Calculer la fft du signal
          fp_plan_fft = fftw_plan_dft_r2c_1d ( ul_nb_donnees , ( double * ) p_d_donnees , fc_sortie_fft , FFTW_ESTIMATE ) ;
          fftw_execute ( fp_plan_fft ) ;
          fftw_destroy_plan ( fp_plan_fft ) ;
 
          // Convoluer par le filtre Butterworth d'ordre 4, applique dans le sens direct et retrograde pour supprimer la phase ( Hr4 * H4 = Gr4 * G4 = (G4)^2 )
          if ( d_frequence_coupure_ph <= 0 ) // Filtre passe bas uniquement
            {
              for ( i = 0 ; i < ul_dim_fft ; i++ )
                {
                  // w=2.pi.f ; wc = (w*Fe)/(fc*N) ; G4^2 = 1 / ( 1 + wc^(2*ordre) )
                  d_gain_pb = 1.0 / ( 1.0 + pow ( ( 2 * PI * i * d_frequence_echantillonnage ) / ( d_frequence_coupure_pb * ul_nb_donnees ) , 2 * ordre ) ) ;
                  fc_sortie_fft [ i ] [ 0 ] = fc_sortie_fft [ i ] [ 0 ] * d_gain_pb ;
                  fc_sortie_fft [ i ] [ 1 ] = fc_sortie_fft [ i ] [ 1 ] * d_gain_pb ;
                }
            }
          else
            {
              if ( d_frequence_coupure_pb <= 0 ) // Filtre passe haut uniquement
                {
                  for ( i = 0 ; i < ul_dim_fft ; i++ )
                    {
                      // w=2.pi.f ; wc = (w*Fe)/(fc*N) ; G4^2 = 1 / ( 1 + wc^(2*ordre) )
                      d_gain_ph = 1 - ( 1.0 / ( 1.0 + pow ( ( 2 * PI * i * d_frequence_echantillonnage ) / ( d_frequence_coupure_ph * ul_nb_donnees ) , 2 * ordre ) ) ) ;
                      fc_sortie_fft [ i ] [ 0 ] = fc_sortie_fft [ i ] [ 0 ] * d_gain_ph ;
                      fc_sortie_fft [ i ] [ 1 ] = fc_sortie_fft [ i ] [ 1 ] * d_gain_ph ;
                    }
                }
              else
                {
                  // Filtre passe bas bande
                  for ( i = 0 ; i < ul_dim_fft ; i++ )
                    {
                      // w=2.pi.f ; wc = (w*Fe)/(fc*N) ; G4^2 = 1 / ( 1 + wc^(2*ordre) )
                      d_gain_ph = 1 - ( 1.0 / ( 1.0 + pow ( ( 2 * PI * i * d_frequence_echantillonnage ) / ( d_frequence_coupure_ph * ul_nb_donnees ) , 2 * ordre ) ) ) ;
                      d_gain_pb = 1.0 / ( 1.0 + pow ( ( 2 * PI * i * d_frequence_echantillonnage ) / ( d_frequence_coupure_pb * ul_nb_donnees ) , 2 * ordre ) ) ;
                      fc_sortie_fft [ i ] [ 0 ] = fc_sortie_fft [ i ] [ 0 ] * d_gain_ph * d_gain_pb ;
                      fc_sortie_fft [ i ] [ 1 ] = fc_sortie_fft [ i ] [ 1 ] * d_gain_ph * d_gain_pb ;
                    }
                }
            }
 
          // Calculer l'ifft du signal
          fp_plan_fft = fftw_plan_dft_c2r_1d ( ul_nb_donnees , fc_sortie_fft , p_d_donnees_filtrees , FFTW_ESTIMATE ) ;
          fftw_execute ( fp_plan_fft ) ;
          fftw_destroy_plan ( fp_plan_fft ) ;
 
          // Prise en compte du facteur d'echelle
          for ( i = 0 ; i < ul_dim_fft ; i++ )
            p_d_donnees_filtrees [ i ] /= ul_nb_donnees / 2 ;
 
          fftw_free ( fc_sortie_fft ) ;
        }
    }
}

PS1 : Les données ne sont filtrer sur la FFT que jusqu'à N/2+1, puisque FFTW n'utilise que la moitié de la FFT ... si ma mémoire est bonne.
PS2 : Cela fait suite à ma première discussion plus "théorique" ici.

[Aleph69]

Bonjour,

ce n'est par forcément la cause de ton bogue, mais sur ce type de calcul tu prends le risque que d_gain_pb soit nul (division par un trop grand nombre) :

Code :

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

1
2
3
4
5
6
7
8
 
              for ( i = 0 ; i < ul_dim_fft ; i++ )
                {
                  // w=2.pi.f ; wc = (w*Fe)/(fc*N) ; G4^2 = 1 / ( 1 + wc^(2*ordre) )
                  d_gain_pb = 1.0 / ( 1.0 + pow ( ( 2 * PI * i * d_frequence_echantillonnage ) / ( d_frequence_coupure_pb * ul_nb_donnees ) , 2 * ordre ) ) ;
                  fc_sortie_fft [ i ] [ 0 ] = fc_sortie_fft [ i ] [ 0 ] * d_gain_pb ;
                  fc_sortie_fft [ i ] [ 1 ] = fc_sortie_fft [ i ] [ 1 ] * d_gain_pb ;
                }

J'aurais tendance à passer temporairement par un logarithme :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
 
              for ( i = 0 ; i < ul_dim_fft ; i++ )
                {
                  // w=2.pi.f ; wc = (w*Fe)/(fc*N) ; G4^2 = 1 / ( 1 + wc^(2*ordre) )
                  double tmp = (2* PI*i*d_frequence_echantillonnage) / ( d_frequence_coupure_pb * ul_nb_donnees );
                  log_gain_pb=-log(1.0+pow ( tmp , 2 * ordre ) ) ;
                  double log_fc_sortie_fft_i0 = log_fc_sortie_fft [ i ] [ 0 ] + log_d_gain_pb ;
                  double log_fc_sortie_fft_i1 = log_fc_sortie_fft [ i ] [ 1 ] + d_gain_pb ;
                  fc_sortie_fft[i][0] = exp(log_fc_sortie_fft_i0);
                  fc_sortie_fft[i][1] = exp(log_fc_sortie_fft_i1);
                }

Note : j'ai mis des variables temporaires pour te faciliter la lecture.

Je ne suis pas fan du passage à l'exponentielle sur les deux dernières lignes non plus mais bon.
Selon le traitement que tu fais de ta sortie, tu peux éventuellement t'en passer et retourner les logarithmes.

Je regarde si je vois autre chose.

EDIT : en fait tu peux même faire sauter la puissance en calculant le logarithme de 1.0-1.0/d_gain_pb.

[Aleph69]

J'ai peut-être une bonne piste.

Ici,

unsigned long ul_dim_fft = ( unsigned long ) ( ul_nb_donnees / 2 ) + 1 ;

si ul_nb_donnees n'est pas pair j'ai peur que ul_dim_fft vale zéro.
Tu pourrais remplacer par

unsigned long ul_dim_fft = ( unsigned long ) ( 0.5 * ul_nb_donnees ) + 1 ;

où faire de la division euclidienne avec l'opérateur %.

Remarque : à vérifier mais je crois qu'il faut arrondir "ul_nb_donnees/2" à l'entier inférieur dans FFTW.

A la fin de ton code (rescale), tu as un truc similaire qui traîne :

Code :

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

1
2
3
 
for ( i = 0 ; i < ul_dim_fft ; i++ )
            p_d_donnees_filtrees [ i ] /= ul_nb_donnees / 2 ;

Je pense qu'il y a plus sûr :

Code :

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

1
2
3
 
for ( i = 0 ; i < ul_dim_fft ; i++ )
            p_d_donnees_filtrees [ i ] /= ul_dim_fft-1  ;

[Cyborg]

Bonjour, et merci à toi Aleph69 !!

Désolé, je me suis mal fait comprendre. Bon faut dire que j'ai tapé ce message rapidement hier.

Mon problème ne se situe pas là, mais plutôt sur son application : lorsque j'applique un filtre de Butterworth d'ordre 3 coupant à 40Hz mon signal d'ECG, j'obtiens des résultats franchement différents si j'utilise Scilab ou si j'utilise mon code C.

Il n'y a qu'à voir le résultat sur l'image :
- ECG original en noir ;
- ECG filtré par Scilab (passe-bas à 40Hz)
- ECG filtré avec mon code C.
Ca me fait un passe-haut impeccable (à vue d'oeil, je dirais entre 0.5 et 1Hz).
Le soucis ... c'est que ce n'est absolument pas ce que je lui ait demandé !!

Si tu penses que c'est plus la théorie qui pèche que le code C++, faudrait alors demander à un modo de déplacer le topic dans la section algo. Mais à mes yeux, ça colle avec ce que je cherche à obtenir...

[Aleph69]

Bonjour,

effectivement tu perds beaucoup d'information après filtrage!

Ne met pas de côté trop vite les problèmes numériques.
Si tu divises par des nombres trop grands, tu vas localement obtenir des zéros non souhaités et détruire une partie de l'info de ta transfo de Fourier (en gros tu vas "padder" tes coeffs avec des zéros à partir d'un certain indice i).
Du coup, après transfo inverse, tu récupères un signal dénué d'information.
Cela dit, je ne dis pas du tout que c'est le cas ici.

Après, il y a le problème de compatibilité entre les algos qui indicent les FTT entre 0 et N et celles qui le font entre -N/2 et N/2+1 mais je pense que tu l'as géré.


Tu peux m'envoyer la sortie de scilab et la tienne dans un ou deux fichiers de données?
On ne sait jamais, un truc peut me sauter aux yeux.

EDIT : il faut lire entre "-N/2 et +N/2". Désolé pour la bourde.

[Invité]

Bonjour à tous,

Bon, j'ai un soucis. J'ai essayé de programmer un filtre numérique, basé sur un filtre de Butterworth, en appliquant directement le carré Gain du filtre sur mon signal, en utilisant comme librairie de FFT : FFTW.

Mon soucis ... C'est que ça ne fait pas ce que je souhaite. J'ai un signal cardiaque de fréquence environ 60 bpm (au début en tout cas), soit un signal périodique de fréquence 1Hz. Hors, lorsque j'applique mon filtre avec une frequence de coupure à 8Hz, il ne me reste déjà plus rien de mon signal cardiaque !!

Déjà quelle est la frequence d'échantillonage de ton signal ?

Pourquoi veut tu passer par une FFT juste pour faire un filtre ?
Un filtre FIR ou IIR (si la phase est pas trop importante) c'est beaucoup plus simple qu'une FFT.

Et si tu veux vraiment passer par une FFT, utilise plutôt une fonction un peu plus smooth que le carré (genre avec une partie polynomiale en pente douce) ca sera plus exact surtout si tu fais du recoupement.

Pour designer un filtre ce logiciel est très bien : http://www.tobybear.de/p_filterexp.html

[Cyborg]

Ah ah, merci pour le conseil, c'était loin d'être inutile.
J'ai planché dans mon code, affiché les calculs intermédiaires et me suis rendu compte du soucis lié à la fonction de gain.
J'ai repris les fonctions de bases, refait les calculs et me suis simplement rendu compte que le 2*PI n'avait rien à faire dans mon calcul.
(w/wc) = (2.PI.f)/(2.PI.fc) = (f/fc)
Comme seul nous importe le rapport de fréquences pour le filtre de Butterworth, les 2.PI permettant de passer des fréquences angulaires aux fréquences classiques s'annulent.

Cela m'a permit de vérifier d'ailleurs qu'aucun élément n'est mis à zéro. Merci pour l'attention que tu m'as accordé, Aleph69.

Déjà quelle est la frequence d'échantillonage de ton signal ?

500Hz.

Pourquoi veut tu passer par une FFT juste pour faire un filtre ?

Plus la dimension des données est importante, et plus il est rentable de travailler dans le domaine fréquentiel. Enfin que je sache. Vu que mon nombre max de données est environ 58 000 000 (2kHz x 8h d'acquisition), je pense que le calcul est plus rentable ainsi.
Mais j'suis pas un crac des filtres non plus, j'suis preneur de tout conseil!!

Et si tu veux vraiment passer par une FFT, utilise plutôt une fonction un peu plus smooth que le carré (genre avec une partie polynomiale en pente douce) ca sera plus exact surtout si tu fais du recoupement.

Et voilà, j'suis encore dépassé côté théorie.
J'utilise le carré de la fonction de Gain car cela correspond à l'application du filtre dans le sens directe, puis dans le sens rétrograde. En procédant ainsi, la phase s'annule, tout simplement. Et le gain finale correspond au produit des gains des filtres, donc dans le cas présent au carré de la fonction de gain du filtre de Butterworth.

Pour designer un filtre ce logiciel est très bien : http://www.tobybear.de/p_filterexp.html

C'est bien la première fois qu'on me propose un logiciel payant pour remplacer quelque chose je code moi-même...

[Invité]

OK j'ai lu le code et apparemment tu veux implémenter un filtre FIR assez long. Donc la FFT ça se tient, j'ai rien dit.

[Cyborg]

Y'a pas de mal !!! J'suis pas à l'abri d'une erreur non plus.

Note : Le filtre de Butterworth est un filtre à Réponse Impulsionnelle Infinie.