Implémentation sous Scilab

// lifting_D4.sce
// ce programme implemente une transformation en ondelettes par le
/// lifting scheme de Sweldens pour les ondelettes Daubechies(4) (D4)
/// 
// Globalement par rapport a la transformee "classique" on
// procede d'abord au sous echantillonnage avant d'appliquer le filtre
// on separe coefficient pairs et impair (Split) 


clear
xdel(winsid());
stacksize(60000000);


/// Le signal a tester est un signal ECG de frequence
//d'echantillonnage de 500Hz

mtlb_load('donnees_rlp_sci.mat');   // d1,d2,d3,VR,VL,VF,V3,V4,V5,V6
sig=d2; // on choisit la derivation d2

N=4096; // Nb de points a considerer

//////-----------------------------------------------------------------------------------------

   
     
  h0=(1+sqrt(3))/(4*sqrt(2));
  h1=(3+sqrt(3))/(4*sqrt(2));
  h2=(3-sqrt(3))/(4*sqrt(2));
  h3=(1-sqrt(3))/(4*sqrt(2));
     

//////------------------------------------------------------------------------------------------


sig=round(sig(1:N));/////////////////// ATTENTION on arrondit au depart pour travailler sur des entiers
sigold=sig;



     
//--------------------------------------------------------------------------------
//-  Notations et donnees
//-
//- 
//--------------------------------------------------------------------------------

// on va prendre l'ondelette de Daubechies D4 ... on utilise la decomposition polyphase (Valens)
// P est la matrice polyphase
// gam designe les coefficients d'ondelettes
// lam designe les coefficients d'echelle
// Pred est la matrice permettant de faire la prediction
// Upda est la mtrice permettant de faire la mise a jour 
// Mnor est la matrice de normalisation
//
// xe est le vecteur de composantes paires
// xo est le vecteur de composantes impaires
// X est le vecteur xe xo


Wsig=[];

reso=3; // nombre de resolutions

K=(sqrt(3)+1)/sqrt(2);

Mnor=[1/K 0; 0 K];  

P1=[1 K-K^2 ; 0 1]; 
P2=[1 0 ; -1/K 1 ];
P3=[1 K-1 ; 0 1 ];
P4=[1 0 ; 1 1];
 


 

P1i=[1  (-K+K^2) ; 0 1]; 
P2i=[1 0;  1/K 1 ];
P3i=[1  -K+1; 0 1 ];
P4i=[1  0; -1 1];
 


Pred=[1 0; sqrt(3) 1];
Upda=[1 -sqrt(3)/4; 0 1];    /// Attention Updazm(1,2) fait intervenir echantillon precedent
Updazm=[0 -((sqrt(3)-2)/4); 0 0]; 
Pred2=[1 0; 0 1];                               /// Attention Pred2zp(2,1) fait intervenir echantillon suivant
Pred2zp=[0 0; -1 0]; 


Predi=[1 0; -sqrt(3) 1];
Updai=[1 sqrt(3)/4  ; 0 1];   /// Attention Updaizm(1,2) fait intervenir echantillon precedent
Updaizm=[0 ((sqrt(3)-2)/4)  ; 0 0]; 
Pred2i=[1 0; 0 1];                                /// Attention Predi2zp(2,1) fait intervenir echantillon suivant
Pred2izp=[0 0; 1 0]; 

//-----------------------------------------------------------
//--- transformation duale (\tilde{P}) ----------------------
//-----------------------------------------------------------


for iter=1:1:reso

  xe=sig(2:2:$);
  xo=sig(1:2:$);

  X=[xe ; xo];


   Px=Pred(1,1)*X(1,:)+(Pred(1,2)*X(2,:));    // d   (details) passe-haut= coeff d'ondelettes 
   Py=Pred(2,1)*X(1,:)+(Pred(2,2)*X(2,:));  // s  (smooth)= passe-bas= coeff d'echelle 
  
   Updax=Upda(1,1)*Px+(Upda(1,2)*Py+Updazm(1,2)*[0 Py(1:1:$-1) ]);  // d
   Upday=Upda(2,1)*Px+(Upda(2,2)*Py);   //s

   Pxx=Pred2(1,1)*Updax+(Pred2(1,2)*Upday);   //d 
   Pyy=(Pred2(2,1)*Updax+Pred2zp(2,1)*[Updax(2:1:$) 0])+(Pred2(2,2)*Upday);  //s
  
   R4x=P4(1,1)*Pxx+(P4(1,2)*Pyy); //d 
   R4y=P4(2,1)*Pxx+(P4(2,2)*Pyy); //s
  
   R3x=P3(1,1)*R4x+(P3(1,2)*R4y); //d 
   R3y=P3(2,1)*R4x+(P3(2,2)*R4y); //s

   R2x=P2(1,1)*R3x+(P2(1,2)*R3y); //d 
   R2y=(P2(2,1)*R3x)+(P2(2,2)*R3y); //s

   R1x=P1(1,1)*R2x+(P1(1,2)*R2y); //d 
   R1y=P1(2,1)*R2x+(P1(2,2)*R2y); //s


  Resu=[R1x ; R1y];

  gam=Resu(1,:); //d
  lam=Resu(2,:); //s
 

  //--------------- Affichage ------------------------------------------

  xset('window',max(winsid()+1));
  subplot(211),plot2d([[1:length(gam)]',  [1:length(lam)]'] ,[gam', lam'],[2, 3],'181',"gama (ondelette) @lambda (echelle)");
  xtitle("coefficients d ondelettes et d echelle "," temps "," amplitude ");
  subplot(212),plot2d([[2:length(gam)]'] ,[gam(2:1:$)'],[2],'181',"coefficients d ondelettes");
  xtitle("coefficients d ondelettes "," temps "," amplitude ");
  Wsig=[gam Wsig]; // stockage des coefficients d'ondelettes'
  sig=lam;
  
  
  ///-- aparte ---------------------------------------------
  /// autre facon de l'ecrire "in-place" (plus compacte, moins matriciellement) ' 
  
  s=xo+sqrt(3)*xe;
  d=xe-sqrt(3)/4*s-(sqrt(3)-2)/4*[0 s(1:1:$-1)];
  s=s-[d(2:1:$) 0];
  d=(sqrt(3)+1)/sqrt(2)*d;
  s=(sqrt(3)-1)/sqrt(2)*s;

  //--------------- Affichage ------------------------------------------

  //xset("window",max(winsid()+1)),
  //subplot(211),plot2d([[1:length(s)]'],[s'],[3],'081')
  //subplot(212),plot2d([ [1:length(d)]'],[ d'],[2],'081')

   ///-- fin de l' aparte ---------------------------------------------'
 
  /// comparaison avec les  banc de filtres 
  /// wt paquetage --- uvi_wave ---
  ///   wcoef=wt(round(sigold),dH,dG,iter,0,0);  
  ///   xset("window",max(winsid()+1)),
  ///  plot2d([1:length(wcoef)]',wcoef,[2],'081')

end  // fin de l'iteration duale 
 


 WSsig=[lam Wsig];      //stockage des coefficients d'echelle et d'ondelettes

//--------------- Affichage ------------------------------------------

 xset('window',max(winsid()+1));
 plot2d([[1:length(WSsig)]'] ,[WSsig'],[3],'181',"coefficient d ondelettes et d echelle");
 xtitle("coefficient d echelle et d ondelettes", string(reso)+ " resolutions","amplitude");

 xset("window",max(winsid()+1)),
 plot2d([[1:length(lam)]' [1:length(gam)]'],[lam', gam'],[3, 2],'081')


//-----------------------------------------------------------------
//------ transformation primaire (P) ------------------------------
//-----------------------------------------------------------------
  
 
dgam=gam;  // initialisation pour le in-place
slam=lam;

 for iter=1:1:reso
   
   
  rx=zeros(1,length(Resu));   

  R1xi=P1i(1,1)*Resu(1,:)+(P1i(1,2)*Resu(2,:));
  R1yi=P1i(2,1)*Resu(1,:)+(P1i(2,2)*Resu(2,:));

  R2xi=P2i(1,1)*R1xi+(P2i(1,2)*R1yi);
  R2yi=P2i(2,1)*R1xi+(P2i(2,2)*R1yi);

  R3xi=P3i(1,1)*R2xi+(P3i(1,2)*R2yi);
  R3yi=P3i(2,1)*R2xi+(P3i(2,2)*R2yi);

  R4xi=P4i(1,1)*R3xi+(P4i(1,2)*R3yi);
  R4yi=P4i(2,1)*R3xi+(P4i(2,2)*R3yi);

  Pxxi=Pred2i(1,1)*R4xi+(Pred2i(1,2)*R4yi);   
  Pyyi=(Pred2i(2,1)*R4xi+Pred2izp(2,1)*[R4xi(2:1:$) 0])+(Pred2i(2,2)*R4yi);  
  
  Updaxi=Updai(1,1)*Pxxi+(Updai(1,2)*Pyyi+Updaizm(1,2)*[0 Pyyi(1:1:$-1)]);  
  Updayi=Updai(2,1)*Pxxi+(Updai(2,2)*Pyyi);   
   
  Pxi=Predi(1,1)*Updaxi+(Predi(1,2)*Updayi);    
  Pyi=Predi(2,1)*Updaxi+(Predi(2,2)*Updayi); 
  
  Orig=[Pxi ; Pyi]; 

  rxe=Orig(1,:);
  rxo=Orig(2,:);

  rx(1:2:length(Resu))=(rxo);
  rx(2:2:length(Resu))=(rxe);

  if iter<reso
     gam=WSsig(length(rx)+1:2*length(rx));
     lam=rx;
  end

   ///  -- aparte --------------------------------------------------------------
   /// autre facon de l'ecrire "in-place" (moins matriciellement) ' 
     rrx=zeros(1,length(Resu));   
     d=dgam;
     s=slam; 
     
     s=(sqrt(3)+1)/sqrt(2)*s;
     d=(sqrt(3)-1)/sqrt(2)*d;
     s=s+[d(2:1:$) 0];
     rrxe=d+sqrt(3)/4*s+(sqrt(3)-2)/4*[0 s(1:1:$-1)];
     rrxo=s-sqrt(3)*rrxe;
     rrx(1:2:length(Resu))=(rrxo);
     rrx(2:2:length(Resu))=(rrxe);
      
     xset('window',max(winsid()+1));
     plot2d([[1:length(rrx)]'] ,[rrx'],[2],'181',"signal reconstruit in place");
     xtitle("signal reconstruit in place ","temps","tension");
   
     if iter<reso
      dgam=WSsig(length(rrx)+1:2*length(rrx));
      slam=rrx;
     end
   
   ///  -- fin de l aparte --------------------------------------------------------------


  Resu=[gam ; lam];
  
  
  

 end

//--------------- Affichage ------------------------------------------

xset('window',max(winsid()+1));
subplot(211),plot2d([[1:length(rx)]',  [1:length(sigold)]'] ,[rx', sigold'],[2, 3],'181',"signal reconstruit@signal original");
xtitle("signal reconstruit et signal original","temps","tension");
subplot(212),plot2d([[1:length(sigold)]'] ,[(rx-sigold)'],[2],'181',"erreur entre les courbes");
xtitle("difference entre le signal reconstruit et le signal original","temps","ecart");

// fin du programme

On présente encore quelques illustrations du programme : le signal analysé est un signal ECG (fréquence d'échantillonnage 500 Hertz). Le résultat de l'analyse (i.e. coefficients d'ondelettes et coefficients d'échelle) sont présentés sur la Figure 9. Les infimes erreurs de reconstruction de la Figure 10 est due à la présence des rationnels $\sqrt{3}$ et $\sqrt{2}$ dans les calculs.

Figure 9: Coefficients d'ondelettes et d'échelle de la résolution 1

Figure 10: Signal original, signal reconstruit et erreur (différence entre les deux signaux

Figure 11: Coefficients d'échelle et d'ondelettes sur 3 niveaux de résolutions

Un exemple de coefficients déchelle et d'ondelettes obtenus pour 3 niveaux de résolution est présenté sur la Figure 11.

(On peut également, comme pour l'ondelette de Haar, effectuer avec l'ondelette de Daubechies la transformation d'entiers en coefficients d'ondelettes entier.)

Voilà, That's all Folks !!!