def problem4():
	# read in tada.wav
	rate, tada = wavfile.read('tada.wav')
	
	# upon inspection, we find that tada.wav is a stereo audio file. 
	# we create stereo white noise that lasts 10 seconds
	L_white = sp.int16(sp.random.randint(-32767,32767,rate*10))
	R_white = sp.int16(sp.random.randint(-32767,32767,rate*10))
	white = sp.zeros((len(L_white),2))
	white[:,0] = L_white
	white[:,1] = R_white
	
	# pad tada signal with zeros
	padded_tada = sp.zeros_like(white)
	padded_tada[:len(tada)] = tada
	ptada = padded_tada
	
	# fourier transforms
	ftada = sp.fft(ptada,axis=0)
	fwhite = sp.fft(white,axis=0)
	
	# inverse transform of convolution
	out = sp.ifft((ftada*fwhite),axis=0)
	
	# prepping output and writing file
	out = sp.real(out)
	scaled = sp.int16(out / sp.absolute(out).max() * 32767)
	wavfile.write('my_tada_conv.wav',rate,scaled)

def ipfft(pft, xs, ys):
	if(xs > 2 * len(pft) + 1 and ys > 2 * len(pft) + 1):
		return fftshift(real(ifft2(fftshift(polar2rect(pft, xs, ys)))))
	t = fftshift(real(ifft2(fftshift(polar2rect(pft, 2*len(pft)+1, 2*len(pft)+1)))))
	tx = len(pft) - xs / 2
	ty = len(pft) - ys / 2
	return t[tx:(tx+xs),ty:(ty+ys)]

def rlsloo_ll1( V, D, Y, lambd):
	"""
	Computes cs and the actual LOO errors for a single value of lambda. (lambd)
	"""
        n = V.shape[0]
        cl = Y.shape[1]
        
        inner  = 1/(D + lambd)
	inner = inner.conj()
        VtY = sp.dot(V.T, Y)
	VtY = VtY.conj()

        # Because of signs of D are flipped (scipy.linalg.eig returns
        # flipped signs for complex part of the eigenvalues)
        in_dot = sp.ones((n,1)) * inner
        ViD = V * in_dot
        cs = sp.dot(ViD, VtY)
        dGi = sp.sum(ViD*V, axis = 1)
        # -- till here works fine
        #check matrix dimensions
        looerrs = cs.ravel()/sp.real(dGi.ravel())
	looerrs = sp.real(looerrs)
        cs = sp.real(cs.transpose())

        return cs.ravel(), looerrs

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 add_shape(self, f):
		#Create shape
		S = Shape(f, self.R, self.SHAPE_R)
	
		#Add to shape list
		S.shape_num = len(self.shape_list)
		self.shape_list.append(S)
		
		row = []
		for k in range(len(self.shape_list)):
			T = self.shape_list[k]
			ift = real(ipfft(pft_mult(pft_rotate(S.pft, 2.*pi/6.), T.pft), 2*self.SHAPE_R+1,2*self.SHAPE_R+1))
			Spad = imrotate(cpad(S.indicator, array([2*self.SHAPE_R+1,2*self.SHAPE_R+1])), 360./6.)
			Tpad = cpad(T.indicator, array([2*self.SHAPE_R+1,2*self.SHAPE_R+1]))
			pind = real(fftconvolve(Spad, Tpad, mode='same'))
			imshow(pind)
			imshow(ift)
			obst = to_ind(pind, 0.001)
			imshow(obst)
			cutoff = best_cutoff(ift, obst, S.radius + T.radius)
			print cutoff
			imshow(to_ind(ift, cutoff))
			row.append(cutoff * self.tarea)
		self.cutoff_matrix.append(row)

		return S

def root_locus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True, 
               PrintGain=True):
    """Calculate the root locus by finding the roots of 1+k*TF(s)
    where TF is self.num(s)/self.den(s) and each k is an element
    of kvect.

    Parameters
    ----------
    sys : linsys
        Linear input/output systems (SISO only, for now)
    kvect : gain_range (default = None)
        List of gains to use in computing diagram
    Plot : boolean (default = True)
        If True, plot magnitude and phase
    PrintGain: boolean (default = True)
        If True, report mouse clicks when close to the root-locus branches,
        calculate gain, damping and print
    Return values
    -------------
    rlist : list of computed root locations
    """

    # Convert numerator and denominator to polynomials if they aren't
    (nump, denp) = _systopoly1d(sys);

    # Compute out the loci
    mymat = _RLFindRoots(sys, kvect)
    mymat = _RLSortRoots(sys, mymat)

    # Create the plot
    if (Plot):
        f = pylab.figure()
        if PrintGain:
            cid = f.canvas.mpl_connect(
                'button_release_event', partial(_RLFeedbackClicks, sys=sys))
        ax = pylab.axes();

        # plot open loop poles
        poles = array(denp.r)
        ax.plot(real(poles), imag(poles), 'x')

        # plot open loop zeros
        zeros = array(nump.r)
        if zeros.any():
            ax.plot(real(zeros), imag(zeros), 'o')

        # Now plot the loci
        for col in mymat.T:
            ax.plot(real(col), imag(col), plotstr)

        # Set up plot axes and labels
        if xlim:
            ax.set_xlim(xlim)
        if ylim:
            ax.set_ylim(ylim)
        ax.set_xlabel('Real')
        ax.set_ylabel('Imaginary')

    return mymat

def confMap(shape,mapfunc):
	shapemapped = [None]*len(shape)
	for i in range(0,len(shape)):
		shapemapped[i] = mapfunc(shape[i])
	
	plt.scatter(sp.real(shape),sp.imag(shape),color='r')
	plt.scatter(sp.real(shapemapped),sp.imag(shapemapped),color='b')
	plt . show ()

def pltFunction (cavity1,cavity2,cavity3,plotType):
  if (plotType==0):
    return sp.absolute(cavity1[:])**2,sp.absolute(cavity2[:])**2,sp.absolute(cavity3[:])**2
  elif (plotType==1):
    return sp.real(cavity1[:]),sp.real(cavity2[:]),sp.real(cavity3[:])
  elif (plotType==2):
    return sp.imag(cavity1[:]),sp.imag(cavity2[:]),sp.imag(cavity3[:])
  else: 
    return cavity1, cavity2, cavity3

    def monitor_limitedFrames_MP(self, i):
        if self.centroid==True:
            fit_func = correlateFrames.cent
        else:
            fit_func = correlateFrames.fitCorrGaussian_v2
        corr_func = correlateFrames.corr
        
        xdim_fft, ydim_fft = self.fft_target.shape
        #i = 0

        #print "Dry Run? ", self.dry_run
        
        if (i<self.frames) and (not self.quit):
            self.active=True
            time.sleep(self.monitorDelay)
            self.times.append(time.clock())
            if self.fr_change>0 and i>(self.unactive_until + self.initialFrames-1):
                ind = int(scipy.floor(i/self.fr_change))
                self.conversion = self.conversions[ind]

            self.stage_zs.append(self.piezo.getPosition(3)-self.zoffset)
            target_gs_ht, defocus_ht_diff, target_xval, target_yval = self.getImage_doCorrs_inThread()
            
            self.target_signal.append(target_gs_ht)
            self.defocus_signal.append(defocus_ht_diff)

            latestSig = defocus_ht_diff / target_gs_ht

            self.xdrift.append(scipy.real(target_xval) - self.gaussFitRegion)
            self.ydrift.append(scipy.real(target_yval) - self.gaussFitRegion)


            if True:
                if i>(self.unactive_until + self.initialFrames-1):
                    if i==self.unactive_until + self.initialFrames:
                        self.initial_ds = scipy.mean(self.defocus_signal[self.unactive_until:])
                        self.initial_z = scipy.mean(self.stage_zs)
                        self.initial_ts = scipy.mean(self.target_signal[self.unactive_until:])
                        self.sig0 = self.initial_ds/self.initial_ts
                        print self.initial_z
                    if i>(self.unactive_until+self.initialFrames):
                        if self.use_multiplane and self.newplane>0:
                            mp = int(i)/int(self.newplane)
                        else:
                            mp = 0
                        if not self.dry_run:
                            self.react_z(self.initial_ds, self.initial_ts, self.initial_z, toprint=False, multiplane=mp)
                            self.react_xy(rolling_av=self.rollingAvXY, toprint=False)
                    initialSigFound = True
                else:
                    initialSigFound = False

            if initialSigFound:
                return self.stage_zs[-1], self.xdrift[-1], self.ydrift[-1], self.sig0, latestSig
            else:
                return self.stage_zs[-1], self.xdrift[-1], self.ydrift[-1], -1, latestSig
        return 0,0,0,0,0

def rceps(x): 
	y = sp.real(ifft(sp.log(sp.absolute(fft(x)))))
	n = len(x) 
	if (n%2) == 1:
		ym = np.hstack((y[0], 2*y[1:n/2], np.zeros(n/2-1)))
	else:
		ym = np.hstack((y[0], 2*y[1:n/2], y[n/2+1], np.zeros(n/2-1)))
	ym = sp.real(ifft(sp.exp(fft(ym)))) 
	return (y, ym)

def sphericalPot(x,y,shift=0,radius=1,scale=1):
    from scipy import real,sqrt
    size_x = x.max()
    size_y = y.max()
    Vbottom = 0
    x = x-size_x/2
    left_sphere = real(sqrt(radius**2-(x**2+(y-shift)**2)))*heaviside(y-shift)+real(sqrt(radius**2-x**2))*heaviside(-(y-shift))
    right_sphere = real(sqrt(radius**2-(x**2+(y-size_y+shift)**2)))*heaviside(-(y-size_y+shift))+real(sqrt(radius**2-x**2))*heaviside((y-size_y+shift))
    V = Vbottom +scale*(left_sphere+right_sphere)
    return V

def RenormalizeFactor(excfile,gsfile,channel=None,Nsamp=1,O=None,q=None):
    if not type(excfile)==list:
        excfile=[excfile]
    if not type(gsfile)==list:
        gsfile=[gsfile]
    exat=GetAttr(excfile[0])
    gsat=GetAttr(gsfile[0])
    L=exat['L']
    if q==None:
        q=sc.array([exat['qx'],exat['qy']])
    if 'phasex' in exat.keys():
        shift=sc.array([exat['phasex']/2.0,exat['phasey']/2.0])
    else:
        shift=sc.array([exat['phase_shift_x']/2.0,exat['phase_shift_y']/2.0])
    phi=exat['phi']
    neel=exat['neel']
    qx,qy,Sq=GetSq(gsfile)
    kx,ky=sf.fermisea(L,L,shift)
    qidx=ml.find((qx==q[0])*(qy==q[1]))
    if O==None:
        _,O,_,_=GetEigSys(excfile,Nsamp)
    pk=None
    sqq=None
    if channel==None:
        channel=exat['channel']
    if channel=='trans':
        pk=sc.squeeze(sf.phiktrans(kx,ky,q[0]/L,q[1]/L,[phi,neel]))
        sqq=sc.real(0.5*(Sq[0,1,qidx]+Sq[0,2,qidx]))
    elif channel=='long':
        pkup=sc.squeeze(sf.phiklong(kx,ky,q[0]/L,q[1]/L,1,[phi,neel]))
        pkdo=sc.squeeze(sf.phiklong(kx,ky,q[0]/L,q[1]/L,-1,[phi,neel]))
        if (q[0]/L==0.5 and q[1]/L==0.5) or (q[0]/L==0 and q[1]/L==0):
            pk=sc.zeros(2*sc.shape(pkup)[0]+1,complex)
        else:
            pk=sc.zeros(2*sc.shape(pkup)[0],complex)
        pk[0:2*sc.shape(pkup)[0]:2]=pkup
        pk[1:2*sc.shape(pkdo)[0]:2]=pkdo
        if (qx[0]/L==0.5 and q[1]/L==0.5) or (q[0]/L==0 and q[1]/L==0):
            if neel==0:
                pk[-1]=0
            else:
                pk[-1]=sum(neel/sf.omega(kx,ky,[phi,neel]))
        sqq=Sq[0,0,qidx]
    else:
        raise(InputFileError('In file \''+excfile+'\', channel=\''+str(channel)+'\'. Should be \'trans\' or \'long\''))
    sqe=sc.einsum('i,jik,k->j',sc.conj(pk),O,pk)
    out=sc.zeros(Nsamp)
    for n in range(Nsamp):
        if abs(sqq)<1e-6 or abs(sqe[n])<1e-6:
            warnings.warn('Probably ill-defined renormalization, returns 1 for sample {0} out of {1}'.format(n,Nsamp),UserWarning)
            out[n]=1
        else:
            out[n]=sc.real(sqq/sqe[n])
    return out

	def grad(self, A, B, q):
		#Compute relative transformation
		pr = norm(q.x)
		
		if(pr >= A.radius + B.radius ):
			return array([0., 0., 0., 0.])
		
		#Load shape parameters
		fa  = A.pft
		fb  = B.pft
		da  = A.pdft
		db  = B.pdft
		ea 	= A.energy
		eb 	= B.energy
		
		#Estimate cutoff threshold
		cutoff = self.__get_cutoff(A, B)
		
		#Compute coordinate coefficients
		m   = 2.j * pi / (2. * self.SHAPE_R + 1) * pr
		phi = atan2(q.x[1], q.x[0])
		
		#Set up initial sums
		s_0	 = real(fa[0][0] * fb[0][0])
		s_x  = 0.
		s_y  = 0.
		s_ta = 0.
		s_tb = 0.
		
		for r in range(1, self.R):
		
			#Compute theta terms
			dtheta = 2. * pi / len(fa[r])
			theta  = arange(len(fa[r])) * dtheta
		
			#Construct multiplier / v
			mult = exp((m * r) * cos(theta + phi)) * r * dtheta
			u 	 = pds.shift(conjugate(fb[r]), q.theta) * mult
			v 	 = fa[r] * u
			
			#Check for early out
			s_0  += sum(real(v))
			if(s_0 + min(ea[r], eb[r]) <= cutoff):
				return array([0.,0.,0.,0.])
				
			#Sum up gradient vectors
			v     = real(1.j * v)
			s_x  -= sum(v * sin(theta + phi) )
			s_y  -= sum(v * cos(theta + phi) )
			s_t  += sum(real(da[r] * u))
		
		if(s_0 <= cutoff):
			return array([0., 0., 0., 0.])
		return array([s_x, s_y, s_ta, s_tb, s_0])

def predictor_step(phi,args,direction):
    assert abs(direction)==1
    
    Q = field_to_quasidensity(phi,args.L)
    Lambda = sp.exp(-args.t)*args.Delta + diags(Q,0)
    rho = sp.real(spsolve(Lambda, args.R-Q))
    
    delta_t = direction*args.epsilon/sp.real(sp.sqrt(sp.sum(rho*Q*rho)))
    delta_phi = phi + delta_t*rho
    
    return [delta_phi, delta_t]

def wigner(psi,xvec,yvec,g=sqrt(2)):
    """Wigner function for a state vector or density matrix 
    at points xvec+i*yvec.
    
    Parameters
    ----------
    state : qobj 
        A state vector or density matrix.
    
    xvec : array_like
        x-coordinates at which to calculate the Wigner function.
    
    yvec : array_like
        y-coordinates at which to calculate the Wigner function.
        
    g : float
        Scaling factor for a = 0.5*g*(x+iy), default g=sqrt(2).
    
    Returns
    --------
    W : array
        Values representing the Wigner function calculated over the specified range [xvec,yvec].
    
    
    """
    if psi.type=='ket' or psi.type=='oper':
        M=prod(psi.shape[0])
    elif psi.type=='bra':
        M=prod(psi.shape[1])
    else:
        raise TypeError('Input state is not a valid operator.')
    X,Y = meshgrid(xvec, yvec)
    amat = 0.5*g*(X + 1.0j*Y)
    wmat=zeros(shape(amat))
    Wlist=array([zeros(shape(amat),dtype=complex) for k in range(M)])
    Wlist[0]=exp(-2.0*abs(amat)**2)/pi
    if psi.type=='ket' or psi.type=='bra':
        psi=ket2dm(psi)
    wmat=real(psi[0,0])*real(Wlist[0])
    for n in range(1,M):
        Wlist[n]=(2.0*amat*Wlist[n-1])/sqrt(n)
        wmat+= 2.0*real(psi[0,n]*Wlist[n])
    for m in range(M-1):
        temp=copy(Wlist[m+1])
        Wlist[m+1]=(2.0*conj(amat)*temp-sqrt(m+1)*Wlist[m])/sqrt(m+1)
        for n in range(m+1,M-1):
            temp2=(2.0*amat*Wlist[n]-sqrt(m+1)*temp)/sqrt(n+1)
            temp=copy(Wlist[n+1])
            Wlist[n+1]=temp2
        wmat+=real(psi[m+1,m+1]*Wlist[m+1])
        for k in range(m+2,M):
            wmat+=2.0*real(psi[m+1,k]*Wlist[k])
    return 0.5*wmat*g**2

def sliceFieldArr(Ein,NX,NY,NZ,SliceX,SliceY,SliceZ,figNum,fieldName):
    """ Visualizing field with slices 
    """
    vxy1=Ein[:,:,SliceZ,0]
    vxy2=Ein[:,:,SliceZ,1]
    vxy3=Ein[:,:,SliceZ,2]
    vyz1=Ein[SliceX,:,:,0]
    vyz2=Ein[SliceX,:,:,1]
    vyz3=Ein[SliceX,:,:,2]
    vxz1=Ein[:,SliceY,:,0]
    vxz2=Ein[:,SliceY,:,1]
    vxz3=Ein[:,SliceY,:,2]
    plt.ion()
    fig=plt.figure(figNum)
    fig.clear()
    plt.subplot(331)
    plt.imshow(sci.real(vxy1))
    plt.title('real($'+fieldName+'_{x}$), xy-plane')
    plt.colorbar()
    plt.subplot(332)
    plt.imshow(sci.real(sci.swapaxes(vxz1,0,1)))
    plt.title('real($'+fieldName+'_{x}$), xz-plane')
    plt.colorbar()
    plt.subplot(333)
    plt.imshow(sci.real(sci.swapaxes(vyz1,0,1)))
    plt.title('real($'+fieldName+'_{x}$), yz-plane')
    plt.colorbar()
    plt.subplot(334)
    plt.imshow(sci.real(vxy2))
    plt.title('real($'+fieldName+'_{y}$), xy-plane')
    plt.colorbar()
    plt.subplot(335)
    plt.imshow(sci.real(sci.swapaxes(vxz2,0,1)))
    plt.title('real($'+fieldName+'_{y}$), xz-plane')
    plt.colorbar()
    plt.subplot(336)
    plt.imshow(sci.real(sci.swapaxes(vyz2,0,1)))
    plt.title('real($'+fieldName+'_{y}$), yz-plane')
    plt.colorbar()
    plt.subplot(337)
    plt.imshow(sci.real(vxy3))
    plt.title('real($'+fieldName+'_{z}$), xy-plane')
    plt.colorbar()
    plt.subplot(338)
    plt.imshow(sci.real(sci.swapaxes(vxz3,0,1)))
    plt.title('real($'+fieldName+'_{z}$), xz-plane')
    plt.colorbar()
    plt.subplot(339)
    plt.imshow(sci.real(sci.swapaxes(vyz3,0,1)))
    plt.title('real($'+fieldName+'_{z}$), yz-plane')
    plt.colorbar()
    
    fig.canvas.draw()
    plt.ioff()
    
    return

def _wigner_laguerre(rho, xvec, yvec, g, parallel):
    """
    Using Laguerre polynomials from scipy to evaluate the Wigner function for
    the density matrices :math:`|m><n|`, :math:`W_{mn}`. The total Wigner
    function is calculated as :math:`W = \sum_{mn} \\rho_{mn} W_{mn}`.
    """

    M = np.prod(rho.shape[0])
    X, Y = meshgrid(xvec, yvec)
    A = 0.5 * g * (X + 1.0j * Y)
    W = zeros(np.shape(A))

    # compute wigner functions for density matrices |m><n| and
    # weight by all the elements in the density matrix
    B = 4 * abs(A) ** 2
    if sp.isspmatrix_csr(rho.data):
        # for compress sparse row matrices
        if parallel:
            iterator = (
                (m, rho, A, B) for m in range(len(rho.data.indptr) - 1))
            W1_out = parfor(_par_wig_eval, iterator)
            W += sum(W1_out)
        else:
            for m in range(len(rho.data.indptr) - 1):
                for jj in range(rho.data.indptr[m], rho.data.indptr[m + 1]):
                    n = rho.data.indices[jj]

                    if m == n:
                        W += real(rho[m, m] * (-1) ** m * genlaguerre(m, 0)(B))

                    elif n > m:
                        W += 2.0 * real(rho[m, n] * (-1) ** m *
                                        (2 * A) ** (n - m) *
                                        sqrt(factorial(m) / factorial(n)) *
                                        genlaguerre(m, n - m)(B))
    else:
        # for dense density matrices
        B = 4 * abs(A) ** 2
        for m in range(M):
            if abs(rho[m, m]) > 0.0:
                W += real(rho[m, m] * (-1) ** m * genlaguerre(m, 0)(B))
            for n in range(m + 1, M):
                if abs(rho[m, n]) > 0.0:
                    W += 2.0 * real(rho[m, n] * (-1) ** m *
                                    (2 * A) ** (n - m) *
                                    sqrt(factorial(m) / factorial(n)) *
                                    genlaguerre(m, n - m)(B))

    return 0.5 * W * g ** 2 * np.exp(-B / 2) / pi

def KramersKronigFFT(ImX_A):
	'''	Hilbert transform used to calculate real part of a function from its imaginary part
	uses piecewise cubic interpolated integral kernel of the Hilbert transform
	use only if len(ImX_A)=2**m-1, uses fft from scipy.fftpack  '''
	X_A = sp.copy(ImX_A)
	N = int(len(X_A))
	## be careful with the data type, orherwise it fails for large N
	if N > 3e6: A = sp.arange(3,N+1,dtype='float64')
	else:       A = sp.arange(3,N+1)  
	X1 = 4.0*sp.log(1.5)
	X2 = 10.0*sp.log(4.0/3.0)-6.0*sp.log(1.5)
	## filling the kernel
	if N > 3e6: Kernel_A = sp.zeros(N-2,dtype='float64')
	else:       Kernel_A = sp.zeros(N-2)
	Kernel_A = (1-A**2)*((A-2)*sp.arctanh(1.0/(1-2*A))+(A+2)*sp.arctanh(1.0/(1+2*A)))\
	+((A**3-6*A**2+11*A-6)*sp.arctanh(1.0/(3-2*A))+(A+3)*(A**2+3*A+2)*sp.arctanh(1.0/(2*A+3)))/3.0
	Kernel_A = sp.concatenate([-sp.flipud(Kernel_A),sp.array([-X2,-X1,0.0,X1,X2]),Kernel_A])/sp.pi
	## zero-padding the functions for fft
	ImXExt_A = sp.concatenate([X_A[int((N-1)/2):],sp.zeros(N+2),X_A[:int((N-1)/2)]])
	KernelExt_A = sp.concatenate([Kernel_A[N:],sp.zeros(1),Kernel_A[:N]])
	## performing the fft
	ftReXExt_A = -fft(ImXExt_A)*fft(KernelExt_A)
	ReXExt_A = sp.real(ifft(ftReXExt_A))
	ReX_A = sp.concatenate([ReXExt_A[int((3*N+3)/2+1):],ReXExt_A[:int((N-1)/2+1)]])
	return ReX_A

 def __init__(self, output='out', input='in', \
              mag=None, phase=None, coh=None, \
              freqlim=[], maglim=[], phaselim=[], \
              averaged='not specified', \
              seedfreq=-1, seedphase=0,
              labels=[], legloc=-1, compin=[]):
     self.output = output
     self.input = input
     if len(compin) > 0:
        if mag is None:
           self.mag = squeeze(colwise(abs(compin)))
        if phase is None:
           self.phase = squeeze(colwise(arctan2(imag(compin),real(compin))*180.0/pi))
     else:
         self.mag = squeeze(mag)
         self.phase = squeeze(phase)
     self.coh = coh
     self.averaged = averaged
     self.seedfreq = seedfreq
     self.seedphase = seedphase
     self.freqlim = freqlim
     self.maglim = maglim
     self.phaselim = phaselim
     self.labels = labels
     self.legloc = legloc

def gabor2d(gsw, gsh, gx0, gy0, wfreq, worient, wphase, shape):
    """ Generate a gabor 2d array
    
    Inputs:
      gsw -- standard deviation of the gaussian envelope (width)
      gsh -- standard deviation of the gaussian envelope (height)
      gx0 -- x indice of center of the gaussian envelope
      gy0 -- y indice of center of the gaussian envelope
      wfreq -- frequency of the 2d wave
      worient -- orientation of the 2d wave
      wphase -- phase of the 2d wave
      shape -- shape tuple (height, width)

    Outputs:
      gabor -- 2d gabor with zero-mean and unit-variance

    """
    
    height, width = shape
    y, x = N.mgrid[0:height, 0:width]
    
    X = x * N.cos(worient) * wfreq
    Y = y * N.sin(worient) * wfreq
	
    env = N.exp( -.5 * ( ((x-gx0)**2./gsw**2.) + ((y-gy0)**2./gsh**2.) ) )
    wave = N.exp( 1j*(2*N.pi*(X+Y) + wphase) )
    gabor = N.real(env * wave)
    
    gabor -= gabor.mean()
    gabor /= fastnorm(gabor)
    
    return gabor