//lifting_Haar.sce // ce programme implemente une transformation en ondelettes par le /// lifting scheme 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; /// -- aparte -------------------------------------------------------------- /// autre facon de l'ecrire "in-place" (plus compacte, moins matriciellement) ' d=xo-xe; s=xo+d/2; s=2/sqrt(2)*s; d=-1/sqrt(2)*d; //--------------- 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 ----------------------------------------------------------- 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) ------------------------------ //----------------------------------------------------------------- dgam=gam; // initialisation pour la reconstruction in place slam=lam; 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 /// -- aparte /// autre facon de l'ecrire "in-place" (moins matriciellement) ' rrx=zeros(1,length(Resu)); d=dgam; s=slam; d=-sqrt(2)*d; s=sqrt(2)*s; rrxo=s-d/2; rrxe=rrxo-d; 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 laparte 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 dans les calculs.
Un exemple de coefficients déchelle et d'ondelettes obtenus pour 3 niveaux de résolution est présenté sur la Figure 5.