//lifting_Haar_int.sce
// ce programme implemente une transformation en ondelettes par le
/// lifting sceme de Sweldens
// 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
sig=round(sig(1:N));/////////////////// ATTENTION on arrondit au depart pour travailler sur des entiers
sigold=sig;
//--------------------------------------------------------------------------------
//- Maintenant le lifting d'ondelettes
//-
//-
//--------------------------------------------------------------------------------
// on va prendre l'ondelette de Haar ... 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 resolution
K=1/sqrt(2);
Mnor=[K 0; 0 -1/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; -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];
PX=Pred*X;
Updax=Upda(1,1)*PX(1,:)+round(Upda(1,2)*PX(2,:));
Upday=Upda(2,1)*PX(1,:)+round(Upda(2,2)*PX(2,:));
R4x=P4(1,1)*Updax+round(P4(1,2)*Upday);
R4y=P4(2,1)*Updax+round(P4(2,2)*Upday);
R3x=P3(1,1)*R4x+round(P3(1,2)*R4y);
R3y=P3(2,1)*R4x+round(P3(2,2)*R4y);
R2x=P2(1,1)*R3x+round(P2(1,2)*R3y);
R2y=round(P2(2,1)*R3x)+round(P2(2,2)*R3y);
R1x=P1(1,1)*R2x+round(P1(1,2)*R2y);
R1y=P1(2,1)*R2x+round(P1(2,2)*R2y);
Resu=[R1x ; R1y];
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));
R1xi=P1i(1,1)*Resu(1,:)+round(P1i(1,2)*Resu(2,:));
R1yi=P1i(2,1)*Resu(1,:)+round(P1i(2,2)*Resu(2,:));
R2xi=P2i(1,1)*R1xi+round(P2i(1,2)*R1yi);
R2yi=round(P2i(2,1)*R1xi)+round(P2i(2,2)*R1yi);
R3xi=P3i(1,1)*R2xi+round(P3i(1,2)*R2yi);
R3yi=P3i(2,1)*R2xi+round(P3i(2,2)*R2yi);
R4xi=P4i(1,1)*R3xi+round(P4i(1,2)*R3yi);
R4yi=P4i(2,1)*R3xi+round(P4i(2,2)*R3yi);
MR=[R4xi ; R4yi];
Updaxi=(Updai(1,1)*MR(1,:))+round(Updai(1,2)*MR(2,:));
Updayi=(Updai(2,1)*MR(1,:))+round(Updai(2,2)*MR(2,:));
Orig=Predi*[Updaxi ; Updayi];
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 a alors obtenu les résultats des Figures 6 et 7. On a pas les infimes d'erreurs de reconstruction précédentes grâce aux arrondis.
On peut remarquer dans le programme précédent que l'on a perdu en partie le côté pratique du formalisme matriciel (à cause des opérations d'arrondis). À noter aussi que l'inconvénient majeur de conserver des entiers est d'avoir augmenté la dynamique du signal de sortie d'un facteur 2 (cf. Figure 7).
