Source code for mne_hfo.simulate

"""
Some parts of code are recoded from package Anderson Brito da Silva's pyhfo.

Reference: (https://github.com/britodasilva/pyhfo)
"""
import numpy as np
from scipy.stats import norm


# ----- Noise types -----
[docs]def simulate_pink_noise(N): """ Create a pink noise (1/f) with N points. Parameters ---------- N: int Number of samples to be returned Returns ------- pink_noise_array: numpy array 1D array of pink noise """ M = N # ensure that the N is even if N % 2: N += 1 x = np.random.randn(N) # generate a white noise X = np.fft.fft(x) # FFT # prepare a vector for 1/f multiplication nPts = int(N / 2 + 1) n = range(1, nPts + 1) n = np.sqrt(n) # multiplicate the left half of the spectrum X[range(nPts)] = X[range(nPts)] / n # prepare a right half of the spectrum - a copy of the left one X[range(nPts, N)] = np.real(X[range(int(N / 2 - 1), 0, -1)]) X[range(nPts, N)] -= 1j * np.imag(X[range(int(N / 2 - 1), 0, -1)]) y = np.fft.ifft(X) # IFFT y = np.real(y) # normalising y -= np.mean(y) y /= np.sqrt(np.mean(y ** 2)) # returning size of N if M % 2 == 1: y = y[:-1] return y
[docs]def simulate_brown_noise(N): """ Create a brown noise (1/f²) with N points. Parameters ---------- N: int Number of samples to be returned Returns ------- brown_noise_array: numpy array 1D array of brown noise """ M = N # ensure that the N is even if N % 2: N += 1 x = np.random.randn(N) # generate a white noise X = np.fft.fft(x) # FFT # prepare a vector for 1/f² multiplication nPts = int(N / 2 + 1) n = range(1, nPts + 1) # multiplicate the left half of the spectrum X[range(nPts)] = X[range(nPts)] / n # prepare a right half of the spectrum - a copy of the left one X[range(nPts, N)] = np.real(X[range(int(N / 2 - 1), 0, -1)]) X[range(nPts, N)] -= 1j * np.imag(X[range(int(N / 2 - 1), 0, -1)]) y = np.fft.ifft(X) # IFFT y = np.real(y) # normalising y -= np.mean(y) y /= np.sqrt(np.mean(y ** 2)) # returning size of N if M % 2 == 1: y = y[:-1] return y
# ----- Artifacts -----
[docs]def simulate_delta(fs=5000, decay_dur=None): """ Delta function with exponential decay. Parameters ---------- fs: int Sampling frequency of the signal (default=5000) decay_dur: float Decay duration before returning to 0 in seconds Returns ------- delta: numpy array 1D numpy array with delta function """ if decay_dur is None: decay_dur = np.random.random() decay_N = int(fs * decay_dur) return_value = 0.001 # This is the value where decay finishes decay_factor = np.log(return_value) / -decay_N t = np.linspace(0, decay_N, decay_N, endpoint=False) decay = np.exp(-t * decay_factor) delta = np.concatenate([[0], decay]) return delta
[docs]def simulate_line_noise(fs=5000, freq=50, numcycles=None): """ Line noise artifact. Parameters ---------- fs: int Sampling frequency of the signal (default=5000) freq: float Line noise frequency (default=50) ncycles: float Number of cycles to create Returns ------- line_noise: numpy array 1D numpy array with line noise """ if numcycles is None: numcycles = np.random.randint(3, 50) dur_samps = int((numcycles / freq) * fs) x = np.arange(dur_samps) y = np.sin(2 * np.pi * freq * x / fs) return y
[docs]def simulate_artifact_spike(fs=5000, dur=None): """ Artifact like spike (sharp, not gaussian). Parameters ---------- fs: int Sampling frequency of the signal (default=5000) dur: float Duration of the event in seconds Returns ------- artifact_spike: numpy array 1D numpy array with artifact spike """ if dur is None: dur = round(np.random.random() / 10, 3) N = int(fs * dur) if not N % 2: # Check if the number is odd - we want to have proper spike N -= 1 y = np.zeros(N) y[:int(N / 2) + 1] = np.linspace(0, 1, int(N / 2) + 1) y[-int(N / 2):] = np.linspace(1, 0, int(N / 2) + 1)[1:] return y
# ----- HFO ----- def _wavelet(numcycles, f, fs): """ Create a wavelet. Parameters ---------- numcycles: int Number of cycles (gaussian window) f: float Central frequency fs: float Signal sampling rate (default=5000) Returns ------- wave: numpy array 1D numpy array with waveform time: numpy array 1D numpy array with the time vector """ N = int((fs * numcycles) / f) time = np.linspace((-numcycles / 2) / float(f), (numcycles / 2) / float(f), N) # time vector std = numcycles / (2 * np.pi * f) # standard deviation wave = np.exp(2 * 1j * np.pi * f * time) wave *= np.exp(-(time ** 2) / (2 * (std ** 2))) # waveform return wave, time
[docs]def simulate_hfo(fs=5000, freq=None, numcycles=None): """ Create a simulated HFO signal. Parameters ---------- fs: float Sampling rate of the signal (default=5000) freq: float Frequency of the artificial HFO (default=None - random frequency between 80 nad 600 Hz) numcycles: int Number of HFO cycles (default=None - cycles between 9 - 15) Returns ------- wave: numpy array 1D numpy array with waveform time: numpy array 1D numpy array with the time vector """ if numcycles is None: numcycles = np.random.randint(9, 15) if freq is None: freq = np.random.randint(80, 600) wave, time = _wavelet(numcycles, freq, fs) return np.real(wave), time
# ----- Spike -----
[docs]def simulate_spike(fs=5000, dur=None): """ Create a simple gaussian spike. Parameters ---------- fs: float Sampling rate (default=5000) dur: float Spike duration in seconds Returns ------- spike numpy array 1D numpy array with a sipke """ if dur is None: dur = round(np.random.random() * 0.5, 2) x = np.linspace(-4, 4, int(fs * dur)) # 4 stds spike_dist = norm.pdf(x, loc=0, scale=1) # Create gaussian shape # Normalize so that the peak is at 1 spike = spike_dist * 1 / max(spike_dist) return spike