Python scipy.fftpack.ifftn() Examples
The following are 27
code examples of scipy.fftpack.ifftn().
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example.
You may also want to check out all available functions/classes of the module
scipy.fftpack
, or try the search function
.
![](https://www.programcreek.com/common/static/images/search.png)
Example #1
Source File: Convolution.py From ClearMap with GNU General Public License v3.0 | 5 votes |
def _irfftn(a, s=None, axes=None): return fft.ifftn(a, shape = s, axes = axes).real;
Example #2
Source File: segmentation_labelling.py From kaggle-heart with MIT License | 5 votes |
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
Example #3
Source File: data.py From kaggle-heart with MIT License | 5 votes |
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
Example #4
Source File: fourier_transform.py From diffsims with GNU General Public License v3.0 | 5 votes |
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")
Example #5
Source File: fourier_transform.py From diffsims with GNU General Public License v3.0 | 5 votes |
def ifftn(a, s=None, axes=None, norm=None, **_): return _ifftn(a, s, axes, norm)
Example #6
Source File: fourier_transform.py From diffsims with GNU General Public License v3.0 | 5 votes |
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
Example #7
Source File: anisotropic.py From rivuletpy with BSD 3-Clause "New" or "Revised" License | 5 votes |
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))
Example #8
Source File: filtertools.py From Splunking-Crime with GNU Affero General Public License v3.0 | 5 votes |
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
Example #9
Source File: numpy_backend.py From kymatio with BSD 3-Clause "New" or "Revised" License | 5 votes |
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
Example #10
Source File: utils.py From Automated-Cardiac-Segmentation-and-Disease-Diagnosis with MIT License | 5 votes |
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
Example #11
Source File: fractals.py From gprMax with GNU General Public License v3.0 | 5 votes |
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
Example #12
Source File: test_basic.py From GraphicDesignPatternByPython with MIT License | 5 votes |
def test_ifftn(self): overwritable = (np.complex128, np.complex64) for dtype in self.dtypes: self._check_nd(ifftn, dtype, overwritable)
Example #13
Source File: test_basic.py From GraphicDesignPatternByPython with MIT License | 5 votes |
def test_invalid_sizes(self): assert_raises(ValueError, ifftn, [[]]) assert_raises(ValueError, ifftn, [[1,1],[2,2]], (4, -3))
Example #14
Source File: test_basic.py From GraphicDesignPatternByPython with MIT License | 5 votes |
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)
Example #15
Source File: test_basic.py From GraphicDesignPatternByPython with MIT License | 5 votes |
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)
Example #16
Source File: test_basic.py From Computable with MIT License | 5 votes |
def test_ifftn(self): overwritable = (np.complex128, np.complex64) for dtype in self.dtypes: self._check_nd(ifftn, dtype, overwritable)
Example #17
Source File: test_basic.py From Computable with MIT License | 5 votes |
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)
Example #18
Source File: test_basic.py From Computable with MIT License | 5 votes |
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)
Example #19
Source File: filtertools.py From vnpy_crypto with MIT License | 5 votes |
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
Example #20
Source File: m_vhartree_pbc.py From pyscf with Apache License 2.0 | 5 votes |
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
Example #21
Source File: fractals.py From gprMax with GNU General Public License v3.0 | 4 votes |
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)
Example #22
Source File: mtm.py From spectrum with BSD 3-Clause "New" or "Revised" License | 4 votes |
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)
Example #23
Source File: convolution.py From lenstronomy with MIT License | 4 votes |
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'.")
Example #24
Source File: dpss.py From wonambi with GNU General Public License v3.0 | 4 votes |
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)
Example #25
Source File: filtertools.py From Splunking-Crime with GNU Affero General Public License v3.0 | 4 votes |
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
Example #26
Source File: AOFFT.py From soapy with GNU General Public License v3.0 | 4 votes |
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
Example #27
Source File: filtertools.py From vnpy_crypto with MIT License | 4 votes |
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