From 78ec406598721c8ceba56179426c05054a91e68e Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 12:37:03 +0200
Subject: [PATCH 01/11] FIX: preserve input intact when normalizing variance

---
 src/meegsim/utils.py | 13 +++++++------
 tests/test_utils.py  |  7 ++++---
 2 files changed, 11 insertions(+), 9 deletions(-)

diff --git a/src/meegsim/utils.py b/src/meegsim/utils.py
index f178456..ac3f4b8 100644
--- a/src/meegsim/utils.py
+++ b/src/meegsim/utils.py
@@ -79,15 +79,16 @@ def normalize_variance(data):
 
     Returns
     -------
-    data: array
+    data_norm: array
         Normalized time series. The variance of each row is equal to 1.
     """
+    # NOTE: make a copy to keep the original waveform intact
+    data_norm = data.copy()
+    if data_norm.ndim == 1:
+        return data_norm / np.std(data_norm)
 
-    if data.ndim == 1:
-        return data / np.std(data)
-
-    data /= np.std(data, axis=-1)[:, np.newaxis]
-    return data
+    data_norm /= np.std(data_norm, axis=-1)[:, np.newaxis]
+    return data_norm
 
 
 def _extract_hemi(src):
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 35c206b..5714f8f 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -96,11 +96,12 @@ def test_combine_stcs_overlap():
 
 def test_normalize_variance():
     data = np.random.randn(10, 1000)
+    data_orig = data.copy()
     normalized = normalize_variance(data)
 
-    # Should not change the shape but should change the norm
-    assert data.shape == normalized.shape
-    assert np.allclose(np.var(normalized, axis=1), 1)
+    assert data.shape == normalized.shape, "Array shape should not be changed"
+    assert np.allclose(np.var(normalized, axis=1), 1), "Expected variance to be 1"
+    assert np.allclose(data_orig, data), "The input waveform should not be modified"
 
 
 def test_get_sfreq():

From bc03989e26c0c59acea19e2133c436b3434ae09e Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 14:08:15 +0200
Subject: [PATCH 02/11] ENH: control amplitude envelope of the coupled waveform

---
 src/meegsim/coupling.py | 109 ++++++++++++++++++++++++++++++++--------
 1 file changed, 88 insertions(+), 21 deletions(-)

diff --git a/src/meegsim/coupling.py b/src/meegsim/coupling.py
index 81ba690..4b24c80 100644
--- a/src/meegsim/coupling.py
+++ b/src/meegsim/coupling.py
@@ -7,13 +7,46 @@
 from scipy.stats import vonmises
 from scipy.signal import butter, filtfilt, hilbert
 
-from meegsim._check import check_numeric
+from meegsim._check import check_numeric, check_option
 from meegsim.snr import get_variance, amplitude_adjustment_factor
 from meegsim.utils import normalize_variance
 from meegsim.waveform import narrowband_oscillation
 
 
-def constant_phase_shift(waveform, sfreq, phase_lag, m=1, n=1, random_state=None):
+def _get_envelope(waveform, envelope, sfreq, fmin=None, fmax=None, random_state=None):
+    check_option(
+        "the amplitude envelope of the coupled waveform", envelope, ["same", "random"]
+    )
+
+    if envelope == "same":
+        return np.abs(waveform)
+
+    if fmin is None or fmax is None:
+        raise ValueError(
+            "Frequency limits are required for generating the envelope of the coupled waveform"
+        )
+    times = np.arange(waveform.size) / sfreq
+    random_waveform = narrowband_oscillation(
+        1, times, fmin=fmin, fmax=fmax, random_state=random_state
+    )
+    random_waveform = hilbert(random_waveform)
+
+    # TODO: here we could also mix original and random envelope with different
+    # values of SNR to achieve smooth control over the resulting envelope correlation
+    return np.abs(random_waveform)
+
+
+def ppc_constant_phase_shift(
+    waveform,
+    sfreq,
+    phase_lag,
+    envelope="random",
+    fmin=None,
+    fmax=None,
+    m=1,
+    n=1,
+    random_state=None,
+):
     """
     Generate a time series that is phase coupled to the input time series with
     a constant phase lag.
@@ -32,21 +65,30 @@ def constant_phase_shift(waveform, sfreq, phase_lag, m=1, n=1, random_state=None
         The input signal to be processed. It can be a real or complex time series.
 
     sfreq : float
-        Sampling frequency of the signal, in Hz. This argument is not used in this
-        function but is accepted for consistency with other coupling methods.
+        Sampling frequency of the signal, in Hz.
 
     phase_lag : float
         Constant phase lag to apply to the waveform in radians.
 
-    m : int, optional
+    envelope : {"same", "random"}
+        Controls the amplitude envelope of the coupled waveform to be either randomly
+        generated or to be the same as the envelope of the input waveform.
+
+    fmin : float, optional
+        Lower cutoff frequency for the oscillation.
+
+    fmax : float, optional
+        Upper cutoff frequency for the oscillation.
+
+    m : float, optional
         Multiplier for the base frequency of the output oscillation, default is 1.
 
-    n : int, optional
+    n : float, optional
         Multiplier for the base frequency of the input oscillation, default is 1.
 
     random_state : None, optional
-        This parameter is accepted for consistency with other coupling functions
-        but not used since no randomness is involved.
+        Random state can be fixed to provide reproducible results. Otherwise, the
+        results may differ between function calls.
 
     Returns
     -------
@@ -55,17 +97,37 @@ def constant_phase_shift(waveform, sfreq, phase_lag, m=1, n=1, random_state=None
     """
     if not np.iscomplexobj(waveform):
         waveform = hilbert(waveform)
-
-    waveform_amp = np.abs(waveform)
     waveform_angle = np.angle(waveform)
-    waveform_coupled = waveform_amp * np.exp(
-        1j * m / n * waveform_angle + 1j * phase_lag
+
+    waveform_amp = _get_envelope(waveform, envelope, sfreq, fmin, fmax, random_state)
+    waveform_coupled = np.real(
+        waveform_amp * np.exp(1j * m / n * waveform_angle + 1j * phase_lag)
+    )
+    if envelope == "same":
+        return normalize_variance(waveform_coupled)
+
+    # NOTE: if the envelope was modified, we filter the result again in the target
+    # frequency range to suppress distortions due to merging amplitude envelope and
+    # phase from different time series
+    b, a = butter(
+        N=2, Wn=np.array([m / n * fmin, m / n * fmax]) / sfreq * 2, btype="bandpass"
     )
-    return normalize_variance(np.real(waveform_coupled))
+    waveform_coupled = filtfilt(b, a, waveform_coupled)
+
+    return normalize_variance(waveform_coupled)
 
 
 def ppc_von_mises(
-    waveform, sfreq, phase_lag, kappa, fmin, fmax, m=1, n=1, random_state=None
+    waveform,
+    sfreq,
+    phase_lag,
+    kappa,
+    fmin,
+    fmax,
+    envelope="random",
+    m=1,
+    n=1,
+    random_state=None,
 ):
     """
     Generate a time series that is phase coupled to the input time series with
@@ -102,6 +164,10 @@ def ppc_von_mises(
     fmax: float
         Upper cutoff frequency of the base frequency harmonic (in Hz).
 
+    envelope : {"same", "random"}
+        Controls the amplitude envelope of the coupled waveform to be either randomly
+        generated or to be the same as the envelope of the input waveform.
+
     m : int, optional
         Multiplier for the base frequency of the output oscillation, default is 1.
 
@@ -121,23 +187,22 @@ def ppc_von_mises(
     if not np.iscomplexobj(waveform):
         waveform = hilbert(waveform)
 
-    waveform_amp = np.abs(waveform)
+    waveform_amp = _get_envelope(waveform, envelope, sfreq, fmin, fmax, random_state)
     waveform_angle = np.angle(waveform)
-    n_samples = len(waveform)
+    n_samples = waveform.size
 
     ph_distr = vonmises.rvs(
         kappa, loc=phase_lag, size=n_samples, random_state=random_state
     )
-    tmp_waveform = np.real(
+    waveform_coupled = np.real(
         waveform_amp * np.exp(1j * m / n * waveform_angle + 1j * ph_distr)
     )
     b, a = butter(
         N=2, Wn=np.array([m / n * fmin, m / n * fmax]) / sfreq * 2, btype="bandpass"
     )
-    tmp_waveform = filtfilt(b, a, tmp_waveform)
-    waveform_coupled = waveform_amp * np.exp(1j * np.angle(hilbert(tmp_waveform)))
+    waveform_coupled = filtfilt(b, a, waveform_coupled)
 
-    return normalize_variance(np.real(waveform_coupled))
+    return normalize_variance(waveform_coupled)
 
 
 def _shifted_copy_with_noise(waveform, sfreq, phase_lag, snr, fmin, fmax, random_state):
@@ -146,7 +211,9 @@ def _shifted_copy_with_noise(waveform, sfreq, phase_lag, snr, fmin, fmax, random
     waveform and (2) mixing it with noise to achieve a desired level of signal-to-noise
     ratio, which determines the resulting phase-phase and amplitude-amplitude coupling.
     """
-    shifted_waveform = constant_phase_shift(waveform, sfreq, phase_lag)
+    shifted_waveform = ppc_constant_phase_shift(
+        waveform, sfreq, phase_lag, envelope="same"
+    )
     signal_var = get_variance(shifted_waveform, sfreq, fmin, fmax, filter=True)
 
     # NOTE: we use another randomly generated narrowband oscillation as noise here

From e27a941e18bc3cf3ce5b068e5374973c474b0c19 Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 14:09:34 +0200
Subject: [PATCH 03/11] TEST: fix and extend tests for ppc_constant_phase_shift

---
 tests/test_coupling.py    | 52 +++++++++++++++++++++++++++++++--------
 tests/test_integration.py |  4 +--
 2 files changed, 44 insertions(+), 12 deletions(-)

diff --git a/tests/test_coupling.py b/tests/test_coupling.py
index 9d7b34d..19bf88e 100644
--- a/tests/test_coupling.py
+++ b/tests/test_coupling.py
@@ -6,7 +6,7 @@
 from scipy.signal import hilbert
 
 from meegsim.coupling import (
-    constant_phase_shift,
+    ppc_constant_phase_shift,
     ppc_von_mises,
     ppc_shifted_copy_with_noise,
     _get_required_snr,
@@ -17,22 +17,23 @@
 from utils.prepare import prepare_sinusoid
 
 
-def prepare_inputs():
+def prepare_inputs(sfreq=250, duration=60):
     n_series = 2
-    fs = 1000
-    times = np.arange(0, 1, 1 / fs)
+    times = np.arange(0, duration, 1 / sfreq)
     return n_series, len(times), times
 
 
 @pytest.mark.parametrize(
     "phase_lag", [np.pi / 4, np.pi / 3, np.pi / 2, np.pi, 2 * np.pi]
 )
-def test_constant_phase_shift(phase_lag):
+def test_ppc_constant_phase_shift_same_envelope(phase_lag):
     # Test with a simple sinusoidal waveform
     _, _, times = prepare_inputs()
     waveform = np.sin(2 * np.pi * 10 * times)
 
-    result = constant_phase_shift(waveform, get_sfreq(times), phase_lag)
+    result = ppc_constant_phase_shift(
+        waveform, get_sfreq(times), phase_lag, envelope="same"
+    )
 
     waveform = hilbert(waveform)
     result = hilbert(result)
@@ -41,20 +42,47 @@ def test_constant_phase_shift(phase_lag):
     plv = np.abs(cplv)
     test_angle = np.angle(cplv)
 
-    assert plv >= 0.99, f"Test failed: plv is smaller than 0.99. plv = {plv}"
+    assert plv >= 0.95, f"Expected PLV to be at least 0.95, got {plv}"
+    assert (
+        (np.abs(test_angle) - phase_lag) <= 0.01
+    ), f"Test failed: angle is different from phase_lag. difference = {np.round((np.abs(test_angle) - phase_lag),2)}"
+
+
+@pytest.mark.parametrize(
+    "phase_lag", [np.pi / 4, np.pi / 3, np.pi / 2, np.pi, 2 * np.pi]
+)
+def test_ppc_constant_phase_shift_random_envelope(phase_lag):
+    # Test with a simple sinusoidal waveform
+    _, _, times = prepare_inputs()
+    waveform = np.sin(2 * np.pi * 10 * times)
+
+    result = ppc_constant_phase_shift(
+        waveform, get_sfreq(times), phase_lag, envelope="random", fmin=9.5, fmax=10.5
+    )
+
+    waveform = hilbert(waveform)
+    result = hilbert(result)
+
+    cplv = compute_plv(waveform, result, m=1, n=1, plv_type="complex")
+    plv = np.abs(cplv)
+    test_angle = np.angle(cplv)
+
+    assert plv >= 0.95, f"Expected PLV to be at least 0.95, got {plv}"
     assert (
         (np.abs(test_angle) - phase_lag) <= 0.01
     ), f"Test failed: angle is different from phase_lag. difference = {np.round((np.abs(test_angle) - phase_lag),2)}"
 
 
 @pytest.mark.parametrize("m, n", [(2, 1), (3, 1), (5 / 2, 1)])
-def test_constant_phase_shift_harmonics(m, n):
+def test_ppc_constant_phase_shift_harmonics_same_envelope(m, n):
     # Test with different m and n harmonics
     _, _, times = prepare_inputs()
     waveform = np.sin(2 * np.pi * 10 * times)
     phase_lag = np.pi / 3
 
-    result = constant_phase_shift(waveform, get_sfreq(times), phase_lag, m=m, n=n)
+    result = ppc_constant_phase_shift(
+        waveform, get_sfreq(times), phase_lag, envelope="same", m=m, n=n
+    )
 
     waveform = hilbert(waveform)
     result = hilbert(result)
@@ -63,7 +91,7 @@ def test_constant_phase_shift_harmonics(m, n):
     plv = np.abs(cplv)
     test_angle = np.angle(cplv)
 
-    assert plv >= 0.9, f"Test failed: plv is smaller than 0.9. plv = {plv}"
+    assert plv >= 0.9, f"Expected PLV to be at least 0.9, got {plv}"
     assert (
         (np.abs(test_angle) - phase_lag) <= 0.1
     ), f"Test failed: angle is different from phase_lag. difference = {np.round((np.abs(test_angle) - phase_lag),2)}"
@@ -120,6 +148,10 @@ def test_ppc_von_mises_harmonics(m, n):
 @pytest.mark.parametrize(
     "coupling_fun,params",
     [
+        (
+            ppc_constant_phase_shift,
+            dict(phase_lag=np.pi / 4, envelope="random", fmin=8, fmax=12),
+        ),
         (ppc_von_mises, dict(kappa=1, phase_lag=np.pi / 4, fmin=8, fmax=12)),
         (
             ppc_shifted_copy_with_noise,
diff --git a/tests/test_integration.py b/tests/test_integration.py
index afca711..592cab0 100644
--- a/tests/test_integration.py
+++ b/tests/test_integration.py
@@ -4,7 +4,7 @@
 
 import numpy as np
 
-from meegsim.coupling import constant_phase_shift, ppc_von_mises
+from meegsim.coupling import ppc_constant_phase_shift, ppc_von_mises
 from meegsim.location import select_random
 from meegsim.simulate import SourceSimulator
 from meegsim.waveform import narrowband_oscillation, white_noise
@@ -62,7 +62,7 @@ def test_builtin_methods():
     )
 
     # Define several edges to test graph traversal and built-in coupling methods
-    sim.set_coupling(("point1", "point2"), method=constant_phase_shift, phase_lag=0)
+    sim.set_coupling(("point1", "point2"), method=ppc_constant_phase_shift, phase_lag=0)
     sim.set_coupling(
         coupling={
             ("point2", "patch3"): dict(kappa=0.1, phase_lag=-np.pi / 6),

From 36940be13422d787b4a9993a7599ab4379ec7681 Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 14:10:38 +0200
Subject: [PATCH 04/11] TEST: fix the integration test

---
 tests/test_integration.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/tests/test_integration.py b/tests/test_integration.py
index 592cab0..c9f831e 100644
--- a/tests/test_integration.py
+++ b/tests/test_integration.py
@@ -62,7 +62,13 @@ def test_builtin_methods():
     )
 
     # Define several edges to test graph traversal and built-in coupling methods
-    sim.set_coupling(("point1", "point2"), method=ppc_constant_phase_shift, phase_lag=0)
+    sim.set_coupling(
+        ("point1", "point2"),
+        method=ppc_constant_phase_shift,
+        phase_lag=0,
+        fmin=8,
+        fmax=12,
+    )
     sim.set_coupling(
         coupling={
             ("point2", "patch3"): dict(kappa=0.1, phase_lag=-np.pi / 6),

From 865ddf23eeffb1de350ff4ec36d8876b0996125e Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 14:13:20 +0200
Subject: [PATCH 05/11] DOC: update the method name in autosummary

---
 docs/api/coupling.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/api/coupling.rst b/docs/api/coupling.rst
index e6cd55b..d111368 100644
--- a/docs/api/coupling.rst
+++ b/docs/api/coupling.rst
@@ -6,6 +6,6 @@ Coupling methods
 .. autosummary::
    :toctree: ../generated/
 
-   ppc_shifted_copy_with_noise
+   ppc_constant_phase_shift
    ppc_von_mises
-   constant_phase_shift
+   ppc_shifted_copy_with_noise

From 382068a005fef2e64ad3eb11cb143a79cf03e247 Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 14:40:25 +0200
Subject: [PATCH 06/11] TEST: rework tests for the von Mises function

---
 src/meegsim/utils.py   |   5 --
 tests/test_coupling.py | 123 +++++++++++++++++++++++++++++++++--------
 2 files changed, 101 insertions(+), 27 deletions(-)

diff --git a/src/meegsim/utils.py b/src/meegsim/utils.py
index ac3f4b8..a46d305 100644
--- a/src/meegsim/utils.py
+++ b/src/meegsim/utils.py
@@ -3,7 +3,6 @@
 import warnings
 
 from mne.io.constants import FIFF
-from scipy.special import i1, i0
 
 
 logger = logging.getLogger("meegsim")
@@ -196,10 +195,6 @@ def unpack_vertices(vertices_lists):
     return unpacked_vertices
 
 
-def theoretical_plv(kappa):
-    return i1(kappa) / i0(kappa)
-
-
 def vertices_to_mne(vertices, src):
     """
     Convert the vertices to the MNE format (list of lists).
diff --git a/tests/test_coupling.py b/tests/test_coupling.py
index 19bf88e..c345e98 100644
--- a/tests/test_coupling.py
+++ b/tests/test_coupling.py
@@ -12,7 +12,7 @@
     _get_required_snr,
     _shifted_copy_with_noise,
 )
-from meegsim.utils import get_sfreq, theoretical_plv
+from meegsim.utils import get_sfreq
 
 from utils.prepare import prepare_sinusoid
 
@@ -32,7 +32,7 @@ def test_ppc_constant_phase_shift_same_envelope(phase_lag):
     waveform = np.sin(2 * np.pi * 10 * times)
 
     result = ppc_constant_phase_shift(
-        waveform, get_sfreq(times), phase_lag, envelope="same"
+        waveform, get_sfreq(times), phase_lag, envelope="same", random_state=1234
     )
 
     waveform = hilbert(waveform)
@@ -42,7 +42,7 @@ def test_ppc_constant_phase_shift_same_envelope(phase_lag):
     plv = np.abs(cplv)
     test_angle = np.angle(cplv)
 
-    assert plv >= 0.95, f"Expected PLV to be at least 0.95, got {plv}"
+    assert plv >= 0.9, f"Expected PLV to be at least 0.9, got {plv}"
     assert (
         (np.abs(test_angle) - phase_lag) <= 0.01
     ), f"Test failed: angle is different from phase_lag. difference = {np.round((np.abs(test_angle) - phase_lag),2)}"
@@ -57,7 +57,13 @@ def test_ppc_constant_phase_shift_random_envelope(phase_lag):
     waveform = np.sin(2 * np.pi * 10 * times)
 
     result = ppc_constant_phase_shift(
-        waveform, get_sfreq(times), phase_lag, envelope="random", fmin=9.5, fmax=10.5
+        waveform,
+        get_sfreq(times),
+        phase_lag,
+        envelope="random",
+        fmin=9.5,
+        fmax=10.5,
+        random_state=1234,
     )
 
     waveform = hilbert(waveform)
@@ -67,7 +73,7 @@ def test_ppc_constant_phase_shift_random_envelope(phase_lag):
     plv = np.abs(cplv)
     test_angle = np.angle(cplv)
 
-    assert plv >= 0.95, f"Expected PLV to be at least 0.95, got {plv}"
+    assert plv >= 0.9, f"Expected PLV to be at least 0.9, got {plv}"
     assert (
         (np.abs(test_angle) - phase_lag) <= 0.01
     ), f"Test failed: angle is different from phase_lag. difference = {np.round((np.abs(test_angle) - phase_lag),2)}"
@@ -81,7 +87,13 @@ def test_ppc_constant_phase_shift_harmonics_same_envelope(m, n):
     phase_lag = np.pi / 3
 
     result = ppc_constant_phase_shift(
-        waveform, get_sfreq(times), phase_lag, envelope="same", m=m, n=n
+        waveform,
+        get_sfreq(times),
+        phase_lag,
+        envelope="same",
+        m=m,
+        n=n,
+        random_state=1234,
     )
 
     waveform = hilbert(waveform)
@@ -97,15 +109,31 @@ def test_ppc_constant_phase_shift_harmonics_same_envelope(m, n):
     ), f"Test failed: angle is different from phase_lag. difference = {np.round((np.abs(test_angle) - phase_lag),2)}"
 
 
-@pytest.mark.parametrize("kappa", [0.001, 0.1, 0.5, 1, 5, 10, 50])
-def test_ppc_von_mises(kappa):
-    # Test kappas that are reliable (more than 0.5)
-    _, _, times = prepare_inputs()
+@pytest.mark.parametrize(
+    "kappa,lo,hi",
+    [
+        (0.001, 0, 0.2),
+        (0.01, 0, 0.2),
+        (0.1, 0.1, 0.4),
+        (0.3, 0.3, 0.7),
+        (0.5, 0.5, 0.9),
+        (1, 0.7, 1),
+        (10, 0.9, 1),
+    ],
+)
+def test_ppc_von_mises_same_envelope_kappa(kappa, lo, hi):
+    _, _, times = prepare_inputs(duration=120)
     waveform = np.sin(2 * np.pi * 10 * times)
-    phase_lag = 0
+    phase_lag = np.pi / 4
 
     result = ppc_von_mises(
-        waveform, get_sfreq(times), phase_lag, kappa=kappa, fmin=8, fmax=12
+        waveform,
+        get_sfreq(times),
+        phase_lag,
+        kappa=kappa,
+        fmin=8,
+        fmax=12,
+        random_state=1234,
     )
 
     waveform = hilbert(waveform)
@@ -113,11 +141,52 @@ def test_ppc_von_mises(kappa):
 
     cplv = compute_plv(waveform, result, m=1, n=1, plv_type="complex")
     plv = np.abs(cplv)
-    plv_theoretical = theoretical_plv(kappa)
+    test_angle = np.abs(np.angle(cplv))
+
+    # NOTE: lower and upper bounds were selected by simulating multiple time series,
+    # this test should prevent large deviations from the expected result due to
+    # errors in the processing
+    assert lo <= plv <= hi, f"Expected PLV to be between {lo} and {hi}, got {plv}"
+    if kappa >= 0.5:
+        assert np.allclose(test_angle, phase_lag, atol=0.1), (
+            f"Expected the actual phase lag ({test_angle:.2f}) to be within "
+            f"0.1 from the desired one ({phase_lag:.2f})"
+        )
+
+
+@pytest.mark.parametrize("kappa,lo,hi", [(0.01, 0, 0.2), (10, 0.8, 1)])
+def test_ppc_von_mises_random_envelope_kappa(kappa, lo, hi):
+    _, _, times = prepare_inputs(duration=120)
+    waveform = np.sin(2 * np.pi * 10 * times)
+    phase_lag = np.pi / 4
 
-    assert (
-        plv >= plv_theoretical
-    ), f"Test failed: plv is smaller than theoretical. plv = {plv}, plv_theoretical = {plv_theoretical}"
+    result = ppc_von_mises(
+        waveform,
+        get_sfreq(times),
+        phase_lag,
+        kappa=kappa,
+        envelope="random",
+        fmin=8,
+        fmax=12,
+        random_state=1234,
+    )
+
+    waveform = hilbert(waveform)
+    result = hilbert(result)
+
+    cplv = compute_plv(waveform, result, m=1, n=1, plv_type="complex")
+    plv = np.abs(cplv)
+    test_angle = np.abs(np.angle(cplv))
+
+    # NOTE: we check only extreme cases here to make sure the random envelope does
+    # not break the result completely. Still, random envelope might decrease PLV a bit,
+    # so the bounds are a bit wider
+    assert lo <= plv <= hi, f"Expected PLV to be between {lo} and {hi}, got {plv}"
+    if kappa >= 0.5:
+        assert np.allclose(np.abs(test_angle), phase_lag, atol=0.1), (
+            f"Expected the actual phase lag ({test_angle:.2f}) to be within "
+            f"0.1 from the desired one ({phase_lag:.2f})"
+        )
 
 
 @pytest.mark.parametrize("m, n", [(2, 1), (3, 1), (5 / 2, 1)])
@@ -125,11 +194,20 @@ def test_ppc_von_mises_harmonics(m, n):
     # Test with different m and n harmonics
     _, _, times = prepare_inputs()
     waveform = np.sin(2 * np.pi * 10 * times)
-    phase_lag = 0
+    phase_lag = np.pi / 4
     kappa = 10
 
     result = ppc_von_mises(
-        waveform, get_sfreq(times), phase_lag, m=m, n=n, kappa=kappa, fmin=8, fmax=12
+        waveform,
+        get_sfreq(times),
+        phase_lag,
+        m=m,
+        n=n,
+        envelope="same",
+        kappa=kappa,
+        fmin=8,
+        fmax=12,
+        random_state=1234,
     )
 
     waveform = hilbert(waveform)
@@ -139,10 +217,11 @@ def test_ppc_von_mises_harmonics(m, n):
     plv = np.abs(cplv)
     test_angle = np.angle(cplv)
 
-    assert plv >= 0.8, f"Test failed: plv is smaller than 0.8. plv = {plv}"
-    assert (
-        (np.abs(test_angle) - phase_lag) <= 0.1
-    ), f"Test failed: angle is different from phase_lag. difference = {np.round((np.abs(test_angle) - phase_lag),2)}"
+    assert 0.7 <= plv, f"Expected PLV to be at least 0.7, got {plv}"
+    assert np.allclose(np.abs(test_angle), phase_lag, atol=0.1), (
+        f"Expected the actual phase lag ({test_angle:.2f}) to be within "
+        f"0.1 from the desired one ({phase_lag:.2f})"
+    )
 
 
 @pytest.mark.parametrize(

From 69c6013e3ecdee0ccdb7c6375038be9386bc566d Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 15:00:51 +0200
Subject: [PATCH 07/11] TEST: tests for _get_envelope

---
 src/meegsim/coupling.py |  8 +++++---
 tests/test_coupling.py  | 28 ++++++++++++++++++++++++++++
 2 files changed, 33 insertions(+), 3 deletions(-)

diff --git a/src/meegsim/coupling.py b/src/meegsim/coupling.py
index 4b24c80..0680858 100644
--- a/src/meegsim/coupling.py
+++ b/src/meegsim/coupling.py
@@ -17,6 +17,8 @@ def _get_envelope(waveform, envelope, sfreq, fmin=None, fmax=None, random_state=
     check_option(
         "the amplitude envelope of the coupled waveform", envelope, ["same", "random"]
     )
+    if not np.iscomplexobj(waveform):
+        waveform = hilbert(waveform)
 
     if envelope == "same":
         return np.abs(waveform)
@@ -40,9 +42,9 @@ def ppc_constant_phase_shift(
     waveform,
     sfreq,
     phase_lag,
-    envelope="random",
     fmin=None,
     fmax=None,
+    envelope="random",
     m=1,
     n=1,
     random_state=None,
@@ -70,7 +72,7 @@ def ppc_constant_phase_shift(
     phase_lag : float
         Constant phase lag to apply to the waveform in radians.
 
-    envelope : {"same", "random"}
+    envelope : str, {"same", "random"}
         Controls the amplitude envelope of the coupled waveform to be either randomly
         generated or to be the same as the envelope of the input waveform.
 
@@ -164,7 +166,7 @@ def ppc_von_mises(
     fmax: float
         Upper cutoff frequency of the base frequency harmonic (in Hz).
 
-    envelope : {"same", "random"}
+    envelope : str, {"same", "random"}
         Controls the amplitude envelope of the coupled waveform to be either randomly
         generated or to be the same as the envelope of the input waveform.
 
diff --git a/tests/test_coupling.py b/tests/test_coupling.py
index c345e98..1bca367 100644
--- a/tests/test_coupling.py
+++ b/tests/test_coupling.py
@@ -9,6 +9,7 @@
     ppc_constant_phase_shift,
     ppc_von_mises,
     ppc_shifted_copy_with_noise,
+    _get_envelope,
     _get_required_snr,
     _shifted_copy_with_noise,
 )
@@ -273,6 +274,33 @@ def test_get_required_snr(coh, expected_snr):
     assert np.isclose(_get_required_snr(coh), expected_snr)
 
 
+def test_get_envelope_same():
+    sfreq = 250
+    waveform = prepare_sinusoid(f=10, sfreq=sfreq, duration=60)
+    waveform_amp = np.abs(hilbert(waveform))
+    envelope = _get_envelope(waveform, envelope="same", sfreq=sfreq)
+    assert np.allclose(waveform_amp, envelope), "Expected envelope to match input"
+
+
+def test_get_envelope_random():
+    sfreq = 250
+    fmin, fmax = 8, 12
+    seed = 1234
+    waveform = prepare_sinusoid(f=10, sfreq=sfreq, duration=60)
+    waveform_amp = np.abs(hilbert(waveform))
+    envelope = _get_envelope(
+        waveform,
+        envelope="random",
+        sfreq=sfreq,
+        fmin=fmin,
+        fmax=fmax,
+        random_state=seed,
+    )
+    assert not np.allclose(
+        waveform_amp, envelope
+    ), "Expected envelope not to match input"
+
+
 @pytest.mark.parametrize(
     "target_coh,tol",
     [

From 4d0e2d89af18210ced8fcf5c3e40d99c5cf48368 Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 15:12:02 +0200
Subject: [PATCH 08/11] Docstrings, comments

---
 src/meegsim/coupling.py | 24 +++++++++++++++---------
 1 file changed, 15 insertions(+), 9 deletions(-)

diff --git a/src/meegsim/coupling.py b/src/meegsim/coupling.py
index 0680858..aaab7f0 100644
--- a/src/meegsim/coupling.py
+++ b/src/meegsim/coupling.py
@@ -77,10 +77,12 @@ def ppc_constant_phase_shift(
         generated or to be the same as the envelope of the input waveform.
 
     fmin : float, optional
-        Lower cutoff frequency for the oscillation.
+        Lower cutoff frequency for the oscillation that gives rise to the random
+        amplitude envelope (only if the ``envelope`` is set to ``"random"``).
 
     fmax : float, optional
-        Upper cutoff frequency for the oscillation.
+        Upper cutoff frequency for the oscillation that gives rise to the random
+        amplitude envelope (only if the ``envelope`` is set to ``"random"``).
 
     m : float, optional
         Multiplier for the base frequency of the output oscillation, default is 1.
@@ -89,8 +91,8 @@ def ppc_constant_phase_shift(
         Multiplier for the base frequency of the input oscillation, default is 1.
 
     random_state : None, optional
-        Random state can be fixed to provide reproducible results. Otherwise, the
-        results may differ between function calls.
+        Random state can be fixed to provide reproducible results if the envelope
+        is generated randomly. Otherwise, the results may differ between function calls.
 
     Returns
     -------
@@ -99,9 +101,9 @@ def ppc_constant_phase_shift(
     """
     if not np.iscomplexobj(waveform):
         waveform = hilbert(waveform)
-    waveform_angle = np.angle(waveform)
 
     waveform_amp = _get_envelope(waveform, envelope, sfreq, fmin, fmax, random_state)
+    waveform_angle = np.angle(waveform)
     waveform_coupled = np.real(
         waveform_amp * np.exp(1j * m / n * waveform_angle + 1j * phase_lag)
     )
@@ -109,8 +111,8 @@ def ppc_constant_phase_shift(
         return normalize_variance(waveform_coupled)
 
     # NOTE: if the envelope was modified, we filter the result again in the target
-    # frequency range to suppress distortions due to merging amplitude envelope and
-    # phase from different time series
+    # frequency range to suppress possible distortions due to merging amplitude
+    # envelope and phase from different time series
     b, a = butter(
         N=2, Wn=np.array([m / n * fmin, m / n * fmax]) / sfreq * 2, btype="bandpass"
     )
@@ -170,10 +172,10 @@ def ppc_von_mises(
         Controls the amplitude envelope of the coupled waveform to be either randomly
         generated or to be the same as the envelope of the input waveform.
 
-    m : int, optional
+    m : float, optional
         Multiplier for the base frequency of the output oscillation, default is 1.
 
-    n : int, optional
+    n : float, optional
         Multiplier for the base frequency of the input oscillation, default is 1.
 
     random_state : None (default) or int
@@ -199,6 +201,10 @@ def ppc_von_mises(
     waveform_coupled = np.real(
         waveform_amp * np.exp(1j * m / n * waveform_angle + 1j * ph_distr)
     )
+
+    # NOTE: we filter the result again in the target frequency range to suppress
+    # possible distortions due to separate adjustment of the phase and amplitude
+    # of the coupled time series
     b, a = butter(
         N=2, Wn=np.array([m / n * fmin, m / n * fmax]) / sfreq * 2, btype="bandpass"
     )

From e3665852097402e7b7dcaf53dac8d01e955e58cc Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 15:15:33 +0200
Subject: [PATCH 09/11] Docstring

---
 src/meegsim/coupling.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/meegsim/coupling.py b/src/meegsim/coupling.py
index aaab7f0..7c7e768 100644
--- a/src/meegsim/coupling.py
+++ b/src/meegsim/coupling.py
@@ -74,7 +74,7 @@ def ppc_constant_phase_shift(
 
     envelope : str, {"same", "random"}
         Controls the amplitude envelope of the coupled waveform to be either randomly
-        generated or to be the same as the envelope of the input waveform.
+        generated (default) or to be the same as the envelope of the input waveform.
 
     fmin : float, optional
         Lower cutoff frequency for the oscillation that gives rise to the random
@@ -170,7 +170,7 @@ def ppc_von_mises(
 
     envelope : str, {"same", "random"}
         Controls the amplitude envelope of the coupled waveform to be either randomly
-        generated or to be the same as the envelope of the input waveform.
+        generated (default) or to be the same as the envelope of the input waveform.
 
     m : float, optional
         Multiplier for the base frequency of the output oscillation, default is 1.

From 426a979981a83a0bd41bdb29a2fca37e11eab611 Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Fri, 27 Jun 2025 15:18:11 +0200
Subject: [PATCH 10/11] MAINT: changelog entry

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4b2f61a..cfe85b0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,6 +19,7 @@ noise ([#58](https://github.com/ctrltz/meegsim/pull/58))
 - A method for setting phase-phase coupling by adding noise to the shifted copy of input waveform ([#71](https://github.com/ctrltz/meegsim/pull/71))
 - Function to convert the sources to mne.Label ([#73](https://github.com/ctrltz/meegsim/pull/73))
 - Quick list-like access to the simulated sources ([#82](https://github.com/ctrltz/meegsim/pull/82))
+- Control over the amplitude envelope of the coupled waveform ([#87](https://github.com/ctrltz/meegsim/pull/87))
 
 ### Changed
 

From a272079b9571a7df226c9ab5fdd686f9b37c903e Mon Sep 17 00:00:00 2001
From: Nikolai Kapralov <4dvlup@gmail.com>
Date: Tue, 22 Jul 2025 16:11:52 +0200
Subject: [PATCH 11/11] Minor edit

---
 src/meegsim/coupling.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/meegsim/coupling.py b/src/meegsim/coupling.py
index 80a35c7..003d1c9 100644
--- a/src/meegsim/coupling.py
+++ b/src/meegsim/coupling.py
@@ -92,7 +92,7 @@ def ppc_constant_phase_shift(
 
     random_state : None, optional
         Random state can be fixed to provide reproducible results if the envelope
-        is generated randomly. Otherwise, the results may differ between function calls.
+        is generated randomly. If not set, the results may differ between function calls.
 
     Returns
     -------