possible improvement to check in noise sandlot calculation

[jobo322]

I did compare the paper of sanplot with the current implementation by AI and there are some differences and possibilities to improve the result. It is something to check in the near future and only open this ticket to keep the result of the revision.

Verdict

The overall approach is appropriate for a global 1D NMR noise estimate:

1. Copy and sort the intensities by decreasing magnitude.
2. Separate positive and negative distributions.
3. Estimate the noise from order statistics.
4. Refine the estimate by excluding points above factorStd × noise.
5. Calculate an SNR and generate SAN-plot data.

This matches the SAN-plot concept in the paper for centered, real-valued spectra with approximately white Gaussian noise. However, the current implementation has several important issues. I would not yet consider the returned value quantitatively reliable, particularly in magnitudeMode.

Critical issue: Rayleigh/magnitude mode is using the wrong quantile

For real spectra, the current calculation is mathematically consistent with the positive half of a Gaussian distribution:

-simpleNormInvRaw(cutOffDist / 2)

This is because a fraction $q$ of the positive points corresponds approximately to a tail probability of $q/2$ in the complete Gaussian distribution.

For magnitude spectra, all values are positive and follow a Rayleigh distribution. If q is the fraction of points above a value in the descending SAN series, the correct Rayleigh relationship is:

$$ P(R \ge r) = q $$

and therefore:

$$ r = \sigma_R \sqrt{-2\ln(q)} $$

The current code uses:

-simpleNormInvMagnitude(cutOffDist / 2)

The utility implements the Rayleigh inverse CDF as -sqrt(-2 * log(1 - p)) in simpleNormInvMagnitude.ts. Consequently, the code currently evaluates the Rayleigh quantile at $q/2$, whereas it should use the CDF probability $1-q$.

For example, at $q=0.5$:

• expected Rayleigh factor: $\sqrt{-2\ln(0.5)} \approx 1.177$
• current factor: $\sqrt{-2\ln(0.75)} \approx 0.758$

This causes the magnitude noise estimate to be substantially overestimated and also affects cutoff selection in determineCutOff().

A shared helper would make the distinction explicit:

function getQuantileFactor(q: number, magnitudeMode: boolean): number {
  return magnitudeMode
    ? -simpleNormInvMagnitude(1 - q)
    : -simpleNormInvRaw(q / 2);
}

Both calculateNoiseLevel() and determineCutOff() should use this helper.

SNR does not follow the paper's definition

The paper defines:

$$ \mathrm{SNR} = \frac{S_{\max}}{2\sigma_N} $$

The function returns:

snr: skyPoint / noiseLevelPositive

in xNoiseSanPlot.ts. Since the raw Gaussian correction produces an estimate of $\sigma_N$, the paper-compatible result should be:

snr: skyPoint / (2 * noiseLevelPositive)

The current result is therefore approximately twice the SNR defined in Equation 1 of the paper. This may be intentional for compatibility with existing behavior, but it should either be changed or clearly documented.

In magnitude mode, the result should not be described as the conventional Gaussian NMR SNR. The paper explicitly notes that the Rayleigh noise parameter has a different meaning.

Options are currently ignored

Several public options are declared but never passed through:

• considerList is not passed to determineCutOff() at xNoiseSanPlot.ts.
• logBaseY is not passed to generateSanPlot() at xNoiseSanPlot.ts.
• fromTo is not passed to generateSanPlot().

Therefore, user-supplied values for these options have no effect.

The cutoff call should include:

determineCutOff(signPositive, {
  magnitudeMode,
  considerList: options.considerList,
});

The SAN plot call should also forward the plotting options.

cutOffFromPositive: false does not do what its documentation says

The documentation says that setting cutOffFromPositive to false calculates the cutoff from the negative distribution. However, the current code sets cutOff to null:

const cutOff =
  options.cutOff ??
  (cutOffFromPositive
    ? determineCutOff(signPositive, { magnitudeMode })
    : null);

Then calculateNoiseLevel() independently determines a cutoff for each sign at xNoiseSanPlot.ts.

If the intention is to calculate one common cutoff from the selected sign, it should be:

const cutoffSign = cutOffFromPositive ? signPositive : signNegative;

const cutOff =
  options.cutOff ??
  determineCutOff(cutoffSign, {
    magnitudeMode,
    considerList: options.considerList,
  });

For magnitude data, signNegative is empty, so cutOffFromPositive: false should probably be rejected or automatically treated as positive mode.

Other correctness issues

Scaling factors below one are ignored

At xNoiseSanPlot.ts, scaling is only applied when:

if (scaleFactor > 1)

The documented behavior says that the input is multiplied by scaleFactor. Therefore 0.5 currently has no effect. It should generally be:

if (scaleFactor !== 1) {
  for (let i = 0; i < input.length; i++) {
    input[i] *= scaleFactor;
  }
}

Empty and degenerate inputs are not handled

For an empty array, skyPoint can become undefined, and the returned SNR becomes NaN. Constant-zero data can produce 0 / 0.

The function should explicitly handle:

• empty input;
• no positive values;
• zero noise;
• fewer points than required for the quantile procedure.

Positive SAN range has an off-by-one issue

The positive range is passed as:

positive: { from: 0, to: lastPositiveValueIndex }

but array.slice(from, to) excludes to, so the last positive point is omitted. The range should probably use:

positive: { from: 0, to: lastPositiveValueIndex + 1 }

Cutoff selection is heuristic, not a direct fit

determineCutOff() chooses the window with the smallest unnormalized variance of local sigma estimates at xNoiseSanPlot.ts. This is a reasonable heuristic, but it can select a flat signal or artifact region instead of the true noise region. A flat region is not necessarily a noise region.

The paper describes matching the SAN plot to the theoretical distribution. The current implementation approximates this by selecting a locally stable region, so it should be validated with simulated data containing:

• pure Gaussian noise with known $\sigma$;
• Gaussian noise plus positive peaks;
• asymmetric positive and negative artifacts;
• baseline offsets;
• different signal densities;
• magnitude data with known Rayleigh scale.

The current test mainly checks one fixture, scaling invariance, masking, and relative positive/negative values in xNoiseSanPlot.test.ts. It does not verify recovery of a known noise standard deviation.

Practical recommendation for 1D NMR

For a conventional phased real spectrum:

xNoiseSanPlot(realSpectrum, {
  magnitudeMode: false,
  fixOffset: true,
});

is conceptually suitable if:

• the spectrum is reasonably baseline corrected;
• the noise is approximately stationary and Gaussian;
• the number of signal points is small relative to the total spectrum;
• strong artifacts or broad unresolved signals are masked.

For magnitude spectra, magnitudeMode: true is conceptually appropriate, but the Rayleigh quantile correction should be fixed before relying on its noise value.

The most important changes are:

1. Correct the Rayleigh quantile calculation.
2. Decide whether snr should follow the paper's factor of $2$.
3. Forward considerList, logBaseY, and fromTo.
4. Add validation for empty, zero-noise, and magnitude-only cases.
5. Add pure Gaussian and Rayleigh simulations with known noise levels.