import numpy as np
from scipy.fftpack import fftn, ifftn
__all__ = [
'fftdeconvolve'
]
[docs]def fftdeconvolve(x, h, eps=1e-8, axis=-1):
"""
Deconvolution by naive division in frequency domain.
Args:
x (np.ndarray): array to deconvolve.
h (np.ndarray): impulse response of deconvolution. Should have the same shape as `x`.
eps (float): threshold for the Fourier coefficients H to prevent division by small number.
axis (int): axis along which to devoncolve.
Returns:
y (np.ndarray): resulting array with same shape as `x`.
"""
x = np.asarray(x)
h = np.asarray(h)
if x.shape[axis] != h.shape[axis]:
raise ValueError('Signals with different lengths. Signals should have same length.')
# Compute FFTs.
X = fftn(x, axes=axis)
H = fftn(h, axes=axis)
if eps > 0:
# Clip H to prevent division by small number.
H_abs = np.abs(H)
H[H_abs < eps] = eps
# Divide in freq domain.
Y = np.divide(X, H)
# Transform back to time domain.
y = ifftn(Y, axes=axis).real
return y