dedebenui/scgenerator
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

change: noise improvements

Browse Source
This commit is contained in:
Benoît Sierro
2023-08-28 11:40:02 +02:00
parent 4ea23bedda
commit b89a76db16
3 changed files with 69 additions and 58 deletions
Download Patch File Download Diff File Expand all files Collapse all files

87
src/scgenerator/noise.py
View File

@@ -1,7 +1,7 @@
from __future__ import annotations
from dataclasses import dataclass, field
from typing import Callable, ClassVar
from typing import Callable, ClassVar, Sequence
import numpy as np
from scipy.integrate import cumulative_trapezoid
@@ -10,35 +10,6 @@ from scipy.interpolate import interp1d
from scgenerator import math
def irfftfreq(freq: np.ndarray, retstep: bool = False):
"""
Given an array of positive only frequency, this returns the corresponding time array centered
around 0 that will be aligned with the `numpy.fft.irfft` of a spectrum aligned with `freq`.
if `retstep` is True, the sample spacing is returned as well
"""
df = freq[1] - freq[0]
nt = (len(freq) - 1) * 2
period = 1 / df
dt = period / nt
t = np.linspace(-(period - dt) / 2, (period - dt) / 2, nt)
if retstep:
return t, dt
else:
return t
def log_power(x):
return 10 * np.log10(np.abs(np.where(x == 0, 1e-7, x)))
def integrated_rin(freq: np.ndarray, psd: np.ndarray) -> float:
"""
given a normalized spectrum, computes the total rms RIN in the provided frequency window
"""
return np.sqrt(cumulative_trapezoid(np.abs(psd)[::-1], -freq[::-1], initial=0)[::-1])
@dataclass
class NoiseMeasurement:
freq: np.ndarray
@@ -73,7 +44,7 @@ class NoiseMeasurement:
@classmethod
def from_time_series(
cls, time: np.ndarray, signal: np.ndarray, window: str = "Hann", num_segments: int = 1
cls, signal: Sequence[float], dt: float, window: str = "Hann", num_segments: int = 1
) -> NoiseMeasurement:
"""
compute a PSD from a time-series measurement.
@@ -93,15 +64,15 @@ class NoiseMeasurement:
number of segments to cut the signal in. This will trade lower frequency information for
better variance of the estimated PSD. The default 1 means no cutting.
"""
signal = np.asanyarray(signal)
signal_segments = segments(signal, num_segments)
n = signal_segments.shape[-1]
window_arr = cls._window_functions[window](n)
window_correction = np.sum(window_arr**2) / n
signal_segments = signal_segments * window_arr
freq = np.fft.rfftfreq(n, time[1] - time[0])
dt = time[1] - time[0]
psd = np.fft.rfft(signal_segments) / np.sqrt(0.5 * n / dt)
freq = np.fft.rfftfreq(n, dt)
psd = np.fft.rfft(signal_segments) * np.sqrt(2 * dt / n)
phase = math.mean_angle(psd)
psd = psd.real**2 + psd.imag**2
return cls(freq, psd.mean(axis=0) / window_correction, phase=phase)
@@ -118,9 +89,19 @@ class NoiseMeasurement:
def psd_dBc(self) -> np.ndarray:
return np.log10(self.psd) * 10
@property
def plottable_linear(self) -> tuple[np.ndarray, np.ndarray]:
if self.freq[0] == 0:
return self.freq[1:], self.psd[1:]
else:
return self.freq, self.psd
@property
def plottable_dBc(self) -> tuple[np.ndarray, np.ndarray]:
return self.freq[1:], self.psd_dBc[1:]
if self.freq[0] == 0:
return self.freq[1:], self.psd_dBc[1:]
else:
return self.freq, self.psd_dBc
def sample_spectrum(self, nt: int, dt: float | None = None) -> tuple[np.ndarray, np.ndarray]:
"""
@@ -175,12 +156,42 @@ class NoiseMeasurement:
return time, signal
def integrated_rin(self) -> np.ndarray:
def root_integrated(self) -> np.ndarray:
"""
returns the integrated RIN as fuction of frequency.
returns the sqrt of the integrated spectrum
The 0th component is the total RIN in the frequency range covered by the measurement
Caution! this may be 0 frequency
"""
return integrated_rin(self.freq, self.psd)
return integrated_noise(self.freq, self.psd)
def irfftfreq(freq: np.ndarray, retstep: bool = False):
"""
Given an array of positive only frequency, this returns the corresponding time array centered
around 0 that will be aligned with the `numpy.fft.irfft` of a spectrum aligned with `freq`.
if `retstep` is True, the sample spacing is returned as well
"""
df = freq[1] - freq[0]
nt = (len(freq) - 1) * 2
period = 1 / df
dt = period / nt
t = np.linspace(-(period - dt) / 2, (period - dt) / 2, nt)
if retstep:
return t, dt
else:
return t
def log_power(x):
return 10 * np.log10(np.abs(np.where(x == 0, 1e-7, x)))
def integrated_noise(freq: np.ndarray, psd: np.ndarray) -> float:
"""
given a normalized spectrum, computes the total rms RIN in the provided frequency window
"""
return np.sqrt(cumulative_trapezoid(np.abs(psd)[::-1], -freq[::-1], initial=0)[::-1])
@NoiseMeasurement.window_function("Square")

37
src/scgenerator/physics/pulse.py
View File

@@ -9,7 +9,6 @@ of shape `([what, ever,] n, nt)` (could be just `(n, nt)` for 2D sc-ordered arra
n is the number of spectra at the same z position and nt is the size of the time/frequency grid
"""
import itertools
import os
from dataclasses import astuple, dataclass
from io import BytesIO
@@ -30,8 +29,6 @@ from scgenerator.defaults import default_plotting
from scgenerator.io import DataFile
from scgenerator.physics import units
c = 299792458.0
hbar = 1.05457148e-34
T = TypeVar("T")
#
@@ -391,6 +388,10 @@ def soliton_length(dispersion_length: float):
return pi / 2 * dispersion_length
def center_of_gravitiy(t: np.ndarray, intensity: np.ndarray):
return np.sum(intensity * t) / np.sum(intensity)
def adjust_custom_field(
input_time: np.ndarray,
input_field: np.ndarray,
@@ -402,8 +403,8 @@ def adjust_custom_field(
peak_power: float = None,
) -> np.ndarray:
"""
When loading a custom input field, use this function to conform the custom data to the simulation
parameters.
When loading a custom input field, use this function to conform the custom data to the
simulation parameters.
Parameters
----------
@@ -534,43 +535,41 @@ def technical_noise(rms_noise: float, noise_correlation: float = -0.4):
return psy, 1 + noise_correlation * (psy - 1)
def shot_noise(w: np.ndarray, T: float, dt: float, additional_noise_factor=1.0):
def shot_noise(w: np.ndarray, dt: float):
"""
Parameters
----------
w_c : array, shape (n,)
angular frequencies centered around 0
T : float
length of the time windows
w : array, shape (n,)
angular frequencies
dt : float
resolution of time grid
additional_noise_factor : float
resulting noise spectrum is multiplied by this number
Returns
-------
out : array, shape (n,)
np.ndarray, shape (n,)
noise field to be added on top of initial field in time domain
"""
dw = w[1] - w[0]
rand_phase = np.random.rand(len(w)) * 2 * pi
A_oppm = np.sqrt(hbar * (np.abs(w)) * T) * np.exp(-1j * rand_phase)
out = ifft(A_oppm / dt * np.sqrt(2 * pi))
return out * additional_noise_factor
oppm_phase = np.exp(-1j * rand_phase)
oppm_amp = np.sqrt(0.5 * units.hbar * np.abs(w) / dw)
return ifft(oppm_amp * oppm_phase / dt * np.sqrt(2 * pi))
def finalize_pulse(
pre_field_0: np.ndarray,
quantum_noise: bool,
w: np.ndarray,
time_window: float,
dt: float,
additional_noise_factor: float,
input_transmission: float,
intensity_noise: float | None,
) -> np.ndarray:
if quantum_noise:
pre_field_0 = pre_field_0 + shot_noise(w, time_window, dt, additional_noise_factor)
pre_field_0 = pre_field_0 + shot_noise(w, dt)
ratio = 1
if intensity_noise is not None:

3
src/scgenerator/physics/units.py
View File

@@ -11,7 +11,8 @@ import numpy as np
from numpy import pi
c = 299792458.0
hbar = 1.05457148e-34
h = 6.62607015e-34
hbar = 0.5 * h / np.pi
NA = 6.02214076e23
R = 8.31446261815324
kB = 1.380649e-23