I have a problem that requires computing convolution between the rows of two 2-D arrays, the convolution being performed in the last axis. The problem is that the grid on which the functions are defined is not necessarily equidistant, since there is large variation between grid points. Thus, I decided to approach the problem by mapping the grid to the interval [-1,1], see code below. While this seems to qualitatively reproduce the output as, for example, compared to a standard FFT convolution in equidistant grid,I am guessing the method has large numerical errors associated. Any good ideas on how to approach the problem along the proposed lines?
Note that interpolating into a very dense grid is NOT what I am looking for, as this needs to be repeated for many rows and memory consumption is expensive. Here is the working code for a minimal example comparing the equidistant grid output with an irregular grid of logarithmically spaced points.
import numpy as np
import matplotlib.pyplot as pyt
def create_logspaced_grid( discretization_parameter, w_lim, number_points ):
result = [0.0]
c = 0
if w_lim < 0:
w_lim *=-1
total_half_points = int( 0.5*(number_points-1) )
res_pos = []
for k in range( 0, total_half_points ):
res_pos.append( w_lim*discretization_parameter**( -0.5*k ) )
c += 1
res_neg = -np.flip( res_pos )
result = np.concatenate( ( res_neg, result ) )
result = np.concatenate( (result, res_pos))
result = np.array(result)
result = np.sort( result )
return result
def convolutions_non_equidistant_grid( x_grid, f, g ):
# we pass grid and f,g arguments to the function; the real axis is mapped
# to the [-1,+1] interval with u = tanh( x/D ) as the variable change, where
# D is the biggest scale of the problem.
# The functions f,g must enter as 2D-arrays, where the rows specify a single component
# and the columns the grid points
from scipy.interpolate import interp1d
from scipy.signal import fftconvolve
scale_factor = np.amax(x_grid)*0.9
u = np.tanh(x_grid/scale_factor)
N = x_grid.shape[0] #number of grid points; odd
#create new grid with linearly spaced points; double number of points
un = np.linspace( 0.0, u[-1], N )
un = np.concatenate( ( -np.flip( un[1:] ) , un ) )
ds = np.abs( un[0] - un[1] )
du = ( scale_factor /( 1.0 - un**2 ) )
#1) interpolate functions in the new space
interp_f = interp1d( u, f, kind="linear", axis=-1 )
interp_g = interp1d( u, g, kind="linear", axis=-1 )
f_new = interp_f( un )
g_new = interp_g( un )
#3)Calculate convolution
result = ds*fftconvolve( du*f_new, g_new, mode="same", axes=-1 )
#4) Get back to original grid by interpolating the result
# the result is the function in the original space
interp_res= interp1d( un, result, kind="linear", axis=-1 )
result = interp_res( u )
#compare with dense interpolation
xnew = np.linspace( -scale_factor, scale_factor, 30001 )
y_interp = interp1d( x_grid, f[0,:], kind="linear" )
g_interp = interp1d( x_grid, g[0,:], kind="linear" )
yy = y_interp(xnew)
gg = g_interp(xnew)
res_fft = np.abs( xnew[0] - xnew[1])*fftconvolve( yy, gg, mode="same")
pyt.plot( xnew, np.real(res_fft) )
pyt.plot( xnew, np.imag(res_fft))
pyt.plot( x_grid, np.real(result[0,:]), ls='--', marker='.' )
pyt.plot( x_grid, np.imag(result[0,:]), ls='--', marker='.')
return result
N = 2001
x_lim = 2000.0
x_min = 1e-4
alpha = ( 4/( N - 1 ) )*np.log( x_lim/x_min )
disc_la = np.exp(alpha)
x_grid = create_logspaced_grid( disc_la, x_lim, N )
f = np.array( [ 1.0/( x_grid + 10j ), ( x_grid /( x_grid + 1e-2j ) )] )
g = np.array( [ 1.0/( x_grid + 3j ), ( x_grid /( x_grid + 4.3e-1j ) )] )
convolutions_non_equidistant_grid( x_grid, f, g )