dewarrn1 - 1 year ago 67
Python Question

# Returning a real-valued, phase-scrambled timeseries using Python

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
``````

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.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download