Python scipy.fftpack.ifftn() Examples

def _irfftn(a, s=None, axes=None):
    return fft.ifftn(a, shape = s, axes = axes).real; 

def generate_fft_image(sequence):
    ff1 = fftn(sequence)
    fh = np.absolute(ifftn(ff1[1, :, :]))
    fh[fh < 0.1* np.max(fh)] =0.0
    #image=img_as_ubyte(fh/np.max(fh))
    return fh 

def read_fft_slice(path):
    d = pickle.load(open(path))['data']
    ff1 = fftn(d)
    fh = np.absolute(ifftn(ff1[1, :, :]))
    fh[fh < 0.1 * np.max(fh)] = 0.0
    d = 1. * fh / np.max(fh)
    d = np.expand_dims(d, axis=0)
    return d 

def convolve(arr1, arr2, dx=None, axes=None):
    """
    Performs a centred convolution of input arrays

    Parameters
    ----------
    arr1, arr2 : `numpy.ndarray`
        Arrays to be convolved. If dimensions are not equal then 1s are appended
        to the lower dimensional array. Otherwise, arrays must be broadcastable.
    dx : float > 0, list of float, or `None` , optional
        Grid spacing of input arrays. Output is scaled by
        `dx**max(arr1.ndim, arr2.ndim)`. default=`None` applies no scaling
    axes : tuple of ints or `None`, optional
        Choice of axes to convolve. default=`None` convolves all axes

    """
    if arr2.ndim > arr1.ndim:
        arr1, arr2 = arr2, arr1
        if axes is None:
            axes = range(arr2.ndim)
    arr2 = arr2.reshape(arr2.shape + (1,) * (arr1.ndim - arr2.ndim))

    if dx is None:
        dx = 1
    elif isscalar(dx):
        dx = dx ** (len(axes) if axes is not None else arr1.ndim)
    else:
        dx = prod(dx)

    arr1 = fftn(arr1, axes=axes)
    arr2 = fftn(ifftshift(arr2), axes=axes)
    out = ifftn(arr1 * arr2, axes=axes) * dx
    return require(out, requirements="CA") 

def ifftn(a, s=None, axes=None, norm=None, **_):
        return _ifftn(a, s, axes, norm) 

def plan_ifft(A, n=None, axis=None, norm=None, **_):
        """
        Plans an ifft for repeated use. Parameters are the same as for `pyfftw`'s `ifftn`
        which are, where possible, the same as the `numpy` equivalents.
        Note that some functionality is only possible when using the `pyfftw` backend.

        Parameters
        ----------
        A : `numpy.ndarray`, of dimension `d`
            Array of same shape to be input for the ifft
        n : iterable or `None`, `len(n) == d`, optional
            The output shape of ifft (default=`None` is same as `A.shape`)
        axis : `int`, iterable length `d`, or `None`, optional
            The axis (or axes) to transform (default=`None` is all axes)
        overwrite : `bool`, optional
            Whether the input array can be overwritten during computation
            (default=False)
        planner : {0, 1, 2, 3}, optional
            Amount of effort put into optimising Fourier transform where 0 is low
            and 3 is high (default=`1`).
        threads : `int`, `None`
            Number of threads to use (default=`None` is all threads)
        auto_align_input : `bool`, optional
            If `True` then may re-align input (default=`True`)
        auto_contiguous : `bool`, optional
            If `True` then may re-order input (default=`True`)
        avoid_copy : `bool`, optional
            If `True` then may over-write initial input (default=`False`)
        norm : {None, 'ortho'}, optional
            Indicate whether ifft is normalised (default=`None`)

        Returns
        -------
        plan : function
            Returns the inverse Fourier transform of `B`, `plan() == ifftn(B)`
        B : `numpy.ndarray`, `A.shape`
            Array which should be modified inplace for ifft to be computed. If
            possible, `B is A`.
        """
        return lambda: ifftn(A, n, axis, norm), A 

def bgtensor(img, lsigma, rho=0.2):
    eps = 1e-12
    fimg = fftn(img, overwrite_x=True)

    for s in lsigma:
        jvbuffer = bgkern3(kerlen=math.ceil(s)*6+1, sigma=s, rho=rho)
        jvbuffer = fftn(jvbuffer, shape=fimg.shape, overwrite_x=True) * fimg
        fimg = ifftn(jvbuffer, overwrite_x=True)
        yield hessian3(np.real(fimg)) 

def _pad_nans(x, head=None, tail=None):
    if np.ndim(x) == 1:
        if head is None and tail is None:
            return x
        elif head and tail:
            return np.r_[[np.nan] * head, x, [np.nan] * tail]
        elif tail is None:
            return np.r_[[np.nan] * head, x]
        elif head is None:
            return np.r_[x, [np.nan] * tail]
    elif np.ndim(x) == 2:
        if head is None and tail is None:
            return x
        elif head and tail:
            return np.r_[[[np.nan] * x.shape[1]] * head, x,
                         [[np.nan] * x.shape[1]] * tail]
        elif tail is None:
            return np.r_[[[np.nan] * x.shape[1]] * head, x]
        elif head is None:
            return np.r_[x, [[np.nan] * x.shape[1]] * tail]
    else:
        raise ValueError("Nan-padding for ndim > 2 not implemented")

#original changes and examples in sandbox.tsa.try_var_convolve

# don't do these imports, here just for copied fftconvolve
#get rid of these imports
#from scipy.fftpack import fft, ifft, ifftshift, fft2, ifft2, fftn, \
#     ifftn, fftfreq
#from numpy import product,array 

def fft(x, direction='C2C', inverse=False):
    """
        Interface with torch FFT routines for 2D signals.

        Example
        -------
        x = torch.randn(128, 32, 32, 2)
        x_fft = fft(x, inverse=True)

        Parameters
        ----------
        input : tensor
            complex input for the FFT
        direction : string
            'C2R' for complex to real, 'C2C' for complex to complex
        inverse : bool
            True for computing the inverse FFT.
            NB : if direction is equal to 'C2R', then an error is raised.

    """
    if direction == 'C2R':
        if not inverse:
            raise RuntimeError('C2R mode can only be done with an inverse FFT.')

    if direction == 'C2R':
        output = np.real(ifftn(x, axes=(-3, -2, -1)))
    elif direction == 'C2C':
        if inverse:
            output = ifftn(x, axes=(-3, -2, -1))
        else:
            output = fftn(x, axes=(-3, -2, -1))

    return output 

def read_fft_volume(data4D, harmonic=1):
    zslices = data4D.shape[2]
    tframes = data4D.shape[3]
    data3d_fft = np.empty((data4D.shape[:2]+(0,)))
    for slice in range(zslices):
        ff1 = fftn([data4D[:,:,slice, t] for t in range(tframes)])
        fh = np.absolute(ifftn(ff1[harmonic, :, :]))
        fh[fh < 0.1 * np.max(fh)] = 0.0
        image = 1. * fh / np.max(fh)
        # plt.imshow(image, cmap = 'gray')
        # plt.show()
        image = np.expand_dims(image, axis=2)
        data3d_fft = np.append(data3d_fft, image, axis=2)
    return data3d_fft 

def generate_fractal_surface(self, G):
        """Generate a 2D array with a fractal distribution.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        if self.xs == self.xf:
            surfacedims = (self.ny, self.nz)
        elif self.ys == self.yf:
            surfacedims = (self.nx, self.nz)
        elif self.zs == self.zf:
            surfacedims = (self.nx, self.ny)

        self.fractalsurface = np.zeros(surfacedims, dtype=complextype)

        # Positional vector at centre of array, scaled by weighting
        v1 = np.array([self.weighting[0] * (surfacedims[0]) / 2, self.weighting[1] * (surfacedims[1]) / 2])

        # 2D array of random numbers to be convolved with the fractal function
        R = np.random.RandomState(self.seed)
        A = R.randn(surfacedims[0], surfacedims[1])

        # 2D FFT
        A = fftpack.fftn(A)
        # Shift the zero frequency component to the centre of the array
        A = fftpack.fftshift(A)

        # Generate fractal
        generate_fractal2D(surfacedims[0], surfacedims[1], G.nthreads, self.b, self.weighting, v1, A, self.fractalsurface)

        # Shift the zero frequency component to start of the array
        self.fractalsurface = fftpack.ifftshift(self.fractalsurface)
        # Take the real part (numerical errors can give rise to an imaginary part) of the IFFT
        self.fractalsurface = np.real(fftpack.ifftn(self.fractalsurface))
        # Scale the fractal volume according to requested range
        fractalmin = np.amin(self.fractalsurface)
        fractalmax = np.amax(self.fractalsurface)
        fractalrange = fractalmax - fractalmin
        self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) \
            + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin 

def test_ifftn(self):
        overwritable = (np.complex128, np.complex64)
        for dtype in self.dtypes:
            self._check_nd(ifftn, dtype, overwritable) 

def test_invalid_sizes(self):
        assert_raises(ValueError, ifftn, [[]])
        assert_raises(ValueError, ifftn, [[1,1],[2,2]], (4, -3)) 

def test_random_complex(self):
        for size in [1,2,51,32,64,92]:
            x = random([size,size]) + 1j*random([size,size])
            assert_array_almost_equal_nulp(ifftn(fftn(x)),x,self.maxnlp)
            assert_array_almost_equal_nulp(fftn(ifftn(x)),x,self.maxnlp) 

def test_definition(self):
        x = np.array([[1,2,3],[4,5,6],[7,8,9]], dtype=self.dtype)
        y = ifftn(x)
        assert_equal(y.dtype, self.cdtype)
        assert_array_almost_equal_nulp(y,direct_idftn(x),self.maxnlp)
        x = random((20,26))
        assert_array_almost_equal_nulp(ifftn(x),direct_idftn(x),self.maxnlp)
        x = random((5,4,3,20))
        assert_array_almost_equal_nulp(ifftn(x),direct_idftn(x),self.maxnlp) 

def test_ifftn(self):
        overwritable = (np.complex128, np.complex64)
        for dtype in self.dtypes:
            self._check_nd(ifftn, dtype, overwritable) 

def test_random_complex(self):
        for size in [1,2,51,32,64,92]:
            x = random([size,size]) + 1j*random([size,size])
            assert_array_almost_equal_nulp(ifftn(fftn(x)),x,self.maxnlp)
            assert_array_almost_equal_nulp(fftn(ifftn(x)),x,self.maxnlp) 

def test_definition(self):
        x = np.array([[1,2,3],[4,5,6],[7,8,9]], dtype=self.dtype)
        y = ifftn(x)
        assert_(y.dtype == self.cdtype)
        assert_array_almost_equal_nulp(y,direct_idftn(x),self.maxnlp)
        x = random((20,26))
        assert_array_almost_equal_nulp(ifftn(x),direct_idftn(x),self.maxnlp)
        x = random((5,4,3,20))
        assert_array_almost_equal_nulp(ifftn(x),direct_idftn(x),self.maxnlp) 

def _pad_nans(x, head=None, tail=None):
    if np.ndim(x) == 1:
        if head is None and tail is None:
            return x
        elif head and tail:
            return np.r_[[np.nan] * head, x, [np.nan] * tail]
        elif tail is None:
            return np.r_[[np.nan] * head, x]
        elif head is None:
            return np.r_[x, [np.nan] * tail]
    elif np.ndim(x) == 2:
        if head is None and tail is None:
            return x
        elif head and tail:
            return np.r_[[[np.nan] * x.shape[1]] * head, x,
                         [[np.nan] * x.shape[1]] * tail]
        elif tail is None:
            return np.r_[[[np.nan] * x.shape[1]] * head, x]
        elif head is None:
            return np.r_[x, [[np.nan] * x.shape[1]] * tail]
    else:
        raise ValueError("Nan-padding for ndim > 2 not implemented")

#original changes and examples in sandbox.tsa.try_var_convolve

# don't do these imports, here just for copied fftconvolve
#get rid of these imports
#from scipy.fftpack import fft, ifft, ifftshift, fft2, ifft2, fftn, \
#     ifftn, fftfreq
#from numpy import product,array 

def vhartree_pbc(self, dens, **kw): 
  """  Compute Hartree potential for the density given in an equidistant grid  """
  from scipy.fftpack import fftn, ifftn
  sh = self.mesh3d.shape
  dens = dens.reshape(sh)
  vh = fftn(dens)    
  umom = self.ucell_mom()
  ii = [np.array([i-sh[j] if i>sh[j]//2 else i for i in range(sh[j])]) for j in range(3)]
  gg = [np.array([umom[j]*i for i in ii[j]]) for j in range(3)]
  vh = apply_inv_G2(vh, gg[0], gg[1], gg[2])
  vh = ifftn(vh).real*(4*np.pi)
  return vh 

def generate_fractal_volume(self, G):
        """Generate a 3D volume with a fractal distribution.

        Args:
            G (class): Grid class instance - holds essential parameters describing the model.
        """

        # Scale filter according to size of fractal volume
        if self.nx == 1:
            filterscaling = np.amin(np.array([self.ny, self.nz])) / np.array([self.ny, self.nz])
            filterscaling = np.insert(filterscaling, 0, 1)
        elif self.ny == 1:
            filterscaling = np.amin(np.array([self.nx, self.nz])) / np.array([self.nx, self.nz])
            filterscaling = np.insert(filterscaling, 1, 1)
        elif self.nz == 1:
            filterscaling = np.amin(np.array([self.nx, self.ny])) / np.array([self.nx, self.ny])
            filterscaling = np.insert(filterscaling, 2, 1)
        else:
            filterscaling = np.amin(np.array([self.nx, self.ny, self.nz])) / np.array([self.nx, self.ny, self.nz])

        # Adjust weighting to account for filter scaling
        self.weighting = np.multiply(self.weighting, filterscaling)

        self.fractalvolume = np.zeros((self.nx, self.ny, self.nz), dtype=complextype)

        # Positional vector at centre of array, scaled by weighting
        v1 = np.array([self.weighting[0] * self.nx / 2, self.weighting[1] * self.ny / 2, self.weighting[2] * self.nz / 2])

        # 3D array of random numbers to be convolved with the fractal function
        R = np.random.RandomState(self.seed)
        A = R.randn(self.nx, self.ny, self.nz)

        # 3D FFT
        A = fftpack.fftn(A)
        # Shift the zero frequency component to the centre of the array
        A = fftpack.fftshift(A)

        # Generate fractal
        generate_fractal3D(self.nx, self.ny, self.nz, G.nthreads, self.b, self.weighting, v1, A, self.fractalvolume)

        # Shift the zero frequency component to the start of the array
        self.fractalvolume = fftpack.ifftshift(self.fractalvolume)
        # Take the real part (numerical errors can give rise to an imaginary part) of the IFFT
        self.fractalvolume = np.real(fftpack.ifftn(self.fractalvolume))
        # Bin fractal values
        bins = np.linspace(np.amin(self.fractalvolume), np.amax(self.fractalvolume), self.nbins)
        for j in range(self.ny):
            for k in range(self.nz):
                self.fractalvolume[:, j, k] = np.digitize(self.fractalvolume[:, j, k], bins, right=True) 

def _fftconvolve(in1, in2, mode="full", axis=None):
    """ Convolve two N-dimensional arrays using FFT. See convolve.

    This is a fix of scipy.signal.fftconvolve, adding an axis argument and
    importing locally the stuff only needed for this function

    """
    #Locally import stuff only required for this:
    from scipy.fftpack import fftn, fft, ifftn, ifft
    from scipy.signal.signaltools import _centered
    from numpy import array, product


    s1 = array(in1.shape)
    s2 = array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complexfloating) or
                      np.issubdtype(in2.dtype, np.complexfloating))

    if axis is None:
        size = s1+s2-1
        fslice = tuple([slice(0, int(sz)) for sz in size])
    else:
        equal_shapes = s1==s2
        # allow equal_shapes[axis] to be False
        equal_shapes[axis] = True
        assert equal_shapes.all(), 'Shape mismatch on non-convolving axes'
        size = s1[axis]+s2[axis]-1
        fslice = [slice(l) for l in s1]
        fslice[axis] = slice(0, int(size))
        fslice = tuple(fslice)

    # Always use 2**n-sized FFT
    fsize = (2**np.ceil(np.log2(size))).astype(np.int64)
    if axis is None:
        IN1 = fftn(in1,fsize)
        IN1 *= fftn(in2,fsize)
        ret = ifftn(IN1)[fslice].copy()
    else:
        IN1 = fft(in1,fsize,axis=axis)
        IN1 *= fft(in2,fsize,axis=axis)
        ret = ifft(IN1,axis=axis)[fslice].copy()
    if not complex_result:
        del IN1
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if product(s1,axis=0) > product(s2,axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret,osize)
    elif mode == "valid":
        return _centered(ret,abs(s2-s1)+1) 

def _static_fft(self, image, mode='same'):
        """
        scipy fft convolution with saved static fft kernel

        :param image: 2d numpy array to be convolved
        :return:
        """
        in1 = image
        in1 = np.asarray(in1)
        if self._pre_computed is False:
            self._s1, self._s2, self._complex_result, self._shape, self._fshape, self._fslice, self._sp2 = self._static_pre_compute(image)
            self._pre_computed = True
        s1, s2, complex_result, shape, fshape, fslice, sp2 = self._s1, self._s2, self._complex_result, self._shape, self._fshape, self._fslice, self._sp2
        #if in1.ndim == in2.ndim == 0:  # scalar inputs
        #    return in1 * in2
        #elif not in1.ndim == in2.ndim:
        #    raise ValueError("in1 and in2 should have the same dimensionality")
        #elif in1.size == 0 or in2.size == 0:  # empty arrays
        #    return np.array([])


        # Check that input sizes are compatible with 'valid' mode
        #if _inputs_swap_needed(mode, s1, s2):
            # Convolution is commutative; order doesn't have any effect on output
            # only applicable for 'valid' mode
        #    in1, s1, in2, s2 = in2, s2, in1, s1

        # Pre-1.9 NumPy FFT routines are not threadsafe.  For older NumPys, make
        # sure we only call rfftn/irfftn from one thread at a time.
        if not complex_result and (_rfft_mt_safe or _rfft_lock.acquire(False)):
            try:
                sp1 = np.fft.rfftn(in1, fshape)
                ret = (np.fft.irfftn(sp1 * sp2, fshape)[fslice].copy())
            finally:
                if not _rfft_mt_safe:
                    _rfft_lock.release()
        else:
            # If we're here, it's either because we need a complex result, or we
            # failed to acquire _rfft_lock (meaning rfftn isn't threadsafe and
            # is already in use by another thread).  In either case, use the
            # (threadsafe but slower) SciPy complex-FFT routines instead.
            sp1 = fftpack.fftn(in1, fshape)
            ret = fftpack.ifftn(sp1 * sp2)[fslice].copy()
            if not complex_result:
                ret = ret.real

        if mode == "full":
            return ret
        elif mode == "same":
            return _centered(ret, s1)
        elif mode == "valid":
            return _centered(ret, s1 - s2 + 1)
        else:
            raise ValueError("Acceptable mode flags are 'valid',"
                             " 'same', or 'full'.") 

def fftconvolve(in1, in2, mode="full", axis=None):
    """ Convolve two N-dimensional arrays using FFT. See convolve.

    This is a fix of scipy.signal.fftconvolve, adding an axis argument and
    importing locally the stuff only needed for this function

    """
    s1 = np.array(in1.shape)
    s2 = np.array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complexfloating) or
                      np.issubdtype(in2.dtype, np.complexfloating))

    if axis is None:
        size = s1 + s2 - 1
        fslice = tuple([slice(0, int(sz)) for sz in size])
    else:
        equal_shapes = s1 == s2
        # allow equal_shapes[axis] to be False
        equal_shapes[axis] = True
        assert equal_shapes.all(), 'Shape mismatch on non-convolving axes'
        size = s1[axis] + s2[axis] - 1
        fslice = [slice(l) for l in s1]
        fslice[axis] = slice(0, int(size))
        fslice = tuple(fslice)

    # Always use 2**n-sized FFT
    fsize = 2 ** int(np.ceil(np.log2(size)))
    if axis is None:
        IN1 = fftpack.fftn(in1, fsize)
        IN1 *= fftpack.fftn(in2, fsize)
        ret = fftpack.ifftn(IN1)[fslice].copy()
    else:
        IN1 = fftpack.fft(in1, fsize, axis=axis)
        IN1 *= fftpack.fft(in2, fsize, axis=axis)
        ret = fftpack.ifft(IN1, axis=axis)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if np.product(s1, axis=0) > np.product(s2, axis=0):
            osize = s1
        else:
            osize = s2
        return signaltools._centered(ret, osize)
    elif mode == "valid":
        return signaltools._centered(ret, abs(s2 - s1) + 1) 

def fftconvolveinv(in1, in2, mode="full"):
    """Convolve two N-dimensional arrays using FFT. See convolve.

    copied from scipy.signal.signaltools, but here used to try out inverse filter
    doesn't work or I can't get it to work

    2010-10-23:
    looks ok to me for 1d,
    from results below with padded data array (fftp)
    but it doesn't work for multidimensional inverse filter (fftn)
    original signal.fftconvolve also uses fftn

    """
    s1 = np.array(in1.shape)
    s2 = np.array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1+s2-1

    # Always use 2**n-sized FFT
    fsize = 2**np.ceil(np.log2(size))
    IN1 = fft.fftn(in1,fsize)
    #IN1 *= fftn(in2,fsize) #JP: this looks like the only change I made
    IN1 /= fft.fftn(in2,fsize)  # use inverse filter
    # note the inverse is elementwise not matrix inverse
    # is this correct, NO  doesn't seem to work for VARMA
    fslice = tuple([slice(0, int(sz)) for sz in size])
    ret = fft.ifftn(IN1)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if np.product(s1,axis=0) > np.product(s2,axis=0):
            osize = s1
        else:
            osize = s2
        return trim_centered(ret,osize)
    elif mode == "valid":
        return trim_centered(ret,abs(s2-s1)+1)

#code duplication with fftconvolveinv 

def fft(self, data=None):
        """
        Perform the fft of `data`.

        Parameters:
            data (ndarray, optional): The data to transform. Optional as sometimes it can be faster to access ``inputData`` directly, though if and only if the data will be c-contiguous.

        Returns:
            ndarray: The transformed data
        """
        if self.FFTMODE=="pyfftw":

            if data is not None:
                self.inputData[:] = data
            if self.direction=="FORWARD":
                return self.fftwPlan()
            elif self.direction=="BACKWARD":
                return self.fftwPlan()

        elif self.FFTMODE=="gpu":

            if data is not None:
                self.inputData[:] = data

            inputData_dev = self.reikna_thread.to_device(self.inputData)

            self.reikna_ft_c(self.outputData_dev, inputData_dev, self.inverse)

            self.outputData = self.reikna_thread.from_device(
                                        self.outputData_dev)

            return  self.outputData



        elif self.FFTMODE=="scipy":
            if data is not None:
                 self.inputData = data
            if self.direction=="FORWARD":
                fft = fftpack.fftn(self.inputData,
                                shape=self.size, axes=self.axes)
            elif self.direction=="BACKWARD":
                fft = fftpack.ifftn(self.inputData,
                                shape=self.size,axes=self.axes)
            return fft 

def fftconvolveinv(in1, in2, mode="full"):
    """Convolve two N-dimensional arrays using FFT. See convolve.

    copied from scipy.signal.signaltools, but here used to try out inverse filter
    doesn't work or I can't get it to work

    2010-10-23:
    looks ok to me for 1d,
    from results below with padded data array (fftp)
    but it doesn't work for multidimensional inverse filter (fftn)
    original signal.fftconvolve also uses fftn

    """
    s1 = np.array(in1.shape)
    s2 = np.array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1+s2-1

    # Always use 2**n-sized FFT
    fsize = 2**np.ceil(np.log2(size))
    IN1 = fft.fftn(in1,fsize)
    #IN1 *= fftn(in2,fsize) #JP: this looks like the only change I made
    IN1 /= fft.fftn(in2,fsize)  # use inverse filter
    # note the inverse is elementwise not matrix inverse
    # is this correct, NO  doesn't seem to work for VARMA
    fslice = tuple([slice(0, int(sz)) for sz in size])
    ret = fft.ifftn(IN1)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if np.product(s1,axis=0) > np.product(s2,axis=0):
            osize = s1
        else:
            osize = s2
        return trim_centered(ret,osize)
    elif mode == "valid":
        return trim_centered(ret,abs(s2-s1)+1)

#code duplication with fftconvolveinv