Python signal.hilbert方法代碼示例

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def InstantaneousPhase(syn, obs, nt, dt, eps=0.05):
    """ Instantaneous phase (from Bozdag et al. 2011, eq 27)
    """
    r = _np.real(_analytic(syn))
    i = _np.imag(_analytic(syn))
    phi_syn = _np.arctan2(i, r)

    r = _np.real(_analytic(obs))
    i = _np.imag(_analytic(obs))
    phi_obs = _np.arctan2(i, r)

    phi_rsd = phi_syn - phi_obs
    esyn = abs(_analytic(syn))
    emax = max(esyn**2.)

    wadj = phi_rsd*_np.imag(_analytic(syn))/(esyn**2. + eps*emax) + \
        _np.imag(_analytic(phi_rsd * syn/(esyn**2. + eps*emax)))

    return wadj 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def InstantaneousPhase2(syn, obs, nt, dt, eps=0.):
    esyn = abs(_analytic(syn))
    eobs = abs(_analytic(obs))

    esyn1 = esyn + eps*max(esyn)
    eobs1 = eobs + eps*max(eobs)
    esyn3 = esyn**3 + eps*max(esyn**3)

    diff1 = syn/(esyn1) - obs/(eobs1)
    diff2 = _hilbert(syn)/esyn1 - _hilbert(obs)/eobs1

    part1 = diff1*_hilbert(syn)**2/esyn3 - diff2*syn*_hilbert(syn)/esyn3
    part2 = diff1*syn*_hilbert(syn)/esyn3 - diff2*syn**2/esyn3

    wadj = part1 + _hilbert(part2)
    return wadj


# Migration 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def SynchronisationIndex(MF1, MF2):
  """
  Computes the synchronisation index between two populations given firing data
  MF1 and MF2.
  """

  if len(MF1) != len(MF2):
    raise Exception("Both time series must have the same length")

  # Centre time series on zero
  MF1 = MF1 - np.mean(MF1)
  MF2 = MF2 - np.mean(MF2)

  # Calculate phase using Hilbert transform
  phase1 = np.angle(hilbert(MF1))
  phase2 = np.angle(hilbert(MF2))

  # Calculate synchronisation index
  phi = np.abs((np.exp(1j*phase1) + np.exp(1j*phase2)) / 2.0)

  return np.mean(phi) 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def snr(data,sampling_rate):
    """
    Signal to noise ratio of N cross-correlations.

    Follows method of Clarke et. al, 2011. Measures SNR at each point.

    """	
    data = np.array(data)
    N,t = data.shape
    data_mean = np.mean(data,axis=0)

    # calculate noise and envelope functions
    sigma = np.mean(data**2,axis=0) - (data_mean)**2
    sigma = np.sqrt(sigma/(N-1.))
    s = np.abs(data_mean + 1j*scipy.signal.hilbert(data_mean))

    # smooth with 10 second sliding cosine window 
    # half window length is 5s, i.e. 5 * sampling rate
    sigma = smooth(sigma,half_win=int(sampling_rate*5))
    s = smooth(s,half_win=int(sampling_rate*5))

    return np.real(s/sigma) 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def phaseanalysis(firstharmonic, displayplots=False):
    print('entering phaseanalysis')
    analytic_signal = hilbert(firstharmonic)
    amplitude_envelope = np.abs(analytic_signal)
    instantaneous_phase = np.angle(analytic_signal)
    if displayplots:
        print('making plots')
        fig = pl.figure()
        ax1 = fig.add_subplot(311)
        ax1.set_title('Analytic signal')
        X = np.linspace(0.0, 1.0, num=len(firstharmonic))
        pl.plot(X, analytic_signal.real, 'k', X, analytic_signal.imag, 'r')
        ax2 = fig.add_subplot(312)
        ax2.set_title('Phase')
        pl.plot(X, instantaneous_phase, 'g')
        ax3 = fig.add_subplot(313)
        ax3.set_title('Amplitude')
        pl.plot(X, amplitude_envelope, 'b')
        pl.show()
        pl.savefig('phaseanalysistest.jpg')
    instantaneous_phase = np.unwrap(instantaneous_phase)
    return instantaneous_phase, amplitude_envelope 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def pli(data, channels_num, window_length):
    try:
        from scipy.signal import hilbert
        if data.shape[0] != channels_num:
            data = data.T
        data_hil = hilbert(data)
        # if data_hil.shape != (channels_num, window_length):
        #     raise Exception('PLI: Wrong dimentions!')
        m = np.zeros((channels_num, channels_num))
        for i in range(channels_num):
            for j in range(channels_num):
                if i < j:
                    m[i, j] = abs(np.mean(np.sign(np.imag(data_hil[i] / data_hil[j]))))
                    # m[i, j] = abs(np.mean(np.sign(np.imag(data_hil[:, i] / data_hil[:, j]))))
        return m + m.T
    except:
        print(traceback.format_exc())
        return None 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def _detect_ready_tone(w, fs):
    # get envelope of DC free signal and envelope of BP signal around freq of interest
    h = np.abs(signal.hilbert(w - np.median(w)))
    fh = np.abs(signal.hilbert(dsp.bp(w, si=1 / fs, b=FTONE * np.array([0.9, 0.95, 1.15, 1.1]))))
    dtect = _running_mean(fh / (h + 1e-3), int(fs * 0.1)) > 0.8
    return np.where(np.diff(dtect.astype(int)) == 1)[0]
    # tone = np.sin(2 * np.pi * FTONE * np.arange(0, fs * 0.1) / fs)
    # tone = tone / np.sum(tone ** 2)
    # xc = np.abs(signal.hilbert(signal.correlate(w - np.mean(w), tone))) 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def hilbert(w):
    return np.imag(analytic(w)) 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def Envelope(syn, obs, nt, dt, eps=0.05):
    """ Envelope difference (Yuan et al 2015, eq 9)
    """
    esyn = abs(_analytic(syn))
    eobs = abs(_analytic(obs))
    ersd = esyn-eobs
    return np.sqrt(np.sum(ersd*ersd*dt)) 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def InstantaneousPhase(syn, obs, nt, dt, eps=0.05):
    """ Instantaneous phase from Bozdag et al. 2011
    """
    r = np.real(_analytic(syn))
    i = np.imag(_analytic(syn))
    phi_syn = np.arctan2(i, r)

    r = np.real(_analytic(obs))
    i = np.imag(_analytic(obs))
    phi_obs = np.arctan2(i, r)

    phi_rsd = phi_syn - phi_obs
    return np.sqrt(np.sum(phi_rsd*phi_rsd*dt)) 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def Envelope2(syn, obs, nt, dt, eps=0.):
    """ Envelope amplitude ratio (Yuan et al 2015, eq B-1)
    """
    esyn = abs(_analytic(syn))
    eobs = abs(_analytic(obs))
    raise NotImplementedError 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def Envelope3(syn, obs, nt, dt, eps=0.):
    """ Envelope cross-correlation lag (Yuan et al 2015, eqs B-4)
    """
    esyn = abs(_analytic(syn))
    eobs = abs(_analytic(obs))
    return Traveltime(esyn, eobs, nt, dt) 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def Envelope(syn, obs, nt, dt, eps=0.05):
    """ Envelope difference (Yuan et al 2015, eq 16)
    """
    esyn = abs(_analytic(syn))
    eobs = abs(_analytic(obs))
    etmp = (esyn - eobs)/(esyn + eps*esyn.max())
    wadj = etmp*syn - _np.imag(_analytic(etmp*_np.imag(_analytic(syn))))
    return wadj 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def Envelope3(syn, obs, nt, dt, eps=0.):
    """ Envelope cross-correlation lag (Yuan et al 2015, eqs B-2, B-5)
    """
    esyn = abs(_analytic(syn))
    eobs = abs(_analytic(obs))

    erat = _np.zeros(nt)
    erat[1:-1] = (esyn[2:] - esyn[0:-2])/(2.*dt)
    erat[1:-1] /= esyn[1:-1]
    erat *= misfit.Envelope3(syn, obs, nt, dt)

    wadj = -erat*syn + _hilbert(erat*_hilbert(esyn))
    return wadj 

# 需要導入模塊: from scipy import signal [as 別名]
# 或者: from scipy.signal import hilbert [as 別名]
def robust_hilbert(sig, increase_n=False):
    """Compute the Hilbert transform, ignoring any boundaries that are NaN.

    Parameters
    ----------
    sig : 1d array
        Time series.
    increase_n : bool, optional, default: False
        If True, zero pad the signal to length the next power of 2 for the Hilbert transform.
        This is because ``scipy.signal.hilbert`` can be very slow for some signal lengths.

    Returns
    -------
    sig_hilb : 1d array
        The Hilbert transform of the input signal.

    Examples
    --------
    Compute a Hilbert transform of a signal, using zero padding:

    >>> from neurodsp.sim import sim_combined
    >>> sig = sim_combined(n_seconds=10, fs=500,
    ...                    components={'sim_powerlaw': {}, 'sim_oscillation': {'freq': 10}})
    >>> sig_hilb = robust_hilbert(sig, increase_n=True)
    """

    # Extract the signal that is not nan
    sig_nonan, sig_nans = remove_nans(sig)

    # Compute Hilbert transform of signal without nans
    if increase_n:
        sig_len = len(sig_nonan)
        n_components = 2**(int(np.log2(sig_len)) + 1)
        sig_hilb_nonan = hilbert(sig_nonan, n_components)[:sig_len]
    else:
        sig_hilb_nonan = hilbert(sig_nonan)

    # Fill in output hilbert with nans on edges
    sig_hilb = restore_nans(sig_hilb_nonan, sig_nans, dtype=complex)

    return sig_hilb