Surrogate Data Technique

The surrogate data technique has been developed to investigate the phase coherency of waves/structures observed.

reference:
https://royalsocietypublishing.org/rsta/article-abstract/366/1864/447/50838/Experimental-evidence-of-phase-coherence-of?redirectedFrom=fulltext
https://journals.aps.org/pre/abstract/10.1103/PhysRevE.75.046401
https://link.springer.com/article/10.1023/A:1025506124402

%config InlineBackend.figure_format='retina'
import math as mt
import matplotlib.pyplot as plt
import numpy as np
import scipy.fft as sf
import datetime
import bisect
import matplotlib.patches as patches
import matplotlib.dates as mdates
from IPython.display import display,Math
import matplotlib.animation as animation

plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['animation.html'] = 'jshtml'
plt.rcParams['animation.embed_limit'] = 256  # Set limit to 256MB

f1=open('19951008_fx_hi_18595.txt','r')
#f1=open('19951008_fx_hi_18099.txt','r')

geb =np.fromfile(f1, dtype='float64',sep=' ').reshape(-1,11)

strt =0
width=16384
end  =strt+width

bx=geb[:,5]
by=geb[:,6]
bz=geb[:,7]
bb=np.sqrt(bx*bx+by*by+bz*bz)

#start = datetime.datetime(1995, 10, 8,  22, 30, 0)
#stop  = datetime.datetime(1995, 10, 8,  23, 30, 0)
#start_ind= bisect.bisect_left(gemag['Epoch'], start)
#stop_ind = bisect.bisect_left(gemag['Epoch'], stop)
#print(start_ind,stop_ind)
plt.figure(figsize=(80,4))
plt.plot( bx[strt:end],'r-')
plt.plot( by[strt:end],'g-')
plt.plot( bz[strt:end],'b-')
plt.plot( bb[strt:end],'k-')
plt.plot(-bb[strt:end],'k-')
plt.fill_between(np.arange(width),-bb[strt:end],bb[strt:end],alpha=0.2,color='grey')
plt.ylim(-15,15)
plt.show()

The key of this technique is to prepare two extreme datasets with fully-randomised and fully-correlated from the original dataset, and we compare the original data with the others and calculate the coherency of the original data.

We first decompose data into the frequency and phase ranges using Fourier transformations, and then only change the phase data to a randomized set and correlated set ($\phi=0$).
For each dataset, we now calculate the path length defined as, $$L(\tau)=\sum_t|x(t+\tau)-x(t)|,$$ where $\tau$ is the time lag. Finally, we calculate the following equation to obtain the phase coherency index $C_\phi$, $$C_\phi(\tau)=\frac{L_{PRS}(\tau)-L_{ORG}(\tau)}{L_{PRS}(\tau)-L_{PCS}(\tau)},$$ where ORG, PRS, and PCS are the original, randomised, and correlated data set, respectively.

nt=np.size(geb[strt:end,0])#np.size(geb['Epoch'][:])
dt=1/16    #[s]
t=dt*np.arange(nt)
lt=dt*nt
ww=sf.fftshift(sf.fftfreq(nt,dt))

### tapered function
wdw=0.5*(np.tanh((t-lt/20)/20)-np.tanh((t-19*lt/20)/20))#np.r_[np.zeros(nt//8),np.ones(6*nt//8),np.zeros(nt//8)]

plt.plot(wdw)
plt.title('tapered function')
plt.show()

dat=bz[strt:end]-bz[strt:end].mean()
dat=dat*wdw

### original data ###
forg=sf.fft(dat)
amp =np.abs(forg)
phi =np.angle(forg)
#####################

### phase randomized surrogate ###
#phir=phi+np.pi#np.random.permutation(phi)
#for it in range(100): 
#    phirr=phir+np.random.permutation(phir)
#phir=phir/100-np.pi
phir=np.random.permutation(phi)#np.pi*(1-2*np.random.rand(nt))
#phir=np.pi*(1-2*np.random.rand(nt))
phir[0]=0; phir[nt:nt//2:-1]=-phir[1:nt//2] #symmetric
fprs=amp*(np.cos(phir)+1j*np.sin(phir))
datr=np.real(sf.ifft(fprs))
##################################

### phase correlated surrogate ###
phic=np.zeros(nt)
fpcs=amp*(np.cos(phic)+1j*np.sin(phic))
datc=np.real(sf.ifft(fpcs))
##################################

### plot figure ###
fig, ax = plt.subplots(nrows=3,ncols=3, sharex='col', figsize=(18, 8))
### original data ###
ax[0,0].plot(t,dat)
ax[0,1].plot(ww,sf.fftshift(abs(forg)))
ax[0,2].plot(ww,sf.fftshift(phi),'.')
ax[0,0].set_ylabel(r'$B_{ORG}$ [nT]')
ax[0,1].set_ylabel(r'power')
ax[0,2].set_ylabel(r'phase [rad]')
ax[0,1].set_xscale('log')
ax[0,1].set_yscale('log')
ax[0,2].set_xscale('log')
ax[0,2].set_ylim(-np.pi,np.pi)
#####################

### phase-randomized surrogate ###
ax[1,0].plot(t,datr)
ax[1,1].plot(ww,sf.fftshift(abs(fprs)))
ax[1,2].plot(ww,sf.fftshift(phir),'.')
ax[1,0].set_ylabel(r'$B_{PRS}$ [nT]')
ax[1,1].set_ylabel(r'power')
ax[1,1].set_yscale('log')
ax[1,2].set_ylabel(r'phase [rad]')
ax[1,2].set_ylim(-np.pi,np.pi)
##################################

### phase-correlated surrogate ###
ax[2,0].plot(t,datc)
ax[2,1].plot(ww,sf.fftshift(abs(fpcs)))
ax[2,2].plot(ww,sf.fftshift(phic),'.')
ax[2,0].set_ylabel(r'$B_{PCS}$ [nT]')
ax[2,1].set_ylabel(r'power')
ax[2,1].set_yscale('log')
ax[2,2].set_ylabel(r'phase [rad]')
ax[2,2].set_ylim(-np.pi,np.pi)
##################################
plt.show()

def path_length(nt,dat):
    x=np.arange(nt)
    L=np.zeros(nt)
    for tau in range(nt):
        #jp=np.r_[np.arange(nt-tau)+tau, np.arange(tau)]
        #L[tau]=np.sum(np.abs(dat[jp]-dat[x]))
        L[tau]=np.sum(np.abs(np.roll(dat,tau)-dat[x]))
    return L

### calculate path length for each datasets
Lorg=path_length(nt,dat )
Lprs=path_length(nt,datr)
Lpcs=path_length(nt,datc)

### plot the path lengths
plt.figure(figsize=(8,4))
tt=dt*np.arange(nt)
plt.plot(tt,Lorg,label='original')
plt.plot(tt,Lprs,label='random')
plt.plot(tt,Lpcs,label='correlated')
plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.ylabel('L')
plt.xlabel(r'$\tau \text{ [s]}$')
plt.show()

### caluclate and plot the phase coherence index 
Cphi=(Lprs-Lorg)/(Lprs-Lpcs)
plt.figure(figsize=(8,4))
plt.plot(tt,Cphi)
plt.xlabel(r'$\tau \text{ [s]}$')
plt.ylabel('$C_\phi$')
plt.xscale('log')
plt.ylim(0,1)
plt.show()

/var/folders/hl/2w9c9jdj5jqgbpmm44p6whpr0000gn/T/ipykernel_66958/3510904837.py:29: RuntimeWarning: invalid value encountered in divide
  Cphi=(Lprs-Lorg)/(Lprs-Lpcs)