def _coherence_bavg(fxy, fxx, fyy):
    r"""
    Compute the band-averaged coherency between the spectra of two time series.
    input to this function is in the frequency domain

    Parameters
    ----------

    fxy : float array
         The cross-spectrum of the time series

    fyy,fxx : float array
         The spectra of the signals

    Returns
    -------

    float :
        the band-averaged coherence
    """
    if not np.isrealobj(fxx):
        fxx = np.real(fxx)
    if not np.isrealobj(fyy):
        fyy = np.real(fyy)

    return (np.abs(fxy.sum()) ** 2) / (fxx.sum() * fyy.sum())

def coherence_spec(fxy, fxx, fyy):
    r"""
    Compute the coherence between the spectra of two time series.

    Parameters of this function are in the frequency domain.

    Parameters
    ----------

    fxy : array
         The cross-spectrum of the time series

    fyy, fxx : array
         The spectra of the signals

    Returns
    -------

    float : a frequency-band-dependent measure of the linear association
        between the two time series

    See also
    --------
    :func:`coherence`
    """
    if not np.isrealobj(fxx):
        fxx = np.real(fxx)
    if not np.isrealobj(fyy):
        fyy = np.real(fyy)
    c = np.abs(fxy) ** 2 / (fxx * fyy)
    return c

def fftconv(a, b, axes=(0, 1)):
    """
    Compute a multi-dimensional convolution via the Discrete Fourier Transform.

    Parameters
    ----------
    a : array_like
      Input array
    b : array_like
      Input array
    axes : sequence of ints, optional (default (0, 1))
      Axes on which to perform convolution

    Returns
    -------
    ab : ndarray
      Convolution of input arrays, a and b, along specified axes
    """

    if np.isrealobj(a) and np.isrealobj(b):
        fft = rfftn
        ifft = irfftn
    else:
        fft = fftn
        ifft = ifftn
    dims = np.maximum([a.shape[i] for i in axes], [b.shape[i] for i in axes])
    af = fft(a, dims, axes)
    bf = fft(b, dims, axes)
    return ifft(af * bf, dims, axes)

 def same_upto_phase_factor(self, *others, **kwargs):
     if sanity_checking_enabled:
         if not np.isrealobj(self) or any(not np.isrealobj(o) for o in others):
             raise NotImplementedError("phase factor detection complex Tensors is not yet implemented")
     cutoff = kwargs.pop('cutoff', Tensor.same_tensor_cutoff)
     for other in others:
         if (abs(other - self) > cutoff).any() and (abs(other + self) > cutoff).any():
             return False
     return True

    def evaluate(self, ind, **kwargs):
        """
        Note that math functions used in the solutions are imported from either
        utilities.fitness.math_functions or called from numpy.

        :param ind: An individual to be evaluated.
        :param kwargs: An optional parameter for problems with training/test
        data. Specifies the distribution (i.e. training or test) upon which
        evaluation is to be performed.
        :return: The fitness of the evaluated individual.
        """

        dist = kwargs.get('dist', 'training')

        if dist == "training":
            # Set training datasets.
            x = self.training_in
            y = self.training_exp

        elif dist == "test":
            # Set test datasets.
            x = self.test_in
            y = self.test_exp

        else:
            raise ValueError("Unknown dist: " + dist)

        if params['OPTIMIZE_CONSTANTS']:
            # if we are training, then optimize the constants by
            # gradient descent and save the resulting phenotype
            # string as ind.phenotype_with_c0123 (eg x[0] +
            # c[0] * x[1]**c[1]) and values for constants as
            # ind.opt_consts (eg (0.5, 0.7). Later, when testing,
            # use the saved string and constants to evaluate.
            if dist == "training":
                return optimize_constants(x, y, ind)

            else:
                # this string has been created during training
                phen = ind.phenotype_consec_consts
                c = ind.opt_consts
                # phen will refer to x (ie test_in), and possibly to c
                yhat = eval(phen)
                assert np.isrealobj(yhat)

                # let's always call the error function with the
                # true values first, the estimate second
                return params['ERROR_METRIC'](y, yhat)

        else:
            # phenotype won't refer to C
            yhat = eval(ind.phenotype)
            assert np.isrealobj(yhat)

            # let's always call the error function with the true
            # values first, the estimate second
            return params['ERROR_METRIC'](y, yhat)

    def eigenConcat(omega, Q, AB, BB, k):
        """
        Find the eigen update of a matrix [A, B]'[A B] where  A'A = V diag(s) V*
        and AB = A*B, BB = B*B. Q is the set of eigenvectors of A*A and s is the
        vector of eigenvalues. 
        """
        #logging.debug("< eigenConcat >")
        Parameter.checkInt(k, 0, omega.shape[0])
        if not numpy.isrealobj(omega) or not numpy.isrealobj(Q):
            raise ValueError("Eigenvalues and eigenvectors must be real")
        if not numpy.isrealobj(AB) or not numpy.isrealobj(BB):
            raise ValueError("AB and BB must be real")
        if omega.ndim != 1:
            raise ValueError("omega must be 1-d array")
        if omega.shape[0] != Q.shape[1]:
            raise ValueError("Must have same number of eigenvalues and eigenvectors")
        if Q.shape[0] != AB.shape[0]:
            raise ValueError("Q must have the same number of rows as AB")
        if AB.shape[1] != BB.shape[0] or  BB.shape[0]!=BB.shape[1]:
            raise ValueError("AB must have the same number of cols/rows as BB")

        #Check Q is orthogonal
        if __debug__:
            Parameter.checkOrthogonal(Q, tol=EigenUpdater.tol, softCheck=True, arrayInfo = "input Q in eigenConcat()")

        m = Q.shape[0]
        p = BB.shape[0]

        inds = numpy.flipud(numpy.argsort(numpy.abs(omega)))
        Q = Q[:, inds[0:k]]
        omega = omega[inds[0:k]]
        Omega = numpy.diag(omega)

        QAB = Q.conj().T.dot(AB)

        F = numpy.c_[numpy.r_[Omega, QAB.conj().T], numpy.r_[QAB, BB]]
        D = numpy.c_[numpy.r_[Q, numpy.zeros((p, Q.shape[1]))], numpy.r_[numpy.zeros((m, p)), numpy.eye(p)]]

        pi, H = scipy.linalg.eigh(F)

        inds = numpy.flipud(numpy.argsort(numpy.abs(pi)))
        inds = inds[numpy.abs(pi)>EigenUpdater.tol]

        H = H[:, inds[0:k]]
        pi = pi[inds[0:k]]

        V = numpy.dot(D, H)

        #logging.debug("</ eigenConcat >")
        return pi, V

 def __setitem__(self, name, value):
     if not isinstance(value, ndarray):
         value   = array(value)
     if not isrealobj(value):                
         self.__matlab_object.handle.PutFullMatrix(name, self.__name_space, value.real, value.imag)                
     else:
         self.__matlab_object.handle.PutWorkspaceData(name, self.__name_space, value)

def irfft(x, n=None, axis=-1, overwrite_x=0):
    """ irfft(x, n=None, axis=-1, overwrite_x=0) -> y

    Return inverse discrete Fourier transform of real sequence x.
    The contents of x is interpreted as the output of rfft(..)
    function.

    The returned real array contains
      [y(0),y(1),...,y(n-1)]
    where for n is even
      y(j) = 1/n (sum[k=1..n/2-1] (x[2*k-1]+sqrt(-1)*x[2*k])
                                   * exp(sqrt(-1)*j*k* 2*pi/n)
                  + c.c. + x[0] + (-1)**(j) x[n-1])
    and for n is odd
      y(j) = 1/n (sum[k=1..(n-1)/2] (x[2*k-1]+sqrt(-1)*x[2*k])
                                   * exp(sqrt(-1)*j*k* 2*pi/n)
                  + c.c. + x[0])
    c.c. denotes complex conjugate of preceeding expression.

    Optional input: see rfft.__doc__
    """
    tmp = asarray(x)
    if not numpy.isrealobj(tmp):
        raise TypeError,"1st argument must be real sequence"
    if istype(tmp, numpy.float32):
        work_function = fftpack.rfft
    else:
        work_function = fftpack.drfft
    return _raw_fft(tmp,n,axis,-1,overwrite_x,work_function)

def app_luinv_to_spmat(alu_solve, Z):
    """ compute A.-1*Z  where A comes factored

    and with a solve routine for possibly complex Z

    Parameters
    ----------
    alu_solve : callable f(v)
        returning a matrix inverse applied to `v`
    Z : (N,K) ndarray, real or complex
        the inverse is to be applied to

    Returns
    -------
    , : (N,K) ndarray
        matrix inverse applied to ndarray

    """

    Z.tocsc()
    # ainvz = np.zeros(Z.shape)
    ainvzl = []  # to allow for complex values
    for ccol in range(Z.shape[1]):
        zcol = Z[:, ccol].toarray().flatten()
        if np.isrealobj(zcol):
            ainvzl.append(alu_solve(zcol))
        else:
            ainvzl.append(alu_solve(zcol.real) +
                          1j*alu_solve(zcol.imag))

    return np.asarray(ainvzl).T

def rfft(x, n=None, axis=-1, overwrite_x=False,
        planner_effort=None, threads=None,
        auto_align_input=True, auto_contiguous=True):
    '''Perform a 1D real FFT.

    The first three arguments are as per :func:`scipy.fftpack.rfft`;
    the rest of the arguments are documented
    in the :ref:`additional argument docs<interfaces_additional_args>`.
    '''
    if not numpy.isrealobj(x):
        raise TypeError('Input array must be real to maintain '
                'compatibility with scipy.fftpack.rfft.')

    x = numpy.asanyarray(x)
    planner_effort = _default_effort(planner_effort)
    threads = _default_threads(threads)

    complex_output = numpy_fft.rfft(x, n, axis, None, overwrite_x,
            planner_effort, threads, auto_align_input, auto_contiguous)

    output_shape = list(x.shape)
    if n is not None:
        output_shape[axis] = n

    return _complex_to_rfft_output(complex_output, output_shape, axis)

def irfft(
    x,
    n=None,
    axis=-1,
    overwrite_x=False,
    planner_effort="FFTW_MEASURE",
    threads=1,
    auto_align_input=True,
    auto_contiguous=True,
):
    """Perform a 1D real inverse FFT.
    
    The first three arguments are as per :func:`scipy.fftpack.irfft`; 
    the rest of the arguments are documented 
    in the :ref:`additional argument docs<interfaces_additional_args>`.
    """
    if not numpy.isrealobj(x):
        raise TypeError("Input array must be real to maintain " "compatibility with scipy.fftpack.irfft.")

    x = numpy.asanyarray(x)

    if n is None:
        n = x.shape[axis]

    complex_input = _irfft_input_to_complex(x, axis)

    return numpy_fft.irfft(
        complex_input, n, axis, overwrite_x, planner_effort, threads, auto_align_input, auto_contiguous
    )

    def __init__(
        self,
        shape=None,
        dtype=None,
        data=None,
        shape_out=None,
        axes=None,
        normalize="rescale",
        stream=None,
    ):
        if not(__have_cufft__) or not(__have_cufft__):
            raise ImportError("Please install pycuda and scikit-cuda to use the CUDA back-end")

        super(CUFFT, self).__init__(
            shape=shape,
            dtype=dtype,
            data=data,
            shape_out=shape_out,
            axes=axes,
            normalize=normalize,
        )
        self.cufft_stream = stream
        self.backend = "cufft"

        self.configure_batched_transform()
        self.allocate_arrays()
        self.real_transform = np.isrealobj(self.data_in)
        self.compute_forward_plan()
        self.compute_inverse_plan()
        self.refs = {
            "data_in": self.data_in,
            "data_out": self.data_out,
        }
        self.configure_normalization()

def pressureResponse(sim, forcex, forcez):

	nz = sim.nz
	nx = sim.nx

	transform = False
	if(np.isrealobj(forcex)):
		forcex = sim.makeSpect(forcex)
		forcez = sim.makeSpect(forcez)
		transform = True

	sim.assertDkx(forcex)
	sim.assertDkz(forcex)

	div = forcex * sim.dkx + forcez * sim.dkz

	pres = div / (sim.dkz**2 + sim.dkx**2)
	pres[0,0,0] = 0
 	solx = -pres * sim.dkx
	solz = -pres * sim.dkz

	if (transform == True):
		solx = np.fft.irfftn(solx)
		solz = np.fft.irfftn(solz)
		pres = np.fft.irfftn(pres)

	return [solx,solz,pres]

def rfft(x, n=None, axis=-1, overwrite_x=0):
    """ rfft(x, n=None, axis=-1, overwrite_x=0) -> y

    Return discrete Fourier transform of real sequence x.

    The returned real arrays contains
      [y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2))]              if n is even
      [y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2)),Im(y(n/2))]   if n is odd
    where
      y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n)
      j = 0..n-1
    Note that y(-j) = y(n-j).

    Optional input:
      n
        Defines the length of the Fourier transform. If n is not
        specified then n=x.shape[axis] is set. If n<x.shape[axis],
        x is truncated. If n>x.shape[axis], x is zero-padded.
      axis
        The transform is applied along the given axis of the input
        array (or the newly constructed array if n argument was used).
      overwrite_x
        If set to true, the contents of x can be destroyed.

    Notes:
      y == rfft(irfft(y)) within numerical accuracy.
    """
    tmp = asarray(x)
    if not numpy.isrealobj(tmp):
        raise TypeError,"1st argument must be real sequence"
    work_function = fftpack.drfft
    return _raw_fft(tmp,n,axis,1,overwrite_x,work_function)

def acorr(x, axis=-1, onesided=False, scale='none'):
    """Compute autocorrelation of x along given axis.

    Parameters
    ----------
        x : array-like
            signal to correlate.
        axis : int
            axis along which autocorrelation is computed.
        onesided: bool, optional
            if True, only returns the right side of the autocorrelation.
        scale: {'none', 'coeff'}
            scaling mode. If 'coeff', the correlation is normalized such as the
            0-lag is equal to 1.

    Notes
    -----
        Use fft for computation: is more efficient than direct computation for
        relatively large n.
    """
    if not np.isrealobj(x):
        raise ValueError("Complex input not supported yet")
    if not scale in ['none', 'coeff']:
        raise ValueError("scale mode %s not understood" % scale)

    maxlag = x.shape[axis]
    nfft = 2 ** nextpow2(2 * maxlag - 1)

    if axis != -1:
        x = np.swapaxes(x, -1, axis)
    a = _acorr_last_axis(x, nfft, maxlag, onesided, scale)
    if axis != -1:
        a = np.swapaxes(a, -1, axis)
    return a

def dctii(x):
    """Compute a Discrete Cosine Transform, type II.

    The DCT type II is defined as:

        \forall u \in 0...N-1, 
        dct(u) = a(u) sum_{i=0}^{N-1}{f(i)cos((i + 0.5)\pi u}

    Where a(0) = sqrt(1/(4N)), a(u) = sqrt(1/(2N)) for u > 0

    Parameters
    ==========
    x : array-like
        input signal

    Returns
    =======
    y : array-like
        DCT-II

    Note
    ====
    Use fft.
    """
    if not np.isrealobj(x):
        raise ValueError("Complex input not supported")
    n = x.size
    y = np.zeros(n * 4, x.dtype)
    y[1:2*n:2] = x
    y[2*n+1::2] = x[-1::-1]
    y = np.real(fft(y))[:n]
    y[0] *= np.sqrt(.25 / n)
    y[1:] *= np.sqrt(.5 / n)
    return y

def testRealEven(x):
    """
    Inputs:
        x (numpy array)= input signal of length M (M is odd)
    Output:
        The function should return a tuple (isRealEven, dftbuffer, X)
        isRealEven (boolean) = True if the input x is real and even,
                                and False otherwise
        dftbuffer (numpy array, possibly complex) = The M point zero
                                phase windowed version of x
        X (numpy array, possibly complex) = The M point DFT
                                            of dftbuffer
    """
    N = len(x)
    hM1 = np.floor((N+1)/2)
    hM2 = np.floor(N/2)
    dftbuffer = np.zeros(N)
    dftbuffer[:hM1] = x[hM2:]
    dftbuffer[-hM2:] = x[:hM2]

    Mx = fft(dftbuffer)
    Mx = Mx[: len(Mx)/2 + 1]
    xfh = x[: np.floor(len(x)/2)]  # First half
    # add np.float() otherwise len(x)/ is int so floor
    xsh = x[np.ceil(np.float(len(x))/2):]  # second half
    xsh = xsh[::-1]
    if np.isrealobj(x):
        if all(xfh == xsh):
            return (True, dftbuffer, Mx)
        else:
            return (False, dftbuffer, Mx)
    else:
        return (False, dftbuffer, Mx)

    def test_numpy_fft(self):
        """
        Test the numpy backend against native fft.
        Results should be exactly the same.
        """
        trinfos = self.param["transform_infos"]
        trdim = self.param["trdim"]
        ndim = len(self.param["size"])
        input_data = self.param["test_data"].data_refs[ndim].astype(trinfos.modes[self.param["mode"]])
        np_fft, np_ifft = self.transforms[trdim][np.isrealobj(input_data)]

        F = FFT(
            template=input_data,
            axes=trinfos.axes[trdim],
            backend="numpy"
        )
        # Test FFT
        res = F.fft(input_data)
        ref = np_fft(input_data)
        self.assertTrue(np.allclose(res, ref))

        # Test IFFT
        res2 = F.ifft(res)
        ref2 = np_ifft(ref)
        self.assertTrue(np.allclose(res2, ref2))

def _maybe_real(A, B, tol=None):
    """
    Return either B or the real part of B, depending on properties of A and B.

    The motivation is that B has been computed as a complicated function of A,
    and B may be perturbed by negligible imaginary components.
    If A is real and B is complex with small imaginary components,
    then return a real copy of B.  The assumption in that case would be that
    the imaginary components of B are numerical artifacts.

    Parameters
    ----------
    A : ndarray
        Input array whose type is to be checked as real vs. complex.
    B : ndarray
        Array to be returned, possibly without its imaginary part.
    tol : float
        Absolute tolerance.

    Returns
    -------
    out : real or complex array
        Either the input array B or only the real part of the input array B.

    """
    # Note that booleans and integers compare as real.
    if np.isrealobj(A) and np.iscomplexobj(B):
        if tol is None:
            tol = {0:feps*1e3, 1:eps*1e6}[_array_precision[B.dtype.char]]
        if np.allclose(B.imag, 0.0, atol=tol):
            B = B.real
    return B

def irfft(x, n=None, axis=-1, overwrite_x=0):
    """ irfft(x, n=None, axis=-1, overwrite_x=0) -> y

    Return inverse discrete Fourier transform of real sequence x.
    The contents of x is interpreted as the output of rfft(..)
    function.

    The returned real array contains
      [y(0),y(1),...,y(n-1)]
    where for n is even
      y(j) = 1/n (sum[k=1..n/2-1] (x[2*k-1]+sqrt(-1)*x[2*k])
                                   * exp(sqrt(-1)*j*k* 2*pi/n)
                  + c.c. + x[0] + (-1)**(j) x[n-1])
    and for n is odd
      y(j) = 1/n (sum[k=1..(n-1)/2] (x[2*k-1]+sqrt(-1)*x[2*k])
                                   * exp(sqrt(-1)*j*k* 2*pi/n)
                  + c.c. + x[0])
    c.c. denotes complex conjugate of preceeding expression.

    Optional input: see rfft.__doc__
    """
    tmp = _asfarray(x)
    if not numpy.isrealobj(tmp):
        raise TypeError,"1st argument must be real sequence"

    try:
        work_function = _DTYPE_TO_RFFT[tmp.dtype]
    except KeyError:
        raise ValueError("type %s is not supported" % tmp.dtype)

    return _raw_fft(tmp,n,axis,-1,overwrite_x,work_function)