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.