# -*- coding: utf-8 -*-
import numpy as np
from ..misc import check_random_state, check_random_state_children
from ..signal import signal_distort, signal_simulate, signal_smooth
[docs]
def rsp_simulate(
duration=10,
length=None,
sampling_rate=1000,
noise=0.01,
respiratory_rate=15,
method="breathmetrics",
random_state=None,
random_state_distort="spawn",
):
"""**Simulate a respiratory signal**
Generate an artificial (synthetic) respiratory signal of a given duration
and rate.
Parameters
----------
duration : int
Desired length of duration (s).
sampling_rate : int
The desired sampling rate (in Hz, i.e., samples/second).
length : int
The desired length of the signal (in samples).
noise : float
Noise level (amplitude of the laplace noise).
respiratory_rate : float
Desired number of breath cycles in one minute.
method : str
The model used to generate the signal. Can be ``"sinusoidal"`` for a simulation based on a
trigonometric sine wave that roughly approximates a single respiratory cycle. If
``"breathmetrics"`` (default), will use an advanced model desbribed by
`Noto, et al. (2018) <https://github.com/zelanolab/breathmetrics>`_.
random_state : None, int, numpy.random.RandomState or numpy.random.Generator
Seed for the random number generator. See for ``misc.check_random_state`` for further information.
random_state_distort : {'legacy', 'spawn'}, None, int, numpy.random.RandomState or numpy.random.Generator
Random state to be used to distort the signal. If ``"legacy"``, use the same random state used to
generate the signal (discouraged as it creates dependent random streams). If ``"spawn"``, spawn
independent children random number generators from the random_state argument. If any of the other types,
generate independent children random number generators from the random_state_distort provided (this
allows generating multiple version of the same signal distorted by different random noise realizations).
See Also
--------
rsp_clean, rsp_findpeaks, signal_rate, rsp_process, rsp_plot
Returns
-------
array
Vector containing the respiratory signal.
Examples
--------
.. ipython:: python
import pandas as pd
import neurokit2 as nk
rsp1 = nk.rsp_simulate(duration=30, method="sinusoidal")
rsp2 = nk.rsp_simulate(duration=30, method="breathmetrics")
@savefig p_rsp_simulate1.png scale=100%
pd.DataFrame({"RSP_Simple": rsp1, "RSP_Complex": rsp2}).plot(subplots=True)
@suppress
plt.close()
References
----------
* Noto, T., Zhou, G., Schuele, S., Templer, J., & Zelano, C. (2018). Automated analysis of
breathing waveforms using BreathMetrics: A respiratory signal processing toolbox. Chemical Senses, 43(8), 583-597.
"""
# Seed the random generator for reproducible results
rng = check_random_state(random_state)
random_state_distort = check_random_state_children(random_state, random_state_distort, n_children=1)
# Generate number of samples automatically if length is unspecified
if length is None:
length = duration * sampling_rate
if method.lower() in ["sinusoidal", "sinus", "simple"]:
rsp = _rsp_simulate_sinusoidal(
duration=duration, sampling_rate=sampling_rate, respiratory_rate=respiratory_rate
)
else:
rsp = _rsp_simulate_breathmetrics(
duration=duration,
sampling_rate=sampling_rate,
respiratory_rate=respiratory_rate,
rng=rng,
)
rsp = rsp[0:length]
# Add random noise
if noise > 0:
rsp = signal_distort(
rsp,
sampling_rate=sampling_rate,
noise_amplitude=noise,
noise_frequency=[5, 10, 100],
noise_shape="laplace",
random_state=random_state_distort[0],
silent=True,
)
return rsp
# =============================================================================
# Simple Sinusoidal Model
# =============================================================================
def _rsp_simulate_sinusoidal(duration=10, sampling_rate=1000, respiratory_rate=15):
"""Generate an artificial (synthetic) respiratory signal by trigonometric sine wave that
roughly approximates a single respiratory cycle."""
# Generate values along the length of the duration
rsp = signal_simulate(
duration=duration,
sampling_rate=sampling_rate,
frequency=respiratory_rate / 60,
amplitude=0.5,
)
return rsp
# =============================================================================
# BreathMetrics Model
# =============================================================================
def _rsp_simulate_breathmetrics_original(
nCycles=100,
sampling_rate=1000,
breathing_rate=0.25,
average_amplitude=0.5,
amplitude_variance=0.1,
phase_variance=0.1,
inhale_pause_percent=0.3,
inhale_pause_avgLength=0.2,
inhale_pauseLength_variance=0.5,
exhale_pause_percent=0.3,
exhale_pause_avgLength=0.2,
exhale_pauseLength_variance=0.5,
pause_amplitude=0.1,
pause_amplitude_variance=0.2,
signal_noise=0.1,
rng=None,
):
"""Simulates a recording of human airflow data by appending individually constructed sin waves and pauses in
sequence. This is translated from the matlab code available `here.
<https://github.com/zelanolab/breathmetrics/blob/master/simulateRespiratoryData.m>`_ by Noto, et al. (2018).
Parameters
----------
nCycles : int or float
number of breathing cycles to simulate.
sampling_rate : int
sampling rate.
breathing_rate : float
average breathing rate.
average_amplitude : float
average amplitude of inhales and exhales.
amplitude_variance: float
variance in respiratory amplitudes.
phase_variance: float
variance in duration of individual breaths.
inhale_pause_percent : float
percent of inhales followed by a pause.
inhale_pause_avgLength : float
average length of inhale pauses.
inhale_pauseLength_variance : float
variance in inhale pause length.
exhale_pause_percent : float
percent of exhales followed by a pause.
exhale_pause_avgLength : float
average length of exhale pauses.
exhale_pauseLength_variance : float
variance in exhale pause length.
pause_amplitude : float
noise amplitude of pauses.
pause_amplitude_variance : float
variance in pause noise.
signal_noise : float
percent of noise saturation in the simulated signal.
Returns
----------
signal
vector containing breathmetrics simulated rsp signal.
"""
# Define additional parameters
sample_phase = sampling_rate / breathing_rate
inhale_pause_phase = np.round(inhale_pause_avgLength * sample_phase).astype(int)
exhale_pause_phase = np.round(exhale_pause_avgLength * sample_phase).astype(int)
# Normalize variance by average breath amplitude
amplitude_variance_normed = average_amplitude * amplitude_variance
amplitudes_with_noise = rng.standard_normal(nCycles) * amplitude_variance_normed + average_amplitude
amplitudes_with_noise[amplitudes_with_noise < 0] = 0
# Normalize phase by average breath length
phase_variance_normed = phase_variance * sample_phase
phases_with_noise = np.round(rng.standard_normal(nCycles) * phase_variance_normed + sample_phase).astype(int)
phases_with_noise[phases_with_noise < 0] = 0
# Normalize pause lengths by phase and variation
inhale_pauseLength_variance_normed = inhale_pause_phase * inhale_pauseLength_variance
inhale_pauseLengths_with_noise = np.round(
rng.standard_normal(nCycles) * inhale_pauseLength_variance_normed + inhale_pause_phase
).astype(int)
inhale_pauseLengths_with_noise[inhale_pauseLengths_with_noise < 0] = 0
exhale_pauseLength_variance_normed = exhale_pause_phase * exhale_pauseLength_variance
exhale_pauseLengths_with_noise = np.round(
rng.standard_normal(nCycles) * exhale_pauseLength_variance_normed + inhale_pause_phase
).astype(int)
# why inhale pause phase?
exhale_pauseLengths_with_noise[exhale_pauseLengths_with_noise < 0] = 0
# Normalize pause amplitudes
pause_amplitude_variance_normed = pause_amplitude * pause_amplitude_variance
# Initialize empty vector to fill with simulated data
simulated_respiration = []
# Initialize parameters to save
inhale_onsets = np.zeros(nCycles)
exhale_onsets = np.zeros(nCycles)
inhale_pause_onsets = np.zeros(nCycles)
exhale_pause_onsets = np.zeros(nCycles)
inhale_lengths = np.zeros(nCycles)
inhale_pauseLengths = np.zeros(nCycles)
exhale_lengths = np.zeros(nCycles)
exhale_pauseLengths = np.zeros(nCycles)
inhale_peaks = np.zeros(nCycles)
exhale_troughs = np.zeros(nCycles)
i = 1
for c in range(nCycles):
# Determine length of inhale pause for this cycle
if rng.uniform() < inhale_pause_percent:
this_inhale_pauseLength = inhale_pauseLengths_with_noise[c]
this_inhale_pause = rng.standard_normal(this_inhale_pauseLength) * pause_amplitude_variance_normed
this_inhale_pause[this_inhale_pause < 0] = 0
else:
this_inhale_pauseLength = 0
this_inhale_pause = []
# Determine length of exhale pause for this cycle
if rng.uniform() < exhale_pause_percent:
this_exhale_pauseLength = exhale_pauseLengths_with_noise[c]
this_exhale_pause = rng.standard_normal(this_exhale_pauseLength) * pause_amplitude_variance_normed
this_exhale_pause[this_exhale_pause < 0] = 0
else:
this_exhale_pauseLength = 0
this_exhale_pause = []
# Determine length of inhale and exhale for this cycle to main
# breathing rate
cycle_length = phases_with_noise[c] - (this_inhale_pauseLength + this_exhale_pauseLength)
# If pauses are longer than the time alloted for this breath, set them
# to 0 so a real breath can be simulated. This will deviate the
# statistics from those initialized but is unavaoidable at the current
# state
if (cycle_length <= 0) or (cycle_length < min(phases_with_noise) / 4):
this_inhale_pauseLength = 0
this_inhale_pause = []
this_exhale_pauseLength = 0
this_exhale_pause = []
cycle_length = phases_with_noise[c] - (
this_inhale_pauseLength + this_exhale_pauseLength
)
# Compute inhale and exhale for this cycle
this_cycle = np.sin(np.linspace(0, 2 * np.pi, cycle_length)) * amplitudes_with_noise[c]
half_cycle = np.round(len(this_cycle) / 2).astype(int)
this_inhale = this_cycle[0:half_cycle]
this_inhale_length = len(this_inhale)
this_exhale = this_cycle[half_cycle:]
this_exhale_length = len(this_exhale)
# Save parameters for checking
inhale_lengths[c] = this_inhale_length
inhale_pauseLengths[c] = this_inhale_pauseLength
exhale_lengths[c] = this_exhale_length
exhale_pauseLengths[c] = this_exhale_pauseLength
inhale_onsets[c] = i
exhale_onsets[c] = i + this_inhale_length + this_inhale_pauseLength
if len(this_inhale_pause) > 0:
inhale_pause_onsets[c] = i + this_inhale_length
else:
inhale_pause_onsets[c] = np.nan
if len(this_exhale_pause) > 0:
exhale_pause_onsets[c] = (
i + this_inhale_length + this_inhale_pauseLength + this_exhale_length
)
else:
exhale_pause_onsets[c] = np.nan
# Compose breath from parameters
this_breath = np.hstack([this_inhale, this_inhale_pause, this_exhale, this_exhale_pause])
# Compute max flow for inhale and exhale for this breath
max_ID = np.argmax(this_breath)
min_ID = np.argmin(this_breath)
inhale_peaks[c] = i + max_ID
exhale_troughs[c] = i + min_ID
# Append breath to simulated resperation vector
simulated_respiration = np.hstack([simulated_respiration, this_breath])
i = i + len(this_breath) - 1
# Smooth signal
simulated_respiration = signal_smooth(
simulated_respiration, kernel="boxzen", size=sampling_rate / 2
)
if signal_noise == 0:
signal_noise = 0.0001
noise_vector = rng.uniform(size=simulated_respiration.shape) * average_amplitude
simulated_respiration = simulated_respiration * (1 - signal_noise) + noise_vector * signal_noise
raw_features = {
"Inhale Onsets": inhale_onsets,
"Exhale Onsets": exhale_onsets,
"Inhale Pause Onsets": inhale_pause_onsets,
"Exhale Pause Onsets": exhale_pause_onsets,
"Inhale Lengths": inhale_lengths / sampling_rate,
"Inhale Pause Lengths": inhale_pauseLengths / sampling_rate,
"Exhale Lengths": exhale_lengths / sampling_rate,
"Exhale Pause Lengths": exhale_pauseLengths / sampling_rate,
"Inhale Peaks": inhale_peaks,
"Exhale Troughs": exhale_troughs,
}
if len(inhale_pauseLengths[inhale_pauseLengths > 0]) > 0:
avg_inhale_pauseLength = np.mean(inhale_pauseLengths[inhale_pauseLengths > 0])
else:
avg_inhale_pauseLength = 0
if len(exhale_pauseLengths[exhale_pauseLengths > 0]) > 0:
avg_exhale_pauseLength = np.mean(exhale_pauseLengths[exhale_pauseLengths > 0])
else:
avg_exhale_pauseLength = 0
estimated_breathing_rate = (1 / np.mean(np.diff(inhale_onsets))) * sampling_rate
feature_stats = {
"Breathing Rate": estimated_breathing_rate,
"Average Inhale Length": np.mean(inhale_lengths / sampling_rate),
"Average Inhale Pause Length": avg_inhale_pauseLength / sampling_rate,
"Average Exhale Length": np.mean(exhale_lengths / sampling_rate),
"Average Exhale Pause Length": avg_exhale_pauseLength / sampling_rate,
}
return simulated_respiration, raw_features, feature_stats
def _rsp_simulate_breathmetrics(duration=10, sampling_rate=1000, respiratory_rate=15, rng=None):
n_cycles = int(respiratory_rate / 60 * duration)
# Loop until it doesn't fail
rsp = False
while rsp is False:
# Generate a longer than necessary signal so it won't be shorter
rsp, _, __ = _rsp_simulate_breathmetrics_original(
nCycles=int(n_cycles * 1.5),
sampling_rate=sampling_rate,
breathing_rate=respiratory_rate / 60,
signal_noise=0,
rng=rng,
)
return rsp
