From 9b331c1b88b2486e939a01f1e5b35094d82b49f4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= <benoit.sierro@bluewin.ch>
Date: Tue, 29 Aug 2023 08:53:40 +0200
Subject: [PATCH] change: moved dB to math

---
 examples/chang2011.py            |  2 +-
 src/scgenerator/math.py          | 30 ++++++++++++++++++++++
 src/scgenerator/noise.py         | 12 ++++-----
 src/scgenerator/physics/units.py | 44 --------------------------------
 4 files changed, 36 insertions(+), 52 deletions(-)

diff --git a/examples/chang2011.py b/examples/chang2011.py
index f4c6cb1..17dda45 100644
--- a/examples/chang2011.py
+++ b/examples/chang2011.py
@@ -39,7 +39,7 @@ quit()
 interp = interp1d(stats["z"], solution, axis=0)
 z = np.linspace(0, params.length, 128)
 plt.imshow(
-    sc.units.to_log(interp(z)),
+    sc.math.to_dB(interp(z)),
     vmin=-50,
     extent=sc.get_extent(sc.units.THz_inv(params.w), z),
     origin="lower",
diff --git a/src/scgenerator/math.py b/src/scgenerator/math.py
index e486f2a..52a949d 100644
--- a/src/scgenerator/math.py
+++ b/src/scgenerator/math.py
@@ -114,6 +114,36 @@ def sigmoid(x):
     return 1 / (np.exp(-x) + 1)
 
 
+def to_dB(arr: np.ndarray, ref=None, axis=None) -> np.ndarray:
+    """
+    converts unitless values in dB
+
+    Parameters
+    ----------
+    arr : np.ndarray
+        array of non-negative values. Any values 0 or below will be mapped to the minimum
+        positive value in dB
+    ref : float, optional
+        reference value corresponding to 0dB (default : max(arr))
+    axis : int | None, optional
+        on which axis to apply the transformation. If `ref` is given, this has no effect
+
+    Returns
+    ----------
+    np.ndarray
+        array in dB
+    """
+    if axis is not None and arr.ndim > 1 and ref is None:
+        return np.apply_along_axis(to_dB, axis, arr)
+
+    if ref is None:
+        ref = np.max(arr)
+    m = arr / ref
+    above_0 = m > 0
+    m = 10 * np.log10(m, out=np.zeros_like(m) - 10 * np.log10(m[above_0].min()), where=above_0)
+    return m
+
+
 def u_nm(n, m):
     """
     returns the mth zero of the Bessel function of order n-1
diff --git a/src/scgenerator/noise.py b/src/scgenerator/noise.py
index b606101..9c06c8c 100644
--- a/src/scgenerator/noise.py
+++ b/src/scgenerator/noise.py
@@ -16,7 +16,6 @@ class NoiseMeasurement:
     psd: np.ndarray
     phase: np.ndarray | None = None
     psd_interp: interp1d = field(init=False)
-    is_uniform: bool = field(default=False, init=False)
     _window_functions: ClassVar[dict[str, Callable[[int], np.ndarray]]] = {}
 
     @classmethod
@@ -78,16 +77,13 @@ class NoiseMeasurement:
         return cls(freq, psd.mean(axis=0) / window_correction, phase=phase)
 
     def __post_init__(self):
-        df = np.diff(self.freq)
-        if df.std() / df.mean() < 1e-12:
-            self.is_uniform = True
         self.psd_interp = interp1d(
             self.freq, self.psd, fill_value=(0, self.psd[-1]), bounds_error=False
         )
 
     @property
     def psd_dBc(self) -> np.ndarray:
-        return np.log10(self.psd) * 10
+        return math.to_dB(self.psd)
 
     @property
     def plottable_linear(self) -> tuple[np.ndarray, np.ndarray]:
@@ -127,7 +123,9 @@ class NoiseMeasurement:
             )
 
         f = np.linspace(0, fmax, nt // 2 + 1)
-        return f, self.psd_interp(f)
+        interp = self.psd_interp(f)
+        interp[0] = 0
+        return f, interp
 
     def time_series(
         self, nt: int | np.ndarray, phase: np.ndarray | None = None, dt: float | None = None
@@ -152,7 +150,7 @@ class NoiseMeasurement:
         time, dt = irfftfreq(freq, True)
 
         amp = np.sqrt(spec) * np.exp(1j * phase)
-        signal = np.fft.irfft(amp) * np.sqrt(0.5 * len(time) / dt)
+        signal = np.fft.irfft(amp) * np.sqrt(0.5 * nt / dt)
 
         return time, signal
 
diff --git a/src/scgenerator/physics/units.py b/src/scgenerator/physics/units.py
index dcf79c7..ffe37c0 100644
--- a/src/scgenerator/physics/units.py
+++ b/src/scgenerator/physics/units.py
@@ -361,50 +361,6 @@ def to_WL(spectrum: np.ndarray, lambda_: np.ndarray) -> np.ndarray:
     return m
 
 
-def to_log(arr, ref=None):
-    """
-    takes the log of each 1D array relative to the max of said array. Useful
-    to plot spectrum evolution, but use to_log2D for spectrograms
-
-    Parameters
-    ----------
-        arr : 1D array of real numbers. >1D array : operation is applied on axis=0
-        ref : reference value corresponding to 0dB (default : max(arr))
-
-    Returns
-    ----------
-        arr array in dB
-    """
-    if arr.ndim > 1:
-        return np.apply_along_axis(to_log, -1, arr, ref)
-    else:
-        if ref is None:
-            ref = np.max(arr)
-        m = arr / ref
-        m = 10 * np.log10(m, out=np.zeros_like(m) - 100, where=m > 0)
-        return m
-
-
-def to_log2D(arr, ref=None):
-    """
-    computes the log of a 2D array
-
-    Parameters
-    ----------
-        arr : 2D array of real numbers
-        ref : reference value corresponding to 0dB
-
-    Returns
-    ----------
-        arr array in dB
-    """
-    if ref is None:
-        ref = np.max(arr)
-    m = arr / ref
-    m = 10 * np.log10(m, out=np.zeros_like(m) - 100, where=m > 0)
-    return m
-
-
 class PlotRange(tuple):
     left: float = property(itemgetter(0))
     right: float = property(itemgetter(1))