#!/usr/bin/python
# vim: set expandtab ts=4 sw=4:
"""
Simulation functions.
Routines:
"""
import logging
import numpy as np
from scipy import signal
from .support import ensure_equal_dims
# Housekeeping for logging
logger = logging.getLogger(__name__)
# Joint Instantaneous Frequency Functions
[docs]def compute_joint_if(freq, amp, phase, sample_rate=128, seconds=2):
"""Compute joint instantaneous frequency from a set of oscillations.
This function implements a signal simulator based on the methods in Fabus
et al (2021) [1]_. freq, amp and phase inputs should be tuples/lists of
user defined values.
Parameters
----------
freq, amp, phase : {tuple, list, np.ndarray}
Frequency, Amplitude and Phase values for each component.
These are lists or tuples containing a single value per component.
sample_rate and seconds must then also be defined.
sample_rate : {None, float}
Sampling frequency of the data used if user defined harmonic values are passed in
seconds : {None, float}
Amount of seconds to generate if user defined harmonic values are passed in
Returns
-------
joint_if : ndarray
Vector containing the joint instantaneous frequency signal
joint_sig : ndarray
Array containing the time-domain signal for each harmonic component
Notes
-----
Example usage - compute joint instantaneous frequency from user defined harmonic values
>>> f = (5, 10, 15)
>>> a = (1, 1/3, 1/9)
>>> p = (0, 0, 0)
>>> joint_if, joint_sig = compute_joint_if(f, a, p, 128, 10)
References
----------
.. [1] Fabus, M., Woolrich, M., Warnaby, C. and Quinn, A., 2021. Understanding
Harmonic Structures Through Instantaneous Frequency. BiorXiv
https://doi.org/10.1101/2021.12.21.473676
"""
time_vect = np.linspace(0, seconds, int(seconds*sample_rate))
# Work with numpy arrays internally
freq = 2*np.pi*np.array(freq)
amp = np.array(amp)
phase = np.array(phase)
ensure_equal_dims([freq, amp, phase], ['freq', 'amp', 'phase'], 'compute_joint_if')
num_comps = freq.shape[0]
num = np.zeros((num_comps, num_comps, time_vect.shape[0]))
denom_sin = np.zeros((num_comps, time_vect.shape[0]))
denom_cos = np.zeros((num_comps, time_vect.shape[0]))
sig = np.zeros((num_comps, time_vect.shape[0]))
for n in range(num_comps):
denom_cos[n, :] = amp[n] * np.cos(freq[n] * time_vect + phase[n])
denom_sin[n, :] = amp[n] * np.sin(freq[n] * time_vect + phase[n])
sig[n, :] = amp[n] * np.cos(freq[n] * time_vect + phase[n])
for m in range(num_comps):
fd = freq[n] - freq[m]
pd = phase[n] - phase[m]
num[n, m, :] = freq[m] * amp[n] * amp[m] * np.cos(fd * time_vect + pd)
joint_if = np.sum(num, axis=(0, 1)) / ((np.sum(denom_cos, axis=0)**2) + (np.sum(denom_sin, axis=0)**2))
joint_if = joint_if / (2*np.pi)
return joint_if, sig
[docs]def abreu2010(f, nonlin_deg, nonlin_phi, sample_rate, seconds):
"""Simulate a non-linear waveform using equation 9 in [1]_.
Parameters
----------
f : float
Fundamental frequency of generated signal
nonlin_deg : float
Degree of non-linearity in generated signal
nonlin_phi : float
Skew in non-linearity of generated signal
sample_rate : float
The sampling frequency of the generated signal
seconds : float
The number of seconds of data to generate
Returns
-------
ndarray
Simulated signal containing non-linear wave
References
----------
.. [1] Abreu, T., Silva, P. A., Sancho, F., & Temperville, A. (2010).
Analytical approximate wave form for asymmetric waves. Coastal Engineering,
57(7), 656–667. https://doi.org/10.1016/j.coastaleng.2010.02.005
"""
time_vect = np.linspace(0, seconds, int(seconds * sample_rate))
factor = np.sqrt(1 - nonlin_deg**2)
num = nonlin_deg * np.sin(nonlin_phi) / 1 + np.sqrt(1 - nonlin_deg**2)
num = num + np.sin(2 * np.pi * f * time_vect)
denom = 1 - nonlin_deg * np.cos(2 * np.pi * f * time_vect + nonlin_phi)
return factor * (num / denom)
[docs]def ar_oscillator(freq, sample_rate, seconds, r=.95, noise_std=None, random_seed=None):
"""Create a simulated oscillation using an autoregressive filter.
A simple filter is defined by direct pole placement and applied to white
noise to generate a random signal with a defined oscillatory peak frequency
that exhibits random variability frequency, amplitude and waveform.
Parameters
----------
freq : float
Peak resonant frequency of the simulated filter.
sample_rate : float
Sampling frequency for the simulation
seconds : float
Number of seconds of data to simulate
r : float (0 < r < 1)
Pole magnitude of simulated autoregressive resonance.
noise_std : float
Scaling of optional noise to add to simulation. Scaling is relative to
standard-deviation of the simulated data.
random_seed : int
Optional random seed generation
Returns
-------
ndarray
A simulated time course.
"""
if random_seed is not None:
np.random.seed(random_seed)
if freq > 0:
freq_rads = (2 * np.pi * freq) / sample_rate
a1 = np.array([1, -2*r*np.cos(freq_rads), (r**2)])
else:
a1 = np.poly(r)
num_samples = int(sample_rate * seconds)
x = signal.filtfilt(1, a1, np.random.randn(1, num_samples)).T
if noise_std is not None:
noise = np.std(x)*noise_std*np.random.randn(1, num_samples).T
x = x + noise
if random_seed is not None:
np.random.seed() # restore defaults
return x