@@ -14,44 +14,55 @@ function gsd(x, y, options){
1414 var options = Object . create ( options || { } ) ;
1515 if ( options . minMaxRatio === undefined ) options . minMaxRatio = 0.00025 ;
1616 if ( options . broadRatio === undefined ) options . broadRatio = 0.00 ;
17- if ( options . noiseLevel === undefined ) options . noiseLevel = 0 ;
17+ if ( options . noiseLevel === undefined ) options . noiseLevel = undefined ;
18+ if ( options . noiseFactor === undefined ) options . noiseFactor = 3 ;
1819 if ( options . maxCriteria === undefined ) options . maxCriteria = true ;
1920 if ( options . smoothY === undefined ) options . smoothY = true ;
2021 if ( options . realTopDetection === undefined ) options . realTopDetection = false ;
2122
2223 var sgOptions = extend ( { } , sgDefOptions , options . sgOptions ) ;
2324
24- //Transform y to use the standard algorithm.
25- var yCorrection = { m : 1 , b : 0 } ;
26- if ( ! options . maxCriteria || options . noiseLevel > 0 ) {
27- y = [ ] . concat ( y ) ;
28- if ( ! options . maxCriteria ) {
29- yCorrection = { m : - 1 , b : stats . array . max ( y ) } ;
30- for ( var i = 0 ; i < y . length ; i ++ ) {
31- y [ i ] = - y [ i ] + yCorrection . b ;
25+ //console.log(JSON.stringify(stats.array.minMax(y)));
26+ if ( options . noiseLevel === undefined ) {
27+ //We have to know if x is equally spaced
28+ var maxDx = 0 , minDx = Number . MAX_VALUE , tmp ;
29+ for ( var i = 0 ; i < x . length - 1 ; i ++ ) {
30+ var tmp = Math . abs ( x [ i + 1 ] - x [ i ] ) ;
31+ if ( tmp < minDx ) {
32+ minDx = tmp ;
3233 }
33- options . noiseLevel = - options . noiseLevel + yCorrection . b ;
34- }
35- if ( options . noiseLevel > 0 ) {
36- for ( var i = 0 ; i < y . length ; i ++ ) {
37- if ( Math . abs ( y [ i ] ) < options . noiseLevel ) {
38- y [ i ] = 0 ;
39- }
34+ if ( tmp > maxDx ) {
35+ maxDx = tmp ;
4036 }
4137 }
42- }
4338
44- //We have to know if x is equally spaced
45- var maxDx = 0 , minDx = Number . MAX_VALUE , tmp ;
46- for ( var i = 0 ; i < x . length - 1 ; i ++ ) {
47- var tmp = Math . abs ( x [ i + 1 ] - x [ i ] ) ;
48- if ( tmp < minDx ) {
49- minDx = tmp ;
39+ if ( ( maxDx - minDx ) / maxDx < 0.05 ) {
40+
41+ options . noiseLevel = getNoiseLevel ( y ) ;
42+ //console.log(options.noiseLevel+" "+stats.array.median(y));
5043 }
51- if ( tmp > maxDx ) {
52- maxDx = tmp ;
44+ else {
45+ options . noiseLevel = 0 ;
46+ }
47+ }
48+ //console.log("options.noiseLevel "+options.noiseLevel);
49+ y = [ ] . concat ( y ) ;
50+ var yCorrection = { m :1 , b :options . noiseLevel } ;
51+ if ( ! options . maxCriteria ) {
52+ yCorrection . m = - 1 ;
53+ yCorrection . b *= - 1 ;
54+ }
55+
56+ for ( var i = 0 ; i < y . length ; i ++ ) {
57+ y [ i ] = yCorrection . m * y [ i ] - yCorrection . b ;
58+ }
59+
60+ for ( var i = 0 ; i < y . length ; i ++ ) {
61+ if ( y [ i ] < 0 ) {
62+ y [ i ] = 0 ;
5363 }
5464 }
65+
5566 //If the max difference between delta x is less than 5%, then, we can assume it to be equally spaced variable
5667 var Y = y ;
5768 if ( ( maxDx - minDx ) / maxDx < 0.05 ) {
@@ -66,7 +77,7 @@ function gsd(x, y, options){
6677 var dY = SG ( y , x , { windowSize :sgOptions . windowSize , polynomial :sgOptions . polynomial , derivative :1 } ) ;
6778 var ddY = SG ( y , x , { windowSize :sgOptions . windowSize , polynomial :sgOptions . polynomial , derivative :2 } ) ;
6879 }
69-
80+
7081 var X = x ;
7182 var dx = x [ 1 ] - x [ 0 ] ;
7283 var maxDdy = 0 ;
@@ -122,9 +133,6 @@ function gsd(x, y, options){
122133 }
123134 }
124135 }
125- if ( options . realTopDetection ) {
126- realTopDetection ( minddY , X , Y ) ;
127- }
128136 //
129137 //console.log(intervalL.length+" "+minddY.length+" "+broadMask.length);
130138 var signals = [ ] ;
@@ -156,6 +164,7 @@ function gsd(x, y, options){
156164 //console.log(height);
157165 if ( Math . abs ( Y [ minddY [ j ] ] ) > options . minMaxRatio * maxY ) {
158166 signals . push ( {
167+ i :minddY [ j ] ,
159168 x : frequency ,
160169 y : ( Y [ minddY [ j ] ] - yCorrection . b ) / yCorrection . m ,
161170 width :Math . abs ( intervalR [ possible ] - intervalL [ possible ] ) , //widthCorrection
@@ -165,6 +174,16 @@ function gsd(x, y, options){
165174 }
166175 }
167176
177+
178+ if ( options . realTopDetection ) {
179+ realTopDetection ( signals , X , Y ) ;
180+ }
181+
182+ //Correct the values to fit the original spectra data
183+ for ( var j = 0 ; j < signals . length ; j ++ ) {
184+ signals [ j ] . base = options . noiseLevel ;
185+ }
186+
168187 signals . sort ( function ( a , b ) {
169188 return a . x - b . x ;
170189 } ) ;
@@ -173,14 +192,35 @@ function gsd(x, y, options){
173192
174193}
175194
195+ function getNoiseLevel ( y ) {
196+ var mean = 0 , stddev = 0 ;
197+ var length = y . length , i = 0 ;
198+ for ( i = 0 ; i < length ; i ++ ) {
199+ mean += y [ i ] ;
200+ }
201+ mean /= length ;
202+ var averageDeviations = new Array ( length ) ;
203+ for ( i = 0 ; i < length ; i ++ )
204+ averageDeviations [ i ] = Math . abs ( y [ i ] - mean ) ;
205+ averageDeviations . sort ( ) ;
206+ if ( length % 2 == 1 ) {
207+ stddev = averageDeviations [ ( length - 1 ) / 2 ] / 0.6745 ;
208+ } else {
209+ stddev = 0.5 * ( averageDeviations [ length / 2 ] + averageDeviations [ length / 2 - 1 ] ) / 0.6745 ;
210+ }
211+
212+ return stddev ;
213+ }
214+
176215function realTopDetection ( peakList , x , y ) {
177216 //console.log(peakList);
178217 //console.log(x);
179218 //console.log(y);
180219 var listP = [ ] ;
181220 var alpha , beta , gamma , p , currentPoint ;
182221 for ( var j = 0 ; j < peakList . length ; j ++ ) {
183- currentPoint = peakList [ j ] ; //peakList[j][2];
222+ currentPoint = peakList [ j ] . i ; //peakList[j][2];
223+ var tmp = currentPoint ;
184224 //The detected peak could be moved 1 or 2 unit to left or right.
185225 if ( y [ currentPoint - 1 ] >= y [ currentPoint - 2 ]
186226 && y [ currentPoint - 1 ] >= y [ currentPoint ] ) {
@@ -211,10 +251,12 @@ function realTopDetection(peakList, x, y){
211251 beta = 20 * Math . log10 ( y [ currentPoint ] ) ;
212252 gamma = 20 * Math . log10 ( y [ currentPoint + 1 ] ) ;
213253 p = 0.5 * ( alpha - gamma ) / ( alpha - 2 * beta + gamma ) ;
214-
215- x [ peakList [ j ] ] = x [ currentPoint ] + ( x [ currentPoint ] - x [ currentPoint - 1 ] ) * p ;
216- y [ peakList [ j ] ] = y [ currentPoint ] - 0.25 * ( y [ currentPoint - 1 ]
217- - [ currentPoint + 1 ] ) * p ; //signal.peaks[j].intensity);
254+ //console.log("p: "+p);
255+ //console.log(x[currentPoint]+" "+tmp+" "+currentPoint);
256+ peakList [ j ] . x = x [ currentPoint ] + ( x [ currentPoint ] - x [ currentPoint - 1 ] ) * p ;
257+ peakList [ j ] . y = y [ currentPoint ] - 0.25 * ( y [ currentPoint - 1 ]
258+ - y [ currentPoint + 1 ] ) * p ; //signal.peaks[j].intensity);
259+ //console.log(y[tmp]+" "+peakList[j].y);
218260 }
219261 }
220262}
0 commit comments