Phase randomization.

    This function randomizes the spectral phase of the input data along the
    last dimension.

    Parameters
    ----------
    data : array
        Input array.

    Returns
    -------
    out : array
        Array of same shape as data.

    Notes
    -----
    The algorithm randomizes the phase component of the input's complex Fourier
    transform.

    Examples
    --------
    .. plot::
        :include-source:

        from pylab import *
        from scot.datatools import randomize_phase
        np.random.seed(1234)
        s = np.sin(np.linspace(0,10*np.pi,1000))
        x = np.vstack([s, np.sign(s)])
        y = randomize_phase(x)
        subplot(2,1,1)
        title('Phase randomization of sine wave and rectangular function')
        plot(x.T + [1.5, -1.5]), axis([0,1000,-3,3])
        subplot(2,1,2)
        plot(y.T + [1.5, -1.5]), axis([0,1000,-3,3])
        plt.show()

def randomize_phase(data, random_state=None):
    """Phase randomization.

    This function randomizes the spectral phase of the input data along the
    last dimension.

    Parameters
    ----------
    data : array
        Input array.

    Returns
    -------
    out : array
        Array of same shape as data.

    Notes
    -----
    The algorithm randomizes the phase component of the input's complex Fourier
    transform.

    Examples
    --------
    .. plot::
        :include-source:

        from pylab import *
        from scot.datatools import randomize_phase
        np.random.seed(1234)
        s = np.sin(np.linspace(0,10*np.pi,1000))
        x = np.vstack([s, np.sign(s)])
        y = randomize_phase(x)
        subplot(2,1,1)
        title('Phase randomization of sine wave and rectangular function')
        plot(x.T + [1.5, -1.5]), axis([0,1000,-3,3])
        subplot(2,1,2)
        plot(y.T + [1.5, -1.5]), axis([0,1000,-3,3])
        plt.show()
    """
    rng = check_random_state(random_state)
    data = np.asarray(data)
    data_freq = np.fft.rfft(data)
    data_freq = np.abs(data_freq) * np.exp(1j*rng.random_sample(data_freq.shape)*2*np.pi)
    return np.fft.irfft(data_freq, data.shape[-1])