def invert(self, estimate):
     """Invert the estimate to produce slopes.
     
     Parameters
     ----------
     estimate : array_like
         Phase estimate to invert.
     
     Returns
     -------
     xs : array_like
         Estimate of the x slopes.
     ys : array_like
         Estimate of the y slopes.
     
     
     """
     if self.manage_tt:
         estimate, ttx, tty = remove_tiptilt(self.ap, estimate)
     
     est_ft = fftpack.fftn(estimate) / 2.0
     
     xs_ft = self.gx * est_ft
     ys_ft = self.gy * est_ft
     
     xs = np.real(fftpack.ifftn(xs_ft))
     ys = np.real(fftpack.ifftn(ys_ft))
     
     if self.manage_tt and not self.suppress_tt:
         xs += ttx
         ys += tty
     
     return (xs, ys)

	def SaveData (self, fname, verbose = True):

		if (verbose):
			print ("  Saving measurement to %s ... " % fname)

		start = time.clock()

		f = h5py.File(fname, 'w')

		f['data_r'] = np.squeeze(np.real(self.data).transpose())
		f['data_i'] = np.squeeze(np.imag(self.data).transpose())

		if (self.noise != 0):
			f['noise_r'] = np.squeeze(np.real(self.noise).transpose())
			f['noise_i'] = np.squeeze(np.imag(self.noise).transpose())
		
		if (self.acs != 0):
			f['acs_r'] = np.squeeze(np.real(self.acs).transpose())
			f['acs_i'] = np.squeeze(np.imag(self.acs).transpose())
		
		if (self.sync.any() != 0):
                        f['sync']  = self.sync.transpose()

		f.close()

		if (verbose):
			print '    ... saved in %(time).1f s.\n' % {"time": time.clock()-start}

		return

def FFT_Correlation(x,y):
    """
    FFT-based correlation, much faster than numpy autocorr.
    x and y are row-based vectors of arbitrary lengths.
    This is a vectorized implementation of O(N*log(N)) flops.
    """

    lengthx = x.shape[0]
    lengthy = y.shape[0]

    x = np.reshape(x,(1,lengthx))
    y = np.reshape(y,(1,lengthy))

    length = np.array([lengthx, lengthy]).min()
    
    x = x[:length]
    y = y[:length]
    
    fftx = fft(x, 2 * length - 1, axis=1) #pad with zeros
    ffty = fft(y, 2 * length - 1, axis=1)

    corr_xy = fft.ifft(fftx * np.conjugate(ffty), axis=1)
    corr_xy = np.real(fft.fftshift(corr_xy, axes=1)) #should be no imaginary part

    corr_yx = fft.ifft(ffty * np.conjugate(fftx), axis=1)
    corr_yx = np.real(fft.fftshift(corr_yx, axes=1))

    corr = 0.5 * (corr_xy[:,length:] + corr_yx[:,length:]) / range(1,length)[::-1]
    return np.reshape(corr,corr.shape[1])

def update_all_baseline_plots(i, fig, crawler, lines, norm_cross=False, forward=True):

    if forward:
        try:
            crawler.forward()
        except EOFError as err:
            print err
            raw_input("End of File. Press enter to quit.")
            sys.exit()

    burst = crawler

    for k in range(len(BASELINES)):
        if k < 4:
            #autos
            lines[k].set_data(FREQS, 10*np.log10(burst.autos[BASELINES[k]]))
            #overlays
            lines[-(k+1)].set_data(FREQS,10*np.log10(burst.autos[BASELINES[k]]))

        elif norm_cross:

			norm_val = np.array(burst.cross[BASELINES[k]])/np.sqrt(np.array(burst.autos[BASELINES[k][0]*2])*np.array(burst.autos[BASELINES[k][1]*2]))
			lines[k]['real'].set_data(FREQS, np.real(norm_val))
			lines[k]['imag'].set_data(FREQS, np.imag(norm_val))
        else:
			lines[k].set_data(FREQS, 10*np.log10(np.abs(np.real(burst.cross[BASELINES[k]]))))



    return lines

    def _exportDataToText(self, file):
        """Saves textual data to a file

        """
        Nt = self.data.shape[0]
        N = self.data.shape[1]
        # all elements [real] + (all elements - diagonal) [imaginary] + time
        Na = N + 1 + N*(N-1) 

        out = numpy.zeros((Nt, Na), dtype=numpy.float64)   
        
        for i in range(Nt):
            #print("%%%%")
            # time
            out[i,0] = self.TimeAxis.data[i]
            #print(0)
            # populations
            for j in range(N):
               out[i,j+1] = numpy.real(self.data[i,j,j])
               #print(j+1)
            # coherences
            l = 0
            for j in range(N):
                for k in range(j+1,N):
                    out[i,N+1+l] = numpy.real(self.data[i,j,k])
                    #print(N+1+l)
                    l += 1
                    out[i,N+1+l] = numpy.imag(self.data[i,j,k])
                    #print(N+1+l)
                    l += 1
                    
        numpy.savetxt(file, out)

    def eta_fil(self, x, V_app, apprx=(0, 0, 0, 0)):
        m_eff = self.m_r * const.electron_mass

        mpmath.mp.dps = 20
        x0 = Symbol('x0')  # eta_fil
        x1 = Symbol('x1')  # eta_ac
        x2 = Symbol('x2')  # eta_hop
        x3 = Symbol('x3')  # V_tunnel

        f0 = const.Boltzmann * self.T / (1 - self.alpha) / const.elementary_charge / self.z * \
             ln(self.A_fil/self.A_ac*(exp(- self.alpha * const.elementary_charge * self.z / const.Boltzmann / self.T * x0) - 1) + 1) - x1# eta_ac = f(eta_fil) x1 = f(x0)
        f1 = x*2*const.Boltzmann*self.T/self.a/self.z/const.elementary_charge*\
             asinh(self.j_0et/self.j_0hop*(exp(- self.alpha * const.elementary_charge * self.z / const.Boltzmann / self.T * x0) - 1)) - x2# eta_hop = f(eta_fil)
        f2 = x1 - x0 + x2 - x3

        f3 = -V_app + ((self.C * 3 * sqrt(2 * m_eff * ((4+x3/2)*const.elementary_charge)) / 2 / x * (const.elementary_charge / const.Planck)**2 * \
             exp(- 4 * const.pi * x / const.Planck * sqrt(2 * m_eff * ((4+x3/2)*const.elementary_charge))) * self.A_fil*x3)
                       + (self.j_0et*self.A_fil*(exp(-self.alpha*const.elementary_charge*self.z*x0/const.Boltzmann/self.T) - 1))) * (self.R_el + self.R_S + self.rho_fil*(self.L - x) / self.A_fil) \
             + x3

        eta_fil, eta_ac, eta_hop, V_tunnel = nsolve((f0, f1, f2, f3), [x0, x1, x2, x3], apprx)
        eta_fil = np.real(np.complex128(eta_fil))
        eta_ac = np.real(np.complex128(eta_ac))
        eta_hop = np.real(np.complex128(eta_hop))
        V_tunnel = np.real(np.complex128(V_tunnel))
        current = ((self.C * 3 * sqrt(2 * m_eff * ((4+V_tunnel)*const.elementary_charge)) / 2 / x * (const.elementary_charge / const.Planck)**2 * \
            exp(- 4 * const.pi * x / const.Planck * sqrt(2 * m_eff * ((4+V_tunnel)*const.elementary_charge))) * self.A_fil*V_tunnel)
                       + (self.j_0et*self.A_fil*(exp(-self.alpha*const.elementary_charge*self.z*eta_fil/const.Boltzmann/self.T) - 1)))
        print(eta_fil, eta_ac, eta_hop, V_tunnel)
        # print(eta_ac - eta_fil + eta_hop - V_tunnel)
        return eta_fil, eta_ac, eta_hop, V_tunnel, current

def window_fn_matrix(Q,N,num_remov=None,save_tag=None,lms=None):
    Q = n.matrix(Q); N = n.matrix(N)
    Ninv = uf.pseudo_inverse(N,num_remov=None) # XXX want to remove dynamically
    #print Ninv 
    info = n.dot(Q.H,n.dot(Ninv,Q))
    M = uf.pseudo_inverse(info,num_remov=num_remov)
    W = n.dot(M,info)

    if save_tag!=None:
        foo = W[0,:]
        foo = n.real(n.array(foo))
        foo.shape = (foo.shape[1]),
        print foo.shape
        p.scatter(lms[:,0],foo,c=lms[:,1],cmap=mpl.cm.PiYG,s=50)
        p.xlabel('l (color is m)')
        p.ylabel('W_0,lm')
        p.title('First Row of Window Function Matrix')
        p.colorbar()
        p.savefig('{0}/{1}_W.pdf'.format(fig_loc,save_tag))
        p.clf()

        print 'W ',W.shape
        p.imshow(n.real(W))
        p.title('Window Function Matrix')
        p.colorbar()
        p.savefig('{0}/{1}_W_im.pdf'.format(fig_loc,save_tag))
        p.clf()


    return W

def _test_random(test_case,inarr,outarr,tol):
    tc = test_case
    # Test that applying a transform and its inverse to reasonably long, random
    # input gives back the (appropriately scaled) input. We must allow for numerical
    # error, and it seems more reliable to check using normwise error (than elementwise).
    #
    # First test IFFT(FFT(random))
    # The numpy randn(n) provides an array of n numbers drawn from standard normal
    if dtype(inarr).kind == 'c':
        inarr._data[:] = randn(len(inarr)) +1j*randn(len(inarr))
        # If we're going to do a HC2R transform we must worry about DC/Nyquist imaginary
        if dtype(outarr).kind == 'f':
            inarr._data[0] = real(inarr[0])
            if (len(outarr)%2)==0:
                inarr._data[len(inarr)-1] = real(inarr[len(inarr)-1])
    else:
        inarr._data[:] = randn(len(inarr))
    incopy = type(inarr)(inarr)
    outarr.clear()
    # An FFT followed by IFFT gives Array scaled by len(Array), but for
    # Time/FrequencySeries there should be no scaling.
    if type(inarr) == pycbc.types.Array:
        incopy *= len(inarr)
    with tc.context:
        pycbc.fft.fft(inarr, outarr)
        pycbc.fft.ifft(outarr, inarr)
        emsg="IFFT(FFT(random)) did not reproduce original array to within tolerance {0}".format(tol)
        if isinstance(incopy,ts) or isinstance(incopy,fs):
            tc.assertTrue(incopy.almost_equal_norm(inarr,tol=tol,dtol=tol),
                          msg=emsg)
        else:
            tc.assertTrue(incopy.almost_equal_norm(inarr,tol=tol),
                          msg=emsg)
    # Perform arithmetic on outarr and inarr to pull them off of the GPU:
    outarr *= 1.0
    inarr *= 1.0
    # Now the same for FFT(IFFT(random))
    if dtype(outarr).kind == 'c':
        outarr._data[:] = randn(len(outarr))+1j*randn(len(outarr))
        # If we're going to do a HC2R transform we must worry about DC/Nyquist imaginary
        if dtype(inarr).kind == 'f':
            outarr._data[0] = real(outarr[0])
            if (len(inarr)%2)==0:
                outarr._data[len(outarr)-1] = real(outarr[len(outarr)-1])
    else:
        outarr._data[:] = randn(len(outarr))
    inarr.clear()
    outcopy = type(outarr)(outarr)
    if type(outarr) == pycbc.types.Array:
        outcopy *= len(inarr)
    with tc.context:
        pycbc.fft.ifft(outarr, inarr)
        pycbc.fft.fft(inarr, outarr)
        emsg="FFT(IFFT(random)) did not reproduce original array to within tolerance {0}".format(tol)
        if isinstance(outcopy,ts) or isinstance(outcopy,fs):
            tc.assertTrue(outcopy.almost_equal_norm(outarr,tol=tol,dtol=tol),
                          msg=emsg)
        else:
            tc.assertTrue(outcopy.almost_equal_norm(outarr,tol=tol),
                          msg=emsg)

    def icd(self, G1, G2):
        """Incomplete Cholesky decomposition
        """
        
        # remove mean. avoid standard calculation N0 = eye(N)-1/N*ones(N);
        G1 = G1 - numpy.array(numpy.mean(G1, 0), ndmin=2, copy=False)
        G2 = G2 - numpy.array(numpy.mean(G2, 0), ndmin=2, copy=False)

        R, D = self.method(G1, G2, self.reg)
        
        #solve generalized eigenvalues problem
        betas, alphas = scipy.linalg.eig(R,D)
        ind = numpy.argsort(numpy.real(betas))
        max_ind = ind[-1]
        alpha = alphas[:, max_ind]
        alpha = alpha/numpy.linalg.norm(alpha)
        beta = numpy.real(betas[max_ind])
        
        N1 = G1.shape[1]
        alpha1 = alpha[:N1]
        alpha2 = alpha[N1:]
        
        y1 = dot(G1, alpha1)
        y2 = dot(G2, alpha2)

        self.alpha1 = alpha1
        self.alpha2 = alpha2
        
        return (y1, y2, beta)                

def closest_point_on_the_parabola(a, px, py):
    thingy1 = 2. * a * py
    thingy2 = np.sqrt(-3 + 0j)
    thingy3 = 2 ** (1 / 3.)
    thingy4 = 2 ** (2 / 3.)
    thingy = (-108. * a ** 4 * px + np.sqrt(11664. * a ** 8 * px ** 2 - 864. * a ** 6 * (-1 + thingy1) ** 3 + 0j)) ** (
    1 / 3.)
    Aone = (thingy3 * (-1. + thingy1))
    Atwo = thingy
    Athree = thingy
    Afour = (6. * thingy3 * a ** 2)
    Bone = ((1. + thingy2) * (-1. + thingy1))
    Btwo = (thingy4 * thingy)
    Bthree = ((1. - thingy2) * thingy)
    Bfour = (12. * thingy3 * a ** 2)
    Cone = (1. - thingy2) * (-1 + thingy1)
    Ctwo = thingy4 * thingy
    Cthree = (1. + thingy2) * thingy
    Cfour = 12. * thingy3 * a ** 2

    A = -np.real(Aone / Atwo + Athree / Afour)
    B = np.real(Bone / Btwo + Bthree / Bfour)
    C = np.real(Cone / Ctwo + Cthree / Cfour)

    solns = [A, B, C]
    solns_temp = []
    for soln in solns:
        solns_temp.append(np.abs(soln - px))

    val, idx = min((val, idx) for (idx, val) in enumerate(solns_temp))
    return solns[idx]

    def simulate(self, steps, time, collect_dynamic = False):
        """!
        @brief Performs static simulation of oscillatory network.
        
        @param[in] steps (uint): Number simulation steps.
        @param[in] time (double): Time of simulation.
        @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
        
        @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' is True, than return dynamic for the whole simulation time,
                 otherwise returns only last values (last step of simulation) of output dynamic.
        
        @see simulate()
        @see simulate_dynamic()
        
        """
        
        dynamic_amplitude, dynamic_time = ([], []) if collect_dynamic is False else ([self.__amplitude], [0]);
        
        step = time / steps;
        int_step = step / 10.0;
        
        for t in numpy.arange(step, time + step, step):
            self.__amplitude = self.__calculate(t, step, int_step);
            
            if collect_dynamic is True:
                dynamic_amplitude.append([ numpy.real(amplitude)[0] for amplitude in self.__amplitude ]);
                dynamic_time.append(t);
        
        if collect_dynamic is False:
            dynamic_amplitude.append([ numpy.real(amplitude)[0] for amplitude in self.__amplitude ]);
            dynamic_time.append(time);

        output_sync_dynamic = fsync_dynamic(dynamic_amplitude, dynamic_time);
        return output_sync_dynamic;

def instantaneous_frequency(data, fs, fk):
    """
    Instantaneous frequency of a signal.

    Computes the instantaneous frequency of the given data which can be
    windowed or not. The instantaneous frequency is determined by the time
    derivative of the analytic signal of the input data.

    :type data: :class:`~numpy.ndarray`
    :param data: Data to determine instantaneous frequency of.
    :param fs: Sampling frequency.
    :param fk: Coefficients for calculating time derivatives
        (calculated via central difference).
    :return: **omega[, domega]** - Instantaneous frequency of input data, Time
        derivative of instantaneous frequency (windowed only).
    """
    x = envelope(data)
    if len(x[0].shape) > 1:
        omega = np.zeros(x[0].shape[0], dtype=np.float64)
        i = 0
        for row in x[0]:
            f = np.real(row)
            h = np.imag(row)
            # faster alternative to calculate f_add
            f_add = np.hstack(([f[0]] * (np.size(fk) // 2), f, [f[np.size(f) - 1]] * (np.size(fk) // 2)))
            fd = signal.lfilter(fk, 1, f_add)
            # correct start and end values of time derivative
            fd = fd[np.size(fk) - 1 : np.size(fd)]
            # faster alternative to calculate h_add
            h_add = np.hstack(([h[0]] * (np.size(fk) // 2), h, [h[np.size(h) - 1]] * (np.size(fk) // 2)))
            hd = signal.lfilter(fk, 1, h_add)
            # correct start and end values of time derivative
            hd = hd[np.size(fk) - 1 : np.size(hd)]
            omega_win = abs(((f * hd - fd * h) / (f * f + h * h)) * fs / 2 / np.pi)
            omega[i] = np.median(omega_win)
            i = i + 1
        # faster alternative to calculate omega_add
        omega_add = np.hstack(
            ([omega[0]] * (np.size(fk) // 2), omega, [omega[np.size(omega) - 1]] * (np.size(fk) // 2))
        )
        domega = signal.lfilter(fk, 1, omega_add)
        # correct start and end values of time derivative
        domega = domega[np.size(fk) - 1 : np.size(domega)]
        return omega, domega
    else:
        omega = np.zeros(np.size(x[0]), dtype=np.float64)
        f = np.real(x[0])
        h = np.imag(x[0])
        # faster alternative to calculate f_add
        f_add = np.hstack(([f[0]] * (np.size(fk) // 2), f, [f[np.size(f) - 1]] * (np.size(fk) // 2)))
        fd = signal.lfilter(fk, 1, f_add)
        # correct start and end values of time derivative
        fd = fd[np.size(fk) - 1 : np.size(fd)]
        # faster alternative to calculate h_add
        h_add = np.hstack(([h[0]] * (np.size(fk) // 2), h, [h[np.size(h) - 1]] * (np.size(fk) // 2)))
        hd = signal.lfilter(fk, 1, h_add)
        # correct start and end values of time derivative
        hd = hd[np.size(fk) - 1 : np.size(hd)]
        omega = abs(((f * hd - fd * h) / (f * f + h * h)) * fs / 2 / np.pi)
        return omega

def getArgumentVariable(_ComplexVariable):

	#Debug
	'''
	print('l 31 Numscipier')
	print('_ComplexVariable is ')
	print(_ComplexVariable)
	print('')
	'''

	#import
	import numpy as np

	#return
	return 2.*np.arctan(
	np.imag(_ComplexVariable)/(
	        np.sqrt(
	            np.imag(
	                _ComplexVariable
	            )**2+np.real(
	                _ComplexVariable
	            )**2)+np.real(
	                _ComplexVariable
	            )
	    )
	);

def test_morlet():
    """Test morlet with and without zero mean"""
    Wz = morlet(1000, [10], 2., zero_mean=True)
    W = morlet(1000, [10], 2., zero_mean=False)

    assert_true(np.abs(np.mean(np.real(Wz[0]))) < 1e-5)
    assert_true(np.abs(np.mean(np.real(W[0]))) > 1e-3)

def _select_function(sort, typ):
    if typ in ['F','D']:
        if callable(sort):
            # assume the user knows what they're doing
            sfunction = sort
        elif sort == 'lhp':
            sfunction = lambda x,y: (np.real(x/y) < 0.0)
        elif sort == 'rhp':
            sfunction = lambda x,y: (np.real(x/y) >= 0.0)
        elif sort == 'iuc':
            sfunction = lambda x,y: (abs(x/y) <= 1.0)
        elif sort == 'ouc':
            sfunction = lambda x,y: (abs(x/y) > 1.0)
        else:
            raise ValueError("sort parameter must be None, a callable, or "
                "one of ('lhp','rhp','iuc','ouc')")
    elif typ in ['f','d']:
        if callable(sort):
            # assume the user knows what they're doing
            sfunction = sort
        elif sort == 'lhp':
            sfunction = lambda x,y,z: (np.real((x+y*1j)/z) < 0.0)
        elif sort == 'rhp':
            sfunction = lambda x,y,z: (np.real((x+y*1j)/z) >= 0.0)
        elif sort == 'iuc':
            sfunction = lambda x,y,z: (abs((x+y*1j)/z) <= 1.0)
        elif sort == 'ouc':
            sfunction = lambda x,y,z: (abs((x+y*1j)/z) > 1.0)
        else:
            raise ValueError("sort parameter must be None, a callable, or "
                "one of ('lhp','rhp','iuc','ouc')")
    else:  # to avoid an error later
        raise ValueError("dtype %s not understood" % typ)
    return sfunction

def writeToDatFile(data, file, field="E", punits="mm", funits="V/m", pscale=1.0, fscale=1.0):
    """
    Write to field map to DAT data file.
    """

    nx = data.nx
    ny = data.ny
    nz = data.nz

    file.write("       x [{0}]         y [{0}]         z [{0}]     {1}xRe [{2}]     {1}yRe [{2}]     {1}zRe [{2}]     {1}xIm [{2}]     {1}yIm [{2}]     {1}zIm [{2}]  \r\n".format(punits, field, funits))
    file.write("------------------------------------------------------------------------------------------------------------------------------------------\r\n")

    for x in xrange(nx):
        px = pscale * data.px[x]
        for y in xrange(ny):
            py = pscale * data.py[y]
            for z in xrange(nz):
                pz = pscale * data.pz[z]
                fxre = fscale * numpy.real(data.fx[x,y,z])
                fyre = fscale * numpy.real(data.fy[x,y,z])
                fzre = fscale * numpy.real(data.fz[x,y,z])
                fxim = fscale * numpy.imag(data.fx[x,y,z])
                fyim = fscale * numpy.imag(data.fy[x,y,z])
                fzim = fscale * numpy.imag(data.fz[x,y,z])
                file.write("%13.1f  %13.1f  %13.1f  %13.6g  %13.6g  %13.6g  %13.6g  %13.6g  %13.6g\r\n" % (px, py, pz, fxre, fyre, fzre, fxim, fyim, fzim))

def p_sector5(MSA=[], C="seq", klist=[0, 1, 2, 3], Niter=20000, r=0.0001):
    if MSA.__class__ != msa:
        print "The function needs a msa as input"
        return False
    if C not in ["seq", "pos"]:
        print "gotta pick seq or pos"
        return False
    print "Calculating the SCA matrices"
    Cp, Cs = SCA5(MSA)

    print "Computing eigenmodes with numpy/LAPACK"
    if C == "seq":
        lbd, ev = numpy.linalg.eig(Cs)
    elif C == "pos":
        lbd, ev = numpy.linalg.eig(Cp)

        # ordering
    print "Sorting eigenmodes"
    idx = lbd.argsort()
    lbdsort = numpy.real(lbd[idx[::-1]])
    evsort = numpy.real(ev[:, idx[::-1]])

    print "top 10 eigenvalues"
    print lbdsort[:10]

    print "Calculating Independent Component Analysis"
    W, c = ICA(evsort, klist, Niter, r)

    return lbdsort, evsort, W, numpy.dot(W, evsort[:, klist].T).T

    def compress_data(X, k):
        Xtemp = X.dot(X.T)
        if len(Xtemp) == 0:
            return None
        e_vals, e_vecs = np.linalg.eig(Xtemp)

        e_vals = np.maximum(0.0, np.real(e_vals))
        e_vecs = np.real(e_vecs)

        idx = np.argsort(e_vals)[::-1]

        e_vals = e_vals[idx]
        e_vecs = e_vecs[:, idx]

        # Truncate to k dimensions
        if k < len(e_vals):
            e_vals = e_vals[:k]
            e_vecs = e_vecs[:, :k]

        # Normalize by the leading singular value of X
        Z = np.sqrt(e_vals.max())

        if Z > 0:
            e_vecs = e_vecs / Z

        return e_vecs.T.dot(X)

def _default_response_frequencies(A, n):
    """Compute a reasonable set of frequency points for bode plot.

    This function is used by `bode` to compute the frequency points (in rad/s)
    when the `w` argument to the function is None.

    Parameters
    ----------
    A : ndarray
        The system matrix, which is square.
    n : int
        The number of time samples to generate.

    Returns
    -------
    w : ndarray
        The 1-D array of length `n` of frequency samples (in rad/s) at which
        the response is to be computed.
    """
    vals = linalg.eigvals(A)
    # Remove poles at 0 because they don't help us determine an interesting
    # frequency range. (And if we pass a 0 to log10() below we will crash.)
    poles = [pole for pole in vals if pole != 0]
    # If there are no non-zero poles, just hardcode something.
    if len(poles) == 0:
        minpole = 1
        maxpole = 1
    else:
        minpole = min(abs(real(poles)))
        maxpole = max(abs(real(poles)))
    # A reasonable frequency range is two orders of magnitude before the
    # minimum pole (slowest) and two orders of magnitude after the maximum pole
    # (fastest).
    w = numpy.logspace(numpy.log10(minpole) - 2, numpy.log10(maxpole) + 2, n)
    return w

def filters_bank(M, N, J, L=8):
    filters = {}
    filters['psi'] = []

    offset_unpad = 0
    for j in range(J):
        for theta in range(L):
            psi = {}
            psi['j'] = j
            psi['theta'] = theta
            psi_signal = morlet_2d(M, N, 0.8 * 2**j, (int(L - L / 2 - 1) - theta) * np.pi / L, 3.0 / 4.0 * np.pi / 2**j,offset=offset_unpad)  # The 5 is here just to match the LUA implementation :)
            psi_signal_fourier = fft.fft2(psi_signal)
            for res in range(j + 1):
                psi_signal_fourier_res = crop_freq(psi_signal_fourier, res)
                psi[res] = torch.FloatTensor(np.stack((np.real(psi_signal_fourier_res), np.imag(psi_signal_fourier_res)), axis=2))
                # Normalization to avoid doing it with the FFT!
                psi[res].div_(M * N // 2**(2 * j))
            filters['psi'].append(psi)

    filters['phi'] = {}
    phi_signal = gabor_2d(M, N, 0.8 * 2**(J - 1), 0, 0, offset=offset_unpad)
    phi_signal_fourier = fft.fft2(phi_signal)
    filters['phi']['j'] = J
    for res in range(J):
        phi_signal_fourier_res = crop_freq(phi_signal_fourier, res)
        filters['phi'][res] = torch.FloatTensor(np.stack((np.real(phi_signal_fourier_res), np.imag(phi_signal_fourier_res)), axis=2))
        filters['phi'][res].div_(M * N // 2 ** (2 * J))

    return filters