Linear-phase filterbank for radial filter design

examples/em32.txt

@@ -0,0 +1,32 @@
+0.000000000000000000e+00 1.204277183876087287e+00 4.200000000000000261e-02
+5.585053606381854552e-01 1.570796326794896558e+00 4.200000000000000261e-02
+0.000000000000000000e+00 1.937315469713705829e+00 4.200000000000000261e-02
+5.724679946541400888e+00 1.570796326794896558e+00 4.200000000000000261e-02
+0.000000000000000000e+00 5.585053606381854552e-01 4.200000000000000261e-02
+7.853981633974482790e-01 9.599310885968812546e-01 4.200000000000000261e-02
+1.204277183876087287e+00 1.570796326794896558e+00 4.200000000000000261e-02
+7.853981633974482790e-01 2.181661564992912083e+00 4.200000000000000261e-02
+0.000000000000000000e+00 2.583087292951607772e+00 4.200000000000000261e-02
+5.497787143782137953e+00 2.181661564992912083e+00 4.200000000000000261e-02
+5.078908123303499167e+00 1.570796326794896558e+00 4.200000000000000261e-02
+5.497787143782137953e+00 9.599310885968812546e-01 4.200000000000000261e-02
+1.588249619314839878e+00 3.665191429188092154e-01 4.200000000000000261e-02
+1.570796326794896558e+00 1.012290966156711214e+00 4.200000000000000261e-02
+1.570796326794896558e+00 2.111848394913138804e+00 4.200000000000000261e-02
+1.553343034274953238e+00 2.775073510670984067e+00 4.200000000000000261e-02
+3.141592653589793116e+00 1.204277183876087287e+00 4.200000000000000261e-02
+3.700098014227978460e+00 1.570796326794896558e+00 4.200000000000000261e-02
+3.141592653589793116e+00 1.937315469713705829e+00 4.200000000000000261e-02
+2.583087292951607772e+00 1.570796326794896558e+00 4.200000000000000261e-02
+3.141592653589793116e+00 5.585053606381854552e-01 4.200000000000000261e-02
+3.926990816987241395e+00 9.599310885968812546e-01 4.200000000000000261e-02
+4.345869837465880181e+00 1.570796326794896558e+00 4.200000000000000261e-02
+3.926990816987241395e+00 2.181661564992912083e+00 4.200000000000000261e-02
+3.141592653589793116e+00 2.583087292951607772e+00 4.200000000000000261e-02
+2.356194490192344837e+00 2.181661564992912083e+00 4.200000000000000261e-02
+1.937315469713705829e+00 1.570796326794896558e+00 4.200000000000000261e-02
+2.356194490192344837e+00 9.599310885968812546e-01 4.200000000000000261e-02
+4.694935687864746576e+00 3.665191429188092154e-01 4.200000000000000261e-02
+4.712388980384689674e+00 1.012290966156711214e+00 4.200000000000000261e-02
+4.712388980384689674e+00 2.129301687433081902e+00 4.200000000000000261e-02
+4.729842272904632772e+00 2.775073510670984067e+00 4.200000000000000261e-02

examples/linph-filterbank-butterworth-fd.py

@@ -0,0 +1,80 @@
+"""Linear-phase filterbank.
+
+- The target magnitude responses for the filter-bank is designed
+  by using zero-phase Butterworth responses.
+  (not to confused with typical (minphase) Butterworth filters)
+
+    Reference
+    ---------
+    Franz Zotter, "A linear-phase filter-bank approach to process
+    rigid spherical microphone array recordings", in Proc. IcETRAN,
+    Palic, Serbia, 2018. (see Fig. 7)
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from micarray.util import maxre_sph, modal_norm, db
+from micarray.modal.radial \
+    import crossover_frequencies, spherical_hn2, tf_linph_filterbank
+
+c = 343
+R = 0.049
+N = 4
+max_boost = 30
+f_xo = crossover_frequencies(N, R, max_boost, modal_weight=maxre_sph)
+
+fmin, fmax, numf = 10, 20000, 2000
+f = np.logspace(np.log10(fmin), np.log10(fmax), num=numf)
+kr = 2 * np.pi * f / c * R
+
+H_fbank = tf_linph_filterbank(f_xo, f, type='butter')
+H_tot = np.sum(H_fbank, axis=0)
+
+# Prototpye radial filters
+H_proto = np.stack([
+        1j**(-n-1) * (kr)**2 * spherical_hn2(n, kr, derivative=True)
+        for n in range(N+1)])
+
+# Equalized radial filters
+H_radial = np.zeros_like(H_proto)
+for i, Hi in enumerate(H_fbank):
+    ai = maxre_sph(i)
+    ai *= 1 / modal_norm(ai)
+    for n, Hn in enumerate(H_proto[:i+1]):
+        H_radial[n] += Hi * ai[n]
+H_radial *= modal_norm(maxre_sph(N))
+
+
+# Plots
+# Filter-bank
+fig, ax = plt.subplots()
+for i, Hi in enumerate(H_fbank):
+    ax.semilogx(f, db(Hi), lw=3, label='${}$'.format(i), alpha=0.5)
+for fx in f_xo:
+    ax.semilogx(fx, 0, 'kv')
+    ax.text(fx, 0, '{:0.1f} Hz'.format(fx), rotation=30,
+            horizontalalignment='left', verticalalignment='bottom')
+ax.semilogx(f, db(H_tot), 'k:', label='Sum')
+ax.set_xlim(fmin, fmax)
+ax.set_ylim(-100, 12)
+ax.set_xscale('log')
+ax.grid(True)
+ax.set_xlabel('frequency in Hz')
+ax.set_ylabel('magnitude in dB')
+ax.legend(title='subband')
+plt.savefig('./linph-filterbank-fd.png', bbox_inches='tight')
+
+# Normalized radial filters
+fig, ax = plt.subplots()
+for n, (Hr, Hp) in enumerate(zip(H_radial, H_proto)):
+    ax.semilogx(f, db(Hp), c='k', ls=':')
+    ax.semilogx(f, db(Hr * Hp), lw=3, label='${}$'.format(n), alpha=0.5)
+ax.hlines(max_boost, xmin=fmin, xmax=fmax, colors='C3', linestyle='--',
+          label='max. boost')
+ax.set_xlim(fmin, fmax)
+ax.set_ylim(-23, 33)
+ax.set_xscale('log')
+ax.grid(True)
+ax.set_xlabel('frequency in Hz')
+ax.set_ylabel('magnitude in dB')
+ax.legend(title='order', loc='lower right', ncol=2)
+plt.savefig('./linph-filterbank-butterworth-fd.png', bbox_inches='tight')

examples/linph-filterbank-butterworth-td.py

@@ -0,0 +1,89 @@
+"""Impulse responses of linear-phase filterbank and radial filters.
+
+    Reference
+    ---------
+    Franz Zotter, "A linear-phase filter-bank approach to process
+    rigid spherical microphone array recordings", in Proc. IcETRAN,
+    Palic, Serbia, 2018.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from micarray.util import maxre_sph, modal_norm
+from micarray.modal.radial \
+    import crossover_frequencies, spherical_hn2, tf_linph_filterbank
+
+
+c = 343
+fs = 44100
+Nfft = 2048
+
+R = 0.049
+N = 5
+max_boost = 30
+f_xo = crossover_frequencies(N, R, max_boost, modal_weight=maxre_sph)
+
+fmin, fmax, numf = 0, fs/2, int(Nfft/2)+1
+f = np.linspace(fmin, fmax, num=numf)
+f[0] = 0.5 * f[1]
+kr = 2 * np.pi * f / c * R
+H_fbank = tf_linph_filterbank(f_xo, f, type='butter')
+
+# Prototype radial filters
+H_proto = np.stack([1j**(-n-1) * (kr)**2
+                    * spherical_hn2(n, kr, derivative=True)
+                    for n in range(N+1)])
+
+H_radial = np.zeros_like(H_proto)
+for i, Hi in enumerate(H_fbank):
+    ai = maxre_sph(i)
+    ai *= 1 / modal_norm(ai)
+    for n, Hr in enumerate(H_proto[:i+1]):
+        H_radial[n] += Hr * Hi * ai[n]
+H_radial *= modal_norm(maxre_sph(N))
+
+# inverse DFT
+h_fbank = np.stack([np.fft.irfft(Hi, n=Nfft, axis=-1) for Hi in H_fbank])
+h_radial = np.stack([np.fft.irfft(Hi, n=Nfft, axis=-1) for Hi in H_radial])
+h_fbank = np.roll(h_fbank, int(Nfft/2), axis=-1)
+h_radial = np.roll(h_radial, int(Nfft/2), axis=-1)
+
+t = ((np.arange(Nfft) - Nfft/2) / fs) * 1000
+t_R = t - R/c*1000
+
+
+# Plots
+def decorate_subplots(axes, **kwargs):
+    for ax in axes.flat:
+        if ax.is_first_col():
+            ax.set_ylabel('Amplitude in dB')
+        if ax.is_last_row():
+            ax.set_xlabel('Time in ms')
+        ax.grid(True)
+        ax.set(**kwargs)
+
+
+# Impulse responses
+figsize = (8, 10)
+gridspec_kw = {'wspace': 0.1}
+tlim = -2, 2
+ylim = -40, None
+
+# each filter bank
+fig, ax = plt.subplots(figsize=figsize, ncols=2, nrows=3,
+                       sharex=True, sharey=True, gridspec_kw=gridspec_kw)
+for i, (axi, hi) in enumerate(zip(ax.flat[:N+1], h_fbank)):
+    hi *= 1 / np.max(np.abs(hi))
+    axi.plot(t, hi)
+    axi.set_title('Subband #{}'.format(i))
+decorate_subplots(ax, xlim=tlim)
+plt.savefig('./linph-filterbank-td.png')
+
+# each order
+fig, ax = plt.subplots(figsize=figsize, ncols=2, nrows=3,
+                       sharex=True, sharey=True, gridspec_kw=gridspec_kw)
+for i, (axi, hi) in enumerate(zip(ax.flat[:N+1], h_radial)):
+    hi *= 1 / np.max(np.abs(hi))
+    axi.plot(t_R, hi)
+    axi.set_title('Order #{}'.format(i))
+decorate_subplots(ax, xlim=tlim)
+plt.savefig('./radialfilters-td.png')

examples/linph-filterbank-crossover-frequencies.py

@@ -0,0 +1,59 @@
+"""Cross-over frequencies for linear-phase filterbank.
+
+- filter bank design for Ambisonic encoding of rigid sphere signals
+- determine cross-over frequencies based on pre-defined maximum boost
+- exploit small argument approximation of spherical Hankel functions
+
+    Reference
+    ---------
+    Franz Zotter, "A linear-phase filter-bank approach to process
+    rigid spherical microphone array recordings", in Proc. IcETRAN,
+    Palic, Serbia, 2018. (see Fig. 6)
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from micarray.util import maxre_sph, double_factorial, modal_norm, db
+from micarray.modal.radial import crossover_frequencies, spherical_hn2
+
+c = 343
+N = 4
+R = 0.049
+max_boost = 30
+f_xo = crossover_frequencies(N, R, max_boost, modal_weight=maxre_sph)
+band_gain = [1 / modal_norm(maxre_sph(n)) for n in range(N+1)]
+
+fmin, fmax, numf = 10, 20000, 2000
+f = np.logspace(np.log10(fmin), np.log10(fmax), num=numf)
+kr = 2 * np.pi * f / c * R
+
+# Plot
+fig, ax = plt.subplots(figsize=(6, 4))
+for n in range(N+1):
+
+    # Analytic radial filter
+    radfilt = band_gain[n] * (kr)**2 * spherical_hn2(n, kr, derivative=True)
+    ax.semilogx(f, db(radfilt), lw=3, alpha=0.5, zorder=0,
+                label='${}$'.format(n))
+
+    if n > 0:
+        fn = f_xo[n-1]
+        krn = 2 * np.pi * fn / c * R
+
+        # Low-frequency (small argument) approximation
+        lf_approx = band_gain[n] * double_factorial(2*n-1) * (n+1) / (kr)**n
+        gain_at_crossover = \
+            band_gain[n] * (krn)**2 * spherical_hn2(n, krn, derivative=True)
+
+        ax.semilogx(f, db(lf_approx), c='black', ls=':')
+        ax.semilogx(fn, db(gain_at_crossover), 'C0o')
+        ax.text(fn, db(gain_at_crossover), '{:3.1f} Hz'.format(fn),
+                ha='left', va='bottom', rotation=55)
+ax.hlines(max_boost, xmin=fmin, xmax=fmax,
+          colors='C3', linestyle='--', label='max. boost')
+ax.set_xlim(fmin, fmax)
+ax.set_ylim(-10, 90)
+ax.grid(True)
+ax.set_xlabel('frequency in Hz')
+ax.set_ylabel('magnitude in dB')
+ax.legend(title='Order', ncol=2)
+plt.savefig('./crossover-frequencies.png', bbox_inches='tight')

examples/linph-filterbank-modal-beamforming-em32.py

@@ -0,0 +1,134 @@
+"""Fourth-order Ambisonics microphone Eigenmike em32.
+
+- filter bank design for Ambisonic encoding of em32 signals
+
+    Reference
+    ---------
+    Franz Zotter, "A linear-phase filter-bank approach to process
+    rigid spherical microphone array recordings",
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import micarray
+import soundfile as sf
+from scipy.special import sph_harm
+from scipy.signal import unit_impulse, sosfilt, fftconvolve as conv,\
+                         freqz, butter
+from micarray.util import db, tapering_window
+from micarray.modal.radial import crossover_frequencies, sos_radial_filter,\
+                                  tf_equalized_radial_filters
+
+c = 343
+fs = 44100
+Nfft = 4096
+Nfir = 1024
+
+array_order = 4
+azi_m, colat_m, R = np.loadtxt('em32.txt').T
+R = R[0]
+Ynm_m = micarray.modal.angular.sht_matrix(array_order, azi_m, colat_m)
+
+# Incident plane wave captured by mic array
+sf_order = 20
+Nimp = 2048
+imp = unit_impulse(Nimp)
+azi_pw, colat_pw = 0 * np.pi, 0.5 * np.pi
+delay, sos = sos_radial_filter(sf_order, R, fs=fs, setup='rigid')
+sos_irs = np.stack([sosfilt(sos[n], imp) for n in range(sf_order+1)])
+snm = np.column_stack([
+        sos_irs[n] * np.conj(sph_harm(m, n, azi_pw, colat_pw))
+        for n in range(sf_order+1)
+        for m in range(-n, n+1)])
+snm *= 4 * np.pi
+Ynm_s = micarray.modal.angular.sht_matrix(sf_order, azi_m, colat_m)
+s = np.real(np.squeeze(np.matmul(Ynm_s, snm[:, :, np.newaxis])))
+
+# Radial filters
+max_boost = 30
+f_xo = crossover_frequencies(array_order, R, max_boost)
+Nfft = 2048
+f_dft = np.fft.rfftfreq(Nfft, d=1/fs)
+f_dft[0] = 0.1 * f_dft[1]
+H_radial = tf_equalized_radial_filters(array_order, R, f_dft, max_boost,
+                                       type='butter')
+
+# Low-pass filtering
+b, a = butter(2, 20000/fs, btype='low')
+H_lpf = freqz(b, a, worN=Nfft//2+1)[1]
+H_radial *= H_lpf
+
+# Temporal windowing
+h_radial = np.stack([np.fft.irfft(Hi, n=Nfft, axis=-1) for Hi in H_radial])
+h_radial = np.roll(h_radial, int(Nfir/2), axis=-1)[:, :Nfir]
+h_radial *= tapering_window(Nfir, int(Nfir/4))
+pre_delay = -Nfir / 2 / fs
+
+# Save to wave file
+sf.write('em2-radial-filters.wav', h_radial.T, samplerate=fs, subtype='FLOAT')
+
+# beamforming
+bf_order = array_order
+N_angle = 360
+azi_l, colat_l = np.linspace(-np.pi, np.pi, num=N_angle), 0.5 * np.pi
+Ynm_l = micarray.modal.angular.sht_matrix(bf_order, azi_l, colat_l)
+snm_hat = np.squeeze(np.matmul(np.linalg.pinv(Ynm_m), s[:, :, np.newaxis]))
+ynm = np.column_stack([
+        conv(h_radial[n], snm_hat[:, n**2+n+m])
+        for n in range(bf_order+1)
+        for m in range(-n, n+1)])
+y = np.real(np.squeeze(np.matmul(Ynm_l, ynm[:, :, np.newaxis])))
+
+# frequency responses
+fmin, fmax, numf = 20, fs/2, 1000
+f = np.logspace(np.log10(fmin), np.log10(fmax), num=numf, endpoint=True)
+Y = np.column_stack([freqz(yi, 1, worN=f, fs=fs)[1] for yi in y.T])
+
+# critical frequencies
+f_alias = c * array_order / 2 / np.pi / R
+f_sf = c * sf_order / 2 / np.pi / R
+
+
+# plots
+def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg):
+    cax = plt.axes([ax.get_position().xmax + pad, ax.get_position().y0,
+                    width, ax.get_position().height], **kwarg)
+    plt.colorbar(im, cax=cax)
+
+
+im_kw = {'cmap': 'Blues', 'vmin': -60, 'vmax': None}
+phimin, phimax = np.rad2deg(azi_l[0]), np.rad2deg(azi_l[-1])
+phiticks = np.arange(phimin, phimax+90, 90)
+tmin = (delay + pre_delay) * 1000
+tmax = tmin + (Nfir + Nimp - 1)/fs * 1000
+tlim = -1.5, 1.5
+flim = fmin, fmax
+
+# Impulse responses
+fig, ax = plt.subplots(figsize=(4, 4))
+im = ax.imshow(db(y), extent=[phimin, phimax, tmin, tmax],
+               origin='lower', interpolation='none', **im_kw)
+ax.axis('tight')
+ax.set_ylim(tlim)
+ax.set_xticks(phiticks)
+ax.set_xlabel('azimuth in deg')
+ax.set_ylabel('time in ms')
+add_cbar(ax, im, xlabel='dB')
+plt.savefig('./em32-td.png', bbox_inches='tight')
+
+# Transfer functions
+fig, ax = plt.subplots(figsize=(4, 4))
+phi_deg = np.rad2deg(azi_l)
+im = ax.pcolormesh(phi_deg, f, db(Y), **im_kw)
+ax.plot(np.zeros_like(f_xo), f_xo, 'k+')
+[plt.text(0, fi, '{:0.1f} Hz'.format(fi), va='bottom', ha='left', rotation=30)
+ for fi in f_xo]
+ax.hlines(f_alias, phimin, phimax, color='r', linestyle='--')
+ax.hlines(f_sf, phimin, phimax, color='k', linestyle='--')
+ax.axis('tight')
+ax.set_xticks(phiticks)
+ax.set_ylim(flim)
+ax.set_yscale('log')
+ax.set_xlabel('azimuth in deg')
+ax.set_ylabel('frequency in Hz')
+add_cbar(ax, im, xlabel='dB')
+plt.savefig('./em32-fd.png', bbox_inches='tight')