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

Uploaded

function  [xout2,xout3,no_pic,signe_d1]=mz_2_savgol2(xin)
    
    // xin: une matrice de 2 colonnes 
    // 1° colonne = valeurs de m/z successives 
    // 2° colonne = intensités de signal 
 
    // xout: un vecteur-ligne ou matrice de 3 colonnes
    // dans l'ordre: m/z_min m/z_savgol m/z_max  
    
    // xout2: une matrice de 2 colonnes 
    // 1° colonne = valeurs de m/z remplacées par leur barycentre 
    // 2° colonne = intensités de signal 
    
    // xout3: une matrice de 2 colonnes 
    // 1° colonne: m/z (uniques)
    // 2° colonne: somme des intensités de signal 
 
 
//    // marche bien avec l'exemple ci-dessous, 9 ou 18 premières lignes  
//    xin_1=[
//    775.08677  804.23633;
//    775.09492   19743.041;
//    775.10307   62518.859;
//    775.11123   100056.77;
//    775.11938   108296.29;
//    775.12753   83583.969;
//    775.13568   45296.258;
//    775.14384   15782.141;
//    775.15199   2796.1245;
//    775.16014   3257.4946;
//    775.1683    17404.266;
//    775.17645   41285.391;
//    775.1846    61480.969;
//    775.19276   65116.719;
//    775.20091   50944.094;
//    775.20906   29804.268;
//    775.21722   12828.898;
//    775.22537   2491.0898];
//    
// ou sur: 
//   355.6344    8. 
//   355.64264   8. 
//   355.65088   10.
//   355.65909   7. 
//   355.66733   4. 
//   355.67557   11.
//   355.68381   40.

// ou cet exemple 
//
//   x=[102.12401   0.    ;   
//   102.12428   0.       ;
//   102.12457   0.       ;
//   102.12485   0. ;
//   102.12513   0.  ;     
//   102.1254    0.   ;    
//   102.12568   18844.752;
//   102.12596   70650.805;
//   102.12623   130295.48;
//   102.12651   187543.03;
//   102.12679   521590.84;
//   102.12707   1696004. ;
//   102.12734   3549963.2;
//   102.12762   5541599. ;
//   102.1279    5764036. ;
//   102.12817   3820761.2;
//   102.12845   1933240. ;
//   102.12873   544016.81;
//   102.129     120472.32;
//   102.12928   50397.285;
//   102.12956   20005.551;
//   102.12983   0.       ;
//   102.1301    0. ;
//   102.13038   0.  ;     
//   102.13065   0.   ;    
//   102.13093   0. ];


// 
    // choix du traitement selon la taille de n
    
    n=size(xin,1);
    no_pic=zeros(n,1);
    
    if n<3 then 
        all_pics=[];    // on ne traite pas 
       
        
    elseif n==3 then
        all_pics=zeros(1,3);
        all_pics(1)=xin(1,1);
        all_pics(3)=xin(3,1);
        all_pics(2)=xin(2,1);
   
        
    elseif n==4 then 
        all_pics=zeros(1,3);
        all_pics(1)=xin(1,1);  // m/z min
        all_pics(3)=xin(4,1);  // m/z max 
        x1=xin(2,1);
        x2=xin(3,1);
        y1=xin(2,2);  
        y2=xin(3,2);   
        x0=(x2*y1-x1*y2)/(y1-y2);
        all_pics(2)=x0;
   
        
    elseif n>4 then 
        // calcul des dérivées 1°
        d1=savgol_xl(xin(:,2)',5,1); 
        // Obtention des signes: 
        signe_d1=sign(d1.d);
        
        // recherche des valeurs de signal nulles //8dec20
        ind_null=find(xin(:,2)==0);
        signe_d1(ind_null)=0;
        
  //  pause    
        // correction des valeurs nulles 
   
   // 8dec20     
//        // correction des bords de Savgol:
//        signe_d1(1)=1;
//        signe_d1(2)=1;
//        signe_d1($-1)=-1;
//        signe_d1($)=-1;
      
        // suppression le 31 mai 18 --------------------------  -> pas de lissage si pics intercallés, ex: ex: -1 -1 1 -1 -1

        // identification des pics par des n°: 1, 2, 3, etc.  27juin17
        no_pic=zeros(n,1);
        no_pic(1)=1;
        for i=2:n;
           if signe_d1(i-1)<0 & signe_d1(i)>=0 then
               no_pic(i)=no_pic(i-1)+1;
           else
               no_pic(i)=no_pic(i-1);
           end
        end
  
        
        // identification des pics
        all_pics=[];
        // pour chaque pic
        i=1;
        pic=zeros(1,3);
        
        while i<=n;                
            if pic(1)==0 then 
                pic(1)=xin(i,1);
            end
            if i<n & signe_d1(i+1)< signe_d1(i)  then   // un sommet 
                x1=xin(i,1);
                x2=xin(i+1,1);
                y1=d1.d(i);
                y2=d1.d(i+1);
                x0=(x2*y1-x1*y2)/(y1-y2);             // calcul précis du m/z 
                pic(2)=x0;
            elseif i<n & signe_d1(i+1)> signe_d1(i) & pic(2)==0 then
                pic=zeros(1,3);                         // 3oct18 fin de la boucle sans rien faire 

            elseif i<n & signe_d1(i+1)> signe_d1(i) & pic(2)~=0 then   // une vallée 
                pic(3)=xin(i,1);
                if all_pics==[] then 
                    all_pics=pic;
                else
                    all_pics=[all_pics;pic];
                end 
                pic=zeros(1,3);   
            elseif i==n & pic(1)~=0 & pic(2)~=0 then 
                pic(3)=xin(i,1);
                if all_pics==[] then 
                    all_pics=pic;
                else
                    all_pics=[all_pics;pic];
                end 
            end
            i=i+1;
        end
        
        if size(all_pics,1)== max(no_pic) then // 3oct18
            break
        end
    end

    // sortie du résultat
    xout=all_pics;
    
    xout2=xin;   // 27juin17
    n_out=size(xout,1);
    if no_pic($)~=0 & no_pic($) ~=[] then
        n_pics=no_pic($); 
        //if n_out ~=n_pics then                //3oct18
            //disp('ligne 154')
            //pause
            //error('n_out different de n_pics')
        //end 
        for i=1:n_out;
            liste_pics=find(no_pic==i);
            xout2(liste_pics,1)=xout(i,2);
        end
    end
    
    // somme des signaux 9mai18
    xout2_unique=unique(xout2(:,1));
    n_unique=max(size(xout2_unique));
    xout3=zeros(n_unique,2);
    xout3(:,1)=xout2_unique;
    for i=1:n_unique;
        tri=find(xout2(:,1)==xout3(i,1));
        sum_signals=sum(xout2(tri,:),'r');
        xout3(i,2)=sum_signals(2);
    end
    
    
    // tri de xout3 (enlever les valeurs de signal nulles)  7dec20
    tri2=find(xout3(:,2))~=0;
    xout3=xout3(tri2,:);  
   
    
endfunction