def METHOD_NAME(signal, method="IAAFT", random_state=None, **kwargs):
    """**Create Signal Surrogates**
    Generate a surrogate version of a signal. Different methods are available, such as:
    * **random**: Performs a random permutation of the signal value. This way, the signal
      distribution is unaffected and the serial correlations are cancelled, yielding a whitened
      signal with an distribution identical to that of the original.
    * **IAAFT**: Returns an Iterative Amplitude Adjusted Fourier Transform (IAAFT) surrogate.
      It is a phase randomized, amplitude adjusted surrogates that have the same power spectrum
      (to a very high accuracy) and distribution as the original data, using an iterative scheme.
    Parameters
    ----------
    signal : Union[list, np.array, pd.Series]
        The signal (i.e., a time series) in the form of a vector of values.
    method : str
        Can be ``"random"`` or ``"IAAFT"``.
    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.
    **kwargs
        Other keywords arguments, such as ``max_iter`` (by default 1000).
    Returns
    -------
    surrogate : array
        Surrogate signal.
    Examples
    --------
    Create surrogates using different methods.
    .. ipython:: python
      import neurokit2 as nk
      import matplotlib.pyplot as plt
      signal = nk.signal_simulate(duration = 1, frequency = [3, 5], noise = 0.1)
      surrogate_iaaft = nk.signal_surrogate(signal, method = "IAAFT")
      surrogate_random = nk.signal_surrogate(signal, method = "random")
      @savefig p_signal_surrogate1.png scale=100%
      plt.plot(surrogate_random, label = "Random Surrogate")
      plt.plot(surrogate_iaaft, label = "IAAFT Surrogate")
      plt.plot(signal, label = "Original")
      plt.legend()
      @suppress
      plt.close()
    As we can see, the signal pattern is destroyed by random surrogates, but not in the IAAFT one.
    And their distributions are identical:
    .. ipython:: python
      @savefig p_signal_surrogate2.png scale=100%
      plt.plot(*nk.density(signal), label = "Original")
      plt.plot(*nk.density(surrogate_iaaft), label = "IAAFT Surrogate")
      plt.plot(*nk.density(surrogate_random), label = "Random Surrogate")
      plt.legend()
      @suppress
      plt.close()
    However, the power spectrum of the IAAFT surrogate is preserved.
    .. ipython:: python
      f = nk.signal_psd(signal, max_frequency=20)
      f["IAAFT"] = nk.signal_psd(surrogate_iaaft, max_frequency=20)["Power"]
      f["Random"] = nk.signal_psd(surrogate_random, max_frequency=20)["Power"]
      @savefig p_signal_surrogate3.png scale=100%
      f.plot("Frequency", ["Power", "IAAFT", "Random"])
      @suppress
      plt.close()
    References
    ----------
    * Schreiber, T., & Schmitz, A. (1996). Improved surrogate data for nonlinearity tests. Physical
      review letters, 77(4), 635.
    """
    # TODO: when discrete signal is detected, run surrogate of markov chains
    # https://github.com/Frederic-vW/eeg_microstates/blob/eeg_microstates3.py#L861
    # Or markov_simulate()
    # Seed the random generator for reproducible results
    rng = check_random_state(random_state)
    method = method.lower()
    if method == "random":
        surrogate = rng.permutation(signal)
    elif method == "iaaft":
        surrogate, _, _ = _signal_surrogate_iaaft(signal, rng=rng, **kwargs)
    return surrogate