causality_tools.py

import numpy
import scipy.signal # as signal # for periodogram
import matplotlib
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft, fftfreq, fftshift # faster than numpy fft

# basic little tools 
# (pigs data has low resolution, and many identical points...)
def count_non_duplicates(x, verbosity=0):
    d1=numpy.count_nonzero(x)
    d2=numpy.count_nonzero(numpy.diff(numpy.sort(x)))
    if verbosity>0: print("%d -> %d" %(d1,d2+1))
    return d2+1

def get_non_duplicates(x):
    xx=numpy.unique(x)
    print(xx.shape, "non duplicates, vs", count_non_duplicates(x))
    return xx

def get_resolution(x, verbosity=0):
    tmp=numpy.diff(numpy.sort(numpy.unique(x)))
    print("min %f, max %f" %(numpy.min(tmp),numpy.max(tmp)))
    tmp=tmp[numpy.where(tmp>(numpy.max(tmp)/tmp.size))]
    dx    = numpy.min(tmp)
    return dx
    
def get_range(x, verbosity=0):
    dx    = get_resolution(x)
    range = (numpy.max(x)-numpy.min(x))/dx
    if verbosity>0: print("range:", range)
    return range


# PSD and searching for time-scales
def plot_PSD(signals, fs, channel_set, channel_names, N_fft=4096, subplot_to_use=-1, extra_label=""):
#channel_set=range(16)
#channel_set = [channel_ECG-1, channel_ECoG_central_left-1, channel_ECoG_central_right-1, 11, 12]
#channel_set = [0,1,2,3]
    N_pts = signals.shape[1]
    N_mag = 16    # "magnifyig factor" to access lower frequencies
    while ((N_fft*N_mag)>N_pts) and (N_mag>1) :
        N_mag=N_mag//2

    if isinstance(subplot_to_use, matplotlib.axes.Axes):
        P1=subplot_to_use
#        print(P1.lines[-1].get_color())
        bool_old = True
    else:
        Fig=plt.figure(figsize=(8,5))
        P1 =Fig.add_subplot(1,1,1)
        bool_old = False

    j=0
    for i in channel_set:
        sig = signals[i,:]
        my_text = channel_names[i]+extra_label
        extra_params = {}
        if bool_old:
            print("i=%d / %d" %(j, len(P1.lines)))
            my_color=P1.lines[j].get_color()
            extra_params = {'color': my_color, 'linestyle': 'dotted'}
        freq, psd = scipy.signal.welch(sig,fs=fs, window='hanning', nperseg=N_fft)
        h=P1.loglog(freq[N_mag:], psd[N_mag:], label=my_text, **extra_params)
        my_color=h[-1].get_color()

        # low frequencies:
        freq, psd = scipy.signal.welch(sig,fs=fs, window='hanning', nperseg=N_fft*N_mag)
        P1.loglog(freq[:N_mag*N_mag], psd[:N_mag*N_mag], color=my_color, label=None)

        j += 2 # we have added 2 plots per dataset !
#        # low-pass filtered signals:
#        sig = signals_LP[i,:]
#        freq, psd = scipy.signal.welch(sig,fs=fs, window='hanning', nperseg=N_fft)
#        P1.loglog(freq[N_mag:]/tau_LP, psd[N_mag:]*tau_LP, '--', color=my_color, label=None)


    P1.set_xlabel(r"$f$ (Hz)", fontsize=18)
    P1.set_ylabel(r"($V/\sqrt{Hz}$)", fontsize=18)
    P1.legend(fontsize=18,  bbox_to_anchor=(1.05, 1))
    #P1.set_xlim(10, 30)
    print("effective frequency : %.0f Hz" %fs)
    return P1


# low pass filter signals
#
# this improves the dynamics (in terms of nb of non-redondant points)
# and this reduces the signal size (so better for ANN algorithms)
#
# signals      : 2-d signal
# tau_LP       : nb of pts to average
# f_resampling : how many point per set of tau_LP to keep (oversampling)
#
# 2022-10-24 : added parameter f_resampling (old default was indeed 1)
#            : tested OK (see notebook "causality_couples_2022-10-21_tests_decel")
def filter_FIR(signals, tau_LP, f_resampling=1):
    Npts=signals.shape[1] # along time
#    print(Npts, "->", (Npts*f_resampling)//tau_LP, "pts in time")
    signals_LP=numpy.zeros((signals.shape[0], (Npts*f_resampling)//tau_LP), dtype="float")
    shift = tau_LP//f_resampling
    for i in range(signals.shape[0]):
        for j in numpy.arange(f_resampling):
            tmp = signals[i,j*shift:]
            Npts_decimated = tmp.size//tau_LP
#            print(Npts_decimated, "(new Npts)")
            tmp = numpy.reshape(tmp[:(Npts_decimated*tau_LP)], (Npts_decimated,tau_LP))
#            print(tmp.shape, "->", numpy.mean(tmp, axis=1).shape, "vs", signals_LP[i,j::f_resampling].shape)
            signals_LP[i,j:Npts_decimated*f_resampling:f_resampling] = numpy.mean(tmp, axis=1)
#        print("%2d : %3d -> %5d non-duplicates" %(i, count_non_duplicates(signals[i,:]), count_non_duplicates(signals_LP[i,:])))
    
    return signals_LP

# Low-Pass filter the signal(s) "signals"
# with a linear causal filter of order 1, with cut-off frequancy fc
#
# signals may be a tensor. In that case, time is assumed to be the last dimension/axis
# fc is the cut-off frequency (expressed in Hertz)
# fs is the sampling frequency (in Hz)
def filter_FFT_LP(signals, fc, fs=1, type="causal"):
    N=signals.shape[-1]   # nb of pts in time
    print(N)
    
    mode   = numpy.arange(N)    # Fourier mode
    f      = (mode-N//2)/N*fs   # true frequency (Hz)
#    f      =fftfreq(N, 1/fs)    # same as the 2 lines above

    # filter construction:
    if type=="causal":
        filter = 1/(1+1j*f/fc)   # causal LP filter of order 1
    else:
        filter = 1/numpy.sqrt(1+(f/fc)**2)   # non-causal LP filter of order 1

    # shift du filtre:
#    filter_s = numpy.zeros(N, dtype='float')
#    filter_s[N//2+1:] = filter[0:N//2]
#    filter_s[0:N//2]  = filter[N//2+1:]
    filter_s = numpy.roll(filter, (N+1)//2)    # same as the 3 lines above
    if (N%2==0): filter_s[N//2] = 0  # annulation de la frequence N/2 si N pair

    # filtrage
    xf = fft(signals) # by default, axis=-1, so along the last axis, ie, along time
    yf = xf*filter_s
    y  = ifft(yf)
    
    return y.real
 
# find appropriate masks for accelerations and decelerations   
#
# data          is the (foetus,mother) data
# inc_threshold is the threshold for definition of acceleration or deceleration
# use_mother    indicates to use mother HR for definition of accel/decel. (default: True)
#
# returned vector is composed of a mask for acceleration, and a mask for decelerations
def get_accel_decel_masks(data, inc_threshold=0., use_mother=True):
    if (use_mother==True):  
        dHR = numpy.diff(data[1,:])
#        print("mother")
    else:                   
        dHR = numpy.diff(data[0,:])
#        print("foetus")
        
    ind_accel = numpy.where(dHR>inc_threshold)
    ind_decel = numpy.where(dHR<-inc_threshold)
    
    ind_accel = numpy.array(ind_accel)+1    # +1 to have correct (causal) date at which the increment was commputed
    ind_decel = numpy.array(ind_decel)+1

    mask_accel  = numpy.zeros(data.shape[1], dtype='int8')
    mask_decel  = numpy.zeros(data.shape[1], dtype='int8')
    mask_accel[ind_accel] = 1
    mask_decel[ind_decel] = 1

    return numpy.array([mask_accel, mask_decel])