// ce programme implemente une transformation en ondelettes par le
/// lifting sceme de Sweldens
// lifting_Haar_int.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);


//getf('polyfit.sci');
//getf('polyval.sci');

// D'abord on va implementer une transformee de Haar basique



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

N=4096;


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




//--------------------------------------------------------------------------------
//- Maintenant le vrai lifting d'ondelettes qui utilise des entiers
//-
//- 
//--------------------------------------------------------------------------------

// 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