Module ocean_science_utilities.wavespectra.timeseries

Expand source code
import numpy as np
from typing import Literal, Optional, Tuple, Union

from ocean_science_utilities.wavespectra.spectrum import (
    FrequencySpectrum,
    FrequencyDirectionSpectrum,
    WaveSpectrum,
)


def surface_timeseries(
    component: Literal["u", "v", "w", "x", "y", "z"],
    sampling_frequency: float,
    signal_length: int,
    spectrum: WaveSpectrum,
    seed: Optional[int] = None,
) -> Tuple[np.typing.NDArray, np.typing.NDArray]:
    """
    Create a timeseries for from a given power spectral density.

    :param component: Wave component to create a timeseries for: u,v,w,x,y,z.
    :param sampling_frequency: Sampling frequency of output signal in Hertz
    :param signal_length: Length of output signal
    :param spectrum: Input power spectrum
    :param seed: Input seed for the random number generator.
    :return:
    """

    nfft = (int(signal_length) // 2) * 2

    frequencies = np.linspace(0, 0.5 * sampling_frequency, nfft // 2, endpoint=False)

    time = np.linspace(0, nfft / sampling_frequency, nfft, endpoint=False)

    timeseries = nfft * np.fft.irfft(
        create_fourier_amplitudes(component, spectrum, frequencies, seed)
    )

    return time, timeseries


def create_fourier_amplitudes(
    component, spectrum: WaveSpectrum, frequencies, seed=None
):
    spectrum = spectrum.interpolate_frequency(frequencies)

    if isinstance(spectrum, FrequencySpectrum):
        radian_directions: Union[float, np.ndarray] = 0.0
        radian_frequency = spectrum.radian_frequency.values
        area = spectrum.frequency_step.values

    elif isinstance(spectrum, FrequencyDirectionSpectrum):
        radian_directions = spectrum.radian_direction.values[None, :]
        radian_frequency = spectrum.radian_frequency.values[:, None]
        area = (
            spectrum.frequency_step.values[:, None]
            * spectrum.direction_step.values[None, :]
        )
    else:
        raise ValueError("Not a spectrum")

    if component == "w":
        factor: Union[float, np.ndarray] = 1j * radian_frequency

    elif component == "u":
        factor = radian_frequency * np.cos(radian_directions)

    elif component == "v":
        factor = radian_frequency * np.sin(radian_directions)

    elif component == "x":
        factor = -1j * np.cos(radian_directions)

    elif component == "y":
        factor = -1j * np.sin(radian_directions)

    else:
        factor = 1.0

    phases = np.random.default_rng(seed=seed).uniform(
        0, 2 * np.pi, spectrum.spectral_shape()
    )
    amplitudes = (
        np.sqrt(area * spectrum.variance_density / 2) * np.exp(1j * phases) * factor
    )

    if isinstance(spectrum, FrequencyDirectionSpectrum):
        amplitudes = np.sum(amplitudes, axis=-1)

    return amplitudes

Functions

def create_fourier_amplitudes(component, spectrum: WaveSpectrum, frequencies, seed=None)
Expand source code
def create_fourier_amplitudes(
    component, spectrum: WaveSpectrum, frequencies, seed=None
):
    spectrum = spectrum.interpolate_frequency(frequencies)

    if isinstance(spectrum, FrequencySpectrum):
        radian_directions: Union[float, np.ndarray] = 0.0
        radian_frequency = spectrum.radian_frequency.values
        area = spectrum.frequency_step.values

    elif isinstance(spectrum, FrequencyDirectionSpectrum):
        radian_directions = spectrum.radian_direction.values[None, :]
        radian_frequency = spectrum.radian_frequency.values[:, None]
        area = (
            spectrum.frequency_step.values[:, None]
            * spectrum.direction_step.values[None, :]
        )
    else:
        raise ValueError("Not a spectrum")

    if component == "w":
        factor: Union[float, np.ndarray] = 1j * radian_frequency

    elif component == "u":
        factor = radian_frequency * np.cos(radian_directions)

    elif component == "v":
        factor = radian_frequency * np.sin(radian_directions)

    elif component == "x":
        factor = -1j * np.cos(radian_directions)

    elif component == "y":
        factor = -1j * np.sin(radian_directions)

    else:
        factor = 1.0

    phases = np.random.default_rng(seed=seed).uniform(
        0, 2 * np.pi, spectrum.spectral_shape()
    )
    amplitudes = (
        np.sqrt(area * spectrum.variance_density / 2) * np.exp(1j * phases) * factor
    )

    if isinstance(spectrum, FrequencyDirectionSpectrum):
        amplitudes = np.sum(amplitudes, axis=-1)

    return amplitudes
def surface_timeseries(component: Literal['u', 'v', 'w', 'x', 'y', 'z'], sampling_frequency: float, signal_length: int, spectrum: WaveSpectrum, seed: Optional[int] = None) ‑> Tuple[numpy.ndarray[Any, numpy.dtype[+ScalarType]], numpy.ndarray[Any, numpy.dtype[+ScalarType]]]

Create a timeseries for from a given power spectral density.

:param component: Wave component to create a timeseries for: u,v,w,x,y,z. :param sampling_frequency: Sampling frequency of output signal in Hertz :param signal_length: Length of output signal :param spectrum: Input power spectrum :param seed: Input seed for the random number generator. :return:

Expand source code
def surface_timeseries(
    component: Literal["u", "v", "w", "x", "y", "z"],
    sampling_frequency: float,
    signal_length: int,
    spectrum: WaveSpectrum,
    seed: Optional[int] = None,
) -> Tuple[np.typing.NDArray, np.typing.NDArray]:
    """
    Create a timeseries for from a given power spectral density.

    :param component: Wave component to create a timeseries for: u,v,w,x,y,z.
    :param sampling_frequency: Sampling frequency of output signal in Hertz
    :param signal_length: Length of output signal
    :param spectrum: Input power spectrum
    :param seed: Input seed for the random number generator.
    :return:
    """

    nfft = (int(signal_length) // 2) * 2

    frequencies = np.linspace(0, 0.5 * sampling_frequency, nfft // 2, endpoint=False)

    time = np.linspace(0, nfft / sampling_frequency, nfft, endpoint=False)

    timeseries = nfft * np.fft.irfft(
        create_fourier_amplitudes(component, spectrum, frequencies, seed)
    )

    return time, timeseries