From 057c6897ba600804aab3cb0ca6c1340c00064a9a Mon Sep 17 00:00:00 2001 From: Till Rettberg Date: Sun, 11 Dec 2016 21:30:12 +0100 Subject: [PATCH 1/4] Add wfs_prefilter --- sfs/time/drivingfunction.py | 61 +++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/sfs/time/drivingfunction.py b/sfs/time/drivingfunction.py index b19fe76..5c3e8ff 100644 --- a/sfs/time/drivingfunction.py +++ b/sfs/time/drivingfunction.py @@ -197,3 +197,64 @@ def apply_delays(signal, delays, fs=None): for channel, cdelay in enumerate(delays_samples): out[cdelay:cdelay + len(signal), channel] = signal return out, offset_samples / fs + + +def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None): + """Get pre-equalization filter for WFS. + + Rising slope with 3dB/oct ('2.5D') or 6dB/oct ('2D' and '3D'). + Constant magnitude below fl and above fu. + Type 1 linear phase FIR filter of order N. + Simple design via "frequency sampling method". + + Parameters + ---------- + N : int, optional + Filter order, shall be even. For odd N, N+1 is used. + dim : str, optional + Dimensionality, must be '2D', '2.5D' or '3D'. + fl : int, optional + Lower corner frequency in Hertz. + fu : int, optional + Upper corner frequency in Hertz. + (Should be around spatial aliasing limit.) + fs : int, optional + Sampling frequency in Hertz. + c : float, optional + Speed of sound. + + Returns + ------- + h : (N+1,) numpy.ndarray + Filter taps. + delay : float + Pre-delay in seconds. + + """ + N = 2*(int(N + 1)//2) # for odd N, use N+1 instead + if fs is None: + fs = defs.fs + if c is None: + c = defs.c + + numbins = int(N/2 + 1) + delta_f = fs / (2*numbins - 1) + f = np.arange(numbins) * delta_f + if dim == '2D' or dim == '3D': + alpha = 1 + elif dim == '2.5D': + alpha = 0.5 + + desired = np.power(2 * np.pi * f / c, alpha) + low_shelf = np.power(2 * np.pi * fl / c, alpha) + high_shelf = np.power(2 * np.pi * fu / c, alpha) + + l_index = int(np.ceil(fl / delta_f)) + u_index = int(min(np.ceil(fu / delta_f), numbins - 1)) + desired[:l_index] = low_shelf + desired[u_index:] = min(high_shelf, desired[u_index]) + + h = np.fft.ifft(np.concatenate((desired, desired[-1:0:-1]))) + h = np.roll(np.real(h), numbins - 1) + delay = (numbins - 1) / fs + return h, delay From 763ec979dc62229d8a5cca6923b94588e3c5e11b Mon Sep 17 00:00:00 2001 From: Till Rettberg Date: Sun, 11 Dec 2016 21:54:12 +0100 Subject: [PATCH 2/4] Fix wfs_prefilter: normalize to unit energy. --- sfs/time/drivingfunction.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sfs/time/drivingfunction.py b/sfs/time/drivingfunction.py index 5c3e8ff..b6c4101 100644 --- a/sfs/time/drivingfunction.py +++ b/sfs/time/drivingfunction.py @@ -256,5 +256,6 @@ def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None): h = np.fft.ifft(np.concatenate((desired, desired[-1:0:-1]))) h = np.roll(np.real(h), numbins - 1) + h = h / np.sqrt(np.sum(abs(h)**2)) # normalize energy delay = (numbins - 1) / fs return h, delay From 5e17e24d0f6df1837c8140034ff5918c12ae432e Mon Sep 17 00:00:00 2001 From: Till Rettberg Date: Mon, 12 Dec 2016 09:00:46 +0100 Subject: [PATCH 3/4] Clean-up wfs_prefilter. --- sfs/time/drivingfunction.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/sfs/time/drivingfunction.py b/sfs/time/drivingfunction.py index b6c4101..1f1b55e 100644 --- a/sfs/time/drivingfunction.py +++ b/sfs/time/drivingfunction.py @@ -210,7 +210,7 @@ def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None): Parameters ---------- N : int, optional - Filter order, shall be even. For odd N, N+1 is used. + Filter order, shall be even. dim : str, optional Dimensionality, must be '2D', '2.5D' or '3D'. fl : int, optional @@ -231,7 +231,8 @@ def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None): Pre-delay in seconds. """ - N = 2*(int(N + 1)//2) # for odd N, use N+1 instead + if N % 2: + raise ValueError('N must be an even int.') if fs is None: fs = defs.fs if c is None: @@ -244,18 +245,13 @@ def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None): alpha = 1 elif dim == '2.5D': alpha = 0.5 - desired = np.power(2 * np.pi * f / c, alpha) low_shelf = np.power(2 * np.pi * fl / c, alpha) high_shelf = np.power(2 * np.pi * fu / c, alpha) + desired = np.clip(desired, low_shelf, high_shelf) - l_index = int(np.ceil(fl / delta_f)) - u_index = int(min(np.ceil(fu / delta_f), numbins - 1)) - desired[:l_index] = low_shelf - desired[u_index:] = min(high_shelf, desired[u_index]) - - h = np.fft.ifft(np.concatenate((desired, desired[-1:0:-1]))) - h = np.roll(np.real(h), numbins - 1) + h = np.fft.irfft(desired, 2*numbins - 1) + h = np.roll(h, numbins - 1) h = h / np.sqrt(np.sum(abs(h)**2)) # normalize energy delay = (numbins - 1) / fs return h, delay From 47e3cc17954160a2f02b9fc045da453605d94111 Mon Sep 17 00:00:00 2001 From: Till Rettberg Date: Tue, 24 Jan 2017 16:56:10 +0100 Subject: [PATCH 4/4] Add wrapper function for WFS prefiltering. Change default parameters of filter design. --- sfs/time/drivingfunction.py | 73 ++++++++++++++++++++++++++++--------- 1 file changed, 55 insertions(+), 18 deletions(-) diff --git a/sfs/time/drivingfunction.py b/sfs/time/drivingfunction.py index 1f1b55e..25c9f35 100644 --- a/sfs/time/drivingfunction.py +++ b/sfs/time/drivingfunction.py @@ -6,6 +6,7 @@ from __future__ import division import numpy as np from numpy.core.umath_tests import inner1d # element-wise inner product +from scipy.signal import fftconvolve from .. import defs from .. import util @@ -199,20 +200,19 @@ def apply_delays(signal, delays, fs=None): return out, offset_samples / fs -def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None): - """Get pre-equalization filter for WFS. +def wfs_25d_fir_prefilter(signal, N=128, fl=50, fu=1200, fs=None, c=None): + """Apply 2.5D pre-equalization to WFS source signal. - Rising slope with 3dB/oct ('2.5D') or 6dB/oct ('2D' and '3D'). - Constant magnitude below fl and above fu. - Type 1 linear phase FIR filter of order N. - Simple design via "frequency sampling method". + (Type 1 linear phase FIR filter of order N. + Rising slope with 3dB/oct between fl and fu. + Constant magnitude below fl and above fu.) Parameters ---------- + signal : (M,) array_like + Mono source signal, M samples. N : int, optional Filter order, shall be even. - dim : str, optional - Dimensionality, must be '2D', '2.5D' or '3D'. fl : int, optional Lower corner frequency in Hertz. fu : int, optional @@ -223,6 +223,47 @@ def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None): c : float, optional Speed of sound. + Returns + ------- + out : (M+N,) numpy.ndarray + Filtered output signal. + t_offset : float + Simulation point in time offset (seconds). + + """ + if fs is None: + fs = defs.fs + if c is None: + c = defs.c + h, delay = _wfs_prefilter_fir('2.5D', N, fl, fu, fs, c) + out = fftconvolve(signal, h) + return out, -delay + + +def _wfs_prefilter_fir(dim, N, fl, fu, fs, c): + """Create pre-equalization filter for WFS. + + Rising slope with 3dB/oct ('2.5D') or 6dB/oct ('2D' and '3D'). + Constant magnitude below fl and above fu. + Type 1 linear phase FIR filter of order N. + Simple design via "frequency sampling method". + + Parameters + ---------- + dim : str + Dimensionality, must be '2D', '2.5D' or '3D'. + N : int + Filter order, shall be even. + fl : int + Lower corner frequency in Hertz. + fu : int + Upper corner frequency in Hertz. + (Should be around spatial aliasing limit.) + fs : int + Sampling frequency in Hertz. + c : float + Speed of sound. + Returns ------- h : (N+1,) numpy.ndarray @@ -233,14 +274,10 @@ def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None): """ if N % 2: raise ValueError('N must be an even int.') - if fs is None: - fs = defs.fs - if c is None: - c = defs.c - numbins = int(N/2 + 1) - delta_f = fs / (2*numbins - 1) - f = np.arange(numbins) * delta_f + bins = int(N/2 + 1) + delta_f = fs / (2*bins - 1) + f = np.arange(bins) * delta_f if dim == '2D' or dim == '3D': alpha = 1 elif dim == '2.5D': @@ -250,8 +287,8 @@ def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None): high_shelf = np.power(2 * np.pi * fu / c, alpha) desired = np.clip(desired, low_shelf, high_shelf) - h = np.fft.irfft(desired, 2*numbins - 1) - h = np.roll(h, numbins - 1) + h = np.fft.irfft(desired, 2*bins - 1) + h = np.roll(h, bins - 1) h = h / np.sqrt(np.sum(abs(h)**2)) # normalize energy - delay = (numbins - 1) / fs + delay = (bins - 1) / fs return h, delay