Composition de Filtre Digitaux en compile-time [RESOLU]

Composition de Filtre Digitaux en compile-time [RESOLU] - C++ - Programmation

Marsh Posté le 19-05-2004 à 17:48:00    

Voila mon problème :
 
Je devellope une série d'application de calcul assez gourmande en temps.
Parmi ceux ci, il existe une quantité astronomique de Filtre à Réponse  
Impulsionelle Fine (FRIF ou FIRF en english). Pour optimiser ces filtrage,  
dont les coeffcieint sont fixes, j'ai décidé de mettre en place un mécanisme  
basé sur les templates pour accélérer la chose.
 
Le principe est le suivant :
 
Soit un filtre F 1x3 de coeff [ -1,0,1 ] (dérivée horizontale), je
declare une instance d'une classe filter :
 

Code :
  1. filter<1,3,-1,0,1> grad_x;


 
Ensuite, j'applique ce filtre à une image de la maniére suivante :
 

Code :
  1. vector<float> image,filtrage;
  2. filtrage = grad_x(image);


 
Rien de bien fulgurant, l'operateur () de filter<W,H,A1,A2,A3> effectue le filtrage.
 
Maintenant j'aimerais bien ecrire :
 

Code :
  1. filter< 1,3, -1,0,1 > grad_x;
  2. filter< 1,3,  1,2,1 > gaussian_x;
  3. vector<float> image,filtrage;
  5. filtrage = (grad_x*gaussian_x)(image);


 
Qui appliquerait successivement le filtre gaussian_x puis grad_x à image.
l'idée est donc de proposer un operateur * sur les filter pour calculer à la compilation
les nouveaux paramétres du filtre composé.
 
/!\ FORMULES /!\
 
Soit un filtre F de taille N [ F0 ... FN-1], un filtre G de taille M [ G0 .. GM-1 ]
La composée de F et de G , noté : F o G est un filtre H tel que :
 
Taille de H : P = M+N-1
 
Et de Coefficient H0 ... HP tel que :
 
Hi = Somme de j=0 à i-1 de F(j)*G(i-j)*delta(i,j)
 
aevc delta(i,j) = 1 si (i-j) < M et j<N, 0 sinon
 
Concretement, le coefficient est un produit de polynomes, minus un certains nombre  
de coefficient non existent. (cf Transformée en Z d'un filtre pour retrouver la formule)
 
/!\ FIN DES FORMULES /!\
 
Voila une partie de la solution de la composition de filtre.  
 

Code :
  1. #include <iostream>
  2. #include <vector>
  4. // initialiseur par operator,
  5. template<class U,int N> class ct_list_initializer
  6. {
  7.     public:
  8.     ct_list_initializer( U* data ) : data_(data) {}
  9.    
  10.     ct_list_initializer<U,N-1> operator,( const U& value )
  11.     {
  12.         *data_ = value;
  13.         return ct_list_initializer<U,N-1>(data_+1);
  14.     }
  16.     private:
  17.     U* data_;
  18. };
  20. // Identifiant d'orientation des filtres
  21. struct vertical {};
  22. struct horizontal {};
  24. // Classe de filtre
  25. template<typename ORIENTATION, typename TYPE, int SIZE>
  26. class filter
  27. {
  28.     public:
  29.    
  30.     typedef filter<ORIENTATION, TYPE, SIZE> self_t;
  31.     typedef TYPE value_type;
  32.    
  33.     filter() {}
  34.    
  35.     TYPE  operator[](int i) const { return coeff[i]; }
  36.     TYPE& operator[](int i)       { return coeff[i]; }
  38.     const TYPE*   begin() const { return &(coeff[0]); }
  39.     TYPE*   begin()  { return &(coeff[0]); }
  40.    
  41.    ct_list_initializer<TYPE,SIZE-2> operator=( const TYPE& val )
  42.     {
  43.         coeff[0] = val;
  44.         return ct_list_initializer<TYPE,SIZE-2>(coeff+1);
  45.     }
  46.    
  47.     private:
  48.     TYPE coeff[SIZE];
  49. };
  51. // Delta de Kroenecker pour (I,J)
  52. template<int I,int J,int M,int N> struct delta
  53. {
  54.     enum { k = ( ( (I-J) < N) && (J-1 < M) ? 1 : 0 ) };
  55. };
  57. // Boucle Interne (sur J)
  58. template<class T,int M, int N,int I,int J> struct inner_loop
  59. {
  60.     enum { go = (J-1) != 0 };
  61.     static void apply( T* result, const T* f, const T* g )
  62.     {
  63.         result[I-1] += f[J-1]*g[I-J]*delta<I,J,M,N>::k;
  64.         inner_loop<T,M,N,I,(go ? (J-1) : 0) >::apply(result,f,g);
  65.     }
  66. };
  68. template<class T,int M, int N,int I> struct inner_loop<T,M,N,I,0>
  69. {
  70.     static void apply( T* result, const T* f, const T* g ) {}
  71. };
  73. // Boucle Externe (sur I)
  74. template<class T,int M, int N,int I> struct outer_loop
  75. {
  76.     enum { go = (I-1) != 0 };
  77.     static void apply( T* result, const T* f, const T* g )
  78.     {
  79.         result[I-1] = T(0); 
  80.         inner_loop<T,M,N,I,I>::apply(result,f,g);   
  81.         outer_loop<T,M,N, (go ? (I-1) : 0) >::apply(result,f,g);
  82.     }
  83. };
  85. template<class T,int M, int N> struct outer_loop<T,M,N,0>
  86. {
  87.     static void apply( T* result, const T* f, const T* g ) {}
  88. };
  90. // operateur * de la classe filter
  91. // Composition de deux filtres
  92. template<class T,class O, int M, int N>
  93. filter<O,T,M+N-1> operator*( const filter<O,T,M>& f, const filter<O,T,N>& g )
  94. {
  95.     filter<O,T,M+N-1> result;
  96.     outer_loop<T,M,N,M+N-1>::apply(result.begin(),f.begin(),g.begin());
  97.     return result;
  98. }
  100. // Test
  101. int main (int argc, char * const argv[])
  102. {           
  103.     filter<horizontal,double,4> f;
  104.     f = 1,2,2,1;
  106.     filter<horizontal,double,3> g;
  107.     g = -1,0,1;
  108.    
  109.     filter<horizontal,double,6> h;
  110.    
  111.     h = g*f;
  112.     for(size_t i=0;i<6;i++) std::cout << h[i] << " ";
  113.     std::cout << std::endl;
  114.    
  115.     return 0;
  116. }


 
Que reste il à faire :
 
-> gérer l'application du filtre à un dataset.
-> gérer les vertical*horizontal.
 
Question : qui connait une formule pour décomposer des filtres séparables ?


Message édité par Joel F le 28-05-2004 à 11:02:52
Reply

Marsh Posté le 19-05-2004 à 17:48:00   

Reply

Marsh Posté le 19-05-2004 à 18:23:25    

Je ne vais absolumment pas pouvoir t'aider étant donné que je ne suis pas assez bon pour ça, mais ça m'intéresse beaucoup car j'ai utilisé ce genre de filtre, et que ça m'intéresserait beaucoup de pouvoir les calculer à la volée, car moi bien sûr je les calculais un par un, puis je faisais le produit de polynôme tout ça avant avec Matlab ..

Reply

Marsh Posté le 19-05-2004 à 18:37:56    

pourquoi t'utilise pas un tableau ? en optimisant bien, il aura disparu au final
 
je comprends pas ou tu les donnes tes A et B .. ton operator* il est louche ... il est binaire et en fonction membre ...  
 
question con : tu connais les typedef ?

Reply

Marsh Posté le 19-05-2004 à 18:42:57    

Taz a écrit :

pourquoi t'utilise pas un tableau ? en optimisant bien, il aura disparu au final


Non, car les coeff ne doivent pas stationner en mémoire , y a du SIMD derriere tout ca et la moindre lecture mémoire me tue les perfs.
 

Taz a écrit :


je comprends pas ou tu les donnes tes A et B .. ton operator* il est louche ... il est binaire et en fonction membre ...  


A0,B0 ... sotn les parametres template des filtres f1 et f2.
et filter<M+N-1,....> est le type de retour.
 

Taz a écrit :


question con : tu connais les typedef ?


 
oui, cela va t il m'aidé ?

Reply

Marsh Posté le 19-05-2004 à 19:03:15    

1) tu devrais faire un test avec un tableau, tu sera surpris du travail du compilateur. je t'assure, c'est très impressionnant : ne ferme pas la porte. bien sur il faut mettre en place un parcours du tableau, ce qui peut être difficile ...
2) le typedef à quoi il sert ? ben à pas écrire
filter<M+N-1 ,  
                A0*B0*d<0,0,M,N>::k,
                A1*B0*d<1,0,M,N>::k + A0*B1*d<0,1,M,N>::k,
                A2*B0*d<2,0,M,N>::k + A1*B1*d<1,1,M,N>::k + A0*B2*d<0,2,M,N>::k,
                A3*B0*d<3,0,M,N>::k + A2*B1*d<2,1,M,N>::k + A1*B2*d<1,2,M,N>::k + A0*B3*d<0,3,M,N>::k,
                A4*B0*d<4,0,M,N>::k + A3*B1*d<3,1,M,N>::k + A2*B2*d<2,2,M,N>::k + A1*B3*d<1,3,M,N>::k + A0*B4*d<0,4,M,N>::k,
                A5*B0*d<5,0,M,N>::k + A4*B1*d<4,1,M,N>::k + A3*B2*d<3,2,M,N>::k + A2*B3*d<2,3,M,N>::k + A1*B4*d<1,4,M,N>::k + A0*B5*d<0,5,M,N>::k,
                // ... etc ...  
 
 
plus d'une fois :D
 
sinon je vois pas trop quoi faire ...tu as déjà déroulé à la main ... peut être avec des classes intermédiaires tu peux laisser le compilateur s'en charger ?  
 
 
sinon tu vas en écrire combien de fois des choses comme ça ? si ça marche, ça vaut vraiment le coup de simplifier ça ?

Reply

Marsh Posté le 19-05-2004 à 19:09:34    

les filtres, j'en ai plus de 80 dans des boucles critiques.  
En scalaire, je gagne déja +110% en vitesse, et en AltiVec, j'evite une PERTE de 95% (en bref je reste aux alentours d'une accél de x8 au lieu de chutait à x1.2) ...
 
AltiVec et la mémoire c'est le mal je t'assure.
 
Mais :
 
1/ Comment ferais tu avec un tabelau, j'ai peut etre oublié une astcue.
 
2/ Ou placer le typedef ??
 
3/ Classes intermediaires : je pourrais ecrire un template recursif qui calcule les Ci récursivement ?


Message édité par Joel F le 19-05-2004 à 19:12:38
Reply

Marsh Posté le 19-05-2004 à 19:36:50    

oui le calculer récursivement, c'est bien. et si tu as un tableau + un paramètre template entier de comptage, tu te balades sur ton tableau et ça roule bien avec ça :D

Reply

Marsh Posté le 19-05-2004 à 19:40:28    

j'ai du mal à voir l'astuche avec le tableau o_O

Reply

Marsh Posté le 19-05-2004 à 19:53:37    

Joel F a écrit :


C'est lourd, c'est moche et inchiable à ecrire ....
 
Y aurait il un moyen de masquer cette liste infâme de parametres ou de simplifier la méthode.


 
un générateur de code?

Reply

Marsh Posté le 19-05-2004 à 20:01:51    

enfait pour les tableaux, j'ai parlé peut être un peu vite ..; j'en était resté sur l'exemple de , ou il allait me chercher valeur comme il faut, mais là sur un tableau, il est incapable
 

Code :
  1. template<typename T, size_t N>
  2. struct StaticAdder
  3. {
  4.   static const T & sum(const T *x, const T &init = T(0))
  5.   {
  6.     return StaticAdder<T, N-1>::sum(x + 1, init + *x);
  7.   }
  8. };
  11. template<typename T>
  12. struct StaticAdder<T, 0>
  13. {
  14.   static const T & sum(const T *, const T &s)
  15.   {
  16.     return s;
  17.   }
  18. };
  21. template<typename T, size_t N>
  22. T static_sum(const T (&x)[N])
  23. {
  24.   BOOST_STATIC_ASSERT(N > 0);
  25.   return StaticAdder<T, N>::sum(x);
  26. }

me génère N-1 addition inlinée :/

Reply

Marsh Posté le 19-05-2004 à 20:01:51   

Reply

Marsh Posté le 19-05-2004 à 20:02:53    

Joel F a écrit :

j'ai du mal à voir l'astuche avec le tableau o_O

bah je pensais qu'il serait capable de travailler correctement : tu as un tableau initialisé littéralement (statiquement), je pensais qu'il aurait était capable d'aller chercher les valeurs ... mais non

Reply

Marsh Posté le 27-05-2004 à 15:26:52    

Up !
 
j'ai reflechi à la chose :
 
La classe filtre :

Code :
  1. template<typename ORIENTATION, typename TYPE, int SIZE>
  2. class filter
  3. {
  4.     public:
  5.    
  6.     typedef filter<ORIENTATION, TYPE, SIZE> self_t;
  7.     typedef TYPE value_type;
  8.    
  9.     filter() {}
  10.    
  11.     TYPE  operator[](int i) const { return coeff[i]; }
  12.     TYPE& operator[](int i)       { return coeff[i]; }
  14.     private:
  15.     TYPE coeff[SIZE];
  16. };


 
ORIENTATION est une classe qui définit le sens (vertical ou horizontal du filtre).
TYPE est le type de données du filtre (float, int etc ...)
SIZE est la taille du filtre.
 
Le filter se rempli via [] ou via , par une surcharge non présente ici.
Je définis * pour les filter de meme orientation par :
 

Code :
  1. template<class T,class O, int M, int N>
  2. filter<O,T,M+N-1> operator*( const filter<O,T,M>& f, const filter<O,T,N>& g )
  3. {
  4.     filter<O,T,M+N-1> result;
  5.    
  6.     for(size_t i=0;i<M+N-1;i++)
  7.     {
  8.         result[i] = 0;
  9.         for(size_t j=0;j<i;j++)
  10.         {
  11.             result[i] += f[j]*g[i-j]*( ((i-j)<N) && (j<M) && (i<N)  ? 1 : 0 );
  12.         }
  13.     }
  14.     return result;
  15. }


 
ensuite operator() fait ma tambouille SIMD.
 
Bref il me ne me reste qu'a metaprogrammer le produit des tableaux :p

Reply

Marsh Posté le 27-05-2004 à 15:47:41    

bah voilà :D c'est bien plus propre comme ça

Reply

Marsh Posté le 27-05-2004 à 15:55:11    

Taz a écrit :

bah voilà :D c'est bien plus propre comme ça


 
N'est-ce pas ^^
Bon je pense que je peut prendre exemple sur ton StaticAdder ?

Reply

Marsh Posté le 27-05-2004 à 16:00:59    

ben comme je t'ai dit, je suis un peu déçu, parce que le résultat est disponible à la compilation, mais au lieu d'y être calculé, j'ai juste déroullée la boucle ...

Reply

Marsh Posté le 27-05-2004 à 16:55:47    

c'est déja ca. Je vais benchmarker le tout.

Reply

Marsh Posté le 27-05-2004 à 17:18:46    

ok, tu me tiens au courant s'il te plaît, ce que tu fais es toujours très intéressant. si tu déroules par bloc, y a peut être moyen de mieux placer tes instructions SIMD, non ?

Reply

Marsh Posté le 27-05-2004 à 17:24:38    

Taz a écrit :

ok, tu me tiens au courant s'il te plaît, ce que tu fais es toujours très intéressant. si tu déroules par bloc, y a peut être moyen de mieux placer tes instructions SIMD, non ?


 
Je te montre ce qui se passe avec AltiVec.
 
On va se placer dans le cas simple du filtre 1x3 horizontal sur une image de flottants.
 
Je génére les vecteurs C1,C2,C3 qui contiennent le coefficient du filtre :
C1 = [ 1 1 1 1 ]
C2 = [ 2 2 2 2 ]
C3 = [ 1 1 1 1 ]
 
Ensuite pour chaque vecteur V et W contigues (ie bloc de 4 éléments), j'effectue l'operation suivante ,je note V = [ a b c d ] et W = [ e f g h ]:
 
V1 = (V,W) >> 1  [ b c d e ]
V2 = (V,W) >> 2  [ c d e f ]
 
T = C1*V+C2*V1;
R = C3*V2+T;
 
J'ecris ensuite R dans le tableau résultat.
 
L'astuce ici est de précalculer els C1 ... Cn pour ne faire qu'une seule boucle. On peut aussi dérouler le calcul de V1,V et de T pour traiter plus de vecteurs simultanement.

Reply

Marsh Posté le 27-05-2004 à 17:47:04    

Si tu veut je te file mon rapport de DEA si tu veut jeter un oeil :p

Reply

Marsh Posté le 27-05-2004 à 17:47:43    

ça serait cool ça (sauf la partie math :D)

Reply

Marsh Posté le 28-05-2004 à 10:22:07    

Bon voila une partie de la solution de la composition de filtre.  
 

Code :
  1. #include <iostream>
  2. #include <vector>
  4. // initialiseur par operator,
  5. template<class U,int N> class ct_list_initializer
  6. {
  7.     public:
  8.     ct_list_initializer( U* data ) : data_(data) {}
  9.    
  10.     ct_list_initializer<U,N-1> operator,( const U& value )
  11.     {
  12.         *data_ = value;
  13.         return ct_list_initializer<U,N-1>(data_+1);
  14.     }
  16.     private:
  17.     U* data_;
  18. };
  20. // Identifiant d'orientation des filtres
  21. struct vertical {};
  22. struct horizontal {};
  24. // Classe de filtre
  25. template<typename ORIENTATION, typename TYPE, int SIZE>
  26. class filter
  27. {
  28.     public:
  29.    
  30.     typedef filter<ORIENTATION, TYPE, SIZE> self_t;
  31.     typedef TYPE value_type;
  32.    
  33.     filter() {}
  34.    
  35.     TYPE  operator[](int i) const { return coeff[i]; }
  36.     TYPE& operator[](int i)       { return coeff[i]; }
  38.     const TYPE*   begin() const { return &(coeff[0]); }
  39.     TYPE*   begin()  { return &(coeff[0]); }
  40.    
  41.    ct_list_initializer<TYPE,SIZE-2> operator=( const TYPE& val )
  42.     {
  43.         coeff[0] = val;
  44.         return ct_list_initializer<TYPE,SIZE-2>(coeff+1);
  45.     }
  46.    
  47.     private:
  48.     TYPE coeff[SIZE];
  49. };
  51. // Delta de Kroenecker pour (I,J)
  52. template<int I,int J,int M,int N> struct delta
  53. {
  54.     enum { k = ( ( (I-J) < N) && (J-1 < M) ? 1 : 0 ) };
  55. };
  57. // Boucle Interne (sur J)
  58. template<class T,int M, int N,int I,int J> struct inner_loop
  59. {
  60.     enum { go = (J-1) != 0 };
  61.     static void apply( T* result, const T* f, const T* g )
  62.     {
  63.         result[I-1] += f[J-1]*g[I-J]*delta<I,J,M,N>::k;
  64.         inner_loop<T,M,N,I,(go ? (J-1) : 0) >::apply(result,f,g);
  65.     }
  66. };
  68. template<class T,int M, int N,int I> struct inner_loop<T,M,N,I,0>
  69. {
  70.     static void apply( T* result, const T* f, const T* g ) {}
  71. };
  73. // Boucle Externe (sur I)
  74. template<class T,int M, int N,int I> struct outer_loop
  75. {
  76.     enum { go = (I-1) != 0 };
  77.     static void apply( T* result, const T* f, const T* g )
  78.     {
  79.         result[I-1] = T(0); 
  80.         inner_loop<T,M,N,I,I>::apply(result,f,g);   
  81.         outer_loop<T,M,N, (go ? (I-1) : 0) >::apply(result,f,g);
  82.     }
  83. };
  85. template<class T,int M, int N> struct outer_loop<T,M,N,0>
  86. {
  87.     static void apply( T* result, const T* f, const T* g ) {}
  88. };
  90. // operateur * de la classe filter
  91. // Composition de deux filtres
  92. template<class T,class O, int M, int N>
  93. filter<O,T,M+N-1> operator*( const filter<O,T,M>& f, const filter<O,T,N>& g )
  94. {
  95.     filter<O,T,M+N-1> result;
  96.     outer_loop<T,M,N,M+N-1>::apply(result.begin(),f.begin(),g.begin());
  97.     return result;
  98. }
  100. // Test
  101. int main (int argc, char * const argv[])
  102. {           
  103.     filter<horizontal,double,4> f;
  104.     f = 1,2,2,1;
  106.     filter<horizontal,double,3> g;
  107.     g = -1,0,1;
  108.    
  109.     filter<horizontal,double,6> h;
  110.    
  111.     h = g*f;
  112.     for(size_t i=0;i<6;i++) std::cout << h[i] << " ";
  113.     std::cout << std::endl;
  114.    
  115.     return 0;
  116. }


 
Que reste il à faire :
 
-> gérer l'application du filtre à un dataset.
-> gérer les vertical*horizontal.
 
Question : qui connait une formule pour décomposer des filtres séparables ?


Message édité par Joel F le 28-05-2004 à 10:58:08
Reply

Marsh Posté le 28-05-2004 à 10:37:52    

(T)0; -> T(0) :D

Reply

Marsh Posté le 28-05-2004 à 10:46:17    

:jap: ok :p

Reply

Marsh Posté le 28-05-2004 à 10:55:01    

fais gaffe il te manque quelques & par ci par là. niveau const, c'est un poil perfectible.
 
tip of the day : la définition d'une fonction membre dans la déclaration de la classe est implicitement inline (ARM 9.3.2)

Reply

Marsh Posté le 28-05-2004 à 10:55:50    

Taz a écrit :

fais gaffe il te manque quelques & par ci par là. niveau const, c'est un poil perfectible.
 
tip of the day : la définition d'une fonction membre dans la déclaration de la classe est implicitement inline (ARM 9.3.2)


 
Hmmm ou ça les & ??
 
Ok pr le tips :p

Reply

Marsh Posté le 28-05-2004 à 11:23:15    

cai bon là

Reply

Marsh Posté le 31-05-2004 à 17:35:41    

joelF :  
Bonjour,  
 
tu utilises quoi, comme librairie FFT, pour le filtrage ?


Message édité par el muchacho le 31-05-2004 à 17:36:27
Reply

Marsh Posté le 31-05-2004 à 20:21:54    

el muchacho a écrit :

joelF :  
Bonjour,  
 
tu utilises quoi, comme librairie FFT, pour le filtrage ?


 
Rien, tout à la main optimisé AltiVec pour Power PC G5

Reply

Marsh Posté le 31-05-2004 à 20:35:54    

A tout hasard, connais-tu FFTW ?
http://www.fftw.org/speed/g5-2GHz/
http://www.fftw.org/accuracy/g5-2GHz/


Message édité par el muchacho le 31-05-2004 à 20:48:19
Reply

Marsh Posté le 31-05-2004 à 21:44:06    

oui mais je n'utilise pas QUE de la FFT.
Mon but estd e fournir une bibliothéque de haut niveau pour AltiVec pas forcement focaliser sur des problematiques de FFT.

Reply

Marsh Posté le    

Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed