diff --git a/src/scgenerator/noise.py b/src/scgenerator/noise.py index c9762e1..b606101 100644 --- a/src/scgenerator/noise.py +++ b/src/scgenerator/noise.py @@ -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") diff --git a/src/scgenerator/physics/pulse.py b/src/scgenerator/physics/pulse.py index 7ca6b51..c68ef62 100644 --- a/src/scgenerator/physics/pulse.py +++ b/src/scgenerator/physics/pulse.py @@ -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: diff --git a/src/scgenerator/physics/units.py b/src/scgenerator/physics/units.py index c4424ac..dcf79c7 100644 --- a/src/scgenerator/physics/units.py +++ b/src/scgenerator/physics/units.py @@ -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