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:

Take a timeseries.
Measure the power and phase in the frequency domain using FFT.
Scramble the existing phase, randomly reassigning phase to frequency.
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

A property of Fourier transform: a real valued time domain signal has a conjugate symmetric frequency domain signal. My suggestion: scramble a half of a frequency domain, make another half be conjugate of it. — Jeon

Elliot Elliot · Accepted Answer · 2016-09-17T08:19:59

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

Returning a real-valued, phase-scrambled timeseries

2 Answers

Solution (may not be optimal)