Returning a real-valued, phase-scrambled timeseries

Question

I'm trying to implement "phase scrambling" of a timeseries in Python using Scipy/Numpy. Briefly, I'd like to:

1. Take a timeseries.
2. Measure the power and phase in the frequency domain using FFT.
3. Scramble the existing phase, randomly reassigning phase to frequency.
4. Return a real-valued (i.e., non-complex) timeseries with scrambled phase using IFFT such that the power spectrum of the timeseries remains constant but the points of the timeseries differ from the original.

I have a script that superficially seems to work (see plots), but I suspect that I'm missing something important. In particular, my returned phase-scrambled timeseries has complex-valued entries instead of real-valued entries, and I'm not sure what to do with that. If any signal-processing folks could weigh in and educate me, I'd greatly appreciate it.

Here's a sample script suitable for Jupyter Notebook:

%matplotlib inline
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft

def phaseScrambleTS(ts):
    """Returns a TS: original TS power is preserved; TS phase is shuffled."""
    fs = fft(ts)
    pow_fs = np.abs(fs) ** 2.
    phase_fs = np.angle(fs)
    phase_fsr = phase_fs.copy()
    np.random.shuffle(phase_fsr)
    fsrp = np.sqrt(pow_fs) * (np.cos(phase_fsr) + 1j * np.sin(phase_fsr))
    tsr = ifft(fsrp)
    return tsr

ts = np.array([0.02, -1.04, 2.50, 2.21, 1.37, -0.05, 0.06, -0.22, -0.48, -0.31, 0.15, 0.99, 0.39, 0.65, 1.13, 0.77, 1.16, 1.35, 0.92, 1.42, 1.58, 1.33, 0.73, 0.98, 0.66, 0.13, -0.19, 2.05, 1.95, 1.25, 1.37, 0.85, 0.00, 1.37, 2.17, 0.69, 1.38, 0.49, 0.52, 0.62, 1.74, 0.67, 0.61, 1.03, 0.38, 0.64, 0.83, 1.16, 1.10, 1.30, 1.98, 0.92, 1.36, -1.49, -0.80, -0.08, 0.01, -0.04, -0.07, -0.20, 0.82, -0.26, 0.83, 0.09, -0.54, -0.45, 0.82, -0.53, -0.88, -0.54, -0.30, 0.52, 0.54, -0.57, 0.73, -0.04, 0.34, 0.59, -0.67, -0.25, -0.44, 0.07, -1.00, -1.88, -2.55, -0.08, -1.13, -0.94, -0.09, -2.08, -1.56, 0.25, -1.87, 0.52, -0.51, -1.42, -0.80, -1.96, -1.42, -1.27, -1.08, -1.79, -0.73, -2.70, -1.14, -1.71, -0.75, -0.78, -1.87, -0.88, -2.15, -1.92, -2.17, -0.98, -1.52, -1.92], dtype=np.float)

N = ts.shape[0]
TR = 2.
x = np.linspace(0.0, N*TR, N)
plt.plot(x, ts)
plt.ylabel('% Sig. Change')
plt.xlabel('Time')
plt.title('RSFC: Time domain')
plt.show()

ts_ps = phaseScrambleTS(ts)
plt.plot(x, ts, x, ts_ps)
plt.ylabel('% Sig. Change')
plt.xlabel('Time')
plt.title('RSFC, Orig. vs. Phase-Scrambled: Time domain')
plt.show()

fs = fft(ts)
fs_ps = fft(ts_ps)
xf = np.linspace(0.0, 1.0/(2.0*TR), N/2)
plt.plot(xf, 2./N * np.abs(fs[0:N/2]), 'b--', xf, 2./N * np.abs(fs_ps[0:N/2]), 'g:')
plt.grid()
plt.ylabel('Amplitude')
plt.xlabel('Freq.')
plt.title('RSFC, Orig. vs. Phase-Scrambled: Freq. domain, Amp.')
plt.show()

Edit: Following up on one of the solutions below, I generalized from the even case to the odd case as follows. I believe that the conditional for detecting non-negligible imaginary components is now unnecessary, but I'll leave it in for posterity.

def phaseScrambleTS(ts):
    """Returns a TS: original TS power is preserved; TS phase is shuffled."""
    fs = fft(ts)
    pow_fs = np.abs(fs) ** 2.
    phase_fs = np.angle(fs)
    phase_fsr = phase_fs.copy()
    if len(ts) % 2 == 0:
        phase_fsr_lh = phase_fsr[1:len(phase_fsr)/2]
    else:
        phase_fsr_lh = phase_fsr[1:len(phase_fsr)/2 + 1]
    np.random.shuffle(phase_fsr_lh)
    if len(ts) % 2 == 0:
        phase_fsr_rh = -phase_fsr_lh[::-1]
        phase_fsr = np.concatenate((np.array((phase_fsr[0],)), phase_fsr_lh,
                                    np.array((phase_fsr[len(phase_fsr)/2],)),
                                    phase_fsr_rh))
    else:
        phase_fsr_rh = -phase_fsr_lh[::-1]
        phase_fsr = np.concatenate((np.array((phase_fsr[0],)), phase_fsr_lh, phase_fsr_rh))
    fsrp = np.sqrt(pow_fs) * (np.cos(phase_fsr) + 1j * np.sin(phase_fsr))
    tsrp = ifft(fsrp)
    if not np.allclose(tsrp.imag, np.zeros(tsrp.shape)):
        max_imag = (np.abs(tsrp.imag)).max()
        imag_str = '\nNOTE: a non-negligible imaginary component was discarded.\n\tMax: {}'
        print imag_str.format(max_imag)
    return tsrp.real

Answer 1

Scipy has rfft and irfft routines to do Fourier transforms of real signals (the r is for real). The function below returns a properly scrambled signal. I should note that the routine as written will only work for even length signals, but fixing it for odd length signals shouldn't be too hard, just check the rfft documentation

from scipy.fftpack import rfft, irfft

def phaseScrambleTS(ts):
    """Returns a TS: original TS power is preserved; TS phase is shuffled."""
    fs = rfft(ts)
    # rfft returns real and imaginary components in adjacent elements of a real array
    pow_fs = fs[1:-1:2]**2 + fs[2::2]**2
    phase_fs = np.arctan2(fs[2::2], fs[1:-1:2])
    phase_fsr = phase_fs.copy()
    np.random.shuffle(phase_fsr)
    # use broadcasting and ravel to interleave the real and imaginary components. 
    # The first and last elements in the fourier array don't have any phase information, and thus don't change
    fsrp = np.sqrt(pow_fs[:, np.newaxis]) * np.c_[np.cos(phase_fsr), np.sin(phase_fsr)]
    fsrp = np.r_[fs[0], fsrp.ravel(), fs[-1]]
    tsr = irfft(fsrp)
    return tsr

Answer 2

A property of Fourier transform: a real valued time domain signal has a conjugate symmetric frequency domain signal. See 110 of Functional relationships of Fourier transform.

Solution (may not be optimal)

Edit a code below in phaseScrambleTs():

np.random.shuffle(phase_fsr)

to:

phase_fsr_lh = phase_fsr[1:len(phase_fsr)/2]
np.random.shuffle(phase_fsr_lh)
phase_fsr_rh = -phase_fsr_lh[::-1]
phase_fsr = np.append(phase_fsr[0],
                      np.append(phase_fsr_lh,
                                np.append(phase_fsr[len(phase_fsr)/2],
                                          phase_fsr_rh)))

The left half of frequency domain takes phase_fsr[1:len(phase_fsr)/2]. Note that the starting index is 1, not 0. And shuffle it. The right half of frequency domain is determined as a sign-inverted of phase_fsr_lh in reverse order (by definition of conjugate). And append all of them, the first element, the left half, the centered element, and the right half.