def test_glm():
    """
    Test PAC function: GLM
    1. Confirm consistency of output with example data
    2. Confirm consistency of output with example data using iir filter
    3. Confirm PAC=1 when expected
    4. Confirm PAC=0 when expected
    """
    # Load data
    data = np.load(os.path.dirname(pacpy.__file__) + '/tests/exampledata.npy')
    assert np.allclose(
        glm(data, data, (13, 30), (80, 200)), 0.03191, atol=10 ** -5)
    assert np.allclose(
        glm(data, data, (13, 30), (80, 200), filterfn=butterf), 0.03476, atol=10 ** -5)

    # Test that the GLM function outputs close to 0 and 1 when expected
    lo, hi = genPAC1(glm_bias=True)
    assert glm(lo, hi, (4, 6), (90, 110)) > 0.99

    lo, hi = genPAC0()
    assert glm(lo, hi, (4, 6), (90, 110)) < 0.01
    
    # Test that Filterfn = False works as expected
    datalo = firf(data, (13,30))
    datahi = firf(data, (80,200))
    pha = np.angle(hilbert(datalo))
    amp = np.abs(hilbert(datahi))
    assert np.allclose(
        glm(pha, amp, (13, 30), (80, 200), filterfn=False),
        glm(data, data, (13, 30), (80, 200)), atol=10 ** -5)

def Envelope(wsyn, wobs, nt, dt, eps=0.05):
    # envelope difference
    esyn = abs(hilbert(wsyn))
    eobs = abs(hilbert(wobs))
    etmp = (esyn - eobs)/(esyn + eps*esyn.max())
    wadj = etmp*wsyn - np.imag(hilbert(etmp*np.imag(hilbert(wsyn))))
    return wadj

def test_ozkurt():
    """
    Test PAC function: Ozkurt
    1. Confirm consistency of output with example data
    2. Confirm consistency of output with example data using iir filter
    3. Confirm PAC=1 when expected
    4. Confirm PAC=0 when expected
    """
    # Load data
    data = np.load(os.path.dirname(pacpy.__file__) + '/tests/exampledata.npy')
    assert np.allclose(
        ozkurt(data, data, (13, 30), (80, 200)), 0.07548, atol=10 ** -5)
    assert np.allclose(
        ozkurt(data, data, (13, 30), (80, 200), filterfn=butterf), 0.07555, atol=10 ** -5)

    # Test that the Ozkurt PAC function outputs close to 0 and 1 when expected
    lo, hi = genPAC1(phabias=.2, fhi=300)
    hif = firf(hi, (100, 400))
    amp = np.abs(hilbert(hif))
    weight = (np.sqrt(len(amp)) * np.sqrt(np.sum(amp ** 2))) / np.sum(amp)
    assert ozkurt(lo, hi, (4, 6), (100, 400)) * weight > 0.99

    lo, hi = genPAC0()
    assert ozkurt(lo, hi, (4, 6), (90, 110)) < 0.001
    
    # Test that Filterfn = False works as expected
    datalo = firf(data, (13,30))
    datahi = firf(data, (80,200))
    pha = np.angle(hilbert(datalo))
    amp = np.abs(hilbert(datahi))
    assert np.allclose(
        ozkurt(pha, amp, (13, 30), (80, 200), filterfn=False),
        ozkurt(data, data, (13, 30), (80, 200)), atol=10 ** -5)

def compute_envelope_matrix_theo(dsyn, obsd, synt, dt, win_idx, taper):
    """
    Compute envelope measurements matrix H and misfit vector G theoretically
    as stated in Appendix in Qinya's paper.
    Attension: not used!!! BUG inside!!!
    """
    istart = win_idx[0]
    iend = win_idx[1]

    syn_array = synt.copy()
    obs_array = obsd.copy()

    syn_analytic = hilbert(taper * syn_array[istart:iend])
    syn_hilbert = np.imag(syn_analytic)
    syn_env = np.abs(syn_analytic)
    dsyn_hilbert = np.imag(hilbert(dsyn))
    env_derivss = \
        syn_env ** (-0.5) * (syn_array[istart:iend] * dsyn +
                             syn_hilbert * dsyn_hilbert)

    A1 = np.dot(env_derivss, env_derivss.transpose()) * dt
    b1 = np.sum(
        (np.abs(hilbert(taper * obs_array[istart:iend])) -
         np.abs(hilbert(taper * syn_array[istart:iend]))) *
        env_derivss * dt, axis=1)
    return A1, b1

def phase_amplitude_coupling(fser, gser, lag=0):
    """
    Compute the product of two time series for calculation of phase-amplitude
    coupling. That is, if the Hilbert transform of fser is A_f exp(\phi_f),
    then the phase-amplitude product of fser and gser is 
    f * g = A_f exp(\phi_g). Lag is the offset of f's amplitude from g's
    phase: A(t) exp(\phi(t - lag)). This is useful for examining asymptotic
    statistics in large-lag regime, where biases in the joint distribution
    of A_f and \phi_g are not due to phase-amplitude coupling, but to the 
    individual signals (cf. Canolty et al., Sciece, 2006, SI).
    Function returns a series of the same length, but first and last lag 
    elements are NaN.
    """
    fh = ssig.hilbert(fser)
    gh = ssig.hilbert(gser)

    Af = np.abs(fh)
    ephig = gh / np.abs(gh)

    if lag == 0:
        pac = Af * ephig
    else:
        pac = Af * np.roll(ephig, lag)
        pac[:lag] = np.nan
        pac[-lag:] = np.nan

    return pac

def test_plv():
    """
    Test PAC function: PLV.
    1. Confirm consistency of output with example data
    2. Confirm consistency of output with example data using firfls filter
    3. Confirm PAC=1 when expected
    4. Confirm PAC=0 when expected
    """
    # Load data
    data = np.load(os.path.dirname(pacpy.__file__) + '/tests/exampledata.npy')
    assert np.allclose(
        plv(data, data, (13, 30), (80, 200)), 0.25114, atol=10 ** -5)
    assert np.allclose(
        plv(data, data, (13, 30), (80, 200), filterfn=firfls), 0.24715, atol=10 ** -5)

    # Test that the PLV function outputs close to 0 and 1 when expected
    lo, hi = genPAC1()
    assert plv(lo, hi, (4, 6), (90, 110)) > 0.99

    lo, hi = genPAC0()
    assert plv(lo, hi, (4, 6), (90, 110)) < 0.01

    # Test that Filterfn = False works as expected
    datalo = firf(data, (13, 30))
    datahi = firf(data, (80, 200))
    datahiamp = np.abs(hilbert(datahi))
    datahiamplo = firf(datahiamp, (13, 30))
    pha1 = np.angle(hilbert(datalo))
    pha2 = np.angle(hilbert(datahiamplo))
    pha1, pha2 = _trim_edges(pha1, pha2)
    assert np.allclose(
        plv(pha1, pha2, (13, 30), (80, 200), filterfn=False),
        plv(data, data, (13, 30), (80, 200)), atol=10 ** -5)

def test_mi_canolty():
    """
    Test PAC function: Canolty MI
    1. Confirm consistency of output with example data
    2. Confirm consistency of output with example data using iir filter
    3. Confirm PAC=1 when expected
    4. Confirm PAC=0 when expected
    """
    # Load data
    data = np.load(os.path.dirname(pacpy.__file__) + '/tests/exampledata.npy')
    assert np.allclose(
        mi_canolty(data, data, (13, 30), (80, 200)), 1.10063, atol=10 ** -5)
    assert np.allclose(mi_canolty(
        data, data, (13, 30), (80, 200), filterfn=butterf), 1.14300, atol=10 ** -5)

    # Test that the Canolty MI function outputs close to 0 and 1 when expected
    lo, hi = genPAC1(phabias=.2, fhi=300)
    hif = firf(hi, (100, 400))
    amp = np.abs(hilbert(hif))
    assert mi_canolty(lo, hi, (4, 6), (100, 400)) / np.mean(amp) > 0.99

    lo, hi = genPAC0()
    assert mi_canolty(lo, hi, (4, 6), (90, 110)) < 0.001
    
    # Test that Filterfn = False works as expected
    datalo = firf(data, (13,30))
    datahi = firf(data, (80,200))
    pha = np.angle(hilbert(datalo))
    amp = np.abs(hilbert(datahi))
    assert np.allclose(
        mi_canolty(pha, amp, (13, 30), (80, 200), filterfn=False),
        mi_canolty(data, data, (13, 30), (80, 200)), atol=10 ** -5)

def pa_series(lo, hi, f_lo, f_hi, fs=1000, filterfn=None, filter_kwargs=None):
    """
    Calculate the phase and amplitude time series

    Parameters
    ----------
    lo : array-like, 1d
        The low frequency time-series to use as the phase component
    hi : array-like, 1d
        The high frequency time-series to use as the amplitude component
    f_lo : (low, high), Hz
        The low frequency filtering range
    f_hi : (low, high), Hz
        The low frequency filtering range
    fs : float
        The sampling rate (default = 1000Hz)
    filterfn : function
        The filtering function, `filterfn(x, f_range, filter_kwargs)`
    filter_kwargs : dict
        Keyword parameters to pass to `filterfn(.)`

    Returns
    -------
    pha : array-like, 1d
        Time series of phase
    amp : array-like, 1d
        Time series of amplitude
        
    Usage
    -----
    >>> import numpy as np
    >>> from scipy.signal import hilbert
    >>> from pacpy.pac import pa_series
    >>> t = np.arange(0, 10, .001) # Define time array
    >>> lo = np.sin(t * 2 * np.pi * 6) # Create low frequency carrier
    >>> hi = np.sin(t * 2 * np.pi * 100) # Create modulated oscillation
    >>> hi[np.angle(hilbert(lo)) > -np.pi*.5] = 0 # Clip to 1/4 of cycle
    >>> pha, amp = pa_series(lo, hi, (4,8), (80,150))
    >>> print pha
    [-1.57079633 -1.53192376 -1.49301802 ..., -1.64840672 -1.6095709 -1.57079634]
    """

    # Arg check
    _x_sanity(lo, hi)
    _range_sanity(f_lo, f_hi)

    # Filter setup
    if filterfn is None:
        filterfn = firf
        filter_kwargs = {}

    # Filter
    xlo = filterfn(lo, f_lo, fs, **filter_kwargs)
    xhi = filterfn(hi, f_hi, fs, **filter_kwargs)

    # Calculate phase time series and amplitude time series
    pha = np.angle(hilbert(xlo))
    amp = np.abs(hilbert(xhi))

    return pha, amp

def envelope(data):
    """
    Envelope of a signal.

    Computes the envelope of the given data which can be windowed or
    not. The envelope is determined by the absolute value of the analytic
    signal of the given data.

    If data are windowed the analytic signal and the envelope of each
    window is returned.

    :type data: :class:`~numpy.ndarray`
    :param data: Data to make envelope of.
    :return: **A_cpx, A_abs** - Analytic signal of input data, Envelope of
        input data.
    """
    nfft = util.next_pow_2(data.shape[-1])
    a_cpx = np.zeros((data.shape), dtype=np.complex64)
    a_abs = np.zeros((data.shape), dtype=np.float64)
    if len(data.shape) > 1:
        i = 0
        for row in data:
            a_cpx[i, :] = signal.hilbert(row, nfft)
            a_abs[i, :] = abs(signal.hilbert(row, nfft))
            i = i + 1
    else:
        a_cpx = signal.hilbert(data, nfft)
        a_abs = abs(signal.hilbert(data, nfft))
    return a_cpx, a_abs

def envelope_function(x, y, z, ttime, dt):
    '''Compute envelope function based on hilbert transform

    Parameters
    ----------
    x: numpy.array
      time series

    y: numpy.array
      time series     

    z: numpy.array
      time series     

    ttime: float
      total time for integral

    dt: float
      time interval of time-series

    Returns
    -------
    (f, out): (numpy.array, numpy.array)
      frequency vector for plotting
      fourier amplitude spectrum
    '''
    sim_ind = np.floor(ttime/dt)
    analytic_x = np.absolute(sig.hilbert(x))
    analytic_y = np.absolute(sig.hilbert(y))
    analytic_z = np.absolute(sig.hilbert(z))
    return (np.sum(analytic_x[:sim_ind])*dt,
            np.sum(analytic_y[:sim_ind])*dt,
            np.sum(analytic_z[:sim_ind])*dt)

def ediff(wsyn, wobs, nt, dt, eps=0.05):
    # envelope difference
    esyn = abs(_signal.hilbert(wsyn))
    eobs = abs(_signal.hilbert(wobs))
    etmp = (esyn - eobs)/(esyn + eps*esyn.max())
    wadj = etmp*wsyn - _np.imag(_signal.hilbert(etmp*_np.imag(_signal.hilbert(wsyn))))
    return wadj

 def _normalize(self):
     SP1 = np.fft.fft(hilbert(self.s1), axis=0)
     SP2 = np.fft.fft(hilbert(self.s2), axis=0)
     indmin = 1 + int(np.round(self.fmin * (self.ts.shape[0] - 2)))
     indmax = 1 + int(np.round(self.fmax * (self.ts.shape[0] - 2)))
     sp1_ana = SP1[(indmin - 1):indmax]
     sp2_ana = SP2[(indmin - 1):indmax]
     return sp1_ana, sp2_ana

    def _process(self):
        print(" - - - In Processing : Hight={}, Width={}".format(self.yarp_image.height(), self.yarp_image.width()))
        
        # Make this into a numpy array
        self._convert()

        beamedAudio = beamformer(self.matrix, rms=False)
        numSamps = beamedAudio.shape[2]

        showEnv = True

        if showEnv:
            
            oneBeam = np.zeros((1, 1, numSamps), dtype=np.float32)
            oneBeam[0][0] = beamedAudio[7][21]

            print(" - - - - Begin Hilbert.")
            analytic_signal = hilbert(oneBeam)
            amplitude_envelope = np.abs(analytic_signal)
            band_passed_amp = bandPass(amplitude_envelope, 5.0, numSamps, 48000)
            
            amp_oneBeam = oneBeam[0][0].copy()
            amp_oneBeam[amp_oneBeam < 0.0] = 0.0

            # Plot.
            self.fig.clear()

            ax0 = self.fig.add_subplot(311)
            ax0.set_ylim(-1, 1)
            ax0.plot(oneBeam[0][0])
            ax0.plot(amp_oneBeam)


            ax1 = self.fig.add_subplot(312)
            ax1.set_ylim(-1, 1)
            ax1.plot(analytic_signal[0][0])
            ax1.plot(amplitude_envelope[0][0])
            
            ax2 = self.fig.add_subplot(313)
            ax2.set_ylim(-1, 1)
            ax2.plot(band_passed_amp[0][0])


        else:
            print(" - - - - Begin Hilbert.")
            analytic_signal    = hilbert(beamedAudio)
            amplitude_envelope = np.abs(analytic_signal)

            print(" - - - - Begin Band Pass.")
            band_passed_amp = bandPass(amplitude_envelope, 5.0, numSamps, 48000)

            reduced_band_pass = np.sqrt(np.sum(band_passed_amp**2, axis=2)) / numSamps

            self.fig.clear()
            plt.imshow(reduced_band_pass)


        plt.pause(0.005)

def InstantaneousPhase(wsyn, wobs, nt, dt, eps=0.05):
    # instantaneous phase 
    r = np.real(hilbert(wsyn))
    i = np.imag(hilbert(wsyn))
    phi_syn = np.arctan2(i,r)

    r = np.real(hilbert(wobs))
    i = np.imag(hilbert(wobs))
    phi_obs = np.arctan2(i,r)

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

def symmetrydemo():
    a=sin(linspace(-5*pi,5*pi,10000))
    b=a+2
    c=a-0.5
    ah,bh,ch=hilbert(a),hilbert(b),hilbert(c)
    ph_a,ph_b,ph_c=unwrap(angle(ah)),unwrap(angle(bh)),unwrap(angle(ch))
    omega_a=diff(ph_a)
    omega_b=diff(ph_b)
    omega_c=diff(ph_c)
    subplot(211),plot(ph_a),plot(ph_b),plot(ph_c)
    subplot(212),plot(omega_a),plot(omega_b),plot(omega_c)
    grid()
    show()
    return a,b,c

def _filter_ph_am(xph, xam, f_ph, f_am, sfreq, filterfn=None, kws_filt=None):
    """Aux function for phase/amplitude filtering for one pair of channels"""
    from pacpy.pac import _range_sanity
    from scipy.signal import hilbert
    filterfn = band_pass_filter if filterfn is None else filterfn
    kws_filt = {} if kws_filt is None else kws_filt

    # Filter the two signals + hilbert/phase
    _range_sanity(f_ph, f_am)
    xph = filterfn(xph, sfreq, *f_ph)
    xam = filterfn(xam, sfreq, *f_am)

    xph = np.angle(hilbert(xph))
    xam = np.abs(hilbert(xam))
    return xph, xam

def square_envelope(corr_o,corr_s,g_speed,
    window_params):
    success = False
    env_s = corr_s.data**2 + np.imag(hilbert(corr_s.data))**2
    env_o = corr_o.data**2 + np.imag(hilbert(corr_o.data))**2
    d_env_1 =  2. * corr_s.data 
    d_env_2 =  (2. * np.imag(hilbert(corr_s.data)))

    u1 = (env_s - env_o) * d_env_1
    u2 = np.imag(hilbert((env_s - env_o) * d_env_2))

    adjt_src = u1 - u2
    
    success = True
    return adjt_src, success

 def __init__(self,hfoObj):
     #signal = sig.detrend(hfoObj.waveform[hfoObj.start_idx:hfoObj.end_idx,0]) # detrending
     fs = hfoObj.sample_rate        
    
     signal = sig.detrend(hfoObj.waveform[3*fs/4:5*fs/4,0])
     PhaseFreqVector= np.arange(1,31,1)
     AmpFreqVector= np.arange(30,990,5)
     PhaseFreq_BandWidth=1
     AmpFreq_BandWidth=10
     Comodulogram=np.zeros((PhaseFreqVector.shape[0],AmpFreqVector.shape[0]))
     nbin=18
     position=np.zeros(nbin)
     winsize = 2*np.pi/nbin
     for j in range(nbin):
         position[j] = -np.pi+j*winsize;
     PHASES = np.zeros((PhaseFreqVector.shape[0],signal.shape[0]))
     for idx,Pf1 in enumerate(PhaseFreqVector):
         print Pf1,
         Pf2 = Pf1 + PhaseFreq_BandWidth
         if signal.shape[0] > 18*np.fix(fs/Pf1):
             b = sig.firwin(3*np.fix(fs/Pf1),[Pf1,Pf2],pass_zero=False,window=('kaiser',0.5),nyq=fs/2)
         else:
             b = sig.firwin(signal.shape[0]/6,[Pf1,Pf2],pass_zero=False,window=('kaiser',0.5),nyq=fs/2)
         PhaseFreq = sig.filtfilt(b,np.array([1]),signal)
         Phase=np.angle(sig.hilbert(PhaseFreq))
         PHASES[idx,:]=Phase;
     print    
     for idx1,Af1 in enumerate(AmpFreqVector):
         print Af1,
         Af2 = Af1 + AmpFreq_BandWidth
         if signal.shape[0] > 18*np.fix(fs/Af1):
             b = sig.firwin(3*np.fix(fs/Af1),[Af1,Af2],pass_zero=False,window=('kaiser',0.5),nyq=fs/2)
         else:
             b = sig.firwin(np.fix(signal.shape[0]/6),[Af1,Af2],pass_zero=False,window=('kaiser',0.5),nyq=fs/2)
         AmpFreq = sig.filtfilt(b,np.array([1]),signal)
         Amp=np.abs(sig.hilbert(AmpFreq))
         for idx2,Pf1 in enumerate(PhaseFreqVector):
             Phase = PHASES[idx2]
             MeanAmp = np.zeros(nbin)
             for j in range(nbin):
                 bol1 = Phase < position[j]+winsize
                 bol2 = Phase >= position[j]
                 I = np.nonzero(bol1 & bol2)[0]
                 MeanAmp[j]=np.mean(Amp[I])
             #MI=(np.log(nbin)-(-np.sum((MeanAmp/np.sum(MeanAmp))*np.log((MeanAmp/np.sum(MeanAmp))))))/np.log(nbin)
             MI =np.log(nbin)-(stat.entropy(MeanAmp)/np.log(nbin))
             Comodulogram[idx2,idx1]=MI;
     plt.contourf(PhaseFreqVector+PhaseFreq_BandWidth/2,AmpFreqVector+AmpFreq_BandWidth/2,Comodulogram.T,100)

def amplitude_envelope(signal, fs):
    """Instantaneous amplitude of tone.
    
    .. seealso:: :func:`scipy.signal.hilbert`
    
    """
    return np.abs(hilbert(signal))

def qea(im):
    H = ss.hilbert(im,axis = 2)
    H = im+1j*H
    ia = np.abs(H)
    ip = np.angle(H)

    h1col = H[1:-1,:,:]
    h0col = H[:-2,:,:]
    h2col = H[2:,:,:]
    ifColSign = np.sign(np.real((h0col-h2col)/(2j*h1col)))
    ifCol = np.arccos((h2col+h0col)/(2*h1col))
    ifCol = (np.abs(ifCol)*ifColSign)/np.pi/2

    ifCol = np.pad(ifCol,((1,1),(0,0),(0,0)), mode='reflect')
    
    h0row = H[:,:-2,:]
    h1row = H[:,1:-1,:]
    h2row = H[:,2:,:]
    #ifxSign = np.sign(np.real((h2x-h0x)/(2j*h1x)))
    ifRow = np.arccos((h2row+h0row)/(2*h1row))
    ifRow = (np.abs(ifRow))/np.pi/2

    ifRow = np.pad(ifRow,((0,0),(1,1),(0,0)), mode='reflect')

    h0time = H[:,:,:-2]
    h1time = H[:,:,1:-1]
    h2time = H[:,:,2:]
    #ifxSign = np.sign(np.real((h2x-h0x)/(2j*h1x)))
    ifTime = np.arccos((h2time+h0time)/(2*h1time))
    ifTime = (np.abs(ifTime))/np.pi/2

    ifTime = np.pad(ifTime,((0,0),(0,0),(1,1)), mode='reflect')
    
    return(ia,ip,ifRow,ifCol,ifTime)