multiplet-analysis - npm Package Compare versions

Comparing version 0.1.0 to 0.2.0

History.md

@@ -0,1 +1,14 @@
+# [0.2.0](https://github.com/cheminfo/multiplet-analysis/compare/v0.1.0...v0.2.0) (2020-07-27)
+
+
+### Bug Fixes
+
+* correctly take the next power of two ([beb9e6d](https://github.com/cheminfo/multiplet-analysis/commit/beb9e6d36cba977e80b39e8df7f403226e8f5234))
+* minor points to pass tests ([8cc92b3](https://github.com/cheminfo/multiplet-analysis/commit/8cc92b39770a32817b8a744e11620cb6bd188d76))
+* remove obsolete test ([63ebe4d](https://github.com/cheminfo/multiplet-analysis/commit/63ebe4dcc367c315bc51d26976c5e6f2d65d3f4d))
+* set code back for clean log of messages associtated to commit ([f651b2f](https://github.com/cheminfo/multiplet-analysis/commit/f651b2ff8d8a389344a7a2940988cdc19676945f))
+* smaller steps in fast mode / fast mode also for debug option ([351bfdd](https://github.com/cheminfo/multiplet-analysis/commit/351bfdd2a89bdd19e260583b998880ccef828887))
+
+
+
 # 0.1.0 (2020-03-29)
@@ -2,0 +15,0 @@

lib/index.js

@@ -58,6 +58,94 @@ 'use strict';
 
+function decofast1(yi, jStar, sign, nbLines, addspace = 0) {
+  let y1 = new Float64Array(yi.length + addspace);
+  for (let scan = 0; scan < addspace; scan++) {
+    y1[scan] = 0;
+  }
+  for (let scan = 0; scan < yi.length; scan++) {
+    y1[scan + addspace] = yi[scan];
+  }
+  if (sign === 1) {
+    if (nbLines === 1) {
+      for (let scan = 0; scan < y1.length - jStar; scan++) {
+        y1[scan + jStar + addspace] -= y1[scan + addspace];
+      }
+    } else {
+      for (let scan = 0; scan < y1.length - jStar * nbLines; scan++) {
+        for (let line = 0; line < nbLines; line++) {
+          y1[scan + (line + 1) * jStar + addspace] -= y1[scan + addspace];
+        }
+      }
+    }
+  } else {
+    if (nbLines === 1) {
+      for (let scan = 0; scan < y1.length - jStar; scan++) {
+        y1[scan + jStar + addspace] += y1[scan + addspace];
+      }
+    } else {
+      for (let scan = 0; scan < y1.length - jStar * nbLines; scan++) {
+        for (let line = 0; line < nbLines; line++) {
+          y1[scan + (line + 1) * jStar + addspace] += y1[scan + addspace];
+        }
+      }
+    }
+  }
+  return y1;
+}
+function decofast2(yi, jStar, sign, nbLines, addspace = 0) {
+  let y2 = new Float64Array(yi.length + addspace);
+
+  for (let scan = 0; scan < yi.length; scan++) {
+    y2[scan] = yi[scan];
+  }
+  for (let scan = 0; scan < addspace; scan++) {
+    y2[scan + yi.length] = 0;
+  }
+  if (sign === 1) {
+    if (nbLines === 1) {
+      for (
+        let scan = y2.length - 1 - addspace;
+        scan >= jStar * nbLines;
+        scan -= 1
+      ) {
+        y2[scan - jStar] -= y2[scan];
+      }
+    } else {
+      for (
+        let scan = y2.length - 1 - addspace;
+        scan >= jStar * nbLines;
+        scan -= 1
+      ) {
+        for (let line = 0; line < nbLines; line++) {
+          y2[scan - (line + 1) * jStar] -= y2[scan];
+        }
+      }
+    }
+  } else {
+    if (nbLines === 1) {
+      for (
+        let scan = y2.length - 1 - addspace;
+        scan >= jStar * nbLines;
+        scan -= 1
+      ) {
+        y2[scan - jStar] += y2[scan];
+      }
+    } else {
+      for (
+        let scan = y2.length - 1 - addspace;
+        scan >= jStar * nbLines;
+        scan -= 1
+      ) {
+        for (let line = 0; line < nbLines; line++) {
+          y2[scan - (line + 1) * jStar] += y2[scan];
+        }
+      }
+    }
+  }
+  return y2;
+}
+
 /**
  *
  * @param {number} [yi]
- * @param {number} [JStar]
+ * @param {number} [jStar]
  * @param {number} [sign=1]  ++ multiplet :1 +- multiplet : -1
@@ -70,3 +158,3 @@ * @param {number} [dir=1] from left to right : 1 -1 from right to left, 0: sum of both
   yi,
-  JStar,
+  jStar,
   sign = 1,
@@ -77,44 +165,52 @@ dir = 1,
 ) {
-  let nbLines = parseInt(2 * multiplicity); // 1 for doublet (spin 1/2) 2, for spin 1, etc... never tested...
-  let y1 = new Array(yi.length);
-  let y2 = new Array(yi.length);
+  let nbLines = parseInt(2 * multiplicity, 10); // 1 for doublet (spin 1/2) 2, for spin 1, etc... never tested...
+  let y1;
+  let y2;
+  if (dir > -1) {
+    y1 = decofast1(yi, jStar, sign, nbLines, 0);
 
-  if (dir > -1) {
-    for (let scan = 0; scan < y1.length; scan++) y1[scan] = yi[scan];
-    for (let scan = 0; scan < y1.length - JStar * nbLines; scan++) {
-      for (let line = 0; line < nbLines; line++) {
-        y1[scan + (line + 1) * JStar] -= sign * y1[scan];
-      }
-    }
     if (dir > 0.5) {
-      return y1.slice(0, y1.length - chopTail * JStar * nbLines);
+      return new Float64Array(
+        y1.buffer,
+        0,
+        y1.length - chopTail * jStar * nbLines,
+      );
     }
   }
   if (dir < 1) {
-    for (let scan = 0; scan < y2.length; scan++) y2[scan] = yi[scan];
-    for (let scan = y2.length - 1; scan >= JStar * nbLines; scan -= 1) {
-      for (let line = 0; line < nbLines; line++) {
-        y2[scan - (line + 1) * JStar] -= sign * y2[scan];
-      }
-    }
+    y2 = decofast2(yi, jStar, sign, nbLines, 0);
+
     if (dir < -0.2) {
-      return y2.slice(chopTail * JStar * nbLines, y2.length);
+      return new Float64Array(
+        y2.buffer,
+        chopTail * jStar * nbLines * 8,
+        y2.length - chopTail * jStar * nbLines,
+      );
     }
   }
+  if (!y1) y1 = new Float64Array(yi.length);
+  if (!y2) y2 = new Float64Array(yi.length);
   if (dir === 0) {
-    for (let scan = 0; scan < y2.length - JStar * nbLines; scan++) {
-      y1[scan] = (y1[scan] + sign * y2[scan + JStar * nbLines]) / 2;
+    for (let scan = 0; scan < y2.length - jStar * nbLines; scan++) {
+      y1[scan] = (y1[scan] + sign * y2[scan + jStar * nbLines]) / 2;
     }
-    return y1.slice(0, y1.length - JStar * nbLines);
+    return new Float64Array(y1.buffer, 0, y1.length - jStar * nbLines);
   }
-  let half = ((y1.length - JStar * nbLines) / 2.0) | 0;
+  let half = ((y1.length - jStar * nbLines) / 2.0) | 0;
   if (dir === 0.1) {
-    for (let scan = half; scan < y2.length - JStar * nbLines; scan++) {
-      y1[scan] = sign * y2[scan + JStar * nbLines];
+    for (let scan = half; scan < y2.length - jStar * nbLines; scan++) {
+      y1[scan] = sign * y2[scan + jStar * nbLines];
     }
-    return y1.slice(0, y1.length - JStar * nbLines);
+    return new Float64Array(y1.buffer, 0, y1.length - jStar * nbLines);
   }
 }
 
-function scalarProduct(y1, y2, sens, incrementForSpeed) {
+function scalarProduct(
+  y1,
+  y2,
+  sens,
+  incrementForSpeed,
+  first = 0,
+  last = y1.length,
+) {
   // sens = 1; crude scalar product
@@ -126,3 +222,3 @@ // sens =-1: flip spectrum first
   if (sens > 0) {
-    for (let index = 0; index < y1.length; index += incrementForSpeed) {
+    for (let index = first; index < last; index += incrementForSpeed) {
       v12 += y1[index] * y2[index];
@@ -134,4 +230,4 @@ v11 += y1[index] * y1[index];
     // Here as if flip left/right array
-    for (let index = 0; index < y1.length; index += incrementForSpeed) {
-      v12 += y1[index] * y2[y1.length - index - 1];
+    for (let index = first; index < last; index += incrementForSpeed) {
+      v12 += y1[index] * y2[last - (index - first) - 1];
       v11 += y1[index] * y1[index];
@@ -146,3 +242,3 @@ v22 += y2[index] * y2[index];
   y, // input vector
-  JStar, // tested value of J in pt
+  jStar, // tested value of J in pt
   sign, // sign of the split 1: ++ -1: +-
@@ -153,6 +249,10 @@ chopTail, // 1: cut tail
 ) {
-  let y1 = [];
-  let y2 = [];
-  y1 = deco(y, JStar, sign, 1, chopTail, multiplicity); // dir left to right
-  y2 = deco(y, JStar, sign, -1, chopTail, multiplicity); // dir right to left
+  //let y1 = deco(y, jStar, sign, 1, chopTail, multiplicity); // dir left to right
+  //let y2 = deco(y, jStar, sign, -1, chopTail, multiplicity); // dir right to left
+
+  let nbLines = parseInt(2 * multiplicity, 10); // 1 for doublet (spin 1/2) 2, for spin 1, etc... never tested...
+  let y1 = decofast1(y, jStar, sign, nbLines, jStar);
+  // console.log(`y1 :: ` + y1);
+  let y2 = decofast2(y, jStar, sign, nbLines, jStar);
+  // console.log(`y2 :: ` + y2);
   return sign * scalarProduct(y1, y2, 1, incrementForSpeed);
@@ -162,36 +262,21 @@ }
 function measureSymShift(y) {
-  let ScalarProductReference;
-  let ScalarProductNewValue;
+  let scalarProductReference;
+  let scalarProductNewValue;
   let movedBy = 0;
-  ScalarProductReference = scalarProduct(y, y, -1, 1);
+  scalarProductReference = scalarProduct(y, y, -1, 1);
   // search left...
   for (let i = 1; i < y.length / 2; i++) {
-    ScalarProductNewValue = scalarProduct(
-      y.slice(i, y.length),
-      y.slice(i, y.length),
-      -1,
-      1,
-    );
-    if (ScalarProductNewValue > ScalarProductReference) {
-      //  console.log(`${ScalarProductNewValue} > ${ScalarProductReference} size: ${y.length}`);
-      ScalarProductReference = ScalarProductNewValue;
+    scalarProductNewValue = scalarProduct(y, y, -1, 1, i, y.length);
+    if (scalarProductNewValue > scalarProductReference) {
+      scalarProductReference = scalarProductNewValue;
       movedBy = i;
     }
   }
-  if (movedBy === 0) {
-    for (let i = 1; i < y.length / 2; i++) {
-      ScalarProductNewValue = scalarProduct(
-        y.slice(0, y.length - i),
-        y.slice(0, y.length - i),
-        -1,
-        1,
-      );
-      if (ScalarProductNewValue > ScalarProductReference) {
-        //  console.log(`${ScalarProductNewValue} > ${ScalarProductReference} size: ${y.length}`);
-        ScalarProductReference = ScalarProductNewValue;
-        movedBy = -i;
-      }
+  for (let i = 1; i < y.length / 2; i++) {
+    scalarProductNewValue = scalarProduct(y, y, -1, 1, 0, y.length - i);
+    if (scalarProductNewValue > scalarProductReference) {
+      scalarProductReference = scalarProductNewValue;
+      movedBy = -i;
     }
   }
-  //console.log('OK  ' + movedBy);
   return movedBy;
@@ -216,4 +301,4 @@ }
 ) {
-  let sca = Array(numberOfPointOutput);
-  let spe = Array(numberOfPointOutput);
+  let scale = new Float64Array(numberOfPointOutput);
+  let spectrum = new Float64Array(numberOfPointOutput);
 
@@ -224,9 +309,9 @@ let scaIncrement = (x[x.length - 1] - x[0]) / (x.length - 1); // delta one pt
   let trueWidth = scaIncrement * x.length;
-  scaIncrement = trueWidth / numberOfPointOutput;
+  let scaIncrementInterp = trueWidth / numberOfPointOutput;
   //scaPt += scaIncrement / 2; // move from limit side to middle of first pt
   // set scale
 
-  for (let loop = 0; loop < numberOfPointOutput; loop++) {
-    sca[loop] = scaPt;
-    scaPt += scaIncrement;
+  for (let i = 0; i < numberOfPointOutput; i++) {
+    scale[i] = scaPt;
+    scaPt += scaIncrementInterp;
   }
@@ -236,7 +321,5 @@
   // prepare input for fft apply fftshift
-  let nextPowerTwoInput = Math.pow(
-    2,
-    Math.round(Math.log(y.length - 1) / Math.log(2.0) + 0.5),
-  );
-  //nextPowerTwoAn = Math.pow(2, Math.ceil(Math.log(y.length) / Math.log(2)));
+  let nextPowerTwoInput = 2 ** Math.ceil(Math.log2(y.length));
+  let nextPowerTwoOut = 2 ** Math.ceil(Math.log2(numberOfPointOutput));
+
   let an = [...Array(nextPowerTwoInput)].map(() => Array(2).fill(0)); // n x 2 array
@@ -246,63 +329,51 @@ const halfNumPt = Math.floor(y.length / 2); // may ignore last pt... if odd number
   const shiftMultiplet = nextPowerTwoInput - y.length;
-  for (let loop = 0; loop < halfNumPt; loop++) {
-    an[shiftMultiplet + loop + halfNumPtB][0] = y[loop]; //Re
-    an[shiftMultiplet + loop + halfNumPtB][1] = 0; //Im
+  for (let i = 0; i < halfNumPt; i++) {
+    an[shiftMultiplet + i + halfNumPtB][0] = y[i]; //Re
+    an[shiftMultiplet + i + halfNumPtB][1] = 0; //Im
   }
-  for (let loop = 0; loop < halfNumPtB; loop++) {
-    an[loop][0] = y[loop + halfNumPt]; //Re
-    an[loop][1] = 0; //Im
+  for (let i = 0; i < halfNumPtB; i++) {
+    an[i][0] = y[i + halfNumPt]; //Re
+    an[i][1] = 0; //Im
   }
-
   let timeDomain = fftJs.ifft(an);
-  /*an = fft(out);
-  for (let loop = 0; loop < 2 * halfNumPt; loop++) {
-    an[loop][1] = 0; //* Math.cos((phase / 180) * Math.PI) -
-  }
-  out = ifft(an);*/
-
   timeDomain[0][0] = timeDomain[0][0] / 2; // divide first point by 2 Re
   timeDomain[0][1] = timeDomain[0][1] / 2; // divide first point by 2 Im
   // move to larger array...
-  let timeDomainZeroFilled = [...Array(numberOfPointOutput)].map(() =>
+  let timeDomainZeroFilled = [...Array(nextPowerTwoOut)].map(() =>
     Array(2).fill(0),
   ); // n x 2 array
-  for (let loop = 0; loop < halfNumPt; loop++) {
-    timeDomainZeroFilled[loop][0] = timeDomain[loop][0]; //* Math.cos((phase / 180) * Math.PI) +
-    timeDomainZeroFilled[loop][1] = timeDomain[loop][1]; //* Math.cos((phase / 180) * Math.PI) -
+  for (let i = 0; i < halfNumPt; i++) {
+    timeDomainZeroFilled[i][0] = timeDomain[i][0]; //* Math.cos((phase / 180) * Math.PI) +
+    timeDomainZeroFilled[i][1] = timeDomain[i][1]; //* Math.cos((phase / 180) * Math.PI) -
   }
-
-  for (let loop = halfNumPt; loop < numberOfPointOutput; loop++) {
+  for (let i = halfNumPt; i < nextPowerTwoInput; i++) {
     // zero filling
-    timeDomainZeroFilled[loop][0] = 0;
-    timeDomainZeroFilled[loop][1] = 0;
+    timeDomainZeroFilled[i][0] = 0;
+    timeDomainZeroFilled[i][1] = 0;
   }
+
   let interpolatedSpectrum = fftJs.fft(timeDomainZeroFilled);
   const halfNumPt2 = Math.floor(numberOfPointOutput / 2);
-  let tmp;
+
   // applies phase change
   let phaseRad = ((addPhaseInterpolation + 0.0) / 180.0) * Math.PI; // this is for testing additional phases
   if (phaseRad !== 0.0) {
-    for (let loop = 0; loop < 2 * halfNumPt2; loop++) {
-      tmp =
-        interpolatedSpectrum[loop][0] * Math.cos(phaseRad) +
-        interpolatedSpectrum[loop][1] * Math.sin(phaseRad); // only Re now...
-      interpolatedSpectrum[loop][1] =
-        -interpolatedSpectrum[loop][0] * Math.sin(phaseRad) +
-        interpolatedSpectrum[loop][1] * Math.cos(phaseRad); // only Re now...
-      interpolatedSpectrum[loop][0] = tmp;
+    for (let i = 0; i < 2 * halfNumPt2; i++) {
+      let tmp =
+        interpolatedSpectrum[i][0] * Math.cos(phaseRad) +
+        interpolatedSpectrum[i][1] * Math.sin(phaseRad); // only Re now...
+      interpolatedSpectrum[i][1] =
+        -interpolatedSpectrum[i][0] * Math.sin(phaseRad) +
+        interpolatedSpectrum[i][1] * Math.cos(phaseRad); // only Re now...
+      interpolatedSpectrum[i][0] = tmp;
     }
   }
+
   let returnedPhase = 0;
-
+  let norm;
   if (appliedPhaseCorrectionType > 0) {
-    let localPhaseRad = 0;
-    let vectx;
-    let vecty;
-    let norm;
-    // let sumNorms;
-
-    for (let loo = 1; loo < 100; loo++) {
-      localPhaseRad = 0;
-      vectx = 0;
-      vecty = 0;
+    for (let i = 1; i < 100; i++) {
+      let localPhaseRad = 0;
+      let vectx = 0;
+      let vecty = 0;
       // sumNorms = 0;
@@ -335,12 +406,9 @@ if (appliedPhaseCorrectionType > 0) {
       norm = Math.sqrt(vecty * vecty + vectx * vectx);
-      phaseRad = -(10.0 / (loo * loo)) * Math.sign(vecty);
+      phaseRad = -(10.0 / (i * i)) * Math.sign(vecty);
 
       returnedPhase -= (180.0 * phaseRad) / Math.PI;
 
-      // console.log('localPhase ' + (180.0 * localPhaseRad) / Math.PI);
-      //console.log('returnedPhase ' + returnedPhase);
-
       if (phaseRad !== 0.0) {
         for (let loop = 0; loop < 2 * halfNumPt2; loop++) {
-          tmp =
+          let tmp =
             interpolatedSpectrum[loop][0] * Math.cos(phaseRad) +
@@ -358,7 +426,8 @@ interpolatedSpectrum[loop][1] * Math.sin(phaseRad); // only Re now...
   // applies fftshift
-  for (let loop = 0; loop < halfNumPt2; loop++) {
-    spe[loop] = interpolatedSpectrum[loop + halfNumPt2][0]; // only Re now...
+  let dropPoints = nextPowerTwoOut - numberOfPointOutput;
+  for (let i = 0; i < halfNumPt2; i++) {
+    spectrum[i] = interpolatedSpectrum[halfNumPt2 + dropPoints + i][0]; // only Re now...
   }
-  for (let loop = 0; loop < halfNumPt2; loop++) {
-    spe[loop + halfNumPt2] = interpolatedSpectrum[loop][0];
+  for (let i = 0; i < halfNumPt2; i++) {
+    spectrum[i + halfNumPt2] = interpolatedSpectrum[i][0];
   }
@@ -368,5 +437,6 @@ if (returnedPhase > 360.0) {
   }
+
   return {
-    spectrum: spe,
-    scale: sca,
+    spectrum,
+    scale,
     phaseCorrectionOnMultipletInDeg: returnedPhase,
@@ -388,2 +458,6 @@ };
   let { x = [], y = [] } = data;
+
+  if (!(x instanceof Float64Array)) x = Float64Array.from(x);
+  if (!(y instanceof Float64Array)) y = Float64Array.from(y);
+
   const {
@@ -393,8 +467,9 @@ frequency = 400,
     maxTestedJ = 20,
-    minTestedJ = 0.5,
+    minTestedJ = 1,
     minimalResolution = 0.01,
-    makeShortCutForSpeed = 0,
-    critFoundJ = 0.95,
+    correctVerticalOffset = true,
+    makeShortCutForSpeed = true,
+    critFoundJ = 0.9,
     sign = 1,
-    chopTail = 1,
+    chopTail = true,
     multiplicity = 0.5,
@@ -406,7 +481,6 @@ symmetrizeEachStep = false,
     appliedPhaseCorrectionType = 0,
+    decreasingJvalues = true,
+    jumpUpAfterFoundValue = 2.0,
   } = options;
 
-  let scalProd = [];
-  let JStarArray = [];
-
   let result = {};
@@ -423,14 +497,25 @@ result.j = [];
 
-  let nextPowerTwo = Math.pow(
-    2,
-    Math.round(
-      Math.log((x.length - 1) * factorResolution) / Math.log(2.0) + 0.5,
-    ),
-  );
-  factorResolution = (factorResolution * nextPowerTwo) / x.length;
+  let nextPowerTwo = 2 ** Math.ceil(Math.log2(x.length * factorResolution));
+  let nextPowerTwoInital = 2 ** Math.ceil(Math.log2(x.length));
 
+  let integerFactorResolution = nextPowerTwo / nextPowerTwoInital;
+
   let movedBy;
-  let sca = [];
-  let spe = [];
+  let scale;
+  let spectrum;
   let topPosJ = 0;
+  // adjust vertical offset
+  if (correctVerticalOffset) {
+    let minValue = 1e20;
+    for (let index = 0; index < y.length; index++) {
+      if (minValue > y[index]) {
+        minValue = y[index];
+      }
+    }
+    if (minValue > 0) {
+      for (let index = 0; index < y.length; index++) {
+        y[index] -= minValue;
+      }
+    }
+  }
   if (resolutionHz > minimalResolution) {
@@ -441,16 +526,17 @@ // need increase resolution
       y,
-      nextPowerTwo,
+      integerFactorResolution * y.length,
       addPhaseInterpolation,
       appliedPhaseCorrectionType,
     );
-    if (debug) console.log(`interpolating`);
+    //if (debug) console.log(`interpolating`);
 
-    spe = returned.spectrum;
-    sca = returned.scale;
+    spectrum = returned.spectrum;
+    scale = returned.scale;
     result.phaseCorrectionOnMultipletInDeg =
       returned.phaseCorrectionOnMultipletInDeg;
   } else {
-    sca = x;
-    spe = y;
+    scale = x;
+    spectrum = y;
   }
+
   /// for testing break symmetry before running...
@@ -469,3 +555,4 @@ /*
 
-  resolutionPpm = Math.abs(sca[0] - sca[sca.length - 1]) / (sca.length - 1);
+  resolutionPpm =
+    Math.abs(scale[0] - scale[scale.length - 1]) / (scale.length - 1);
   resolutionHz = resolutionPpm * frequency;
@@ -480,12 +567,9 @@ let maxTestedPt = Math.trunc(maxTestedJ / resolutionHz);
 
-  for (let jStar = 0; jStar < minTestedPt; jStar++) {
-    JStarArray[jStar] = jStar * resolutionHz;
-    scalProd[jStar] = -1;
-  }
   let incrementForSpeed = 1;
   let curIncrementForSpeed;
 
-  if (!debug) {
-    incrementForSpeed = (1 + 0.5 / minimalResolution) | 0; // 1 could be set better (according to line widht ?!)
-  }
+  //if (!debug) {
+  //incrementForSpeed = (1 + 0.5 / minimalResolution) | 0; // not enough.... some peak sneek between points
+  incrementForSpeed = (1 + 0.3 / minimalResolution) | 0; // 1 could be set better (according to line widht ?!)
+  //}
   for (
@@ -496,22 +580,28 @@ let loopoverJvalues = 1;
   ) {
-    let beforeSymSpe = new Array(spe.length);
+    let scalProd = [];
+    let jStarArray = [];
+    for (let jStar = 0; jStar < minTestedPt; jStar++) {
+      jStarArray[jStar] = jStar * resolutionHz;
+      scalProd[jStar] = 0;
+    }
+    let beforeSymSpe = new Float64Array(spectrum.length);
 
     //symmetrize if requested to
     if (symmetrizeEachStep === true) {
-      movedBy = -measureSymShift(spe);
+      movedBy = -measureSymShift(spectrum);
       if (movedBy > 0) {
-        spe = spe.slice(0, spe.length - movedBy);
-        sca = sca.slice(0, sca.length - movedBy);
+        spectrum = spectrum.slice(0, spectrum.length - movedBy);
+        scale = scale.slice(0, scale.length - movedBy);
       }
       if (movedBy < 0) {
-        spe = spe.slice(-movedBy, spe.length);
-        sca = sca.slice(-movedBy, sca.length);
+        spectrum = spectrum.slice(-movedBy, spectrum.length);
+        scale = scale.slice(-movedBy, scale.length);
       }
       if (debug) {
         // save this to plot it as well
-        for (let index = 0; index < spe.length; index++) {
-          beforeSymSpe[index] = spe[index];
+        for (let index = 0; index < spectrum.length; index++) {
+          beforeSymSpe[index] = spectrum[index];
         }
       }
-      spe = symmetrize(spe);
+      spectrum = symmetrize(spectrum);
     }
@@ -521,10 +611,14 @@
     let gotJValue = false;
-    let LimitCoupling = Math.floor(sca.length / 2) - 1; //limit with respect to size of spectrum (which is reducing at each step)
+    let limitCoupling = scale.length - 1; //limit with respect to size of spectrum (which is reducing at each step)
     let critFoundJLow = critFoundJ - 0.3;
-    if (maxTestedPt > LimitCoupling) {
-      maxTestedPt = LimitCoupling;
+    if (maxTestedPt > limitCoupling) {
+      maxTestedPt = limitCoupling;
     }
-    if (loopoverJvalues > 1 && !debug) {
-      if (maxTestedPt > Math.floor(topPosJ + 1.0 / resolutionHz)) {
-        maxTestedPt = Math.floor(topPosJ + 1.0 / resolutionHz);
+    if (loopoverJvalues > 1 && decreasingJvalues) {
+      if (
+        maxTestedPt > Math.floor(topPosJ + jumpUpAfterFoundValue / resolutionHz)
+      ) {
+        maxTestedPt = Math.floor(
+          topPosJ + jumpUpAfterFoundValue / resolutionHz,
+        );
         //console.log(`recuded size to :: ` + maxTestedPt);
@@ -536,9 +630,13 @@ //console.log(`recuded size to :: ` + maxTestedPt + " = " + (maxTestedPt*resolutionHz));
     let jStarFine;
+    for (let jStar = minTestedPt; jStar < maxTestedPt; jStar++) {
+      jStarArray[jStar] = jStar * resolutionHz;
+      scalProd[jStar] = -1;
+    }
     for (
       let jStar = maxTestedPt;
       jStar >= minTestedPt;
-      jStar -= curIncrementForSpeed
+      jStar -= incrementForSpeed
     ) {
       scalProd[jStar] = measureDeco(
-        spe,
+        spectrum,
         jStar,
@@ -550,3 +648,4 @@ sign,
       );
-      JStarArray[jStar] = jStar * resolutionHz;
+
+      jStarArray[jStar] = jStar * resolutionHz;
       if (!gotJValue) {
@@ -557,9 +656,10 @@ if (scalProd[jStar] > topValue) {
         }
-        if (jStar < maxTestedPt) {
+
+        if (jStar < maxTestedPt - 2 * curIncrementForSpeed) {
           if (
             scalProd[jStar] < scalProd[jStar + curIncrementForSpeed] &&
-            topValue > critFoundJLow
+            scalProd[jStar + curIncrementForSpeed] >=
+              scalProd[jStar + 2 * curIncrementForSpeed] &&
+            scalProd[jStar + curIncrementForSpeed] > critFoundJLow
           ) {
-            // here refine...
-            jStarFine = jStar;
             while (curIncrementForSpeed > 1) {
@@ -569,9 +669,15 @@ //console.log(`curIncrementForSpeed:: ` + curIncrementForSpeed);
               curIncrementForSpeed = Math.floor(curIncrementForSpeed / 2); // get smaller and smaller step
+              let froms = topPosJ - 2 * curIncrementForSpeed; // maybe 1 is enough....
+              while (froms < 0) froms += curIncrementForSpeed;
+              let tos = topPosJ + 2 * curIncrementForSpeed; // maybe 1 is enough....
+              while (tos >= maxTestedPt) tos -= curIncrementForSpeed;
+              topValue = -1; // reset because increased precision may make the top lower
               for (
-                jStarFine = topPosJ - 2 * curIncrementForSpeed;
-                jStarFine < topPosJ + 2 * curIncrementForSpeed;
+                jStarFine = froms;
+                jStarFine <= tos;
                 jStarFine += curIncrementForSpeed
               ) {
+                // if (( scalProd[jStarFine] === -1) || true) { // this was a bad idea because precision reduced with earlyer tests
                 scalProd[jStarFine] = measureDeco(
-                  spe,
+                  spectrum,
                   jStarFine,
@@ -589,4 +695,4 @@ sign,
             }
+            curIncrementForSpeed = incrementForSpeed;
             // end refine
-
             if (topValue > critFoundJ) {
@@ -605,3 +711,3 @@ // ugly force value for tests
 
-              if (debug) console.log(`J:: ${topPosJ * resolutionHz}`);
+              //  if (debug) console.log(`J:: ${topPosJ * resolutionHz}`);
 
@@ -613,3 +719,6 @@ result.j.push({
               gotJValue = true;
-              if (makeShortCutForSpeed) break;
+
+              if (makeShortCutForSpeed) {
+                break;
+              }
             }
@@ -627,5 +736,5 @@ }
         appendDebug(
-          sca,
-          spe,
-          JStarArray,
+          scale,
+          spectrum,
+          jStarArray,
           scalProd,
@@ -637,3 +746,10 @@ loopoverJvalues,
       } else {
-        appendDebug(sca, spe, JStarArray, scalProd, loopoverJvalues, result);
+        appendDebug(
+          scale,
+          spectrum,
+          jStarArray,
+          scalProd,
+          loopoverJvalues,
+          result,
+        );
       }
@@ -646,4 +762,4 @@ }
       // apply here the deconvolution for the next step of the recursive process
-      spe = deco(
-        spe,
+      spectrum = deco(
+        spectrum,
         topPosJ,
@@ -657,5 +773,5 @@ sign,
         let remove = 0.5 * topPosJ * (2 * multiplicity);
-        sca = sca.slice(remove, sca.length - remove);
+        scale = scale.slice(remove, scale.length - remove);
       }
-      if (sca.length !== spe.length) {
+      if (scale.length !== spectrum.length) {
         throw Error('sts');
@@ -667,18 +783,4 @@ }
 
-  let maxAmplitudePosition = maxY({ x: sca, y: spe });
-  result.chemShift = sca[maxAmplitudePosition.index];
-
-  // for non-javascript implementation
-  /*
-  let curTop = Number.NEGATIVE_INFINITY;
-  let curChemShift;
-  for (let i = 0; i < spe.length; i++) {
-    if (spe[i] > curTop) {
-      curTop = spe[i];
-      curChemShift = sca[i];
-    }
-  }
-  result.chemShift = curChemShift;*/
-  // end to be commented
-
+  let maxAmplitudePosition = maxY({ x: scale, y: spectrum });
+  result.chemShift = scale[maxAmplitudePosition.index];
   return result;
@@ -685,0 +787,0 @@ }

package.json

 {
   "name": "multiplet-analysis",
-  "version": "0.1.0",
+  "version": "0.2.0",
   "description": "",
@@ -55,5 +55,5 @@ "main": "lib/index.js",
     "fft-js": "0.0.12",
-    "ml-array-xy-max-y": "^1.0.1",
-    "ml-fft": "^1.3.5"
+    "fft.js": "^4.0.3",
+    "ml-array-xy-max-y": "^1.0.1"
   }
 }

README.md

@@ -8,4 +8,6 @@ # multiplet-analysis
 The goal of this project is to be able to determine the multiplicity of
-a NMR signal as well as the coupling constants.
-
+a NMR signal as well as the coupling constants. It is based on a delta-function deconvolution
+developed by D. Jeannerat and G. Bodenhausen
+and [published](https://www.sciencedirect.com/science/article/abs/pii/S1090780799918451?via%3Dihub) in *J. Magn. Reson.* **1999**, 141(1), p133 doi:10.1006/jmre.1999.1845. 
+More info and discussion in doc/README.md 
 The result of the analysis that is an object composed:
@@ -22,12 +24,13 @@
 
-`node -r esm examples/quadruplet.js`
-`node -r esm examples/ddd.js`
-`node -r esm examples/ddd_ABCD.js`
-`node -r esm examples/doublet.js`
+```
+node -r esm examples/quadruplet.js
+node -r esm examples/ddd.js
+node -r esm examples/ddd_ABCD.js
+node -r esm examples/doublet.js
 
-`node -r esm examples/simulate.js` to simulate from a user-defined spin system
-`node -r esm examples/dd-exp.js` to simulate from a user-defined spin system
+node -r esm examples/simulate.js; # to simulate from a user-defined spin system
+node -r esm examples/dd-exp.js; # to simulate from a user-defined spin system
+```
 
 
-
 ## Developement
@@ -67,3 +70,3 @@
 If you want to execute those scripts written as module you need to use `esm` that is installed as a development dependency.
-
+`npm install esm`
 `node -r esm ./examples/web/exampleGenerateAnnotations.js`
@@ -70,0 +73,0 @@

src/deco.js

@@ -0,5 +1,93 @@
+export function decofast1(yi, jStar, sign, nbLines, addspace = 0) {
+  let y1 = new Float64Array(yi.length + addspace);
+  for (let scan = 0; scan < addspace; scan++) {
+    y1[scan] = 0;
+  }
+  for (let scan = 0; scan < yi.length; scan++) {
+    y1[scan + addspace] = yi[scan];
+  }
+  if (sign === 1) {
+    if (nbLines === 1) {
+      for (let scan = 0; scan < y1.length - jStar; scan++) {
+        y1[scan + jStar + addspace] -= y1[scan + addspace];
+      }
+    } else {
+      for (let scan = 0; scan < y1.length - jStar * nbLines; scan++) {
+        for (let line = 0; line < nbLines; line++) {
+          y1[scan + (line + 1) * jStar + addspace] -= y1[scan + addspace];
+        }
+      }
+    }
+  } else {
+    if (nbLines === 1) {
+      for (let scan = 0; scan < y1.length - jStar; scan++) {
+        y1[scan + jStar + addspace] += y1[scan + addspace];
+      }
+    } else {
+      for (let scan = 0; scan < y1.length - jStar * nbLines; scan++) {
+        for (let line = 0; line < nbLines; line++) {
+          y1[scan + (line + 1) * jStar + addspace] += y1[scan + addspace];
+        }
+      }
+    }
+  }
+  return y1;
+}
+export function decofast2(yi, jStar, sign, nbLines, addspace = 0) {
+  let y2 = new Float64Array(yi.length + addspace);
+
+  for (let scan = 0; scan < yi.length; scan++) {
+    y2[scan] = yi[scan];
+  }
+  for (let scan = 0; scan < addspace; scan++) {
+    y2[scan + yi.length] = 0;
+  }
+  if (sign === 1) {
+    if (nbLines === 1) {
+      for (
+        let scan = y2.length - 1 - addspace;
+        scan >= jStar * nbLines;
+        scan -= 1
+      ) {
+        y2[scan - jStar] -= y2[scan];
+      }
+    } else {
+      for (
+        let scan = y2.length - 1 - addspace;
+        scan >= jStar * nbLines;
+        scan -= 1
+      ) {
+        for (let line = 0; line < nbLines; line++) {
+          y2[scan - (line + 1) * jStar] -= y2[scan];
+        }
+      }
+    }
+  } else {
+    if (nbLines === 1) {
+      for (
+        let scan = y2.length - 1 - addspace;
+        scan >= jStar * nbLines;
+        scan -= 1
+      ) {
+        y2[scan - jStar] += y2[scan];
+      }
+    } else {
+      for (
+        let scan = y2.length - 1 - addspace;
+        scan >= jStar * nbLines;
+        scan -= 1
+      ) {
+        for (let line = 0; line < nbLines; line++) {
+          y2[scan - (line + 1) * jStar] += y2[scan];
+        }
+      }
+    }
+  }
+  return y2;
+}
+
 /**
  *
  * @param {number} [yi]
- * @param {number} [JStar]
+ * @param {number} [jStar]
  * @param {number} [sign=1]  ++ multiplet :1 +- multiplet : -1
@@ -12,3 +100,3 @@ * @param {number} [dir=1] from left to right : 1 -1 from right to left, 0: sum of both
   yi,
-  JStar,
+  jStar,
   sign = 1,
@@ -19,41 +107,42 @@ dir = 1,
 ) {
-  let nbLines = parseInt(2 * multiplicity); // 1 for doublet (spin 1/2) 2, for spin 1, etc... never tested...
-  let y1 = new Array(yi.length);
-  let y2 = new Array(yi.length);
+  let nbLines = parseInt(2 * multiplicity, 10); // 1 for doublet (spin 1/2) 2, for spin 1, etc... never tested...
+  let y1;
+  let y2;
+  if (dir > -1) {
+    y1 = decofast1(yi, jStar, sign, nbLines, 0);
 
-  if (dir > -1) {
-    for (let scan = 0; scan < y1.length; scan++) y1[scan] = yi[scan];
-    for (let scan = 0; scan < y1.length - JStar * nbLines; scan++) {
-      for (let line = 0; line < nbLines; line++) {
-        y1[scan + (line + 1) * JStar] -= sign * y1[scan];
-      }
-    }
     if (dir > 0.5) {
-      return y1.slice(0, y1.length - chopTail * JStar * nbLines);
+      return new Float64Array(
+        y1.buffer,
+        0,
+        y1.length - chopTail * jStar * nbLines,
+      );
     }
   }
   if (dir < 1) {
-    for (let scan = 0; scan < y2.length; scan++) y2[scan] = yi[scan];
-    for (let scan = y2.length - 1; scan >= JStar * nbLines; scan -= 1) {
-      for (let line = 0; line < nbLines; line++) {
-        y2[scan - (line + 1) * JStar] -= sign * y2[scan];
-      }
-    }
+    y2 = decofast2(yi, jStar, sign, nbLines, 0);
+
     if (dir < -0.2) {
-      return y2.slice(chopTail * JStar * nbLines, y2.length);
+      return new Float64Array(
+        y2.buffer,
+        chopTail * jStar * nbLines * 8,
+        y2.length - chopTail * jStar * nbLines,
+      );
     }
   }
+  if (!y1) y1 = new Float64Array(yi.length);
+  if (!y2) y2 = new Float64Array(yi.length);
   if (dir === 0) {
-    for (let scan = 0; scan < y2.length - JStar * nbLines; scan++) {
-      y1[scan] = (y1[scan] + sign * y2[scan + JStar * nbLines]) / 2;
+    for (let scan = 0; scan < y2.length - jStar * nbLines; scan++) {
+      y1[scan] = (y1[scan] + sign * y2[scan + jStar * nbLines]) / 2;
     }
-    return y1.slice(0, y1.length - JStar * nbLines);
+    return new Float64Array(y1.buffer, 0, y1.length - jStar * nbLines);
   }
-  let half = ((y1.length - JStar * nbLines) / 2.0) | 0;
+  let half = ((y1.length - jStar * nbLines) / 2.0) | 0;
   if (dir === 0.1) {
-    for (let scan = half; scan < y2.length - JStar * nbLines; scan++) {
-      y1[scan] = sign * y2[scan + JStar * nbLines];
+    for (let scan = half; scan < y2.length - jStar * nbLines; scan++) {
+      y1[scan] = sign * y2[scan + jStar * nbLines];
     }
-    return y1.slice(0, y1.length - JStar * nbLines);
+    return new Float64Array(y1.buffer, 0, y1.length - jStar * nbLines);
   }
 }

src/index.js

@@ -28,2 +28,6 @@ /* eslint-disable no-console */
   let { x = [], y = [] } = data;
+
+  if (!(x instanceof Float64Array)) x = Float64Array.from(x);
+  if (!(y instanceof Float64Array)) y = Float64Array.from(y);
+
   const {
@@ -33,8 +37,9 @@ frequency = 400,
     maxTestedJ = 20,
-    minTestedJ = 0.5,
+    minTestedJ = 1,
     minimalResolution = 0.01,
-    makeShortCutForSpeed = 0,
-    critFoundJ = 0.95,
+    correctVerticalOffset = true,
+    makeShortCutForSpeed = true,
+    critFoundJ = 0.9,
     sign = 1,
-    chopTail = 1,
+    chopTail = true,
     multiplicity = 0.5,
@@ -46,7 +51,6 @@ symmetrizeEachStep = false,
     appliedPhaseCorrectionType = 0,
+    decreasingJvalues = true,
+    jumpUpAfterFoundValue = 2.0,
   } = options;
 
-  let scalProd = [];
-  let JStarArray = [];
-
   let result = {};
@@ -63,14 +67,25 @@ result.j = [];
 
-  let nextPowerTwo = Math.pow(
-    2,
-    Math.round(
-      Math.log((x.length - 1) * factorResolution) / Math.log(2.0) + 0.5,
-    ),
-  );
-  factorResolution = (factorResolution * nextPowerTwo) / x.length;
+  let nextPowerTwo = 2 ** Math.ceil(Math.log2(x.length * factorResolution));
+  let nextPowerTwoInital = 2 ** Math.ceil(Math.log2(x.length));
 
+  let integerFactorResolution = nextPowerTwo / nextPowerTwoInital;
+
   let movedBy;
-  let sca = [];
-  let spe = [];
+  let scale;
+  let spectrum;
   let topPosJ = 0;
+  // adjust vertical offset
+  if (correctVerticalOffset) {
+    let minValue = 1e20;
+    for (let index = 0; index < y.length; index++) {
+      if (minValue > y[index]) {
+        minValue = y[index];
+      }
+    }
+    if (minValue > 0) {
+      for (let index = 0; index < y.length; index++) {
+        y[index] -= minValue;
+      }
+    }
+  }
   if (resolutionHz > minimalResolution) {
@@ -81,16 +96,17 @@ // need increase resolution
       y,
-      nextPowerTwo,
+      integerFactorResolution * y.length,
       addPhaseInterpolation,
       appliedPhaseCorrectionType,
     );
-    if (debug) console.log(`interpolating`);
+    //if (debug) console.log(`interpolating`);
 
-    spe = returned.spectrum;
-    sca = returned.scale;
+    spectrum = returned.spectrum;
+    scale = returned.scale;
     result.phaseCorrectionOnMultipletInDeg =
       returned.phaseCorrectionOnMultipletInDeg;
   } else {
-    sca = x;
-    spe = y;
+    scale = x;
+    spectrum = y;
   }
+
   /// for testing break symmetry before running...
@@ -109,3 +125,4 @@ /*
 
-  resolutionPpm = Math.abs(sca[0] - sca[sca.length - 1]) / (sca.length - 1);
+  resolutionPpm =
+    Math.abs(scale[0] - scale[scale.length - 1]) / (scale.length - 1);
   resolutionHz = resolutionPpm * frequency;
@@ -120,12 +137,9 @@ let maxTestedPt = Math.trunc(maxTestedJ / resolutionHz);
 
-  for (let jStar = 0; jStar < minTestedPt; jStar++) {
-    JStarArray[jStar] = jStar * resolutionHz;
-    scalProd[jStar] = -1;
-  }
   let incrementForSpeed = 1;
   let curIncrementForSpeed;
 
-  if (!debug) {
-    incrementForSpeed = (1 + 0.5 / minimalResolution) | 0; // 1 could be set better (according to line widht ?!)
-  }
+  //if (!debug) {
+  //incrementForSpeed = (1 + 0.5 / minimalResolution) | 0; // not enough.... some peak sneek between points
+  incrementForSpeed = (1 + 0.3 / minimalResolution) | 0; // 1 could be set better (according to line widht ?!)
+  //}
   for (
@@ -136,22 +150,28 @@ let loopoverJvalues = 1;
   ) {
-    let beforeSymSpe = new Array(spe.length);
+    let scalProd = [];
+    let jStarArray = [];
+    for (let jStar = 0; jStar < minTestedPt; jStar++) {
+      jStarArray[jStar] = jStar * resolutionHz;
+      scalProd[jStar] = 0;
+    }
+    let beforeSymSpe = new Float64Array(spectrum.length);
 
     //symmetrize if requested to
     if (symmetrizeEachStep === true) {
-      movedBy = -measureSymShift(spe);
+      movedBy = -measureSymShift(spectrum);
       if (movedBy > 0) {
-        spe = spe.slice(0, spe.length - movedBy);
-        sca = sca.slice(0, sca.length - movedBy);
+        spectrum = spectrum.slice(0, spectrum.length - movedBy);
+        scale = scale.slice(0, scale.length - movedBy);
       }
       if (movedBy < 0) {
-        spe = spe.slice(-movedBy, spe.length);
-        sca = sca.slice(-movedBy, sca.length);
+        spectrum = spectrum.slice(-movedBy, spectrum.length);
+        scale = scale.slice(-movedBy, scale.length);
       }
       if (debug) {
         // save this to plot it as well
-        for (let index = 0; index < spe.length; index++) {
-          beforeSymSpe[index] = spe[index];
+        for (let index = 0; index < spectrum.length; index++) {
+          beforeSymSpe[index] = spectrum[index];
         }
       }
-      spe = symmetrize(spe);
+      spectrum = symmetrize(spectrum);
     }
@@ -161,10 +181,14 @@
     let gotJValue = false;
-    let LimitCoupling = Math.floor(sca.length / 2) - 1; //limit with respect to size of spectrum (which is reducing at each step)
+    let limitCoupling = scale.length - 1; //limit with respect to size of spectrum (which is reducing at each step)
     let critFoundJLow = critFoundJ - 0.3;
-    if (maxTestedPt > LimitCoupling) {
-      maxTestedPt = LimitCoupling;
+    if (maxTestedPt > limitCoupling) {
+      maxTestedPt = limitCoupling;
     }
-    if (loopoverJvalues > 1 && !debug) {
-      if (maxTestedPt > Math.floor(topPosJ + 1.0 / resolutionHz)) {
-        maxTestedPt = Math.floor(topPosJ + 1.0 / resolutionHz);
+    if (loopoverJvalues > 1 && decreasingJvalues) {
+      if (
+        maxTestedPt > Math.floor(topPosJ + jumpUpAfterFoundValue / resolutionHz)
+      ) {
+        maxTestedPt = Math.floor(
+          topPosJ + jumpUpAfterFoundValue / resolutionHz,
+        );
         //console.log(`recuded size to :: ` + maxTestedPt);
@@ -176,9 +200,13 @@ //console.log(`recuded size to :: ` + maxTestedPt + " = " + (maxTestedPt*resolutionHz));
     let jStarFine;
+    for (let jStar = minTestedPt; jStar < maxTestedPt; jStar++) {
+      jStarArray[jStar] = jStar * resolutionHz;
+      scalProd[jStar] = -1;
+    }
     for (
       let jStar = maxTestedPt;
       jStar >= minTestedPt;
-      jStar -= curIncrementForSpeed
+      jStar -= incrementForSpeed
     ) {
       scalProd[jStar] = measureDeco(
-        spe,
+        spectrum,
         jStar,
@@ -190,3 +218,4 @@ sign,
       );
-      JStarArray[jStar] = jStar * resolutionHz;
+
+      jStarArray[jStar] = jStar * resolutionHz;
       if (!gotJValue) {
@@ -197,9 +226,10 @@ if (scalProd[jStar] > topValue) {
         }
-        if (jStar < maxTestedPt) {
+
+        if (jStar < maxTestedPt - 2 * curIncrementForSpeed) {
           if (
             scalProd[jStar] < scalProd[jStar + curIncrementForSpeed] &&
-            topValue > critFoundJLow
+            scalProd[jStar + curIncrementForSpeed] >=
+              scalProd[jStar + 2 * curIncrementForSpeed] &&
+            scalProd[jStar + curIncrementForSpeed] > critFoundJLow
           ) {
-            // here refine...
-            jStarFine = jStar;
             while (curIncrementForSpeed > 1) {
@@ -209,9 +239,15 @@ //console.log(`curIncrementForSpeed:: ` + curIncrementForSpeed);
               curIncrementForSpeed = Math.floor(curIncrementForSpeed / 2); // get smaller and smaller step
+              let froms = topPosJ - 2 * curIncrementForSpeed; // maybe 1 is enough....
+              while (froms < 0) froms += curIncrementForSpeed;
+              let tos = topPosJ + 2 * curIncrementForSpeed; // maybe 1 is enough....
+              while (tos >= maxTestedPt) tos -= curIncrementForSpeed;
+              topValue = -1; // reset because increased precision may make the top lower
               for (
-                jStarFine = topPosJ - 2 * curIncrementForSpeed;
-                jStarFine < topPosJ + 2 * curIncrementForSpeed;
+                jStarFine = froms;
+                jStarFine <= tos;
                 jStarFine += curIncrementForSpeed
               ) {
+                // if (( scalProd[jStarFine] === -1) || true) { // this was a bad idea because precision reduced with earlyer tests
                 scalProd[jStarFine] = measureDeco(
-                  spe,
+                  spectrum,
                   jStarFine,
@@ -229,4 +265,4 @@ sign,
             }
+            curIncrementForSpeed = incrementForSpeed;
             // end refine
-
             if (topValue > critFoundJ) {
@@ -245,3 +281,3 @@ // ugly force value for tests
 
-              if (debug) console.log(`J:: ${topPosJ * resolutionHz}`);
+              //  if (debug) console.log(`J:: ${topPosJ * resolutionHz}`);
 
@@ -253,3 +289,6 @@ result.j.push({
               gotJValue = true;
-              if (makeShortCutForSpeed) break;
+
+              if (makeShortCutForSpeed) {
+                break;
+              }
             }
@@ -267,5 +306,5 @@ }
         appendDebug(
-          sca,
-          spe,
-          JStarArray,
+          scale,
+          spectrum,
+          jStarArray,
           scalProd,
@@ -277,3 +316,10 @@ loopoverJvalues,
       } else {
-        appendDebug(sca, spe, JStarArray, scalProd, loopoverJvalues, result);
+        appendDebug(
+          scale,
+          spectrum,
+          jStarArray,
+          scalProd,
+          loopoverJvalues,
+          result,
+        );
       }
@@ -286,4 +332,4 @@ }
       // apply here the deconvolution for the next step of the recursive process
-      spe = deco(
-        spe,
+      spectrum = deco(
+        spectrum,
         topPosJ,
@@ -297,5 +343,5 @@ sign,
         let remove = 0.5 * topPosJ * (2 * multiplicity);
-        sca = sca.slice(remove, sca.length - remove);
+        scale = scale.slice(remove, scale.length - remove);
       }
-      if (sca.length !== spe.length) {
+      if (scale.length !== spectrum.length) {
         throw Error('sts');
@@ -307,18 +353,4 @@ }
 
-  let maxAmplitudePosition = maxY({ x: sca, y: spe });
-  result.chemShift = sca[maxAmplitudePosition.index];
-
-  // for non-javascript implementation
-  /*
-  let curTop = Number.NEGATIVE_INFINITY;
-  let curChemShift;
-  for (let i = 0; i < spe.length; i++) {
-    if (spe[i] > curTop) {
-      curTop = spe[i];
-      curChemShift = sca[i];
-    }
-  }
-  result.chemShift = curChemShift;*/
-  // end to be commented
-
+  let maxAmplitudePosition = maxY({ x: scale, y: spectrum });
+  result.chemShift = scale[maxAmplitudePosition.index];
   return result;
@@ -325,0 +357,0 @@ }

src/measureDeco.js

@@ -1,2 +0,2 @@
-import { deco } from './deco';
+import { decofast1, decofast2 } from './deco';
 import { scalarProduct } from './scalarProduct';
@@ -6,3 +6,3 @@
   y, // input vector
-  JStar, // tested value of J in pt
+  jStar, // tested value of J in pt
   sign, // sign of the split 1: ++ -1: +-
@@ -13,7 +13,11 @@ chopTail, // 1: cut tail
 ) {
-  let y1 = [];
-  let y2 = [];
-  y1 = deco(y, JStar, sign, 1, chopTail, multiplicity); // dir left to right
-  y2 = deco(y, JStar, sign, -1, chopTail, multiplicity); // dir right to left
+  //let y1 = deco(y, jStar, sign, 1, chopTail, multiplicity); // dir left to right
+  //let y2 = deco(y, jStar, sign, -1, chopTail, multiplicity); // dir right to left
+
+  let nbLines = parseInt(2 * multiplicity, 10); // 1 for doublet (spin 1/2) 2, for spin 1, etc... never tested...
+  let y1 = decofast1(y, jStar, sign, nbLines, jStar);
+  // console.log(`y1 :: ` + y1);
+  let y2 = decofast2(y, jStar, sign, nbLines, jStar);
+  // console.log(`y2 :: ` + y2);
   return sign * scalarProduct(y1, y2, 1, incrementForSpeed);
 }

src/measureSymShift.js

 import { scalarProduct } from './scalarProduct';
 
 export function measureSymShift(y) {
-  let ScalarProductReference;
-  let ScalarProductNewValue;
+  let scalarProductReference;
+  let scalarProductNewValue;
   let movedBy = 0;
-  ScalarProductReference = scalarProduct(y, y, -1, 1);
+  scalarProductReference = scalarProduct(y, y, -1, 1);
   // search left...
   for (let i = 1; i < y.length / 2; i++) {
-    ScalarProductNewValue = scalarProduct(
-      y.slice(i, y.length),
-      y.slice(i, y.length),
-      -1,
-      1,
-    );
-    if (ScalarProductNewValue > ScalarProductReference) {
-      //  console.log(`${ScalarProductNewValue} > ${ScalarProductReference} size: ${y.length}`);
-      ScalarProductReference = ScalarProductNewValue;
+    scalarProductNewValue = scalarProduct(y, y, -1, 1, i, y.length);
+    if (scalarProductNewValue > scalarProductReference) {
+      scalarProductReference = scalarProductNewValue;
       movedBy = i;
-    } else {
-      // console.log('OK+ ' + movedBy + " " + i + ' ' + ScalarProductNewValue);
     }
   }
-  if (movedBy === 0) {
-    for (let i = 1; i < y.length / 2; i++) {
-      ScalarProductNewValue = scalarProduct(
-        y.slice(0, y.length - i),
-        y.slice(0, y.length - i),
-        -1,
-        1,
-      );
-      if (ScalarProductNewValue > ScalarProductReference) {
-        //  console.log(`${ScalarProductNewValue} > ${ScalarProductReference} size: ${y.length}`);
-        ScalarProductReference = ScalarProductNewValue;
-        movedBy = -i;
-      } else {
-        // console.log('OK- ' + movedBy + ' ' + i + ' ' + ScalarProductNewValue);
-      }
+  for (let i = 1; i < y.length / 2; i++) {
+    scalarProductNewValue = scalarProduct(y, y, -1, 1, 0, y.length - i);
+    if (scalarProductNewValue > scalarProductReference) {
+      scalarProductReference = scalarProductNewValue;
+      movedBy = -i;
     }
   }
-  //console.log('OK  ' + movedBy);
   return movedBy;
 }

src/scalarProduct.js

@@ -1,2 +0,9 @@
-export function scalarProduct(y1, y2, sens, incrementForSpeed) {
+export function scalarProduct(
+  y1,
+  y2,
+  sens,
+  incrementForSpeed,
+  first = 0,
+  last = y1.length,
+) {
   // sens = 1; crude scalar product
@@ -8,3 +15,3 @@ // sens =-1: flip spectrum first
   if (sens > 0) {
-    for (let index = 0; index < y1.length; index += incrementForSpeed) {
+    for (let index = first; index < last; index += incrementForSpeed) {
       v12 += y1[index] * y2[index];
@@ -16,4 +23,4 @@ v11 += y1[index] * y1[index];
     // Here as if flip left/right array
-    for (let index = 0; index < y1.length; index += incrementForSpeed) {
-      v12 += y1[index] * y2[y1.length - index - 1];
+    for (let index = first; index < last; index += incrementForSpeed) {
+      v12 += y1[index] * y2[last - (index - first) - 1];
       v11 += y1[index] * y1[index];
@@ -20,0 +27,0 @@ v22 += y2[index] * y2[index];

src/trigInterpolate.js

@@ -10,4 +10,4 @@ import { fft, ifft } from 'fft-js';
 ) {
-  let sca = Array(numberOfPointOutput);
-  let spe = Array(numberOfPointOutput);
+  let scale = new Float64Array(numberOfPointOutput);
+  let spectrum = new Float64Array(numberOfPointOutput);
 
@@ -18,9 +18,9 @@ let scaIncrement = (x[x.length - 1] - x[0]) / (x.length - 1); // delta one pt
   let trueWidth = scaIncrement * x.length;
-  scaIncrement = trueWidth / numberOfPointOutput;
+  let scaIncrementInterp = trueWidth / numberOfPointOutput;
   //scaPt += scaIncrement / 2; // move from limit side to middle of first pt
   // set scale
 
-  for (let loop = 0; loop < numberOfPointOutput; loop++) {
-    sca[loop] = scaPt;
-    scaPt += scaIncrement;
+  for (let i = 0; i < numberOfPointOutput; i++) {
+    scale[i] = scaPt;
+    scaPt += scaIncrementInterp;
   }
@@ -30,7 +30,5 @@
   // prepare input for fft apply fftshift
-  let nextPowerTwoInput = Math.pow(
-    2,
-    Math.round(Math.log(y.length - 1) / Math.log(2.0) + 0.5),
-  );
-  //nextPowerTwoAn = Math.pow(2, Math.ceil(Math.log(y.length) / Math.log(2)));
+  let nextPowerTwoInput = 2 ** Math.ceil(Math.log2(y.length));
+  let nextPowerTwoOut = 2 ** Math.ceil(Math.log2(numberOfPointOutput));
+
   let an = [...Array(nextPowerTwoInput)].map(() => Array(2).fill(0)); // n x 2 array
@@ -40,63 +38,51 @@ const halfNumPt = Math.floor(y.length / 2); // may ignore last pt... if odd number
   const shiftMultiplet = nextPowerTwoInput - y.length;
-  for (let loop = 0; loop < halfNumPt; loop++) {
-    an[shiftMultiplet + loop + halfNumPtB][0] = y[loop]; //Re
-    an[shiftMultiplet + loop + halfNumPtB][1] = 0; //Im
+  for (let i = 0; i < halfNumPt; i++) {
+    an[shiftMultiplet + i + halfNumPtB][0] = y[i]; //Re
+    an[shiftMultiplet + i + halfNumPtB][1] = 0; //Im
   }
-  for (let loop = 0; loop < halfNumPtB; loop++) {
-    an[loop][0] = y[loop + halfNumPt]; //Re
-    an[loop][1] = 0; //Im
+  for (let i = 0; i < halfNumPtB; i++) {
+    an[i][0] = y[i + halfNumPt]; //Re
+    an[i][1] = 0; //Im
   }
-
   let timeDomain = ifft(an);
-  /*an = fft(out);
-  for (let loop = 0; loop < 2 * halfNumPt; loop++) {
-    an[loop][1] = 0; //* Math.cos((phase / 180) * Math.PI) -
-  }
-  out = ifft(an);*/
-
   timeDomain[0][0] = timeDomain[0][0] / 2; // divide first point by 2 Re
   timeDomain[0][1] = timeDomain[0][1] / 2; // divide first point by 2 Im
   // move to larger array...
-  let timeDomainZeroFilled = [...Array(numberOfPointOutput)].map(() =>
+  let timeDomainZeroFilled = [...Array(nextPowerTwoOut)].map(() =>
     Array(2).fill(0),
   ); // n x 2 array
-  for (let loop = 0; loop < halfNumPt; loop++) {
-    timeDomainZeroFilled[loop][0] = timeDomain[loop][0]; //* Math.cos((phase / 180) * Math.PI) +
-    timeDomainZeroFilled[loop][1] = timeDomain[loop][1]; //* Math.cos((phase / 180) * Math.PI) -
+  for (let i = 0; i < halfNumPt; i++) {
+    timeDomainZeroFilled[i][0] = timeDomain[i][0]; //* Math.cos((phase / 180) * Math.PI) +
+    timeDomainZeroFilled[i][1] = timeDomain[i][1]; //* Math.cos((phase / 180) * Math.PI) -
   }
-
-  for (let loop = halfNumPt; loop < numberOfPointOutput; loop++) {
+  for (let i = halfNumPt; i < nextPowerTwoInput; i++) {
     // zero filling
-    timeDomainZeroFilled[loop][0] = 0;
-    timeDomainZeroFilled[loop][1] = 0;
+    timeDomainZeroFilled[i][0] = 0;
+    timeDomainZeroFilled[i][1] = 0;
   }
+
   let interpolatedSpectrum = fft(timeDomainZeroFilled);
   const halfNumPt2 = Math.floor(numberOfPointOutput / 2);
-  let tmp;
+
   // applies phase change
   let phaseRad = ((addPhaseInterpolation + 0.0) / 180.0) * Math.PI; // this is for testing additional phases
   if (phaseRad !== 0.0) {
-    for (let loop = 0; loop < 2 * halfNumPt2; loop++) {
-      tmp =
-        interpolatedSpectrum[loop][0] * Math.cos(phaseRad) +
-        interpolatedSpectrum[loop][1] * Math.sin(phaseRad); // only Re now...
-      interpolatedSpectrum[loop][1] =
-        -interpolatedSpectrum[loop][0] * Math.sin(phaseRad) +
-        interpolatedSpectrum[loop][1] * Math.cos(phaseRad); // only Re now...
-      interpolatedSpectrum[loop][0] = tmp;
+    for (let i = 0; i < 2 * halfNumPt2; i++) {
+      let tmp =
+        interpolatedSpectrum[i][0] * Math.cos(phaseRad) +
+        interpolatedSpectrum[i][1] * Math.sin(phaseRad); // only Re now...
+      interpolatedSpectrum[i][1] =
+        -interpolatedSpectrum[i][0] * Math.sin(phaseRad) +
+        interpolatedSpectrum[i][1] * Math.cos(phaseRad); // only Re now...
+      interpolatedSpectrum[i][0] = tmp;
     }
   }
+
   let returnedPhase = 0;
-
+  let norm;
   if (appliedPhaseCorrectionType > 0) {
-    let localPhaseRad = 0;
-    let vectx;
-    let vecty;
-    let norm;
-    // let sumNorms;
-
-    for (let loo = 1; loo < 100; loo++) {
-      localPhaseRad = 0;
-      vectx = 0;
-      vecty = 0;
+    for (let i = 1; i < 100; i++) {
+      let localPhaseRad = 0;
+      let vectx = 0;
+      let vecty = 0;
       // sumNorms = 0;
@@ -129,12 +115,9 @@ if (appliedPhaseCorrectionType > 0) {
       norm = Math.sqrt(vecty * vecty + vectx * vectx);
-      phaseRad = -(10.0 / (loo * loo)) * Math.sign(vecty);
+      phaseRad = -(10.0 / (i * i)) * Math.sign(vecty);
 
       returnedPhase -= (180.0 * phaseRad) / Math.PI;
 
-      // console.log('localPhase ' + (180.0 * localPhaseRad) / Math.PI);
-      //console.log('returnedPhase ' + returnedPhase);
-
       if (phaseRad !== 0.0) {
         for (let loop = 0; loop < 2 * halfNumPt2; loop++) {
-          tmp =
+          let tmp =
             interpolatedSpectrum[loop][0] * Math.cos(phaseRad) +
@@ -152,7 +135,8 @@ interpolatedSpectrum[loop][1] * Math.sin(phaseRad); // only Re now...
   // applies fftshift
-  for (let loop = 0; loop < halfNumPt2; loop++) {
-    spe[loop] = interpolatedSpectrum[loop + halfNumPt2][0]; // only Re now...
+  let dropPoints = nextPowerTwoOut - numberOfPointOutput;
+  for (let i = 0; i < halfNumPt2; i++) {
+    spectrum[i] = interpolatedSpectrum[halfNumPt2 + dropPoints + i][0]; // only Re now...
   }
-  for (let loop = 0; loop < halfNumPt2; loop++) {
-    spe[loop + halfNumPt2] = interpolatedSpectrum[loop][0];
+  for (let i = 0; i < halfNumPt2; i++) {
+    spectrum[i + halfNumPt2] = interpolatedSpectrum[i][0];
   }
@@ -162,7 +146,8 @@ if (returnedPhase > 360.0) {
   }
+
   return {
-    spectrum: spe,
-    scale: sca,
+    spectrum,
+    scale,
     phaseCorrectionOnMultipletInDeg: returnedPhase,
   };
 }

New alerts

License Policy Violation

License

This package is not allowed per your license policy. Review the package's license to ensure compliance.

Found 1 instance in 1 package

Fixed alerts

License Policy Violation

License

This package is not allowed per your license policy. Review the package's license to ensure compliance.

Found 1 instance in 1 package

Improved metrics

Total package byte prevSize

66390

increased by10.75%

Lines of code

1657

increased by10.69%

Number of lines in readme file

108

increased by2.86%

Worsened metrics

Number of package files

14

decreased by-6.67%

Dependency changes

• + Addedfft.js@^4.0.3

• + Addedfft.js@4.0.4(transitive)

• - Removedml-fft@^1.3.5

• - Removedml-fft@1.3.5(transitive)