def analytic_signal(vi, windwidth, percover, win):
    nvi = len(vi)
    #h = zeros(vi.shape)
    bli, ble, num = fftsegs(windwidth, percover, nvi)
    for ii in range(len(bli)):
        v = vi[bli[ii]:ble[ii]+1]
        nv = len(v)
        if win == 1:
            fv = fft(v * numpy.hamming(nv))
        else:
            fv = fft(v)
        wind = zeros(v.size)
        # zero negative frequencies, double positive frequencies
        if nv % 2 == 0:
            wind[0] = 1  # keep DC
            wind[(nv / 2)] = 1
            wind[1:(nv / 2)] = 2  # double pos. freq

        else:
            wind[0] = 1
            wind[range(1, (nv + 1) / 2)] = 2
        h = ifft(fv * wind)
    for i in range(len(h)):
        h[i] /= numpy.complex(num[i])
    return h

def filt(frq,tol,data,spd):
    ''' this function filters signals given a filtering frequency (frq), tolerance (tol), data, and sampling freqency (spd) '''
    #define frequency tolerance range
    lowcut = (frq-frq*tol)
    highcut = (frq+frq*tol)
    #conduct fft
    ffta = fft.fft(data)
    bp2 = ffta[:]
    fftb = fft.fftfreq(len(bp2))
    #make every amplitude value 0 that is not in the tolerance range of frequency of interest
    #24 adjusts the frequency to cpd
    for i in range(len(fftb)):
        #spd is samples per day (if hourly = 24)
        if (fftb[i]*spd)>highcut or (fftb[i]*spd)<lowcut:
            bp2[i]=0
    #conduct inverse fft to transpose the filtered frequencies back into time series
    crve = fft.ifft(bp2)    #complex number returned
    #convert back to frequency domain
    fta = fft.fft(crve)
    rl = fta.real
    im = fta.imag
    mag = [math.sqrt(rl[i]**2 + im[i]**2) for i in range(len(rl))] # magnitude
    phase = [math.atan2(im[i],rl[i]) for i in range(len(rl))]       # phase
    yfilt = crve.real       #real component of complex number
    zfilt = crve.imag       #imaginary component of complex nunmber
    return yfilt, zfilt, crve, mag, phase

def InnerProd(ser1, ser2, PSD):
  size = Numeric.shape(ser1)[0]
  pdlen = size/2
  nyquistf = 0.5/15.0   #   !!! hardcoded !!!!
  freqs = Numeric.arange(0,pdlen+1,dtype='d') * (nyquistf / pdlen)
  if(Numeric.shape(ser2)[0] != size):
     print "size of time series must be the same"
     sys.exit(1)
  if(Numeric.shape(PSD)[0] != pdlen):
     print "wrong size of psd: ", pdlen, Numeric.shape(PSD)
     sys.exit(1)
  fourier1 = FFT.fft(ser1)
  fourier2 = FFT.fft(ser2)
  prod = Numeric.zeros(pdlen+1, dtype='d')
  prod[0] = 0.0
  prod[1:pdlen] = numpy.multiply(fourier1[1:pdlen],numpy.conjugate(fourier2[1:pdlen])) + numpy.multiply(fourier1[-1:pdlen:-1],numpy.conjugate(fourier2[-1:pdlen:-1]))
  prod[pdlen] = fourier1[pdlen]*fourier2[pdlen]
  Numeric.divide(prod[1:], PSD, prod[1:]) 
  olap0 = 0.0
  for i in xrange(pdlen):
      if (freqs[i] > fLow and freqs[i]<= fHigh):
           olap0 += prod[i]
  olap0 = 2.0*olap0/float(size)
 # olap0 =  2.0*(numpy.sum(prod[1:]))/float(size) #it must be scaled by dt
  return  olap0

def autocorr_fft(signal, axis = -1):
    """Return full autocorrelation along specified axis. Use fft
    for computation."""
    if N.ndim(signal) == 0:
        return signal
    elif signal.ndim == 1:
        n       = signal.shape[0]
        nfft    = int(2 ** nextpow2(2 * n - 1))
        lag     = n - 1
        a       = fft(signal, n = nfft, axis = -1)
        au      = ifft(a * N.conj(a), n = nfft, axis = -1)
        return N.require(N.concatenate((au[-lag:], au[:lag+1])), dtype = signal.dtype)
    elif signal.ndim == 2:
        n       = signal.shape[axis]
        lag     = n - 1
        nfft    = int(2 ** nextpow2(2 * n - 1))
        a       = fft(signal, n = nfft, axis = axis)
        au      = ifft(a * N.conj(a), n = nfft, axis = axis)
        if axis == 0:
            return N.require(N.concatenate( (au[-lag:], au[:lag+1]), axis = axis), \
                    dtype = signal.dtype)
        else:
            return N.require(N.concatenate( (au[:, -lag:], au[:, :lag+1]), 
                        axis = axis), dtype = signal.dtype)
    else:
        raise RuntimeError("rank >2 not supported yet")

def deconvolve(signal, HRFs, SNRs):
	"""Deconvolve signal using Wiener deconvolution."""
	H = fft(HRFs, len(signal), axis=0)
	wiener_filter = np.conj(H) / (H * np.conj(H) + 1 / SNRs**2)
	deconvolved = np.real(ifft(wiener_filter * fft(signal, axis=0), axis=0))

	return deconvolved

def XcorrNormed(x,y):
    """ method useful to compare binary vectors """
    if len(x) != len(y):
        raise ValueError("signals should be the same size")

    X = fft(x)
    Y = fft(y)

    # Compute classic cross-correlation
    classic_xcorr = (ifft(X*(Y.conj()))).real
    maxlag = len(X)/2
    classic_xcorr = concatenate((classic_xcorr[-maxlag:],classic_xcorr[0:maxlag]));
    ind = abs(classic_xcorr).argmax();
#    ind = classic_xcorr.argmax()
    val = classic_xcorr[ind]

    # normalize
    normx = math.sqrt(sum((x)**2))
    normy = math.sqrt(sum((y)**2))
#    print 'Norm of ' , normx*normy , ' for a value found of ' , val
#    val = float(val)/(float(len(nonzero(x)[0]) + len(nonzero(y)[0]))/2)
    if (normx * normy) != 0:
        val = float(val)/(normx * normy)

    return classic_xcorr , (len(X)/2 - ind) , val;

def convolve(A,B,):
    N = A.shape[0]
    fA = fft.fft(A)
    fB = fft.fft(B)
    one = np.ones(N)
    one[1::2] = -1
    return np.real(fft.ifft(fA*fB)/float(N))

def convFFT(x,y):
    nx=len(x)
    xf=np.pad(x,(nx/2,nx/2),mode='constant')
    yf=np.pad(y,(nx/2,nx/2),mode='constant')
    xf=(fft.fft(fft.fftshift(xf)))
    yf=(fft.fft(fft.fftshift(yf)))
    return fft.fftshift(np.real((fft.ifft(xf*yf))))[nx/2:3*nx/2]

def fft16_a4_2(x):
    global mul_count
    tmp = x
    w3, w2, w1 = get_w_fft16_r2()
    mul_count += flops4Muls(w1)
    mul_count += flops4Muls(w2)
    mul_count += flops4Muls(w3)

    for j in range(8):
        tmp[j:j+9:8] = fft(tmp[j:j+9:8])
    tmp *= w1

    for i in range(2):
        for j in range(4):
            tmp[i*8+j:i*8+j+5:4] = fft(tmp[i*8+j:i*8+j+5:4])
    tmp *= w2

    for i in range(4):
        for j in range(2):
            tmp[i*4+j:i*4+j+3:2] = fft(tmp[i*4+j:i*4+j+3:2])
    tmp *= w3

    for i in range(8):
        tmp[i*2:i*2+2] = fft(tmp[i*2:i*2+2])

    return bit_revers(tmp)

def calc_cospectrum(a,b,**kwargs):
    nfft_time = np.shape(a)[0]
    if 'nfft_time' in kwargs.keys():
        nfft_time = kwargs['nfft_time']
    nlon = np.shape(a)[1]
    fa = fft.fft(a,axis=1)
    fb = fft.fft(b,axis=1)
    nomega = nfft_time/2+1
    nk = nlon/2+1 
    cfa = np.real(fa[:,:nk])
    sfa = np.imag(fa[:,:nk])
    cfb = np.real(fb[:,:nk])
    sfb = np.imag(fb[:,:nk])
    pp = np.zeros([nomega,nk])
    pn = np.zeros([nomega,nk])
    for i in range(nk):
        omega, pcacb = signal.csd(cfa[:,i],cfb[:,i],nperseg=nfft_time)
        omega, psasb = signal.csd(sfa[:,i],sfb[:,i],nperseg=nfft_time)
        omega, pcasb = signal.csd(cfa[:,i],sfb[:,i],nperseg=nfft_time)
        omega, psacb = signal.csd(sfa[:,i],cfb[:,i],nperseg=nfft_time)
        pp[:,i] = np.real(pcacb)+np.real(psasb)+np.imag(pcasb)-np.imag(psacb)
        pn[:,i] = np.real(pcacb)+np.real(psasb)-np.imag(pcasb)+np.imag(psacb)
    p_all = np.zeros([nomega*2,nk])
    p_all[:nomega,:] = np.flipud(pn)
    p_all[nomega:,:] = pp
    sigma = 0.25/np.pi*nomega
    x = np.linspace(-nomega/2,nomega/2,nomega)
    gauss = np.exp(-x**2/(2*sigma**2))
    gauss = gauss/np.sum(gauss)
    for i in range(nk):
        p_all[:,i] = np.convolve(p_all[:,i],gauss,mode='same')
    omega_all = np.concatenate((np.flipud(-omega),omega))
    return p_all,omega_all

 def potential(self):
     density=self.rho_den()
     soft_pot=self.softened_pot()
     ft1=fft(density)
     ft2=fft(soft_pot)
     conv_pot=ifft(ft1*ft2)
     return conv_pot

def b(v, F, I, idx):
    x_i = np.dot(v, F[idx][:, 0])
    y_i = np.dot(v, F[idx][:, 1])
    z_i = np.dot(v, F[idx][:, 2])
    vec = [x_i, y_i, z_i]

    return ifft(fft(I[idx] * fft(l(vec))))

def fwgn_model(fm,fs,N):
    N = int(N)
    Nfft = 2**max(3,nextpow2(2*fm/fs*N))
    Nifft = math.ceil(Nfft*fs/(2*fm))

    doppler_coeff = np.sqrt(doppler_filter(fm,Nfft))

    CGI, CGQ = fft(randn(Nfft)), fft(randn(Nfft))
    f_CGI = CGI * doppler_coeff
    f_CGQ = CGQ * doppler_coeff

    del CGI, CGQ, doppler_coeff
    gc.collect()

    tzeros = np.zeros(abs(Nifft-Nfft))
    filtered_CGI = np.hstack((f_CGI[:Nfft//2], tzeros, f_CGI[Nfft//2:]))
    filtered_CGQ = np.hstack((f_CGQ[:Nfft//2], tzeros, f_CGQ[Nfft//2:]))

    del tzeros, f_CGI, f_CGQ
    gc.collect()

    hI, hQ = ifft(filtered_CGI), ifft(filtered_CGQ)

    del filtered_CGI, filtered_CGQ
    gc.collect()

    rayEnvelope = np.abs(np.abs(hI) + 1j * hQ)
    rayRMS = math.sqrt(np.mean(rayEnvelope[:N]**2))

    # h_{I}(t) - jh_{Q}(t)
    # Here we have the phase shift of pi/2 when multiplying the imaginary
    # portion by -1j
    h = (np.real(hI[:N]) - 1j * np.real(hQ[:N]))/rayRMS

    return h

    def magabs_cs():
        #b.line_plot("magabs_cs", a.frequency, a.MagAbs[:, 0])
        #b.line_plot("magabs_cs", a.frequency, a.MagAbs[:, 257])

        if 0:
            myifft=fft.ifft(a.Magcom[:,500])
            myifft[50:-50]=0.0
            #myifft[:20]=0.0
            #myifft[-20:]=0.0
            bg=fft.fft(myifft)
            filt=[]
            for n in range(len(a.yoko)):
                myifft=fft.ifft(a.Magcom[:,n])
                #b.line_plot("ifft", absolute(myifft))
                myifft[50:-50]=0.0
                #myifft[:20]=0.0
                #myifft[-20:]=0.0
                filt.append(absolute(fft.fft(myifft)-bg))
            b.colormesh("filt", a.frequency, a.yoko, filt)
        if 1:
            myifft=fft.ifft(mag[:,500])
            myifft[40:-40]=0.0
            myifft[:20]=0.0
            myifft[-20:]=0.0
            bg=fft.fft(myifft)
            filt=[]
            for n in range(len(yok)):
                myifft=fft.ifft(mag[:,n])
                #b.line_plot("ifft", absolute(myifft))
                myifft[50:-50]=0.0
                #myifft[:20]=0.0
                #myifft[-20:]=0.0
                filt.append(absolute(fft.fft(myifft)))
            b.colormesh("filt", frq, yok, filt)

 def __rmul__(self,other):
     if isinstance(other,HRR):
         x=ifft(fft(self.v)*fft(other.v)).real
         x=x/norm(x)
         return HRR(data=x)
     else:
         return HRR(data=self.v*other)

def ccor(ts1, ts2):
    '''
    Given two standardized time series, computes their cross-correlation using
    fast fourier transform. Assume that the two time series are of the same
    length.

    Parameters
    ----------
    ts1 : TimeSeries
        A standardized time series
    ts2 : TimeSeries
        Another standardized time series

    Returns
    -------
    The two time series' cross-correlation.
    '''
    # calculate fast fourier transform of the two time series
    fft_ts1 = nfft.fft(ts1.valuesseq)
    fft_ts2 = nfft.fft(ts2.valuesseq)

    # print(len(ts1))
    # print(len(ts2))
    # assert len(ts1) == len(ts2)

    # return cross-correlation, i.e. the convolution of the first fft
    # and the conjugate of the second
    return ((1 / (1. * len(ts1))) *
            nfft.ifft(fft_ts1 * np.conjugate(fft_ts2)).real)

    def InvPotentialVorticity(self,field, length=None):

        if length is None:
            length = 2*pi;

        N = shape(field)[0];

        k = array(range(N),dtype=complex128);
        k = concatenate((range(0,N/2),range(-N/2,0)));
        k *= (2*pi)/length;

        [KX, KY] = meshgrid(k,k);

        """ We are trying to solve d_yy(eta) + d_xx(eta) - eta = p
        Therefore, in Fourier space, it will become
        (-(kx^2 + ky^2) - 1)etaHat = pHat
        """
        delsq = -(KX*KX + KY*KY) - 1 ;
        delsq[0,0] = 1;

        tmp = fft(field,axis=0);
        tmp = fft(tmp,axis=1);

        tmp = tmp/delsq;

        tmp = ifft(tmp,axis=1);
        tmp = ifft(tmp,axis=0);

        return tmp.real;

def D2dWT(x, seeds, seedt, OM, om):
    '''
    Discrete two-dimenstional Wavelet Transform
    -   x: input signal (one-dimensional)
    -   ws: spectral wavelet
    -   wt: temporal wavelet
    '''
    N = len(x)
    ws_t = seeds(np.arange(-100, 100)/100.0, OM)
    ws = np.zeros(200)
    '''
    ws[100:] = ws_t[:100]
    ws[:100] = ws_t[100:]
    '''
    ws = ws_t
    T = int(3.0 / om * 44100)
    wt = seedt(np.arange(T)/44100.0, om)
    WS = fft(ws)
    plt.plot(WS); plt.show()
    w = np.ndarray(shape=(T, 200), dtype='complex')
    for t in range(T):
        w[t] = WS * wt[t]
    y = np.zeros(N-T) 
    for tau in range(N-T):
        y[tau] = sum(fft(x[tau:tau+200]) * w[tau])
    return y

def ild_bare(hrir, cf, **kwdargs):
    '''
    ILD computation routine. called by ild that handles multiprocessing,...
    '''
    samplerate = hrir.samplerate

    # perform some checks and special cases
    if (hrir[:,0] == hrir[:,1]).all():
        return np.zeros(len(cf))

    if (abs(hrir[:,0])<= 10e-6).all() or  (abs(hrir[:,1])<=10e-6).all():
        log_debug('Blank hrirs detected, output will be weird')

    if not isinstance(hrir, Sound):
        hrir = Sound(hrir, samplerate = samplerate)

    fb = Gammatone(Repeat(hrir, len(cf)), np.hstack((cf, cf)))
    filtered_hrirset = fb.process()
    
    ilds = []
    for i in range(len(cf)):
        left = filtered_hrirset[:, i]
        right = filtered_hrirset[:, i+len(cf)]
        # This FFT stuff does a correlate(left, right, 'full')
        Lf = fft(np.hstack((left, np.zeros(len(left)))))
        Rf = fft(np.hstack((right[::-1], np.zeros(len(right)))))
        C = ifft(Lf*Rf).real
        ilds.append(np.sqrt(np.amax(C)/sum(right**2)))
    ilds = np.array(ilds)
    return ilds

def STRW(x, om, OM, up):
    M = len(x)
    N = len(x[0])
   
    # 2d convolution:
    y = np.ndarray(shape=(M, N), dtype='complex')
    h = get_wavelet(om, OM, up, M, N)
    
    # temporal convolution:
    for m in range(M):
        X = fft(x[m])
        H = fft(h[m])
        Z = X*H

        z = ifft(Z)
        y[m] = z

    # spectral convolution
    x = y.transpose()
    h = h.transpose()
    y = y.transpose()
    for n in range(N):
        X = fft(x[n])
        H = fft(h[n])
        Z = X*H
        
        z = ifft(Z)
        y[n] = z

    y = y.transpose()

    return y