pysmo.tools.noise

This module provides support for calculating random synthetic noise that matches the naturally observed amplitude spectrum.

See also

Peterson, J., 1993. Observations and modelling of background seismic noise. Open-file report 93-322, U. S. Geological Survey, Albuquerque, New Mexico.

Classes:

Name	Description
NoiseModel	Class to store seismic noise models.

Functions:

Name	Description
peterson	Generate a noise model using Peterson's models as base.
generate_noise	Generate noise from a noise model.

Examples:

>>> from pysmo.tools.noise import generate_noise, peterson
>>> NLNM = peterson(noise_level=0)
>>> delta = 0.05
>>> npts = 5000
>>> low_noise_seismogram = generate_noise(NLMN, npts, delta)

NoiseModel dataclass

Class to store seismic noise models.

Attributes:

Name	Type	Description
psd	ndarray	Power spectral density of ground acceleration [dB].
T	ndarray	Period [seconds].

Source code in pysmo/tools/noise.py

@dataclass(frozen=True)
class NoiseModel:
    """Class to store seismic noise models.

    Attributes:
        psd: Power spectral density of ground acceleration [dB].
        T: Period [seconds].
    """

    psd: np.ndarray = field(default_factory=lambda: np.array([]))
    T: np.ndarray = field(default_factory=lambda: np.array([]))

    def __post_init__(self) -> None:
        if np.size(self.psd) != np.size(self.T):
            raise ValueError(
                "dB and T arrays are not of same size",
                f"({np.size(self.psd)} != {np.size(self.T)}",
            )
        self.psd.flags.writeable = False
        self.T.flags.writeable = False

generate_noise(model, npts, delta=SEISMOGRAM_DEFAULTS.delta, begin_time=SEISMOGRAM_DEFAULTS.begin_time, return_velocity=False, seed=None)

Generate a random seismogram from a noise model. Realistic seismic noise is generated by assigning random phases to a fixed amplitude spectrum prescribed by a noise model.

Parameters:

Name	Type	Description	Default
model	NoiseModel	Noise model used to compute seismic noise.	required
npts	int	Number of points of generated noise	required
delta	float	Sampling interval of generated noise	delta
begin_time	datetime	Seismogram begin time	begin_time
return_velocity	bool	Return velocity instead of acceleration.	False
seed	int | None	set random seed manually (usually for testing purposes).	None

Returns:

Type	Description
MiniSeismogram	Seismogram with random seismic noise as data.

Source code in pysmo/tools/noise.py

def generate_noise(
    model: NoiseModel,
    npts: int,
    delta: float = SEISMOGRAM_DEFAULTS.delta,
    begin_time: datetime = SEISMOGRAM_DEFAULTS.begin_time,
    return_velocity: bool = False,
    seed: int | None = None,
) -> MiniSeismogram:
    """Generate a random seismogram from a noise model. Realistic seismic noise is
    generated by assigning random phases to a fixed amplitude spectrum prescribed by a
    noise model.

    Parameters:
        model: Noise model used to compute seismic noise.
        npts: Number of points of generated noise
        delta: Sampling interval  of generated noise
        begin_time: Seismogram begin time
        return_velocity: Return velocity instead of acceleration.
        seed: set random seed manually (usually for testing purposes).

    Returns:
        Seismogram with random seismic noise as data.
    """
    # Sampling frequency
    Fs = 1 / delta

    # Nyquist frequency
    Fnyq = 0.5 / delta

    # get next power of 2 of the nunmber of points and calculate frequencies from
    # Fs/NPTS to Fnyq (we skip a frequency of 0 for now to avoid dividing by 0)
    NPTS = int(2 ** np.ceil(np.log2(npts)))
    freqs = np.linspace(Fs / NPTS, Fnyq, NPTS - 1)

    # interpolate psd and recreate amplitude spectrum with the first
    # term=0 (i.e. mean=0).
    Pxx = np.interp(1 / freqs, model.T, model.psd)
    spectrum = np.append(np.array([0]), np.sqrt(10 ** (Pxx / 10) * NPTS / delta * 2))

    # phase is randomly generated
    rng = np.random.default_rng(seed)
    phase = (rng.random(NPTS) - 0.5) * np.pi * 2

    # combine amplitude with phase and perform ifft
    NewX = spectrum * (np.cos(phase) + 1j * np.sin(phase))
    acceleration = np.fft.irfft(NewX)

    start = int((NPTS - npts) / 2)
    end = start + npts
    if return_velocity:
        velocity = cumtrapz(acceleration, dx=delta)
        velocity = velocity[start:end]
        return MiniSeismogram(begin_time=begin_time, delta=delta, data=velocity)
    acceleration = acceleration[start:end]
    return MiniSeismogram(begin_time=begin_time, delta=delta, data=acceleration)

peterson(noise_level)

Generate a noise model by interpolating between Peterson's New Low Noise Model (NLNM) and New High Noice Model (NHNM).

Parameters:

Name	Type	Description	Default
noise_level	float	Determines the noise level of the generated noise model. A noise level of 0 returns the NLNM, 1 returns the NHNM, and anything > 0 and < 1 returns an interpolated model that lies between the NLNM and NHNM.	required

Returns:

Type	Description
NoiseModel	Noise model.

Source code in pysmo/tools/noise.py

def peterson(noise_level: float) -> NoiseModel:
    """Generate a noise model by interpolating between Peterson's
    New Low Noise Model (NLNM) and New High Noice Model (NHNM).

    Parameters:
        noise_level: Determines the noise level of the generated noise model.
                     A noise level of 0 returns the NLNM, 1 returns the NHNM,
                     and anything > 0 and < 1 returns an interpolated model
                     that lies between the NLNM and NHNM.

    Returns:
        Noise model.
    """
    # check for valid input
    if not 0 <= noise_level <= 1:
        raise ValueError(
            f"Parameter noise_level={noise_level} is not within 0-1 range."
        )

    # calculate noise model
    if noise_level == 0:
        return NLNM
    elif noise_level == 1:
        return NHNM
    else:
        T_common = np.unique(np.concatenate((NLNM.T, NHNM.T)))
        NLNM_interp = np.interp(T_common, NLNM.T, NLNM.psd)
        NHNM_interp = np.interp(T_common, NHNM.T, NHNM.psd)
        dB = NLNM_interp + (NHNM_interp - NLNM_interp) * noise_level
        return NoiseModel(psd=dB, T=T_common)