返回一个实值的相位加扰时间序列

Posted 2023-03-12

【中文标题】返回一个实值的相位加扰时间序列

【英文标题】：Returning a real-valued, phase-scrambled timeseries

【发布时间】：2017-01-25 08:41:58

【问题描述】：

我正在尝试使用 Scipy/Numpy 在 Python 中实现时间序列的“相位加扰”。简而言之，我想：

进行时间序列。 使用 FFT 在频域中测量功率和相位。 扰乱现有相位，随机将相位重新分配给频率。 使用 IFFT 返回具有加扰相位的实值（即非复数）时间序列，以使时间序列的功率谱保持不变，但时间序列的点与原始点不同。

我有一个表面上似乎可以工作的脚本（见情节），但我怀疑我遗漏了一些重要的东西。特别是，我返回的相位加扰时间序列具有复值条目而不是实值条目，我不知道该怎么做。如果任何信号处理人员可以参与并教育我，我将不胜感激。

这是一个适用于 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()

编辑：根据以下解决方案之一，我将偶数情况概括为奇数情况，如下所示。我相信检测不可忽略的虚部的条件现在是不必要的，但我会把它留给后代。

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

【问题讨论】：

傅里叶变换的一个性质：实值时域信号具有共轭对称频域信号。我的建议：打乱频域的一半，让另一半与之共轭。

【参考方案1】：

Scipy 有 rfft 和 irfft 例程来对真实信号进行傅里叶变换（r 是真实的）。下面的函数返回一个正确加扰的信号。我应该注意，编写的例程仅适用于偶数长度的信号，但修复奇数长度的信号应该不会太难，只需查看rfft 文档

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

【讨论】：

这是一个很好的工作答案，我很感激。我选择另一个作为“正确”，因为我发现它更易于理解并且能够更容易地将其扩展到奇数长度的 TS 案例，但如果我更了解 Numpy，我怀疑你的解决方案会更受欢迎。谢谢！

【参考方案2】：

傅里叶变换的性质：实值时域信号具有共轭对称频域信号。见110 of Functional relationships of Fourier transform。

解决方案（可能不是最佳的）

在phaseScrambleTs()中编辑下面的代码：

np.random.shuffle(phase_fsr)

到：

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

频域的左半部分取phase_fsr[1:len(phase_fsr)/2]。请注意，起始索引是1，而不是0。并洗牌。频域的右半部分被确定为phase_fsr_lh 的符号反转（根据共轭的定义）。并追加所有元素，第一个元素、左半部分、居中元素和右半部分。

【讨论】：

这行得通，它让我走上了解决奇数长度 TS 案例的正确轨道。谢谢！