view src/mz_8_peakpicking_dw.sci @ 2:7e757edfa3b0 draft default tip

Deleted selected files
author jcb-mpl
date Thu, 12 May 2022 12:15:57 +0000
parents c0a2e5e78725
children
line wrap: on
line source

function  [peaks_list,peaks_detail,no_eic]=mz_8_peakpicking_dw(x_tr_mz3_trouves,list_mz,plage_tr,plage_min_tr,dw_max,signal_min,option_plot)
    
    // x_tr_mz3_trouves:    un div (n x p)
    // list_mz:             les n° des m/z à traiter; valeurs entre 1 et p
    // plage tr:            une valeur; par defaut: +/- 0,2 minutes
    // plage_min_tr:        une valeur; par defaut: 0,05
    // dw_max:              une valeur entre 0 et 2; par défaut: 0.5
    // signal_min:          une valeur; par défaut: 100000
    // option_plot:         1=figures; 0= non (par défaut)
    
    if argn(2)<7 then
        option_plot=0;
    end
    
    if argn(2)<6 then
        signal_min=100000;
    end
    
    if argn(2)<5 then
        dw_max=0.5;
    end
    
    if argn(2)<4 then
        plage_min_tr=0.05;
    end
    
    if argn(2)<3 then
        plage_tr=0.2; 
    end
    
    
    no_eic=[];

    n_mz=size(x_tr_mz3_trouves.d,2);
    list_mz_possible=[1:n_mz]';
    list_mz=list_mz_possible(list_mz);  // pour accepter le $
   
    n=max(size(list_mz));
    
    peaks_list.d=[]; 
    peaks_list.i=[];
    peaks_list.v=['TRmin';'TRmax';'TRmedian';'m/z';'signal'];
    
    peaks_detail=list();
    index_peaks_detail=1;
    
    // transformation de plage_tr (des valeurs de temps) en n_tr (un nbre de TR)   
    ntemp=100;  // pour ne pas commencer au début 
    n_tr=0;
    val_tr=0;
    
    while val_tr<=plage_tr;
        n_tr=n_tr+1;
        val_tr=strtod(x_tr_mz3_trouves.i(ntemp+n_tr))-strtod(x_tr_mz3_trouves.i(ntemp));
    end
    
    nbr_total_tr=size(x_tr_mz3_trouves.i,1);             
    
    for i=1:n;
        xi=x_tr_mz3_trouves(:,list_mz(i));
        
        signal_max=max(xi.d)                // 1 seule colonne 
        
        index_selected=0*strtod((xi.i));  // un vecteur de 0 de même taille que xi.i
        
        while signal_max > signal_min;
            // identification de la plage du pic
            index=find(xi.d==signal_max & index_selected==0);
            index_min=max(1,index-n_tr);
            index_max=min(nbr_total_tr,index+n_tr); 
            //ajustement du min pour ne pas retomber sur une zone deja selectionnee
            index_selected_min=index_selected(index_min);
            while index_selected_min>0;
                index_min=index_min+1;
                index_selected_min=index_selected(index_min);
            end 
             //ajustement du max pour ne pas retomber sur une zone deja selectionnee
            index_selected_max=index_selected(index_max);
            while index_selected_max>0;
                index_max=index_max-1;
                index_selected_max=index_selected(index_max);
            end
           
            // options de tri: DW et plage_min_tr
            xi_sel=xi(index_min:index_max,:);
            
            if (strtod(xi.i(index_max))-strtod(xi.i(index_min)))> plage_min_tr then
            
                dw=mz_6bis_durbin_watson(detrending(xi_sel')');       // car si ldb élevée et régulière: DW= faible 
           
                if dw<dw_max then 
           
                    // calcul du signal de la plage 
                    vect_tr=[index_min:index_max]
                    signal=sum(xi.d(vect_tr));  
           
                    // enregistrement du resultat
                    if peaks_list.d==[] then
                        peaks_list.d=[strtod(xi.i(index_min))  strtod(xi.i(index_max)) strtod(xi.i(index)) strtod(xi.v) signal];
                        peaks_list.i=i;
                    else 
                        peaks_list.d=[peaks_list.d; [strtod(xi.i(index_min)) strtod(xi.i(index_max)) strtod(xi.i(index)) strtod(xi.v) signal]];
                        peaks_list.i=[peaks_list.i;i];
                    end
                
                    peaks_detail(index_peaks_detail)=xi([index_min:index_max],:);
                    index_peaks_detail=index_peaks_detail+1;
        
                    // edition de figure (selon option)
                    if option_plot==1 then 
                        figure;
                        curves(xi([index_min:index_max],:)); 
                    end          
            
                    // rajout des n° d'EICs de x_tr_mz3_trouves ou x_tr_mz4
                    if no_eic==[] then
                        no_eic=i;
                    else
                        no_eic=[no_eic;i];
                    end
                end
            
            end 
            // neutralisation de la plage precedente
            index_selected(index_min:index_max)=1;      // on supprime
            xi.d(index_min:index_max)=0;
            
            signal_max=max(xi.d);
        end
        
    end
    
    
    // préparation des sorties
    peaks_list=div(peaks_list);
    
    

    
endfunction