Mercurial > repos > jcb-mpl > eics_centroids
changeset 0:0d2d20ef6e72 draft
Uploaded
| author | jcb-mpl |
|---|---|
| date | Tue, 27 Apr 2021 13:52:21 +0000 |
| parents | |
| children | 96cffbbe245e |
| files | 2018_07_centroids.xml src/mz_1bis_signal_threshold.sci src/mz_2_barycenter_adjust2.sci src/mz_2_savgol2.sci src/mz_3_remove_zeros.sci src/mz_eic5_step1.sci src/mz_hdf5tosci.sci src/savgol_xl.sci test-data/VI2016_AC_4Areduit_1.dat test-data/VI2016_AC_4Areduit_1.h5 test-data/VI2016_AC_4Areduit_1centroide.mat |
| diffstat | 11 files changed, 836 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/2018_07_centroids.xml Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,207 @@ +<tool id="2021_05_centroids" name="MS .h5 or .dat -> MS-centroid" version="0.0.1"> + +<description> </description> + +<requirements> + <requirement type="package" version="6.0.2">scilab</requirement> +</requirements> + + +<stdio> + <exit_code range="1:" level="fatal" /> +</stdio> + + +<command> + <![CDATA[ + if [ -d $__root_dir__/packages/scilab-6.1.0 ]; then $__root_dir__/packages/scilab-6.1.0/bin/scilab-cli -nb -quit -f $* < ${script_file}; else scilab-cli -nb -quit -f $* < ${script_file}; fi + ]]> +</command> + + <configfiles> + <configfile name="script_file"> + <![CDATA[ exec("$__tool_directory__/src/mz_1bis_signal_threshold.sci",-1); ... + exec("$__tool_directory__/src/mz_2_barycenter_adjust2.sci",-1); ... + exec("$__tool_directory__/src/mz_2_savgol2.sci",-1); ... + exec("$__tool_directory__/src/mz_3_remove_zeros.sci",-1); ... + exec("$__tool_directory__/src/mz_eic5_step1.sci",-1); ... + exec("$__tool_directory__/src/mz_hdf5tosci.sci",-1); ... + if ~isdef('pls') then ... + atomsSystemUpdate(); ... + atomsInstall('FACT'); ... + atomsLoad('FACT'); ... + end; ... + lasterror(); ... + type_sm=${type_sm}; ... + format_io=${format_io}; ... + if type_sm==1 then ... // mode profil Orbitrap + if format_io==0 then ... // format Scilab + load("${filein}"); ... // la variable doit s'appeler res + elseif format_io==1 then ... // format HDF5 + res=mz_hdf5tosci("${filein}"); ... + end; ... + resbary=mz_eic5_step1(res,${signalthresh}); ... + elseif type_sm==0 then ... // mode centroide + if format_io==0 then ... // format Scilab = on ne fait rien + load("${filein}"); ... // la variable doit deja s'appeler resbary + elseif format_io==1 then ... // format HDF5 + res=mz_hdf5tosci("${filein}"); ... + resbary=mz_3_remove_zeros(res); ... + end; ... + end; ... + save("${resbary}",'resbary'); ... + if ~isempty(lasterror(%f)) then ... + write(0,lasterror()); ... + end; ]]> +</configfile> + </configfiles> + + + <inputs> + <param name="format_io" type="select" size="5" format="integer" label="MS format" > + <option value="0" selected="true" > Scilab format .dat </option> + <option value="1" > HDF5 format .h5 </option> + </param> + <param name="filein" format="mat,h5,dat" type="data" label="MS data" /> + <param name="signalthresh" size="5" value="0" type="integer" label="Baseline threshold" help="if not applied during the spectrum acquisition" /> + <param name="type_sm" type="select" size="5" format="integer" label="Mode during spectrum acquisition" help="if centroid: 0 are removed; if profile: centroids are calculated" > + <option value="0" selected="true" > Centroid </option> + <option value="1" > Profile mode with Orbitrap </option> + </param> + </inputs> + + + <outputs> + <data name="resbary" format="mat" label="MS-centroid" > + </data> + </outputs> + + + <tests> + + <test> + <param name="format_io" value="0"/> + <param name="filein" value="VI2016_AC_4Areduit_1.dat"/> + <param name="signalthresh" value="0"/> + <param name="type_sm" value="1"/> + <output name="resbary" file="VI2016_AC_4Areduit_1centroide.mat" compare="sim_size" delta="2000" /> + </test> + + <test> + <param name="format_io" value="1"/> + <param name="filein" value="VI2016_AC_4Areduit_1.h5"/> + <param name="signalthresh" value="0"/> + <param name="type_sm" value="1"/> + <output name="resbary" file="VI2016_AC_4Areduit_1centroide.mat" compare="sim_size" delta="2000" /> + </test> + + </tests> + + + <help> + + +**Author** Jean-Claude Boulet (INRAE). + + +--------------------------------------------------- + +==================== +CENTROID CALCULATION +==================== + +----------- +Description +----------- + + +Centroids are calculated according to the reference cited below. + + +------ +Inputs +------ + +**MS format** + +Two formats are possible: HDF5 obtained from Scilab (.dat) or HDF5 obtained from MSconvert then HDFView (.h5) + +The Scilab format contains a structure with two fields: + +- time: the retention times, a vector; + +- mzdata: a list of same length as time; each element of the list is a matrix of two columns representing a MS1 mass spectrum: m/zvalues and associated signals. + +The HDF5 (.h5) file contains at least the following fields: + +- ChomatogramTime: a vector containing the retention times of each mass spectrum. + +- ChromatogramIndex: the indexes of the signals identified in ChomatogramTime; the last value of ChromatogramIndex is the size of ChomatogramTime. + +- ChromatogramIntensity: the sum of intensities of the mass spectra associated to the RTs of ChomatogramTime; ChomatogramTime and ChromatogramIntensity are vectors of same length. + +- SpectrumIndex: a vecteur of integers corresponding to the numbers of signals of each MS in SpectrumMZ et SpectrumIntensity. For example, if the two first values of SpectrumIndex are 1421 and 2887, then the first MS is between 0 à 1421, the second between 1422 and 2887. The last value of spectrumIndex is the number of values in SpectrumMZ and in SpectrumIntensity. + +- SpectrumMz: the m/z values. For each mass spectrum, the first value is the real m/z value, but the following are obtained by the sum of the previous values. For example, if the three first values are: 220.0501, 0.1171 and 0.9856, then the real m/z values to be calculated will be: 220.0501, 220.1672=220.0501+0.1171, 220.2528=220.0501+0.1171+0.9856. + +- SpectrumIntensity: a vector containing all the intensities recorded. + +Even if the MS spectrum has been obtained with the centroid option of MSconvert, this function should be applied to remove all the zeros kept by the MSconcert peak picking option. + + +**MS data** + +The mass spectrum + + +---------- +Parameters +---------- + + +**Baseline threshold** + +Signals under the threshold are dropped. Default value = 0. + +Most of the peaks should be separated by regions with signal=0. + + +**Mode during spectrum acquisition** + +Profile or centroid. + + + +------- +Outputs +------- + +**MS-centroid** + +For each RT, the signals associated to the same m/z value have been replaced by their centroid. + +The output is a structure with two fields: + +- time: retention times; + +- mzdata: m/z values and centroid signals. + + +</help> + + +<citations> + +<citation type="bibtex">@article{sgcentroid, + title={High-resolution mass spectrometry (HRMS): focus on the m/z values estimated by the Savitzky-Golay first derivative}, + author={Boulet, J.C. and Meudec, E. and Vallverdu-Queralt, A. and Cheynier, V.}, + journal={Rapid Communications in Mass Spectrometry}, + year={2020}, + doi={10.1002/rcm.9036} + } +</citation> + +</citations> + + +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_1bis_signal_threshold.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,22 @@ +function res2=mz_1bis_signal_threshold(res,sigthresh); + + res2=res; + + n=max(size(res2.time)); + + for i=1:n; + resi=res2.mzdata(i); + + resi_col2=resi(:,2); + + // application du seuil + resi_col2(resi_col2<sigthresh)=0; + + resi(:,2)=resi_col2; + + res2.mzdata(i)=resi; + end + + + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_2_barycenter_adjust2.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,28 @@ +function xout=mz_2_barycenter_adjust2(xin) + + // remplace les donnees brutes par: + // (1) mz barycentrique + // (2) somme des signaux + + // xin: une structure avec 2 champs: + // xin.time + // xin.mzdata + + // xout: 2 champs, comme xin + + // re-ecrit totalement et simplifie le 9dec20 + + xout.time=xin.time; + xout.mzdata=list(); + n=max(size(xin.time)); + + for i=1:n; + xtemp_in=xin.mzdata(i); + [nul,xtemp_out]=mz_2_savgol2(xtemp_in); + xout.mzdata(i)=xtemp_out; + end + + + + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_2_savgol2.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,221 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_3_remove_zeros.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,18 @@ +function x_out=mz_3_remove_zeros(x_in) + + // x_in: un fichier au format barycentrique avec 2 champs: .time et .mzdata + // si le calcul centroide a été fait par MSconvert, le fichier contient des m/z avec signal=0: on va les enlever + + // initialisation + x_out.time=x_in.time; + x_out.mzdata=list(); + + n=max(size(x_in.time)); + + for i=1:n; + mzdata_i=x_in.mzdata(i); + index=find(mzdata_i(:,2)>0); + x_out.mzdata(i)=mzdata_i(index,:); + end + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_eic5_step1.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,45 @@ +function resbary=mz_eic5_step1(res,sigthresh) + + // cette fonction mobilise les fonctions suivantes: + // exec('mz_1bis_signal_threshold.sci',-1); + // exec('mz_2_barycenter_adjust2.sci',-1); + // exec(‘mz_2_savgol2.sci’,-1); + + // Arguments: + // res: une variable avec 2 champs: time et mzdata + // sigthresh: un seuil pour un seul signal sur un seul TR ; ex: sigthresh=10000 + + // resbary: un fichier barycentrique + + + res=mz_1bis_signal_threshold(res,sigthresh); + + resbary=mz_2_barycenter_adjust2(res); + + +endfunction + + + + + + + + + + + + + + + + + + + + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_hdf5tosci.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,103 @@ +function x_mz=mz_hdf5tosci(file_in,max_scans) + + // entrees: ------------------------- + // file_in: un fichier au format HDF5 + // max_scans: valeur introduite pour accelerer la fonction de test + + // sorties: --------------------------- + // x_mz: une variable avec deux champs: + // x_mz.time + // x_mz.mzdata une liste dont chaque élément est un vecteur de 2 colonnes: m/z et intensités de signal + + // JCB 2 aout18 + + h5close() + + // lecture du fichier HDF5 ---------------------------------- + a=h5open(file_in,'r','stdio'); // que lecture 6jan20:stdio + //pause + spectrum_chromatogram_index=h5read(a.root.ChromatogramIndex); + spectrum_chromatogram_index=spectrum_chromatogram_index(1); // une valeur + spectrum_time=h5read(a.root.ChomatogramTime); + spectrum_intensity=h5read(a.root.SpectrumIntensity); + spectrum_mz=h5read(a.root.SpectrumMZ); + spectrum_index=h5read(a.root.SpectrumIndex); + //spectrum_index=spectrum_index(1:spectrum_chromatogram_index); + h5close(a) + + // reconstruction du spectre ----------------- + if argn(2)<2 then + nbr_scans=max(size(spectrum_index)); + else + nbr_scans=max_scans; + end + + + // 1- extraction du temps + //25juin20 + time=spectrum_time; + if size(time,1)==1 then + time=time'; + end + //pas=(spectrum_time($)-spectrum_time(1))/(nbr_scans-1); + //time=[spectrum_time(1):pas:spectrum_time($)+pas]; // rajout de 1 par securite pour être sur d'arriver a nbr_scans + //time=time(1:nbr_scans); // nbr exact de RT + + mzdata=list(); + + //disp(nbr_scans,'nbr_scans=') + + for i=1:nbr_scans; + if i==1 then + plage_i=[1,double(spectrum_index(i))]; + else + plage_i=[double(spectrum_index(i-1))+1,double(spectrum_index(i))]; + end + //disp(plage_i,'plage_i=') + // 2- extraction et reconstruction des m/z + // dans le fichier HDF5, ne sont enregistrées que les différences de m/z entre 2 mesures + mz_i=spectrum_mz(plage_i(1):plage_i(2)); + + n=size(mz_i,2); + mz_i2=zeros(1,n); + + mz_i2(1)=mz_i(1); + for j=2:n; + mz_i2(j)=mz_i2(j-1)+mz_i(j); + end + + // 3- extraction des signaux + sig_i=spectrum_intensity(plage_i(1):plage_i(2)); + + // regroupement des m/z et signaux + sp_i=[mz_i2;sig_i]; + + mzdata(i)=sp_i'; + + + end + + // verif. que time est un vecteur colonne + if size(time,2)>1 then + time=time'; + end + + + // amelioration de la coherence des donnees + index=time(2:$)-time(1:$-1); + index2=find(index<0); + if index2 ~=[] then + index3=min(index2); + else; + index3=max(size(time)); + end + + x_mz.time=time(1:index3); + mzdata2=list(); + for i=1:min(index3,size(mzdata)); + mzdata2(i)=mzdata(i); + end + x_mz.mzdata=mzdata2; + + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/savgol_xl.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,192 @@ +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; +