From fde2f86954073fb141c83b491ca3c612a3ac617b Mon Sep 17 00:00:00 2001 From: Sascha Spors Date: Tue, 7 Feb 2017 18:53:14 +0100 Subject: [PATCH 01/10] Interpolation for better results --- sfs/time/source.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index 2d66a53d..05ea0bed 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -8,6 +8,8 @@ from __future__ import division import numpy as np +from scipy.interpolate import interp1d +from scipy.signal import resample from .. import util from .. import defs @@ -59,6 +61,16 @@ def point(xs, signal, t, grid, fs=None, c=None): # evaluate g over grid g_amplitude = 1 / (4 * np.pi * r) g_time = r / c - p = np.interp(t - g_time, np.arange(len(signal)) / fs, signal, + + oversampling = 10 + signal = resample(signal, oversampling * len(signal)) + + p = np.interp(t - g_time, np.arange(len(signal)) / (oversampling*fs), signal, left=0, right=0) + + + + #interpolator = interp1d(np.arange(len(signal)), signal, kind='cubic', bounds_error=False, fill_value=0) + #p = interpolator((t - g_time) * fs) + return p * g_amplitude From 5681f5f9789dd558c87268e33d10b66dcc4cc788 Mon Sep 17 00:00:00 2001 From: Sascha Spors Date: Tue, 7 Feb 2017 19:32:48 +0100 Subject: [PATCH 02/10] Used interpolation for calculation --- sfs/time/source.py | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index 05ea0bed..a5c3fa6f 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -14,7 +14,7 @@ from .. import defs -def point(xs, signal, t, grid, fs=None, c=None): +def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear'): r"""Source model for a point source: 3D Green's function. Calculates the scalar sound pressure field for a given point in @@ -61,16 +61,11 @@ def point(xs, signal, t, grid, fs=None, c=None): # evaluate g over grid g_amplitude = 1 / (4 * np.pi * r) g_time = r / c - - oversampling = 10 - signal = resample(signal, oversampling * len(signal)) - - p = np.interp(t - g_time, np.arange(len(signal)) / (oversampling*fs), signal, - left=0, right=0) - - - - #interpolator = interp1d(np.arange(len(signal)), signal, kind='cubic', bounds_error=False, fill_value=0) - #p = interpolator((t - g_time) * fs) - + + interpolator = interp1d(np.arange(len(signal)), signal, + kind=interpolator_kind, bounds_error=False, + fill_value=0) + + p = interpolator((t - g_time) * fs) + return p * g_amplitude From 70dcd888c29c4828b6f77f421ed938a87f4885b2 Mon Sep 17 00:00:00 2001 From: Sascha Spors Date: Tue, 7 Feb 2017 19:33:28 +0100 Subject: [PATCH 03/10] Removed import --- sfs/time/source.py | 1 - 1 file changed, 1 deletion(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index a5c3fa6f..4601ceff 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -9,7 +9,6 @@ from __future__ import division import numpy as np from scipy.interpolate import interp1d -from scipy.signal import resample from .. import util from .. import defs From 3d6bed5b1159d481ffd52842fa648ec0ea7231a2 Mon Sep 17 00:00:00 2001 From: Sascha Spors Date: Wed, 8 Feb 2017 13:20:34 +0100 Subject: [PATCH 04/10] Implemented sinc interpolation --- sfs/time/source.py | 34 ++++++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index 4601ceff..fcd6082f 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -51,6 +51,7 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear'): """ xs = util.asarray_1d(xs) signal = util.asarray_1d(signal) + t = util.asarray_1d(t) grid = util.as_xyz_components(grid) if fs is None: fs = defs.fs @@ -61,10 +62,35 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear'): g_amplitude = 1 / (4 * np.pi * r) g_time = r / c - interpolator = interp1d(np.arange(len(signal)), signal, - kind=interpolator_kind, bounds_error=False, - fill_value=0) + if interpolator_kind is 'sinc': + p = _sinc_interp(signal, np.arange(len(signal)), + np.array((t - g_time) * fs), fs) + else: + interpolator = interp1d(np.arange(len(signal)), signal, + kind=interpolator_kind, bounds_error=False, + fill_value=0) - p = interpolator((t - g_time) * fs) + p = interpolator((t - g_time) * fs) return p * g_amplitude + + +def _sinc_interp(x, s, u, fs): + """ + Interpolates x, sampled at "s" instants + Output y is sampled at "u" instants ("u" for "upsampled") + + from Matlab: + http://phaseportrait.blogspot.com/2008/06/sinc-interpolation-in-matlab.html + """ + + #if len(x) != len(s): + # raise Exception, 'x and s must be the same length' + + # Find the period + T = s[1] - s[0] + + sincM = np.tile(u, (len(s), 1)) - np.tile(s[:, np.newaxis], (1, len(u))) + y = np.dot(x, np.sinc(sincM/T)) + + return y From e1a88fda87547df10ec122254933db719334eaa2 Mon Sep 17 00:00:00 2001 From: Sascha Spors Date: Wed, 8 Feb 2017 13:44:19 +0100 Subject: [PATCH 05/10] Updated documentation --- sfs/time/source.py | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index fcd6082f..f17709c6 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -64,7 +64,7 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear'): if interpolator_kind is 'sinc': p = _sinc_interp(signal, np.arange(len(signal)), - np.array((t - g_time) * fs), fs) + np.array((t - g_time) * fs)) else: interpolator = interp1d(np.arange(len(signal)), signal, kind=interpolator_kind, bounds_error=False, @@ -75,21 +75,29 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear'): return p * g_amplitude -def _sinc_interp(x, s, u, fs): +def _sinc_interp(x, s, u): """ - Interpolates x, sampled at "s" instants - Output y is sampled at "u" instants ("u" for "upsampled") + Ideal sinc interpolation of a signal + adapted from https://gist.github.com/endolith/1297227 - from Matlab: - http://phaseportrait.blogspot.com/2008/06/sinc-interpolation-in-matlab.html - """ + Parameters + ---------- + x : (N,) array_like + Signal to be interpolated. + s : (N,) array_like + Sampling instants of signal. + u : (N,) array_like + Sampling instants after interpolation. - #if len(x) != len(s): - # raise Exception, 'x and s must be the same length' + Returns + ------- + numpy.ndarray + Interpolated signal + """ - # Find the period + # sampling period T = s[1] - s[0] - + # perform sinc interpolation sincM = np.tile(u, (len(s), 1)) - np.tile(s[:, np.newaxis], (1, len(u))) y = np.dot(x, np.sinc(sincM/T)) From 78e7461ab9224944e641314db2bdfce03b058a65 Mon Sep 17 00:00:00 2001 From: Till Rettberg Date: Wed, 8 Feb 2017 14:23:56 +0100 Subject: [PATCH 06/10] Add time signal zeropadding to point source --- sfs/time/source.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index fcd6082f..a14a1434 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -13,7 +13,7 @@ from .. import defs -def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear'): +def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear', zeropad=True): r"""Source model for a point source: 3D Green's function. Calculates the scalar sound pressure field for a given point in @@ -34,6 +34,8 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear'): Sampling frequency in Hertz. c : float, optional Speed of sound. + zeropad : bool, optional + zerodap time signals with 2 samples for spline interpolation. Returns ------- @@ -61,12 +63,16 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear'): # evaluate g over grid g_amplitude = 1 / (4 * np.pi * r) g_time = r / c - + if zeropad: + time_instants = np.arange(-2, len(signal)+2) + signal = np.concatenate([[0, 0], signal, [0, 0]]) + else: + time_instants = np.arange(len(signal)) if interpolator_kind is 'sinc': - p = _sinc_interp(signal, np.arange(len(signal)), + p = _sinc_interp(signal, time_instants, np.array((t - g_time) * fs), fs) else: - interpolator = interp1d(np.arange(len(signal)), signal, + interpolator = interp1d(time_instants, signal, kind=interpolator_kind, bounds_error=False, fill_value=0) @@ -81,7 +87,7 @@ def _sinc_interp(x, s, u, fs): Output y is sampled at "u" instants ("u" for "upsampled") from Matlab: - http://phaseportrait.blogspot.com/2008/06/sinc-interpolation-in-matlab.html + http://phaseportrait.blogspot.com/2008/06/sinc-interpolation-in-matlab.html """ #if len(x) != len(s): From ec587f3d7d04d5eb283d848b333e4508244cde85 Mon Sep 17 00:00:00 2001 From: Till Rettberg Date: Wed, 15 Feb 2017 12:24:48 +0100 Subject: [PATCH 07/10] Modify Sinc Interpolator to cope with desired times in ndarrays --- sfs/time/source.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index c76b3ab4..dc2604cf 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -105,7 +105,9 @@ def _sinc_interp(x, s, u): # sampling period T = s[1] - s[0] # perform sinc interpolation + shape = u.shape + u = u.flatten() sincM = np.tile(u, (len(s), 1)) - np.tile(s[:, np.newaxis], (1, len(u))) y = np.dot(x, np.sinc(sincM/T)) - + y = np.reshape(y, shape) return y From efdfb5d9ce77e110d90b842a7497f3b041aab791 Mon Sep 17 00:00:00 2001 From: Till Rettberg Date: Fri, 17 Feb 2017 16:33:51 +0100 Subject: [PATCH 08/10] Change API for interpolation in source. Reimplement sinc interpolation. --- sfs/time/source.py | 59 ++++++++++------------------------------------ 1 file changed, 13 insertions(+), 46 deletions(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index dc2604cf..9fc7432b 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -13,7 +13,8 @@ from .. import defs -def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear', zeropad=2): +def point(xs, signal, t, grid, fs=None, c=None, interpolator=interp1d, + zeropad=0, **kwargs): r"""Source model for a point source: 3D Green's function. Calculates the scalar sound pressure field for a given point in @@ -34,7 +35,10 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear', zero Sampling frequency in Hertz. c : float, optional Speed of sound. - zeropad : int + interpolator : function, optional + A function which constructs and returns a 1d interpolator. + (Expamples: interp1d, PchipInterpolator from scipy.interpolator) + zeropad : int, optional pre- & post-zeropad time signals with N samples. (For interpolation only.) @@ -64,50 +68,13 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator_kind='linear', zero # evaluate g over grid g_amplitude = 1 / (4 * np.pi * r) g_time = r / c - if zeropad: - time_instants = np.arange(-zeropad, len(signal)+zeropad) - signal = np.concatenate([np.zeros(zeropad), signal, np.zeros(zeropad)]) - else: - time_instants = np.arange(len(signal)) - if interpolator_kind is 'sinc': - p = _sinc_interp(signal, time_instants, - np.array((t - g_time) * fs)) - else: - interpolator = interp1d(time_instants, signal, - kind=interpolator_kind, bounds_error=False, - fill_value=0) - - p = interpolator((t - g_time) * fs) - + time_instants = np.arange(-zeropad, len(signal)+zeropad) + signal = np.concatenate([np.zeros(zeropad), signal, np.zeros(zeropad)]) + p = interpolator(time_instants, signal, **kwargs)((t - g_time) * fs) return p * g_amplitude -def _sinc_interp(x, s, u): - """ - Ideal sinc interpolation of a signal - adapted from https://gist.github.com/endolith/1297227 - - Parameters - ---------- - x : (N,) array_like - Signal to be interpolated. - s : (N,) array_like - Sampling instants of signal. - u : (N,) array_like - Sampling instants after interpolation. - - Returns - ------- - numpy.ndarray - Interpolated signal - """ - - # sampling period - T = s[1] - s[0] - # perform sinc interpolation - shape = u.shape - u = u.flatten() - sincM = np.tile(u, (len(s), 1)) - np.tile(s[:, np.newaxis], (1, len(u))) - y = np.dot(x, np.sinc(sincM/T)) - y = np.reshape(y, shape) - return y +def sincinterp(x, y, **kwargs): + def f(xnew): + return sum([y[i] * np.sinc(xnew - x[i]) for i in range(len(x))]) + return f From a620a96acf24cf0c6107d910ad0f7688840d56d9 Mon Sep 17 00:00:00 2001 From: Till Rettberg Date: Mon, 27 Feb 2017 14:30:20 +0100 Subject: [PATCH 09/10] Change API for interpolation in source. Add wrapper for linear interpolation. --- sfs/time/source.py | 42 ++++++++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index 9fc7432b..85152de4 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -13,8 +13,32 @@ from .. import defs -def point(xs, signal, t, grid, fs=None, c=None, interpolator=interp1d, - zeropad=0, **kwargs): +def linear_interpolator(x, y): + """1d linear interpolator with zero-padding. + + Parameters + ---------- + x : (N,) array_like + Sampling points. + y : (N,) array_like + Values at sampling points. + + Returns + ------- + function + Piecewise linear interpolant. + + """ + x = util.asarray_1d(x) + y = util.asarray_1d(y) + if len(x) < 3: + x = np.concatenate([np.array([min(x)-1]), x, np.array([max(x)+1])]) + y = np.concatenate([np.array([0]), y, np.array([0])]) + return interp1d(x, y, bounds_error=False, fill_value=0) + + +def point(xs, signal, t, grid, fs=None, c=None, + interpolator=linear_interpolator): r"""Source model for a point source: 3D Green's function. Calculates the scalar sound pressure field for a given point in @@ -37,10 +61,7 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator=interp1d, Speed of sound. interpolator : function, optional A function which constructs and returns a 1d interpolator. - (Expamples: interp1d, PchipInterpolator from scipy.interpolator) - zeropad : int, optional - pre- & post-zeropad time signals with N samples. - (For interpolation only.) + see: linear_interpolator, sinc_interpolator Returns ------- @@ -68,13 +89,14 @@ def point(xs, signal, t, grid, fs=None, c=None, interpolator=interp1d, # evaluate g over grid g_amplitude = 1 / (4 * np.pi * r) g_time = r / c - time_instants = np.arange(-zeropad, len(signal)+zeropad) - signal = np.concatenate([np.zeros(zeropad), signal, np.zeros(zeropad)]) - p = interpolator(time_instants, signal, **kwargs)((t - g_time) * fs) + p = interpolator(np.arange(len(signal)), signal)((t - g_time) * fs) return p * g_amplitude -def sincinterp(x, y, **kwargs): +def sinc_interpolator(x, y): + x = util.asarray_1d(x) + y = util.asarray_1d(y) + def f(xnew): return sum([y[i] * np.sinc(xnew - x[i]) for i in range(len(x))]) return f From 7a6977dd8fb081479efa563dca2e7306894bd511 Mon Sep 17 00:00:00 2001 From: Till Rettberg Date: Wed, 1 Nov 2017 16:11:08 +0100 Subject: [PATCH 10/10] Fix linear interpolation: Always use zeropadding --- sfs/time/source.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/sfs/time/source.py b/sfs/time/source.py index ddc559ba..f71e1bbe 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -31,9 +31,8 @@ def linear_interpolator(x, y): """ x = util.asarray_1d(x) y = util.asarray_1d(y) - if len(x) < 3: - x = np.concatenate([np.array([min(x)-1]), x, np.array([max(x)+1])]) - y = np.concatenate([np.array([0]), y, np.array([0])]) + x = np.concatenate([np.array([min(x)-1]), x, np.array([max(x)+1])]) + y = np.concatenate([np.array([0]), y, np.array([0])]) return interp1d(x, y, bounds_error=False, fill_value=0)