view src/savgol_xl.sci @ 0:0d2d20ef6e72 draft

Uploaded

function[res_savgol]=savgol_xl(saisir,window2,d012,degre_polynom)
    //
    // Algorithme basé sur Savitzsky-Golay (SG), NON modifiée cad projection 
    // sur trois axes non orthogonaux entre eux
    //
    // ------------------------------------------------------------------------
    // Entrées: 
    //    saisir:       matrice au format Saisir ou bien matrice classique 2D
    //    plage:        nombre entier positif impair >=3
    //    d012:         dérivée (0, 1 ou 2)
    //    degre_polynom: degré du polynome; par défaut: 2 
    // Sorties:
    //    res_savgol    spectres après application de d012
    //    filter_d012:  un filtre F tel que: X après dérivée d012= X.F
    // ------------------------------------------------------------------------
    //
    // Avertissement: 
    // SG ne gère pas les demi-plages à droite et à gauche. Celles-ci 
    // sont approximées au mieux par régression aux moindres carrés, mais le 
    // résultat est toujours nettement moins bien que lorsque SG est applicable.
    // En conséquence les extrémités droite et gauche sont perturbées,
    // et la perturbation s'amplifie avec le degré de dérivation
    // 
    // Avec les filtres, c'est encore pire pour d1 et d2 
    // ------------------------------------------------------------------------


    //Vérification/transformation en matrice Saisir:  
    saisir=div(saisir);
                 
    //vérification de l'absence de valeurs manquantes
    test=sum(sum(isnan(saisir.d)));
    if(test>0) then
        xdisp('There are ',test,' Nan values in the data table');
        xdisp('It is necessary to resolve the problem before using savgol');
        error(' ');
    end


    // choix du degré du polynome
    if argn(2)==4 then
        degre=degre_polynom;
    else
        degre=2;       // degré 2 
    end

    if degre<2 | degre>5 then
        error('the degree of the polynom should be between 2 and 5')
    end


    //Initialisation:
    [n,p]=size(saisir.d);
    filter_d0.d=zeros(p,p);  //9mars18
    filter_d1.d=zeros(p,p);
    filter_d2.d=zeros(p,p);
    
    rayon=floor((window2-1)/2);
    plage=2*rayon+1;
  
    // verification de la fenêtre          
    if rayon<1 | plage ~= window2 then 
        error('the window should be an odd integer, equal to 5 or more ');
    end

    d0_savgol=saisir.d;
    d1_savgol=ones(n,p);
    d2_savgol=ones(n,p);
    m=zeros(6,plage);


    // ---------------------------------------------------------
    // calculs de Savitzsky-Golay sur le spectre, cotés exceptés 
    // ---------------------------------------------------------

    //Calcul de m  
    vecteur=[-rayon:1:rayon];
    m(1,:)=vecteur.^0;  
    m(2,:)=vecteur.^1;  
    m(3,:)=vecteur.^2;  
    m(4,:)=vecteur.^3;  
    m(5,:)=vecteur.^4;  
    m(6,:)=vecteur.^5;  
    
    // ajustement de m à la dimension du polynome
    m=m(1:degre+1,:);   

    c=inv(m*m')*m; 

    for a=rayon+1:p-rayon;  // -1;  12mars18
        d0_savgol(:,a)=saisir.d(:,a-rayon:a+rayon)*c(1,:)';
        d1_savgol(:,a)=saisir.d(:,a-rayon:a+rayon)*c(2,:)';
        d2_savgol(:,a)=2*saisir.d(:,a-rayon:a+rayon)*c(3,:)';
        //la multiplication par 2 correspond à la dérivée seconde de: a+bx+cx**2
 
        // calcul des filtres 
        filter_d0.d(a-rayon:a+rayon,a)=c(1,:)';
        filter_d1.d(a-rayon:a+rayon,a)=c(2,:)';
        filter_d2.d(a-rayon:a+rayon,a)=2*c(3,:)';
    end;
  

    // ---------------------------------------------------
    // calculs de Savitzsky-Golay sur les cotés du spectre
    // --------------------------------------------------- 


    // 1 - calcul le plus précis  

    //d0 à gauche
    calcul=saisir.d(:,1:2*rayon+1)*c'*m;
    d0_savgol(:,1:rayon)=calcul(:,1:rayon);
    //d0 à droite
    calcul=saisir.d(:,p-2*rayon:p)*c'*m;
    d0_savgol(:,p-rayon:p)=calcul(:,rayon+1:2*rayon+1);
    
    //d1 à gauche
    d1=d0_savgol(:,2:p)-d0_savgol(:,1:p-1);
    d1_savgol(:,1:rayon)=d1(:,1:rayon);
    calcul=d1_savgol(:,1:2*rayon+1)*m'*inv(m*m')*m;
    d1_savgol(:,1:rayon)=calcul(:,1:rayon);
    //d1 à droite
    d1_savgol(:,p-rayon:p)=d1(:,p-rayon-1:p-1);   //perte d'un indice à cause de la dérivée première
    calcul=d1_savgol(:,p-2*rayon:p)*m'*inv(m*m')*m;
    d1_savgol(:,p-rayon:p)=calcul(:,rayon+1:2*rayon+1);

    //d2 à gauche
    d2=d1_savgol(:,2:p-1)-d1_savgol(:,1:p-2);
    d2_savgol(:,1:rayon)=d2(:,1:rayon);
    calcul=d2_savgol(:,1:2*rayon+1)*m'*inv(m*m')*m;
    d2_savgol(:,1:rayon)=calcul(:,1:rayon);
    //d2 à droite
    d2_savgol(:,p-rayon:p)=d2(:,p-rayon-2:p-2);   //perte de deux indices à cause de la dérivée seconde
    calcul=d2_savgol(:,p-2*rayon:p)*m'*inv(m*m')*m;
    d2_savgol(:,p-rayon:p)=calcul(:,rayon+1:2*rayon+1);


//    // 2- calcul sous forme de filtres (+ généralisable mais moins précis)                    
//    for i=1:rayon;
//        mi=zeros(6,2*i+1); 
//        vi=[-i:1:i];
//        mi(1,:)=vi.^0;  
//        mi(2,:)=vi.^1;  
//        mi(3,:)=vi.^2;  
//        mi(4,:)=vi.^3;  
//        mi(5,:)=vi.^4;  
//        mi(6,:)=vi.^5; 
//        mi=mi(1:degre+1,:);  
//        ci=inv(mi*mi')*mi;                                      //dimensions: (degre+1) x (2*rayon+1)
//        filter_d0.d(1:2*i+1,i+1)=           ci(1,:)';           // à gauche
//        filter_d0.d(p-(2*i+1)+1:p,p-i)=     ci(1,:)';           // à droite 
//        filter_d1.d(1:2*i+1,i+1)=           ci(2,:)';           // à gauche
//        filter_d1.d(p-(2*i+1)+1:p,p-i)=     ci(2,:)';           // à droite 
//        filter_d2.d(1:2*i+1,i+1)=           2*ci(3,:)';         // à gauche
//        filter_d2.d(p-(2*i+1)+1:p,p-i)=     2*ci(3,:)';         // à droite 
//        if i==1 then
//            filter_d0.d(1,1)=1;                                 // meilleur lissage pour 1 valeur = elle-même
//            filter_d0.d(p,p)=1;
//            filter_d1.d(1:2,1)=[-1;1];                          // on soustrait le 1° au 2° 
//            filter_d1.d(p-1:p,p)=[-1;1];
//            filter_d2.d(1:2*i+1,1)=         2*ci(3,:)';         // on recopie les valeurs à coté; pas trouvé mieux
//            filter_d2.d(p-(2*i+1)+1:p,p)=   2*ci(3,:)';
//        end
//    end


    //sélection des sorties       
    if d012==0 then
        res_savgol.d=d0_savgol;
  //      filter_d012=filter_d0;
    elseif d012==1 then 
        res_savgol.d=d1_savgol;
 //       filter_d012=filter_d1;
    elseif d012==2 then
        res_savgol.d=d2_savgol;
 //       filter_d012=filter_d2;
    end
    
    
    // div
    res_savgol.i=saisir.i;
    res_savgol.v=saisir.v;
    res_savgol=div(res_savgol);
    
//    filter_d012.i=saisir.v;
//    filter_d012.v=saisir.v;
//    filter_d012=div(filter_d012);
    


endfunction;