Implémentation sous Scilab

//lifting_Haar.sce
// ce programme implemente une transformation en ondelettes par le
/// lifting sceme de Sweldens
/// (a partir de lifting_int3.sce)
// 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);





mtlb_load('donnees_rlp_sci.mat'); // un signal ECG sert de signal test
// d1,d2,d3,VR,VL,VF,V3,V4,V5,V6
sig=d2;

N=4096;

sig=round(sig(1:N)); // on arrondit pour travailler avec des entiers
sigold=sig;




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



// on va prendre l'ondelette de Haar ... on utilise la decomposition polyphase (Valens, Sweldens)
// 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
// fnor est le facteur de normalisation
// xe est le vecteur de composantes paires
// xo est le vecteur de composantes impaires
// X est le vecteur xe xo



Wsig=[]; // coefficients d'ondelettes

reso=2; // nombre de resolution


fnor=1/sqrt(2);
fnorinv=sqrt(2);


Mnor=[2 0; 0 -1];
Mnorinv=inv(Mnor);

Pred=[1 0; -1 1];
Upda=[1 1/2; 0 1];

Predi=[1 0; 1 1];
Updai=[1 -1/2; 0 1];



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

for iter=1:1:reso
  xe=sig(2:2:$);
  xo=sig(1:2:$);
  X=[xe ; xo];

  Resu=fnor*Mnor*Upda*Pred*X; 

  gam=Resu(2,:);
  lam=Resu(1,:);

  //--------------- 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([[1:length(gam)]'] ,[gam'],[2],'181',"coefficients d ondelettes");
  xtitle("coefficients d ondelettes "," temps "," amplitude ");

  Wsig=[gam Wsig]; // stockage des coefficients d'ondelettes
  sig=lam;
  
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 ondelettes et d echelle","temps","tension");
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) ------------------------------
//-----------------------------------------------------------------


 for iter=1:1:reso
   rx=zeros(1,length(Resu));   
   Orig=fnorinv*Predi*Updai*(Mnorinv)*Resu; 

   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

   Resu=[lam ; gam];
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 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 3. L'erreur de reconstruction de la Figure 4 est due à la présence du rationnel $\sqrt{2}$ dans les calculs.

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

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

Figure 5: 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 5.