def algoChannelSelection(left, right):

    ''' Algorithm which automatically selects the channel with dominant vocals from a stereo flamenco recording
    based on spectral band energies as described in section 2-A-I of

    Kroher, N. & Gomez, E. (2016). Automatic Transcription of Flamenco Singing from Polyphonic Music Recordings.
    ACM / IEEE Transactions on Audio, Speech and Language Processing, 24(5), pp. 901-913.

    :param left: samples of the left audio channel in 44.1kHz
    :param right: samples of the right audio channel in 44.1kHz
    :return: index of the dominant vocal channel (0 = left, 1 = right)
    '''

    # PARAMETERS
    fs = 44100 # sample rate
    wSize = 2048 # window size in samples
    hSize = 2048 # hop size in samples
    fftSize = 2048 # FFT size
    freqGuitLow = 80.0 # lower bound for guitar band
    freqGuitHigh = 400.0 # upper bound for guitar band
    freqVocLow = 500.0 # lower bound for vocal band
    freqVocHigh = 6000.0 # higher bound for vocal band

    # INIT
    window = hanning(wSize)
    numFrames = int(math.floor(float(len(left))/float(wSize)))
    # bin indices corresponding to freqeuncy band limits
    indGuitLow = int(round((freqGuitLow/fs)*fftSize))
    indGuitHigh = int(round((freqGuitHigh/fs)*fftSize))
    indVocLow = int(round((freqVocLow/fs)*fftSize))
    indVocHigh = int(round((freqVocHigh/fs)*fftSize))

    # frame-wise computation of the spectral band ratio
    sbrL = []
    sbrR = []
    for i in range(0,numFrames-100):
        frameL = left[i*hSize:i*hSize+wSize]
        specL = fft(frameL*window) / fftSize
        specL = abs(specL * conj(specL))
        guitMag = sum(specL[indGuitLow:indGuitHigh],0)
        vocMag = sum(specL[indVocLow:indVocHigh],0)
        sbrL.append(20*math.log10(vocMag/guitMag))
        frameR = right[i*hSize:i*wSize+wSize]
        specR = fft(frameR*window) / fftSize
        specR = abs(specR * conj(specR))
        guitMag = sum(specR[indGuitLow:indGuitHigh],0)
        vocMag = sum(specR[indVocLow:indVocHigh],0)
        sbrR.append(20*math.log10(vocMag/guitMag))

    # select channel based on mean SBR
    if mean(sbrL)>=mean(sbrR):
        ind = 0
    else:
        ind = 1

    return ind

def GetEigSys(filename,gsfile=None,Nsamp=1,channel=None,wavefile=None,q=None):
    if type(filename)==str:
        filename=[filename]
    hfile=h5py.File(filename[0],'r')
    attr=GetAttr(filename[0])
    if channel==None:
        channel=attr['channel']
    dpath,args=GetStat(filename,Nsamp)
    dat=sc.array(hfile["/rank-1/data-0"])
    hfile.close()
    N=int(sc.shape(dat)[0]/2)
    L=attr['L']
    shift=None
    if 'phasex' in attr.keys():
        shift=[attr['phasex']/2.0,attr['phasey']/2.0]
    else:
        shift=[attr['phase_shift_x']/2.0,attr['phase_shift_y']/2.0]
    H=sc.zeros([Nsamp,N,N],complex)
    O=sc.zeros([Nsamp,N,N],complex)
    E=sc.zeros([Nsamp,N])
    V=sc.zeros([Nsamp,N,N],complex)
    for sample,b in enumerate(args):
        for d in b:
            hfile=h5py.File(dpath[d][0],'r')
            dat=hfile[dpath[d][1]]
            H[sample,:,:]+=dat[0:N,0:2*N:2]+1j*dat[0:N,1:2*N:2]
            O[sample,:,:]+=dat[N:2*N,0:2*N:2]+1j*dat[N:2*N,1:2*N:2]
            hfile.close()
        H[sample,:,:]=0.5*(H[sample,:,:]+sc.conj(H[sample,:,:].T))/len(b)
        O[sample,:,:]=0.5*(O[sample,:,:]+sc.conj(O[sample,:,:].T))/len(b)
    if channel=='groundstate':
        return H
    fs=None
    refstate=sc.zeros(2*L*L)
    refstate[0::2]=1
    if wavefile==None:
        fs=GetFermiSigns(filename[0],refstate,channel=channel)
    else:
        fs=GetFermiSigns(wavefile,refstate,channel=channel)
    for s in range(sc.shape(H)[0]):
        H[s,:,:]=sc.dot(sc.diag(fs),sc.dot(H[s,:,:],sc.diag(fs)))
        O[s,:,:]=sc.dot(sc.diag(fs),sc.dot(O[s,:,:],sc.diag(fs)))
    ren=sc.ones(Nsamp)
    if gsfile!=None:
        ren=RenormalizeFactor(filename,gsfile,Nsamp=1,channel=channel,O=O,q=q)
    print('{0} pair of (H,O) matrices loaded, now diagonalize'.format(sc.shape(H)[0]))
    H=sc.einsum('ijk,i->ijk',H,ren)
    O=sc.einsum('ijk,i->ijk',O,ren)
    for s in range(sc.shape(H)[0]):
        E[s,:],V[s,:,:]=vln.geneigh(sc.squeeze(H[s,:,:]),sc.squeeze(O[s,:,:]))
    print('diagonalization finished')
    return H,O,E,V

def sqwtransamp(V,O,Lx,Ly,q,shift,phi,neel,r=sc.zeros((1,2)),rp=sc.zeros((1,2))):
    """
    Returns Sq[sample,r,rp,n]=<q,r|q,n><q,n|q,rp>
    """
    sqn=sc.zeros(sc.shape(V)[0:2],complex)
    kx,ky=fermisea(Lx,Ly,shift)
    pkrp=sc.zeros((sc.shape(V)[1],sc.shape(rp)[0]),complex)
    pkr=sc.zeros((sc.shape(V)[1],sc.shape(r)[0]),complex)
    pkrp[0:len(kx),:]=phiktrans(kx,ky,q[0],q[1],[phi,neel],rp)
    pkr[0:len(ky),:]=phiktrans(kx,ky,q[0],q[1],[phi,neel],r)
    OV=sc.einsum('ijk,ikl->ijl',O,V)
    rhs=sc.einsum('ijk,jl->ikl',sc.conj(OV),pkrp)
    lhs=sc.einsum('ij,kil->kjl',sc.conj(pkr),OV)
    sqn=sc.einsum('ijk,ikl->ijlk',lhs,rhs)
    return sqn

def SelfEnergy(GF_A,ChiGamma_A):
	''' calculating the dynamical self-energy from the Schwinger-Dyson equation '''
	N = int((len(En_A)-1)/2)
	## zero-padding the arrays
	exFD_A = sp.concatenate([FD_A[N:],sp.zeros(2*N+3),FD_A[:N]])
	exBE_A = sp.concatenate([BE_A[N:],sp.zeros(2*N+3),BE_A[:N]])
	ImGF_A = sp.concatenate([sp.imag(GF_A[N:]),sp.zeros(2*N+3),sp.imag(GF_A[:N])])
	ImCG_A = sp.concatenate([sp.imag(ChiGamma_A[N:]),sp.zeros(2*N+3),sp.imag(ChiGamma_A[:N])])
	## performing the convolution
	ftImSE1_A = -sp.conj(fft(exBE_A*ImCG_A))*fft(ImGF_A)*dE
	ftImSE2_A = -fft(exFD_A*ImGF_A)*sp.conj(fft(ImCG_A))*dE
	ImSE_A = sp.real(ifft(ftImSE1_A+ftImSE2_A))/sp.pi
	ImSE_A = sp.concatenate([ImSE_A[3*N+4:],ImSE_A[:N+1]])
	Sigma_A = KramersKronigFFT(ImSE_A) + 1.0j*ImSE_A
	return Sigma_A

def bondOrientation2sh(atoms,basis,l,neighbs=None,rcut=None,debug=False):
    atoms = array(atoms)
    basis = array(basis)    
    atoms = rectify(atoms,basis)

    if neighbs==None:
        bounds=[[0,basis[0][0]],[0,basis[1][1]],[0,basis[2][2]]]

        if rcut==None:
            rcut = generateRCut(atoms,basis,debug=debug)
            #print "Automatically generating r-cutoff=",rcut

        neighbs = secondShell( neighbors(atoms,bounds,rcut) )

    #sum the spherical harmonic over ever neighbor pair
    a = 4*np.pi / (2*l+1.)
    Ql=list()
    for i,ineighbs in enumerate(neighbs):
        n=len(ineighbs)

        shij = np.vectorize(complex)(zeros(2*l+1)) #spherical harmonic for bond i-j
        for j in ineighbs:
            shij += pairSphereHarms(atoms[i],minImageAtom(atoms[i],atoms[j],basis),l)/n
        shi = a * sum( scipy.real( scipy.multiply(shij,scipy.conj(shij)) ) )
        Ql.append(shi**0.5)
    
    return Ql,rcut

    def project_to_canonical(self, potential):
        r"""
        Project the Hagedorn wavepacket into the canonical basis.

        :param potential: The potential :math:`V` whose eigenvectors :math:`nu_l` are used for the transformation.

        .. note:: This function is expensive and destructive! It modifies the coefficients of the ``self`` instance.
        """
        # No projection for potentials with a single energy level.
        # The canonical and eigenbasis are identical here.
        if potential.get_number_components() == 1:
            return

        potential.calculate_eigenvectors()

        # Basically an ugly hack to overcome some shortcomings of the matrix function
        # and of the data layout.
        def f(q, x, component):
            x = x.reshape((self.quadrature.get_qr().get_number_nodes(),))
            z = potential.evaluate_eigenvectors_at(x)
            (row, col) = component
            return z[col][row,:]

        F = transpose(conj(self.quadrature.build_matrix(self,self,f)))
        c = self.get_coefficient_vector()
        d = dot(F, c)
        self.set_coefficient_vector(d)

    def evaluate_basis_at(self, nodes, component, prefactor=False):
        r"""
        Evaluate the Hagedorn functions :math:`\phi_k` recursively at the given nodes :math:`\gamma`.

        :param nodes: The nodes :math:`\gamma` at which the Hagedorn functions are evaluated.
        :param component: The index :math:`i` of the component whose basis functions :math:`\phi^i_k` we want to evaluate.
        :param prefactor: Whether to include a factor of :math:`\left(\det\left(Q_i\right)\right)^{-\frac{1}{2}}`.
        :return: Returns a twodimensional array :math:`H` where the entry :math:`H[k,i]` is the value
                 of the :math:`k`-th Hagedorn function evaluated at the node :math:`i`.
        """
        H = zeros((self.basis_size[component], nodes.size), dtype=complexfloating)

        (P, Q, S, p, q) = self.parameters[component]
        Qinv = Q**(-1.0)
        Qbar = conj(Q)
        nodes = nodes.reshape((1,nodes.size))

        H[0] = pi**(-0.25)*self.eps**(-0.5) * exp(1.0j/self.eps**2 * (0.5*P*Qinv*(nodes-q)**2 + p*(nodes-q)))
        H[1] = Qinv*sqrt(2.0/self.eps**2) * (nodes-q) * H[0]

        for k in xrange(2, self.basis_size[component]):
            H[k] = Qinv*sqrt(2.0/self.eps**2)*1.0/sqrt(k) * (nodes-q) * H[k-1] - Qinv*Qbar*sqrt((k-1.0)/k) * H[k-2]

        if prefactor is True:
            sqrtQ, self._cont_sqrt_cache[component] = cont_sqrt(Q, reference=self._cont_sqrt_cache[component])
            H = 1.0/sqrtQ*H

        return H

def sc_expand(cmodes,smodes,N):
    """
    Helper function: combines roles of c_expand, s_expand
    """
    from numpy import sqrt,zeros
    from scipy import conj
    Neven = (N+1)%2
    n = N/2
    factor = N - Neven
    
    # 'positive' modes
    pmodes = cmodes.copy()
    pmodes[0] *= sqrt(2)
    pmodes[1:] += smodes
    # Turn smodes into 'negative' modes
    smodes = (cmodes[1:] - smodes)[::-1]

    # Put it all together
    modes = zeros(2*n+1,dtype='complex128')
    modes[:n],modes[n:] = smodes,pmodes
    if bool(Neven):
        modes[0] += conj(modes[-1])
        return modes[:-1]
    else:
        return modes

def sc_collapse(modes,N):
    """
    Helper function: combines roles of c_collapse, s_collapse
    """
    from numpy import sqrt, append
    from scipy import conj

    NEven = (N%2)==0
    Nq = N + (N+1)%2
    modescopy = modes.copy()
    if NEven:
        modescopy[0] /= 2.
        modescopy = append(modescopy,conj(modescopy[0]))

    tempN = Nq/2
    tempN2 = tempN - ((Nq+1)%2)

    #cmodes = modes[:tempN+1][::-1].copy()
    cmodes = modescopy[:tempN+1][::-1]
    smodes = -cmodes[1:].copy()
    
    cmodes[1:tempN2+1] += modescopy[tempN+1:]
    smodes[:tempN2] += modescopy[tempN+1:]
    cmodes[0] *= sqrt(2)

    return [cmodes*1/2.,smodes*1/2.]

    def evaluate_basis_at(self, nodes, component=None, prefactor=False):
        r"""
        Evaluate the Hagedorn functions :math:`\phi_k` recursively at the given nodes :math:`\gamma`.

        :param nodes: The nodes :math:`\gamma` at which the Hagedorn functions are evaluated.
        :param component: Takes the basis size :math:`K_i` of this component :math:`i` as upper bound for :math:`K`.
        :param prefactor: Whether to include a factor of :math:`\left(\det\left(Q\right)\right)^{-\frac{1}{2}}`.
        :return: Returns a twodimensional :math:`K` times #nodes array :math:`H` where the entry :math:`H[k,i]` is
                 the value of the :math:`k`-th Hagedorn function evaluated at the node :math:`i`.
        """
        if component is not None:
            basis_size = self.basis_size[component]
        else:
            # Evaluate up to maximal :math:`K_i` and slice later if necessary
            basis_size = max(self.basis_size)

        H = zeros((basis_size, nodes.size), dtype=complexfloating)

        Qinv = self.Q**(-1.0)
        Qbar = conj(self.Q)
        nodes = nodes.reshape((1,nodes.size))

        H[0] = pi**(-0.25)*self.eps**(-0.5) * exp(1.0j/self.eps**2 * (0.5*self.P*Qinv*(nodes-self.q)**2 + self.p*(nodes-self.q)))
        H[1] = Qinv*sqrt(2.0/self.eps**2) * (nodes-self.q) * H[0]

        for k in xrange(2, basis_size):
            H[k] = Qinv*sqrt(2.0/self.eps**2)*1.0/sqrt(k) * (nodes-self.q) * H[k-1] - Qinv*Qbar*sqrt((k-1.0)/k) * H[k-2]

        if prefactor is True:
            sqrtQ, self._cont_sqrt_cache = cont_sqrt(self.Q, reference=self._cont_sqrt_cache)
            H = 1.0/sqrtQ*H

        return H

def ffacorr(a):
 """Returns the autocorrelation of a. Expects raw data"""
 z=np.zeros(2*len(a))
 z[:len(a)]=a
 fft=sc.fft(z)
 out=sc.ifft(fft*sc.conj(fft))
 return (out[:len(out)/2])

def prob4(filename='saw.wav', new_rate = 11025, outfile='prob4.wav'):
    """Down-samples a given .wav file to a new rate and saves the resulting
    signal as another .wav file.
    
    Parameters
    ----------
    filename : string, optional
        The name of the .wav sound file to be down-sampled.
        Defaults to 'saw.wav'.
    new_rate : integer, optional
        The down-sampled rate. Defaults to 11025.
    outfile : string, optional
        The name of the new file. Defaults to prob4.wav.

    Returns
    -------
    None
    """
    old_rate, in_sig = wavfile.read(filename)
    fin = fftw.fft(sp.float32(in_sig))
    # Use if scipy_fftpack is unavailable
    # fin = sp.fft(sp.float32(in_sig))
    nsiz = sp.floor(in_sig.size * new_rate / old_rate)
    nsizh = sp.floor(nsiz / 2)
    fout = sp.zeros(nsiz) + 0j
    fout[0:nsizh] = fin[0:nsizh]
    fout[nsiz-nsizh+1:] = sp.conj(sp.flipud(fout[1:nsizh]))
    out = sp.real(sp.ifft(fout))
    out = sp.int16(out/sp.absolute(out).max() * 32767)
    plot_signal(filename)
    wavfile.write('prob4.wav',new_rate,out)
    print ""; plot_signal('prob4.wav')

def down_sample(filename, new_rate, outputfile=None):
    """
    Create a down-sampled copy of the provided .wav file.  Unless overridden, the output
        file will be of the form "down_<orginalname>.wav"
        
    Parameters
    ----------
    filename : string
        input .wav file
    new_rate : int
        sample rate of output file
    outputfile : string
        name of output file
    """

    if outputfile is None:
        outputfile = "down_" + filename

    old_rate, in_sig = wavfile.read(filename)
    in_sig = sp.float32(in_sig)
    fin = sp.fft(in_sig)
    nsiz = sp.floor(in_sig.size * new_rate / old_rate)
    nsizh = sp.floor(nsiz / 2)
    fout = sp.zeros(nsiz)
    fout = fout + 0j
    fout[0:nsizh] = fin[0:nsizh]
    fout[nsiz - nsizh + 1 :] = sp.conj(sp.flipud(fout[1:nsizh]))
    out = sp.ifft(fout)
    out = sp.real(out)  # Take the real component of the signal
    out = sp.int16(out / sp.absolute(out).max() * 32767)
    wavfile.write(outputfile, new_rate, out)

def phiktrans(kx,ky,qx,qy,p,r=sc.zeros((1,2))):
    """
    Returns phi[k,r] such that |q,r>=sum_k phi[k,r]|q,k>
    """
    kqx=kx-qx
    kqy=ky-qy
    pk=sc.zeros((sc.shape(kx)[0],sc.shape(r)[0]),complex)
    pke=sc.conj(uk(kx,ky,1,1,p))*uk(kqx,kqy,-1,-1,p)+\
        sc.conj(vk(kx,ky,1,1,p))*vk(kqx,kqy,-1,-1,p)
    pko=sc.conj(vk(kx,ky,1,1,p))*uk(kqx,kqy,-1,-1,p)+\
        sc.conj(uk(kx,ky,1,1,p))*vk(kqx,kqy,-1,-1,p)
    even=1-sc.mod(r[:,0]+r[:,1],2)
    odd=sc.mod(r[:,0]+r[:,1],2)
    ph=sc.exp(-2j*sc.pi*(sc.einsum('i,j->ij',kx,r[:,0])+sc.einsum('i,j->ij',ky,r[:,1])))
    pk=sc.einsum('ij,j,i->ij',ph,even,pke)+sc.einsum('ij,j,i->ij',ph,odd,pko)
    return pk

def coherent(N, alpha, method='operator'):
    """Generates a coherent state with eigenvalue alpha. 
    
    Constructed using displacement operator on vacuum state.
    
    Parameters
    ----------
    N : int 
        Number of Fock states in Hilbert space.    
    alpha : float/complex 
        Eigenvalue of coherent state.
    method : string {'operator', 'analytic'}
        Method for generating coherent state.
    
    Returns
    -------
    state : qobj
        Qobj quantum object for coherent state
    
    Examples
    --------        
    >>> coherent(5,0.25j)
    Quantum object: dims = [[5], [1]], shape = [5, 1], type = ket
    Qobj data = 
    [[  9.69233235e-01+0.j        ]
     [  0.00000000e+00+0.24230831j]
     [ -4.28344935e-02+0.j        ]
     [  0.00000000e+00-0.00618204j]
     [  7.80904967e-04+0.j        ]]
     
    Notes
    -----
    Select method 'operator' (default) or 'analytic'. With the 
    'operator' method, the coherent state is generated by displacing
    the vacuum state using the displacement operator defined in the
    truncated Hilbert space of size 'N'. This method guarantees that the 
    resulting state is normalized. With 'analytic' method the coherent state 
    is generated using the analytical formula for the coherent state 
    coefficients in the Fock basis. THIS METHOD DOES NOT GUARANTEE THAT THE 
    STATE IS NORMALIZED if truncated to a small number of Fock states, 
    but would in that case give more accurate coefficients.
         
    """
    if method == "operator":

        x=basis(N,0)
        a=destroy(N)
        D=(alpha*a.dag()-conj(alpha)*a).expm()
        return D*x

    elif method == "analytic":

        data = np.zeros([N,1],dtype=complex)
        n = arange(N)
        data[:,0] = np.exp(-(abs(alpha)**2)/2.0)*(alpha**(n))/_sqrt_factorial(n)
        return Qobj(data)

    else:
        raise TypeError("The method option can only take values 'operator' or 'analytic'")

def xcorr(t, x, y, zeropad = True):

    tau = t
    sx = len(x)
    sy = len(y)
    if zeropad == True:
        Xn = sp.fft(x, n = len(x)*2)
        Yn = sp.conj(sp.fft(y, n = len(x)*2))
    else:
        Xn = sp.fft(x)
        Yn = sp.conj(sp.fft(y))

    xcor = sp.real(fftpack.fftshift(sp.ifft(Xn*Yn)))
    dt = t[1]-t[0]
    
    tau = sp.linspace(-len(xcor)/2*dt-dt/2,len(xcor)/2*dt-dt/2,len(xcor))
    return tau, xcor

def Renorm(sqsq,O,Lx,Ly,q,shift,p):
    kx,ky=fermisea(Lx,Ly,shift)
    pp=phiktrans(kx,ky,float(q[0])/Lx,float(q[1])/Ly,[p['phi'],p['neel']])
    b=sc.dot(sc.conj(pp),sc.dot(O,pp))
    r=sqsq/b
    if sc.isnan(r):
        r=1
    return r,b

def XIntegralsFFT(GF_A,Bubble_A,Lambda,BubZero):
	''' calculate X integral to susceptibilities using FFT '''
	N = int((len(En_A)-1)/2)
	Kappa_A  = TwoParticleBubble(GF_A,GF_A**2,'eh')
	Bubble_A = TwoParticleBubble(GF_A,GF_A,'eh')
	#print(Kappa_A[N],Bubble_A[N])
	V_A   = 1.0/(1.0+Lambda*Bubble_A)
	KV_A  = Lambda*Kappa_A*V_A**2
	KmV_A = Lambda*sp.flipud(sp.conj(Kappa_A))*V_A**2
	## zero-padding the arrays
	exFD_A  = sp.concatenate([FD_A[N:],sp.zeros(2*N+2),FD_A[:N+1]])
	ImGF_A  = sp.concatenate([sp.imag(GF_A[N:]),sp.zeros(2*N+2),sp.imag(GF_A[:N+1])])
	ImGF2_A = sp.concatenate([sp.imag(GF_A[N:]**2),sp.zeros(2*N+2),sp.imag(GF_A[:N+1]**2)])
	ImV_A   = sp.concatenate([sp.imag(V_A[N:]),sp.zeros(2*N+2),sp.imag(V_A[:N+1])])
	ImKV_A  = sp.concatenate([sp.imag(KV_A[N:]),sp.zeros(2*N+2),sp.imag(KV_A[:N+1])])
	ImKmV_A = sp.concatenate([sp.imag(KmV_A[N:]),sp.zeros(2*N+2),sp.imag(KmV_A[:N+1])])
	## performing the convolution
	ftImX11_A = -sp.conj(fft(exFD_A*ImV_A))*fft(ImGF2_A)*dE
	ftImX12_A =  fft(exFD_A*ImGF2_A)*sp.conj(fft(ImV_A))*dE
	ftImX21_A = -sp.conj(fft(exFD_A*ImKV_A))*fft(ImGF_A)*dE
	ftImX22_A =  fft(exFD_A*ImGF_A)*sp.conj(fft(ImKV_A))*dE
	ftImX31_A = -sp.conj(fft(exFD_A*ImKmV_A))*fft(ImGF_A)*dE
	ftImX32_A =  fft(exFD_A*ImGF_A)*sp.conj(fft(ImKmV_A))*dE
	## inverse transform
	ImX1_A =  sp.real(ifft(ftImX11_A+ftImX12_A))/sp.pi
	ImX2_A =  sp.real(ifft(ftImX21_A+ftImX22_A))/sp.pi
	ImX3_A = -sp.real(ifft(ftImX31_A+ftImX32_A))/sp.pi
	ImX1_A =  sp.concatenate([ImX1_A[3*N+4:],ImX1_A[:N+1]])
	ImX2_A =  sp.concatenate([ImX2_A[3*N+4:],ImX2_A[:N+1]])
	ImX3_A =  sp.concatenate([ImX3_A[3*N+4:],ImX3_A[:N+1]])
	## getting real part from imaginary
	X1_A = KramersKronigFFT(ImX1_A) + 1.0j*ImX1_A + BubZero # constant part !!!
	X2_A = KramersKronigFFT(ImX2_A) + 1.0j*ImX2_A
	X3_A = KramersKronigFFT(ImX3_A) + 1.0j*ImX3_A
	return [X1_A,X2_A,X3_A]

def conj(x):
    """
    Wrapper for conjugate on a CXData object.
    """
    if isinstance(x, CXData):
        l=[]
        for i in xrange(len(x)):
            l.append(sp.conj(x.data[i]))
        return CXData(data=l)
    elif isinstance(x, CXModal):
        l=[]
        for mode in range(len(x.modes)):
            l.append(conj(x.modes[mode]))
        return CXModal(modes=l)
    elif isinstance(x, np.ndarray):
        return sp.conj(x)
    else:
        raise Exception('Unknown data type passed to conj')

def cross_periodogram(x1, x2, Fs, window, scale_by_freq = True, two_sided = False, detrend = True):
    # Remove mean
    if detrend:
        t1 = s.fftpack.fft((x1-s.mean(x1))*window)
        t2 = s.fftpack.fft((x2-s.mean(x2))*window)
    else:
        t1 = s.fftpack.fft(x1*window)
        t2 = s.fftpack.fft(x2*window)
        
    # Return cross spectral density or just spectrum    
    if scale_by_freq:
        csd = t1*s.conj(t2)/s.sum(window**2)/Fs
    else:
        csd = t1*s.conj(t2)/s.sum(window**2)
    if two_sided:
        return csd
    else:
        return csd[:len(x1)/2+1]