Pythonを使って実数値のフェーズスクランブルtimeseriesを返す

Question

私はScipy/Numpyを使ってPythonでtimeseriesの "phase scrambling"を実装しようとしています。簡単に、私はしたい：

1. 私はしたい：timeseries。
2. FFTを使用して周波数領域のパワーと位相を測定します。
3. 既存のフェーズをスクランブルし、位相を周波数にランダムに再割り当てします。
4. IFFTを使用して、時系列のパワースペクトルは一定のままであるが、時系列の点は元のものとは異なるように、スクランブルされた位相を有する実数値の（すなわち、複雑でない）

私は、表面的にはうまくいくようですが（プロット参照）、私は何か重要なものがないと思っています。特に返されたフェーズ・スクランブルされたtimeseriesは、実数値のエントリではなく複素数のエントリを持ちます。それでどうしたらいいのか分かりません。シグナル処理担当者が私の体重を増やして教育することができれば、私は大いに感謝しています。

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

編集：次のように下に、Iが奇数の場合にも、ケースから一般化ソリューションの1にフォローアップ

ここJupyterノートブックに適したサンプルスクリプトです。私は無視できない虚数成分を検出するための条件は不要になったと思いますが、後世のために残しておきます。

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

フーリエ変換の性質は変換：実数値時間領域信号は共役対称の周波数領域の信号を有しています。 110 of Functional relationships of Fourier transformを参照してください。

ソリューション（最適ではないかもしれない）

編集phaseScrambleTs()に以下のコード：周波数領域の左半分はphase_fsr[1:len(phase_fsr)/2]かかり

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

：

np.random.shuffle(phase_fsr) 

します。開始インデックスは1であり、0ではありません。シャッフルしてください。周波数領域の右半分は、逆順（コンジュゲートの定義による）でphase_fsr_lhの符号反転で決定されます。最初の要素、左半分、中央の要素、右半分のすべてを追加します。

Answer 2

Scipyは、実際の信号のフーリエ変換を行うルーチン（は実際のものです）を持つrfftとirfftルーティンを持っています。以下の関数は、適切にスクランブルされたシグナルを返します。実数値：私だけでも長さ信号のために働くだろう書かれていますが、奇数の長さの信号のためにそれを固定してルーチンがあまりにも難しいことではありませんことに注意しなければならない、ちょうどフーリエ変換の性質変換

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