Repository: neuropsychology/NeuroKit (language: python)

Issue:
MNE 1.0: adjustments
```python
import neurokit2 as nk
import mne

raw = mne.io.read_raw_fif(
        mne.datasets.sample.data_path() + "/MEG/sample/sample_audvis_raw.fif", preload=True
    )
x = nk.mne_channel_extract(raw, what="EOG", name="EOG")
x
```

This used to return a vector (I think?), but now it returns a series (at least based on the bug in the test that caught this)

TODO: need to investigate `mne_channel_extract()` and see if changes are to be made based on the most recent MNE update 
MNE 1.0: adjustments
```python
import neurokit2 as nk
import mne

raw = mne.io.read_raw_fif(
        mne.datasets.sample.data_path() + "/MEG/sample/sample_audvis_raw.fif", preload=True
    )
x = nk.mne_channel_extract(raw, what="EOG", name="EOG")
x
```

This used to return a vector (I think?), but now it returns a series (at least based on the bug in the test that caught this)

TODO: need to investigate `mne_channel_extract()` and see if changes are to be made based on the most recent MNE update 


Discussion / hints:

Fix patch:
diff --git a/README.rst b/README.rst
index 2731f7aed0..38ddbbb2e3 100644
--- a/README.rst
+++ b/README.rst
@@ -213,6 +213,8 @@ Physiological Data Preprocessing
 Simulate physiological signals
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
+You can easily simulate artificial ECG (also `12-Lead multichannel ECGs <https://github.com/neuropsychology/NeuroKit/...>`_), PPG, RSP, EDA, and EMG signals to test your scripts and algorithms.
+
 .. code-block:: python
 
     import numpy as np
diff --git a/docs/studies/ecg_generating_12_leads.rst b/docs/studies/ecg_generating_12_leads.rst
new file mode 100644
index 0000000000..7c797187d1
--- /dev/null
+++ b/docs/studies/ecg_generating_12_leads.rst
@@ -0,0 +1,1 @@
+.. mdinclude:: ../../studies/ecg_generating_12_leads/README.md
diff --git a/docs/studies/index.rst b/docs/studies/index.rst
index a1e7470ff4..dc693cd076 100644
--- a/docs/studies/index.rst
+++ b/docs/studies/index.rst
@@ -7,6 +7,7 @@ Contents:
    :maxdepth: 1
 
    ecg_preprocessing
+   ecg_generating_12_leads
    eog_blinktemplate
    hrv_structure
    erp_gam
diff --git a/neurokit2/ecg/__init__.py b/neurokit2/ecg/__init__.py
index e7f3e48278..d964965a86 100644
--- a/neurokit2/ecg/__init__.py
+++ b/neurokit2/ecg/__init__.py
@@ -17,7 +17,6 @@
 from .ecg_segment import ecg_segment
 from .ecg_simulate import ecg_simulate
 
-
 __all__ = [
     "ecg_simulate",
     "ecg_clean",
diff --git a/neurokit2/ecg/ecg_simulate.py b/neurokit2/ecg/ecg_simulate.py
index e68d4b9bf3..c44eb1a0be 100644
--- a/neurokit2/ecg/ecg_simulate.py
+++ b/neurokit2/ecg/ecg_simulate.py
@@ -2,13 +2,22 @@
 import math
 
 import numpy as np
+import pandas as pd
 import scipy
 
 from ..signal import signal_distort, signal_resample
 
 
 def ecg_simulate(
-    duration=10, length=None, sampling_rate=1000, noise=0.01, heart_rate=70, heart_rate_std=1, method="ecgsyn", random_state=None
+    duration=10,
+    length=None,
+    sampling_rate=1000,
+    noise=0.01,
+    heart_rate=70,
+    heart_rate_std=1,
+    method="ecgsyn",
+    random_state=None,
+    **kwargs,
 ):
     """Simulate an ECG/EKG signal.
 
@@ -35,10 +44,12 @@ def ecg_simulate(
         Desired heart rate standard deviation (beats per minute).
     method : str
         The model used to generate the signal. Can be 'simple' for a simulation based on Daubechies
-        wavelets that roughly approximates a single cardiac cycle. If 'ecgsyn' (default), will use an
-        advanced model desbribed `McSharry et al. (2003) <https://physionet.org/content/ecgsyn/>`_.
+        wavelets that roughly approximates a single cardiac cycle. If 'ecgsyn' (default), will use the model desbribed `McSharry et al. (2003) <https://physionet.org/content/ecgsyn/>`_. If
+        'multileads', will return a DataFrame containing 12-leads (see `12-leads ECG simulation <https://neurokit2.readthedocs.io/en/latest/studies/ecg_generating_12_leads.html>`_).
     random_state : int
         Seed for the random number generator.
+    **kwargs
+        Other keywords parameters for ECGSYN algorithm, such as `lfhfratio`, `ti`, `ai`, `bi`.
 
     Returns
     -------
@@ -50,11 +61,20 @@ def ecg_simulate(
     >>> import pandas as pd
     >>> import neurokit2 as nk
     >>>
+    >>> # Simulate single lead ECG
     >>> ecg1 = nk.ecg_simulate(duration=10, method="simple")
     >>> ecg2 = nk.ecg_simulate(duration=10, method="ecgsyn")
     >>> pd.DataFrame({"ECG_Simple": ecg1,
     ...               "ECG_Complex": ecg2}).plot(subplots=True) #doctest: +ELLIPSIS
     array([<AxesSubplot:>, <AxesSubplot:>], dtype=object)
+    >>>
+    >>> # Simulate 12-leads ECG
+    >>> ecg12 = nk.ecg_simulate(duration=10, method="multileads")
+    >>> ecg12[0:10000].plot(subplots=True)
+    array([<AxesSubplot:>, <AxesSubplot:>, <AxesSubplot:>, <AxesSubplot:>,
+           <AxesSubplot:>, <AxesSubplot:>, <AxesSubplot:>, <AxesSubplot:>,
+           <AxesSubplot:>, <AxesSubplot:>, <AxesSubplot:>, <AxesSubplot:>],
+          dtype=object)
 
     See Also
     --------
@@ -79,36 +99,73 @@ def ecg_simulate(
 
     # Run appropriate method
     if method.lower() in ["simple", "daubechies"]:
-        ecg = _ecg_simulate_daubechies(
+        signals = _ecg_simulate_daubechies(
             duration=duration, length=length, sampling_rate=sampling_rate, heart_rate=heart_rate
         )
     else:
         approx_number_beats = int(np.round(duration * (heart_rate / 60)))
-        ecg = _ecg_simulate_ecgsyn(
-            sfecg=sampling_rate,
-            N=approx_number_beats,
-            Anoise=0,
-            hrmean=heart_rate,
-            hrstd=heart_rate_std,
-            lfhfratio=0.5,
-            sfint=sampling_rate,
-            ti=(-70, -15, 0, 15, 100),
-            ai=(1.2, -5, 30, -7.5, 0.75),
-            bi=(0.25, 0.1, 0.1, 0.1, 0.4),
-        )
+        if method.lower() in ["multi", "multilead", "multileads", "multichannel"]:
+            # Gamma, a (12,5) matrix to modify the five waves' amplitudes of 12 leads (P, Q, R, S, T)
+            gamma = np.array(
+                [
+                    [1, 0.1, 1, 1.2, 1],
+                    [2, 0.2, 0.2, 0.2, 3],
+                    [1, -0.1, -0.8, -1.1, 2.5],
+                    [-1, -0.05, -0.8, -0.5, -1.2],
+                    [0.05, 0.05, 1, 1, 1],
+                    [1, -0.05, -0.1, -0.1, 3],
+                    [-0.5, 0.05, 0.2, 0.5, 1],
+                    [0.05, 0.05, 1.3, 2.5, 2],
+                    [1, 0.05, 1, 2, 1],
+                    [1.2, 0.05, 1, 2, 2],
+                    [1.5, 0.1, 0.8, 1, 2],
+                    [1.8, 0.05, 0.5, 0.1, 2],
+                ]
+            )
+
+            signals, results = _ecg_simulate_ecgsyn(
+                sfecg=sampling_rate,
+                N=approx_number_beats,
+                hrmean=heart_rate,
+                hrstd=heart_rate_std,
+                sfint=sampling_rate,
+                gamma=gamma,
+                **kwargs,
+            )
+        else:
+            signals, results = _ecg_simulate_ecgsyn(
+                sfecg=sampling_rate,
+                N=approx_number_beats,
+                hrmean=heart_rate,
+                hrstd=heart_rate_std,
+                sfint=sampling_rate,
+                gamma=np.ones((1, 5)),
+                **kwargs,
+            )
         # Cut to match expected length
-        ecg = ecg[0:length]
+        for i in range(len(signals)):
+            signals[i] = signals[i][0:length]
 
     # Add random noise
     if noise > 0:
-        ecg = signal_distort(
-            ecg,
-            sampling_rate=sampling_rate,
-            noise_amplitude=noise,
-            noise_frequency=[5, 10, 100],
-            noise_shape="laplace",
-            random_state=random_state,
-            silent=True,
+        for i in range(len(signals)):
+            signals[i] = signal_distort(
+                signals[i],
+                sampling_rate=sampling_rate,
+                noise_amplitude=noise,
+                noise_frequency=[5, 10, 100],
+                noise_shape="laplace",
+                random_state=random_state,
+                silent=True,
+            )
+
+    # Format
+    if len(signals) == 1:
+        ecg = signals[0]
+    else:
+        ecg = pd.DataFrame(
+            np.array(signals).T,
+            columns=["I", "II", "III", "aVR", "aVL", "aVF", "V1", "V2", "V3", "V4", "V5", "V6"],
         )
 
     # Reset random seed (so it doesn't affect global)
@@ -143,10 +200,15 @@ def _ecg_simulate_daubechies(duration=10, length=None, sampling_rate=1000, heart
 
     # Resample
     ecg = signal_resample(
-        ecg, sampling_rate=int(len(ecg) / 10), desired_length=length, desired_sampling_rate=sampling_rate
+        ecg,
+        sampling_rate=int(len(ecg) / 10),
+        desired_length=length,
+        desired_sampling_rate=sampling_rate,
     )
 
-    return ecg
+    # Return the signal in a list to match
+    # with the potential multichanel output of ecgsyn
+    return [ecg]
 
 
 # =============================================================================
@@ -155,7 +217,6 @@ def _ecg_simulate_daubechies(duration=10, length=None, sampling_rate=1000, heart
 def _ecg_simulate_ecgsyn(
     sfecg=256,
     N=256,
-    Anoise=0,
     hrmean=60,
     hrstd=1,
     lfhfratio=0.5,
@@ -163,43 +224,52 @@ def _ecg_simulate_ecgsyn(
     ti=(-70, -15, 0, 15, 100),
     ai=(1.2, -5, 30, -7.5, 0.75),
     bi=(0.25, 0.1, 0.1, 0.1, 0.4),
+    gamma=np.ones((1, 5)),
+    **kwargs,
 ):
-    """This function is a python translation of the matlab script by `McSharry & Clifford (2013)
-
-    <https://physionet.org/content/ecgsyn>`_.
+    """
+    This function is a python translation of the matlab script by `McSharry & Clifford (2013) <https://physionet.org/content/ecgsyn>`_.
 
     Parameters
     ----------
-    % Operation uses the following parameters (default values in []s):
-    % sfecg: ECG sampling frequency [256 Hertz]
-    % N: approximate number of heart beats [256]
-    % Anoise: Additive uniformly distributed measurement noise [0 mV]
-    % hrmean: Mean heart rate [60 beats per minute]
-    % hrstd: Standard deviation of heart rate [1 beat per minute]
-    % lfhfratio: LF/HF ratio [0.5]
-    % sfint: Internal sampling frequency [256 Hertz]
-    % Order of extrema: (P Q R S T)
-    % ti = angles of extrema (in degrees)
-    % ai = z-position of extrema
-    % bi = Gaussian width of peaks
+    sfecg:
+        ECG sampling frequency [256 Hertz]
+    N:
+        approximate number of heart beats [256]
+    Anoise:
+        Additive uniformly distributed measurement noise [0 mV]
+    hrmean:
+        Mean heart rate [60 beats per minute]
+    hrstd:
+        Standard deviation of heart rate [1 beat per minute]
+    lfhfratio:
+        LF/HF ratio [0.5]
+    sfint:
+        Internal sampling frequency [256 Hertz]
+    ti
+        angles of extrema (in degrees). Order of extrema is (P Q R S T).
+    ai
+        z-position of extrema.
+    bi
+        Gaussian width of peaks.
 
     Returns
     -------
     array
         Vector containing simulated ecg signal.
 
-#    Examples
-#    --------
-#    >>> import matplotlib.pyplot as plt
-#    >>> import neurokit2 as nk
-#    >>>
-#    >>> s = _ecg_simulate_ecgsynth()
-#    >>> x = np.linspace(0, len(s)-1, len(s))
-#    >>> num_points = 4000
-#    >>>
-#    >>> num_points = min(num_points, len(s))
-#    >>> plt.plot(x[:num_points], s[:num_points]) #doctest: +SKIP
-#    >>> plt.show() #doctest: +SKIP
+    #    Examples
+    #    --------
+    #    >>> import matplotlib.pyplot as plt
+    #    >>> import neurokit2 as nk
+    #    >>>
+    #    >>> s = _ecg_simulate_ecgsynth()
+    #    >>> x = np.linspace(0, len(s)-1, len(s))
+    #    >>> num_points = 4000
+    #    >>>
+    #    >>> num_points = min(num_points, len(s))
+    #    >>> plt.plot(x[:num_points], s[:num_points]) #doctest: +SKIP
+    #    >>> plt.show() #doctest: +SKIP
 
     """
 
@@ -224,7 +294,11 @@ def _ecg_simulate_ecgsyn(
     if q != qd:
         raise ValueError(
             "Internal sampling frequency (sfint) must be an integer multiple of the ECG sampling frequency"
-            " (sfecg). Your current choices are: sfecg = " + str(sfecg) + " and sfint = " + str(sfint) + "."
+            " (sfecg). Your current choices are: sfecg = "
+            + str(sfecg)
+            + " and sfint = "
+            + str(sfint)
+            + "."
         )
 
     # Define frequency parameters for rr process
@@ -266,26 +340,40 @@ def _ecg_simulate_ecgsyn(
     Tspan = [0, (Nt - 1) * dt]
     t_eval = np.linspace(0, (Nt - 1) * dt, Nt)
 
-    # as passing extra arguments to derivative function is not supported yet in solve_ivp
-    # lambda function is used to serve the purpose
-    result = scipy.integrate.solve_ivp(
-        lambda t, x: _ecg_simulate_derivsecgsyn(t, x, rrn, ti, sfint, ai, bi), Tspan, x0, t_eval=t_eval
-    )
-    X0 = result.y
+    # Initialize results containers
+    results = []
+    signals = []
+
+    # Multichannel modification (#625):
+    # --------------------------------------------------
+    # Loop over the twelve leads modifying ai in the loop to generate each lead's data
+    # Because these are all starting at the same position, it may make sense to grab a random segment within the series to simulate random phase and to forget the initial conditions
+
+    for lead in range(len(gamma)):
+        # as passing extra arguments to derivative function is not supported yet in solve_ivp
+        # lambda function is used to serve the purpose
+        result = scipy.integrate.solve_ivp(
+            lambda t, x: _ecg_simulate_derivsecgsyn(t, x, rrn, ti, sfint, gamma[lead] * ai, bi),
+            Tspan,
+            x0,
+            t_eval=t_eval,
+        )
+        results.append(result)  # store results
+        X0 = result.y  # get signal
+
+        # downsample to required sfecg
+        X = X0[:, np.arange(0, X0.shape[1], q).astype(int)]
 
-    # downsample to required sfecg
-    X = X0[:, np.arange(0, X0.shape[1], q).astype(int)]
+        # Scale signal to lie between -0.4 and 1.2 mV
+        z = X[2, :].copy()
+        zmin = np.min(z)
+        zmax = np.max(z)
+        zrange = zmax - zmin
+        z = (z - zmin) * 1.6 / zrange - 0.4
 
-    # Scale signal to lie between -0.4 and 1.2 mV
-    z = X[2, :].copy()
-    zmin = np.min(z)
-    zmax = np.max(z)
-    zrange = zmax - zmin
-    z = (z - zmin) * 1.6 / zrange - 0.4
+        signals.append(z)
 
-    # include additive uniformly distributed measurement noise
-    eta = 2 * np.random.uniform(len(z)) - 1
-    return z + Anoise * eta  # Return signal
+    return signals, results
 
 
 def _ecg_simulate_derivsecgsyn(t, x, rr, ti, sfint, ai, bi):
diff --git a/neurokit2/eeg/eeg_rereference.py b/neurokit2/eeg/eeg_rereference.py
index fade7c38a5..b16c4a2ca6 100644
--- a/neurokit2/eeg/eeg_rereference.py
+++ b/neurokit2/eeg/eeg_rereference.py
@@ -41,23 +41,23 @@ def eeg_rereference(eeg, reference="average", robust=False, **kwargs):
     ---------
     >>> import neurokit2 as nk
     >>>
-    >>> eeg = nk.mne_data("filt-0-40_raw")
+    >>> raw = nk.mne_data("filt-0-40_raw")
+    >>> eeg = raw.get_data()
     >>>
     >>> # Difference between robust average
     >>> avg = nk.eeg_rereference(eeg, 'average', robust=False)
     >>> avg_r = nk.eeg_rereference(eeg, 'average', robust=True)
     >>>
-    >>> nk.signal_plot([avg.get_data()[0, 0:1000],
-    ...                 avg_r.get_data()[0, 0:1000]])
+    >>> nk.signal_plot([avg[0, 0:1000], avg_r[0, 0:1000]])
     >>>
     >>> # Compare the rerefering of an array vs. the MNE object
-    >>> data_mne = eeg.copy().set_eeg_reference('average', verbose=False).get_data()
-    >>> data_nk = nk.eeg_rereference(eeg.get_data(), 'average')
+    >>> avg_mne = raw.copy().set_eeg_reference('average', verbose=False)
+    >>> nk.signal_plot([avg[0, 0:1000], avg_mne.get_data()[0, 0:1000]])
     >>>
     >>> # Difference between average and LAP
-    >>> lap = nk.eeg_rereference(eeg, 'lap')
+    >>> lap = nk.eeg_rereference(raw, 'lap')
     >>>
-    >>> nk.signal_plot([avg.get_data()[0, 0:1000],
+    >>> nk.signal_plot([avg_mne.get_data()[0, 0:1000],
     ...                 lap.get_data()[0, 0:1000]], standardize=True)
 
     References
@@ -99,8 +99,15 @@ def eeg_rereference_mne(eeg, reference="average", robust=False, **kwargs):
 
     eeg = eeg.copy()
     if reference == "average" and robust is True:
-        eeg._data = eeg_rereference_array(eeg._data, reference=reference, robust=robust)
-        eeg.info["custom_ref_applied"] = True
+        # Assigning "custom_ref_applied" to True throws an error with the
+        # latest MNE. If this error goes away in the future, we might able to
+        # restore this feature.
+        # > eeg._data = eeg_rereference_array(eeg._data, reference=reference, robust=robust)
+        # > eeg.info["custom_ref_applied"] = True
+        raise ValueError(
+            "NeuroKit error: eeg_rereference(): 'robust=True' currently not supported for MNE",
+            " objects.",
+        )
     elif reference in ["lap", "csd"]:
         try:
             import mne
diff --git a/neurokit2/eeg/mne_data.py b/neurokit2/eeg/mne_data.py
index e38f787d9b..61d6b101fe 100644
--- a/neurokit2/eeg/mne_data.py
+++ b/neurokit2/eeg/mne_data.py
@@ -39,7 +39,7 @@ def mne_data(what="raw", path=None):
     # Find path of mne data
     if path is None:
         try:
-            path = mne.datasets.sample.data_path()
+            path = str(mne.datasets.sample.data_path())
         except ValueError:
             raise ValueError(
                 "NeuroKit error: the mne sample data folder does not exist. ",
diff --git a/neurokit2/signal/signal_distort.py b/neurokit2/signal/signal_distort.py
index 34abbb0664..feed1f53a5 100644
--- a/neurokit2/signal/signal_distort.py
+++ b/neurokit2/signal/signal_distort.py
@@ -145,9 +145,17 @@ def signal_distort(
 
     distorted = signal + noise
 
+    # Reset random seed (so it doesn't affect global)
+    np.random.seed(None)
+
     return distorted
 
 
+# ===========================================================================
+# Types of Noise
+# ===========================================================================
+
+
 def _signal_linear_drift(signal):
 
     n_samples = len(signal)
@@ -201,12 +209,21 @@ def _signal_distort_artifacts(
 
 
 def _signal_distort_powerline(
-    signal, signal_sd=None, sampling_rate=1000, powerline_frequency=50, powerline_amplitude=0.1, silent=False
+    signal,
+    signal_sd=None,
+    sampling_rate=1000,
+    powerline_frequency=50,
+    powerline_amplitude=0.1,
+    silent=False,
 ):
 
     duration = len(signal) / sampling_rate
     powerline_noise = signal_simulate(
-        duration=duration, sampling_rate=sampling_rate, frequency=powerline_frequency, amplitude=1, silent=silent
+        duration=duration,
+        sampling_rate=sampling_rate,
+        frequency=powerline_frequency,
+        amplitude=1,
+        silent=silent,
     )
 
     if signal_sd is not None:
@@ -226,7 +243,9 @@ def _signal_distort_noise_multifrequency(
     silent=False,
 ):
     base_noise = np.zeros(len(signal))
-    params = listify(noise_amplitude=noise_amplitude, noise_frequency=noise_frequency, noise_shape=noise_shape)
+    params = listify(
+        noise_amplitude=noise_amplitude, noise_frequency=noise_frequency, noise_shape=noise_shape
+    )
 
     for i in range(len(params["noise_amplitude"])):
 
@@ -252,7 +271,12 @@ def _signal_distort_noise_multifrequency(
 
 
 def _signal_distort_noise(
-    n_samples, sampling_rate=1000, noise_frequency=100, noise_amplitude=0.1, noise_shape="laplace", silent=False
+    n_samples,
+    sampling_rate=1000,
+    noise_frequency=100,
+    noise_amplitude=0.1,
+    noise_shape="laplace",
+    silent=False,
 ):
 
     _noise = np.zeros(n_samples)
@@ -267,7 +291,7 @@ def _signal_distort_noise(
                 f" the sampling rate of {sampling_rate} Hz. Please increase "
                 f" sampling rate to {noise_frequency * 10} Hz or choose "
                 f" frequencies smaller than or equal to {nyquist} Hz.",
-                category=NeuroKitWarning
+                category=NeuroKitWarning,
             )
         return _noise
     # Also make sure that at least one period of the frequency can be
@@ -282,7 +306,7 @@ def _signal_distort_noise(
                 f" Please choose noise frequencies larger than "
                 f" {1 / duration} Hz or increase the duration of the "
                 f" signal above {1 / noise_frequency} seconds.",
-                category=NeuroKitWarning
+                category=NeuroKitWarning,
             )
         return _noise
 
@@ -293,7 +317,9 @@ def _signal_distort_noise(
     elif noise_shape == "laplace":
         _noise = np.random.laplace(0, noise_amplitude, noise_duration)
     else:
-        raise ValueError("NeuroKit error: signal_distort(): 'noise_shape' should be one of 'gaussian' or 'laplace'.")
+        raise ValueError(
+            "NeuroKit error: signal_distort(): 'noise_shape' should be one of 'gaussian' or 'laplace'."
+        )
 
     if len(_noise) != n_samples:
         _noise = signal_resample(_noise, desired_length=n_samples, method="interpolation")
diff --git a/studies/ecg_generating_12_leads/3D.png b/studies/ecg_generating_12_leads/3D.png
new file mode 100644
index 0000000000..487e5b7b29
Binary files /dev/null and b/studies/ecg_generating_12_leads/3D.png differ
diff --git a/studies/ecg_generating_12_leads/ECG12_abnormal.png b/studies/ecg_generating_12_leads/ECG12_abnormal.png
new file mode 100644
index 0000000000..00362c47b4
Binary files /dev/null and b/studies/ecg_generating_12_leads/ECG12_abnormal.png differ
diff --git a/studies/ecg_generating_12_leads/ECG12_normal.png b/studies/ecg_generating_12_leads/ECG12_normal.png
new file mode 100644
index 0000000000..cca13ab577
Binary files /dev/null and b/studies/ecg_generating_12_leads/ECG12_normal.png differ
diff --git a/studies/ecg_generating_12_leads/README.md b/studies/ecg_generating_12_leads/README.md
new file mode 100644
index 0000000000..f7377cdc9f
--- /dev/null
+++ b/studies/ecg_generating_12_leads/README.md
@@ -0,0 +1,65 @@
+# Generating Normal and Abnormal 12-leads Synthetic ECG Data
+
+## Normal Multi-lead ECG
+
+We can generate a normal multi-lead ECG signal with the following code:
+
+```python
+import numpy as np
+import neurokit2 as nk
+
+ecg12 = nk.ecg_simulate(duration=10, method="multileads")
+
+# Visualize results
+ecg12[0:10000].plot(subplots=True)
+```
+
+![](../../studies/ecg_generating_12_leads/ECG12_normal.png)
+
+## Abnormal Multi-lead ECG
+
+Abnormal ECG can be generated by modifying the simulation parameters. We will draw random parameters from a distribution.
+
+```python
+# Normal parameters (used by default)
+# ===================================
+# t, the starting position along the circle of each interval in radius
+ti = np.array((-70, -15, 0, 15, 100))
+# a, the amplitude of each spike
+ai = np.array((1.2, -5, 30, -7.5, 0.75))
+# b, the width of each spike
+bi = np.array((0.25, 0.1, 0.1, 0.1, 0.4))
+
+# Add noise
+# ===============
+ti = np.random.normal(ti, np.ones(5) * 3)
+ai = np.random.normal(ai, np.abs(ai / 5))
+bi = np.random.normal(ai, np.abs(bi / 5))
+
+ecg12 = nk.ecg_simulate(duration=10, method="multileads", ti=ti, ai=ai, bi=bi)
+
+# Visualize results
+ecg12[0:10000].plot(subplots=True)
+```
+![](../../studies/ecg_generating_12_leads/ECG12_abnormal.png)
+
+## Explanation
+
+
+For a better understanding of the above parameters, please read the following.
+
+![](../../studies/ecg_generating_12_leads/table.png)
+
+We can see from the above table that each interval in the 3D trajectory can be fixed by 3 parameters:
+
+- `ti`: The starting position $\theta_i/t_i$ along the circle of each interval in radius
+- `ai`: The amplitude of each spike a
+- `bi`: The width of each wave b. By altering these 3 parameters we can change the shape of the 3D trajectory and thus change the waveform of the resulting ECG.
+
+![](../../studies/ecg_generating_12_leads/3D.png)
+
+**Prebuilt Synthetic Data**
+Prebuilt synthetic data can be found at the following link
+https://drive.google.com/drive/folders/1iqyAlyHAvNWdOvjEGn8Y6C8kItsdZ_GC?usp=sharing
+
+
diff --git a/studies/ecg_generating_12_leads/script.py b/studies/ecg_generating_12_leads/script.py
new file mode 100644
index 0000000000..457e154672
--- /dev/null
+++ b/studies/ecg_generating_12_leads/script.py
@@ -0,0 +1,31 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+import neurokit2 as nk
+
+ecg12 = nk.ecg_simulate(duration=10, method="multileads")
+
+# Visualize results
+ecg12[0:10000].plot(subplots=True)
+plt.savefig("ECG12_normal.png")
+
+# Normal parameters (used by default)
+# ===================================
+# t, the starting position along the circle of each interval in radius
+ti = np.array((-70, -15, 0, 15, 100))
+# a, the amplitude of each spike
+ai = np.array((1.2, -5, 30, -7.5, 0.75))
+# b, the width of each spike
+bi = np.array((0.25, 0.1, 0.1, 0.1, 0.4))
+
+# Add noise
+# ===============
+ti = np.random.normal(ti, np.ones(5) * 3)
+ai = np.random.normal(ai, np.abs(ai / 5))
+bi = np.random.normal(bi, np.abs(bi / 5))
+
+ecg12 = nk.ecg_simulate(duration=10, method="multileads", ti=ti, ai=ai, bi=bi)
+
+# Visualize results
+ecg12[0:10000].plot(subplots=True)
+plt.savefig("ECG12_abnormal.png")
diff --git a/studies/ecg_generating_12_leads/table.png b/studies/ecg_generating_12_leads/table.png
new file mode 100644
index 0000000000..81d1e3de72
Binary files /dev/null and b/studies/ecg_generating_12_leads/table.png differ


Test patch:
diff --git a/tests/tests_eog.py b/tests/tests_eog.py
index 462f16b286..6e607a5c37 100644
--- a/tests/tests_eog.py
+++ b/tests/tests_eog.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
-import numpy as np
-import mne
 import matplotlib.pyplot as plt
+import mne
+import numpy as np
 import pytest
 
 import neurokit2 as nk
@@ -15,10 +15,12 @@ def test_eog_clean():
     assert eog_cleaned.size == eog_signal.size
 
     # test with mne.io.Raw
-    raw = mne.io.read_raw_fif(mne.datasets.sample.data_path() + "/MEG/sample/sample_audvis_raw.fif", preload=True)
+    raw = mne.io.read_raw_fif(
+        mne.datasets.sample.data_path() + "/MEG/sample/sample_audvis_raw.fif", preload=True
+    )
     sampling_rate = raw.info["sfreq"]
 
-    eog_channels = nk.mne_channel_extract(raw, what="EOG", name="EOG")
+    eog_channels = nk.mne_channel_extract(raw, what="EOG", name="EOG").values
     eog_cleaned = nk.eog_clean(eog_channels, sampling_rate, method="agarwal2019")
     assert eog_cleaned.size == eog_channels.size
 
@@ -42,11 +44,13 @@ def test_eog_clean():
 
 def test_eog_findpeaks():
 
-    eog_signal = nk.data('eog_100hz')
+    eog_signal = nk.data("eog_100hz")
     eog_cleaned = nk.eog_clean(eog_signal, sampling_rate=100)
 
     # Test with Neurokit
-    nk_peaks = nk.eog_findpeaks(eog_cleaned, sampling_rate=100, method="neurokit", threshold=0.33, show=False)
+    nk_peaks = nk.eog_findpeaks(
+        eog_cleaned, sampling_rate=100, method="neurokit", threshold=0.33, show=False
+    )
     assert nk_peaks.size == 19
 
     # Test with MNE
@@ -95,21 +99,27 @@ def test_eog_plot():
     np.testing.assert_array_equal(fig.axes[0].get_xticks(), fig.axes[1].get_xticks())
     plt.close(fig)
 
-    with pytest.raises(ValueError, match=r'NeuroKit error: eog_plot.*'):
+    with pytest.raises(ValueError, match=r"NeuroKit error: eog_plot.*"):
         nk.eog_plot(None)
 
 
 def test_eog_eventrelated():
 
-    eog = nk.data('eog_200hz')['vEOG']
+    eog = nk.data("eog_200hz")["vEOG"]
     eog_signals, info = nk.eog_process(eog, sampling_rate=200)
-    epochs = nk.epochs_create(eog_signals, events=[5000, 10000, 15000], epochs_start=-0.1, epochs_end=1.9)
+    epochs = nk.epochs_create(
+        eog_signals, events=[5000, 10000, 15000], epochs_start=-0.1, epochs_end=1.9
+    )
     eog_eventrelated = nk.eog_eventrelated(epochs)
 
     # Test rate features
-    assert np.alltrue(np.array(eog_eventrelated["EOG_Rate_Min"]) < np.array(eog_eventrelated["EOG_Rate_Mean"]))
+    assert np.alltrue(
+        np.array(eog_eventrelated["EOG_Rate_Min"]) < np.array(eog_eventrelated["EOG_Rate_Mean"])
+    )
 
-    assert np.alltrue(np.array(eog_eventrelated["EOG_Rate_Mean"]) < np.array(eog_eventrelated["EOG_Rate_Max"]))
+    assert np.alltrue(
+        np.array(eog_eventrelated["EOG_Rate_Mean"]) < np.array(eog_eventrelated["EOG_Rate_Max"])
+    )
 
     # Test blink presence
     assert np.alltrue(np.array(eog_eventrelated["EOG_Blinks_Presence"]) == np.array([1, 0, 0]))
@@ -130,23 +140,23 @@ def test_eog_eventrelated():
 
 def test_eog_intervalrelated():
 
-    eog = nk.data('eog_200hz')['vEOG']
+    eog = nk.data("eog_200hz")["vEOG"]
     eog_signals, info = nk.eog_process(eog, sampling_rate=200)
 
-    columns = ['EOG_Peaks_N', 'EOG_Rate_Mean']
+    columns = ["EOG_Peaks_N", "EOG_Rate_Mean"]
 
     # Test with signal dataframe
     features = nk.eog_intervalrelated(eog_signals)
 
-    assert all(elem in np.array(features.columns.values, dtype=str) for elem
-               in columns)
+    assert all(elem in np.array(features.columns.values, dtype=str) for elem in columns)
     assert features.shape[0] == 1  # Number of rows
 
     # Test with dict
-    columns.append('Label')
-    epochs = nk.epochs_create(eog_signals, events=[5000, 10000, 15000], epochs_start=-0.1, epochs_end=1.9)
+    columns.append("Label")
+    epochs = nk.epochs_create(
+        eog_signals, events=[5000, 10000, 15000], epochs_start=-0.1, epochs_end=1.9
+    )
     epochs_dict = nk.eog_intervalrelated(epochs)
 
-    assert all(elem in columns for elem
-               in np.array(epochs_dict.columns.values, dtype=str))
+    assert all(elem in columns for elem in np.array(epochs_dict.columns.values, dtype=str))
     assert epochs_dict.shape[0] == len(epochs)  # Number of rows
diff --git a/tests/tests_signal.py b/tests/tests_signal.py
index e8bdac9018..2cb0f252a9 100644
--- a/tests/tests_signal.py
+++ b/tests/tests_signal.py
@@ -8,7 +8,6 @@
 
 import neurokit2 as nk
 
-
 # =============================================================================
 # Signal
 # =============================================================================
@@ -16,11 +15,15 @@
 
 def test_signal_simulate():
     # Warning for nyquist criterion
-    with pytest.warns(nk.misc.NeuroKitWarning, match=r"Skipping requested frequency.*cannot be resolved.*"):
+    with pytest.warns(
+        nk.misc.NeuroKitWarning, match=r"Skipping requested frequency.*cannot be resolved.*"
+    ):
         nk.signal_simulate(sampling_rate=100, frequency=11, silent=False)
 
     # Warning for period duration
-    with pytest.warns(nk.misc.NeuroKitWarning, match=r"Skipping requested frequency.*since its period of.*"):
+    with pytest.warns(
+        nk.misc.NeuroKitWarning, match=r"Skipping requested frequency.*since its period of.*"
+    ):
         nk.signal_simulate(duration=1, frequency=0.1, silent=False)
 
 
@@ -38,7 +41,8 @@ def test_signal_smooth():
 def test_signal_smooth_boxcar():
     signal = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=float)
     np.testing.assert_array_almost_equal(
-        nk.signal_smooth(signal, kernel="boxcar", size=3), [(1 + 1 + 2) / 3, 2, 3, 4, 5, 6, 7, 8, 9, (9 + 10 + 10) / 3]
+        nk.signal_smooth(signal, kernel="boxcar", size=3),
+        [(1 + 1 + 2) / 3, 2, 3, 4, 5, 6, 7, 8, 9, (9 + 10 + 10) / 3],
     )
 
 
@@ -59,14 +63,25 @@ def test_signal_resample():
     downsampled_interpolation = nk.signal_resample(
         signal, method="interpolation", sampling_rate=1000, desired_sampling_rate=500
     )
-    downsampled_numpy = nk.signal_resample(signal, method="numpy", sampling_rate=1000, desired_sampling_rate=500)
-    downsampled_pandas = nk.signal_resample(signal, method="pandas", sampling_rate=1000, desired_sampling_rate=500)
-    downsampled_fft = nk.signal_resample(signal, method="FFT", sampling_rate=1000, desired_sampling_rate=500)
-    downsampled_poly = nk.signal_resample(signal, method="poly", sampling_rate=1000, desired_sampling_rate=500)
+    downsampled_numpy = nk.signal_resample(
+        signal, method="numpy", sampling_rate=1000, desired_sampling_rate=500
+    )
+    downsampled_pandas = nk.signal_resample(
+        signal, method="pandas", sampling_rate=1000, desired_sampling_rate=500
+    )
+    downsampled_fft = nk.signal_resample(
+        signal, method="FFT", sampling_rate=1000, desired_sampling_rate=500
+    )
+    downsampled_poly = nk.signal_resample(
+        signal, method="poly", sampling_rate=1000, desired_sampling_rate=500
+    )
 
     # Upsample
     upsampled_interpolation = nk.signal_resample(
-        downsampled_interpolation, method="interpolation", sampling_rate=500, desired_sampling_rate=1000
+        downsampled_interpolation,
+        method="interpolation",
+        sampling_rate=500,
+        desired_sampling_rate=1000,
     )
     upsampled_numpy = nk.signal_resample(
         downsampled_numpy, method="numpy", sampling_rate=500, desired_sampling_rate=1000
@@ -74,8 +89,12 @@ def test_signal_resample():
     upsampled_pandas = nk.signal_resample(
         downsampled_pandas, method="pandas", sampling_rate=500, desired_sampling_rate=1000
     )
-    upsampled_fft = nk.signal_resample(downsampled_fft, method="FFT", sampling_rate=500, desired_sampling_rate=1000)
-    upsampled_poly = nk.signal_resample(downsampled_poly, method="poly", sampling_rate=500, desired_sampling_rate=1000)
+    upsampled_fft = nk.signal_resample(
+        downsampled_fft, method="FFT", sampling_rate=500, desired_sampling_rate=1000
+    )
+    upsampled_poly = nk.signal_resample(
+        downsampled_poly, method="poly", sampling_rate=500, desired_sampling_rate=1000
+    )
 
     # Check
     rez = pd.DataFrame(
@@ -128,7 +147,9 @@ def test_signal_filter():
     powerline = np.sin(2 * np.pi * 50 * (samples / sampling_rate))
 
     signal_corrupted = signal + powerline
-    signal_clean = nk.signal_filter(signal_corrupted, sampling_rate=sampling_rate, method="powerline")
+    signal_clean = nk.signal_filter(
+        signal_corrupted, sampling_rate=sampling_rate, method="powerline"
+    )
 
     # import matplotlib.pyplot as plt
     # figure, (ax0, ax1, ax2) = plt.subplots(nrows=3, ncols=1, sharex=True)
@@ -155,7 +176,9 @@ def test_signal_findpeaks():
     signal1 = np.cos(np.linspace(start=0, stop=30, num=1000))
     info1 = nk.signal_findpeaks(signal1)
 
-    signal2 = np.concatenate([np.arange(0, 20, 0.1), np.arange(17, 30, 0.1), np.arange(30, 10, -0.1)])
+    signal2 = np.concatenate(
+        [np.arange(0, 20, 0.1), np.arange(17, 30, 0.1), np.arange(30, 10, -0.1)]
+    )
     info2 = nk.signal_findpeaks(signal2)
     assert len(info1["Peaks"]) > len(info2["Peaks"])
 
@@ -187,15 +210,23 @@ def test_signal_rate():  # since singal_rate wraps signal_period, the latter is
     duration = 120
     sampling_rate = 1000
     rsp = nk.rsp_simulate(
-        duration=duration, sampling_rate=sampling_rate, respiratory_rate=15, method="sinuosoidal", noise=0
+        duration=duration,
+        sampling_rate=sampling_rate,
+        respiratory_rate=15,
+        method="sinuosoidal",
+        noise=0,
     )
     rsp_cleaned = nk.rsp_clean(rsp, sampling_rate=sampling_rate)
     signals, info = nk.rsp_peaks(rsp_cleaned)
-    rate = nk.signal_rate(signals, sampling_rate=sampling_rate, desired_length=duration * sampling_rate)
+    rate = nk.signal_rate(
+        signals, sampling_rate=sampling_rate, desired_length=duration * sampling_rate
+    )
     assert rate.shape == (signals.shape[0],)
 
     # Test with dictionary.produced from rsp_findpeaks.
-    rate = nk.signal_rate(info, sampling_rate=sampling_rate, desired_length=duration * sampling_rate)
+    rate = nk.signal_rate(
+        info, sampling_rate=sampling_rate, desired_length=duration * sampling_rate
+    )
     assert rate.shape == (duration * sampling_rate,)
 
 
@@ -263,7 +294,9 @@ def test_signal_power():
 
 def test_signal_timefrequency():
 
-    signal = nk.signal_simulate(duration=50, frequency=5) + 2 * nk.signal_simulate(duration=50, frequency=20)
+    signal = nk.signal_simulate(duration=50, frequency=5) + 2 * nk.signal_simulate(
+        duration=50, frequency=20
+    )
 
     # short-time fourier transform
     frequency, time, stft = nk.signal_timefrequency(
@@ -277,7 +310,9 @@ def test_signal_timefrequency():
     assert np.sum(stft[indices_freq5]) < np.sum(stft[indices_freq20])
 
     # wavelet transform
-    frequency, time, cwtm = nk.signal_timefrequency(signal, method="cwt", max_frequency=50, show=False)
+    frequency, time, cwtm = nk.signal_timefrequency(
+        signal, method="cwt", max_frequency=50, show=False
+    )
 
     assert len(frequency) == cwtm.shape[0]
     assert len(time) == cwtm.shape[1]
@@ -286,7 +321,9 @@ def test_signal_timefrequency():
     assert np.sum(cwtm[indices_freq5]) < np.sum(cwtm[indices_freq20])
 
     # wvd
-    frequency, time, wvd = nk.signal_timefrequency(signal, method="wvd", max_frequency=50, show=False)
+    frequency, time, wvd = nk.signal_timefrequency(
+        signal, method="wvd", max_frequency=50, show=False
+    )
     assert len(frequency) == wvd.shape[0]
     assert len(time) == wvd.shape[1]
     indices_freq5 = np.logical_and(frequency > 3, frequency < 7)
@@ -294,7 +331,9 @@ def test_signal_timefrequency():
     assert np.sum(wvd[indices_freq5]) < np.sum(wvd[indices_freq20])
 
     # pwvd
-    frequency, time, pwvd = nk.signal_timefrequency(signal, method="pwvd", max_frequency=50, show=False)
+    frequency, time, pwvd = nk.signal_timefrequency(
+        signal, method="pwvd", max_frequency=50, show=False
+    )
     assert len(frequency) == pwvd.shape[0]
     assert len(time) == pwvd.shape[1]
     indices_freq5 = np.logical_and(frequency > 3, frequency < 7)
@@ -310,14 +349,21 @@ def test_signal_psd(recwarn):
 
     assert list(out.columns) == ["Frequency", "Power"]
 
+
 def test_signal_distort():
     signal = nk.signal_simulate(duration=10, frequency=0.5, sampling_rate=10)
 
     # Warning for nyquist criterion
-    with pytest.warns(nk.misc.NeuroKitWarning, match=r"Skipping requested noise frequency.*cannot be resolved.*"):
+    with pytest.warns(
+        nk.misc.NeuroKitWarning, match=r"Skipping requested noise frequency.*cannot be resolved.*"
+    ):
         nk.signal_distort(signal, sampling_rate=10, noise_amplitude=1, silent=False)
 
     # Warning for period duration
-    with pytest.warns(nk.misc.NeuroKitWarning, match=r"Skipping requested noise frequency.*since its period of.*"):
+    with pytest.warns(
+        nk.misc.NeuroKitWarning, match=r"Skipping requested noise frequency.*since its period of.*"
+    ):
         signal = nk.signal_simulate(duration=1, frequency=1, sampling_rate=10)
         nk.signal_distort(signal, noise_amplitude=1, noise_frequency=0.1, silent=False)
+
+    signal2 = nk.signal_simulate(duration=10, frequency=0.5, sampling_rate=10)