From 54a1fd2eef01a7b5e4d44c025ebd5ca5f2fb2239 Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Thu, 12 Dec 2019 16:52:43 +0100 Subject: [PATCH 01/10] Linear-phase filterbank for radial filter design --- examples/em32.txt | 32 +++++ examples/linph-filterbank-butterworth-fd.py | 80 ++++++++++++ examples/linph-filterbank-butterworth-td.py | 89 +++++++++++++ .../linph-filterbank-crossover-frequencies.py | 59 +++++++++ ...linph-filterbank-modal-beamforming-em32.py | 119 ++++++++++++++++++ .../linph-filterbank-modal-beamforming-fd.py | 101 +++++++++++++++ .../linph-filterbank-modal-beamforming-td.py | 111 ++++++++++++++++ ...inph-filterbank-modal-beamforming-zylia.py | 119 ++++++++++++++++++ .../linph-filterbank-noise-amplification.py | 54 ++++++++ .../linph-filterbank-sidelobe-suppression.py | 47 +++++++ examples/zylia.txt | 19 +++ 11 files changed, 830 insertions(+) create mode 100644 examples/em32.txt create mode 100644 examples/linph-filterbank-butterworth-fd.py create mode 100644 examples/linph-filterbank-butterworth-td.py create mode 100644 examples/linph-filterbank-crossover-frequencies.py create mode 100644 examples/linph-filterbank-modal-beamforming-em32.py create mode 100644 examples/linph-filterbank-modal-beamforming-fd.py create mode 100644 examples/linph-filterbank-modal-beamforming-td.py create mode 100644 examples/linph-filterbank-modal-beamforming-zylia.py create mode 100644 examples/linph-filterbank-noise-amplification.py create mode 100644 examples/linph-filterbank-sidelobe-suppression.py create mode 100644 examples/zylia.txt diff --git a/examples/em32.txt b/examples/em32.txt new file mode 100644 index 0000000..32a4a59 --- /dev/null +++ b/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 diff --git a/examples/linph-filterbank-butterworth-fd.py b/examples/linph-filterbank-butterworth-fd.py new file mode 100644 index 0000000..9dfa116 --- /dev/null +++ b/examples/linph-filterbank-butterworth-fd.py @@ -0,0 +1,80 @@ +"""Linear-phase filterbank. + +- The target magnitude responses fo the filter-bank is designed + by using the zero-phase Butterworth responses. + (not to confused with typical 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) +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(n), alpha=0.5) +for fx in f_xo: + ax.semilogx(fx, 0, 'kv') + ax.text(fx, 0, '{:0.1f}'.format(fx), rotation=30, + horizontalalignment='left', verticalalignment='bottom') +ax.semilogx(f, db(H_tot), 'k:', label='Tot.') +ax.set_xlim(fmin, fmax) +ax.set_ylim(-100, 12) +ax.set_xscale('log') +ax.grid(True) +ax.set_xlabel('Frequency / Hz') +ax.set_ylabel('Magnitude / 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 / Hz') +ax.set_ylabel('Magnitude / dB') +ax.legend(title='order', loc='lower right', ncol=2) +plt.savefig('./linph-filterbank-butterworth-fd.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-butterworth-td.py b/examples/linph-filterbank-butterworth-td.py new file mode 100644 index 0000000..78fa10a --- /dev/null +++ b/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.0429 +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) + +# 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 / dB') + if ax.is_last_row(): + ax.set_xlabel('Time / 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') diff --git a/examples/linph-filterbank-crossover-frequencies.py b/examples/linph-filterbank-crossover-frequencies.py new file mode 100644 index 0000000..562629d --- /dev/null +++ b/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}'.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 / Hz') +ax.set_ylabel('Magnitude / dB') +ax.legend(title='Order', ncol=2) +plt.savefig('./crossover-frequencies.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-modal-beamforming-em32.py b/examples/linph-filterbank-modal-beamforming-em32.py new file mode 100644 index 0000000..192d4dd --- /dev/null +++ b/examples/linph-filterbank-modal-beamforming-em32.py @@ -0,0 +1,119 @@ +"""Fourth-order Ambisonics (Eigenmike em32). + + 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 +from scipy.special import sph_harm +from scipy.signal import unit_impulse, sosfilt, fftconvolve as conv,\ + kaiser, freqz +from micarray.util import db +from micarray.modal.radial import crossover_frequencies, sos_radial_filter,\ + tf_equalized_radial_filters + +c = 343 +fs = 44100 +Nfft = 2048 + +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) +h_radial = np.stack([np.fft.irfft(Hi, n=Nfft, axis=-1) for Hi in H_radial]) +h_radial = np.roll(h_radial, int(Nfft/2), axis=-1) +h_radial *= kaiser(Nfft, beta=8.6) +pre_delay = -Nfft / 2 / fs + +# 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 + (Nfft + 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 / deg') +ax.set_ylabel('Time / 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}'.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 / deg') +ax.set_ylabel('Frequency / Hz') +add_cbar(ax, im, xlabel='dB') +plt.savefig('./em32-fd.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-modal-beamforming-fd.py b/examples/linph-filterbank-modal-beamforming-fd.py new file mode 100644 index 0000000..c01b06d --- /dev/null +++ b/examples/linph-filterbank-modal-beamforming-fd.py @@ -0,0 +1,101 @@ +"""Angular-spectral responses. + +- No spatial sampling + + 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, point_spread, modal_norm, db +from micarray.modal.radial import crossover_frequencies, tf_linph_filterbank + +c = 343 +N = 5 +R = 0.049 +max_boost = 30 + +fmin, fmax, numf = 10, 20000, 500 +f = np.logspace(np.log10(fmin), np.log10(fmax), num=numf) +f_xo = crossover_frequencies(N, R, max_boost, modal_weight=maxre_sph) +H_fbank = tf_linph_filterbank(f_xo, f) + +# Look directions +azimin, azimax, numazi = -np.pi, np.pi, 361 +azi = np.linspace(azimin, azimax, num=numazi, endpoint=True) + +# Beamformer output +Y_band = np.zeros((N+1, numazi, numf), dtype='complex') +Y_order = np.zeros((N+1, numazi, numf), dtype='complex') +for i, Hi in enumerate(H_fbank): + ps = point_spread(i, azi, modal_weight=maxre_sph, equalization='diffuse') + Yi = ps[:, :, np.newaxis] * Hi + Y_order[:i+1, :] += Yi + Y_band[i, :] = np.sum(Yi, axis=0) +normN = modal_norm(maxre_sph(N)) +Y_band *= normN +Y_order *= normN +Y = np.sum(Y_order, axis=0) + + +# 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) + + +def decorate_singleplot(ax, **kwarg): + ax.axis('tight') + ax.set_yscale('log') + ax.set(**kwarg) + + +def decorate_subplots(axes, im, **kwarg): + for axi in axes.flat: + decorate_singleplot(axi, **kwarg) + if axi.is_last_row(): + axi.set_xlabel('Angle / deg') + if axi.is_first_col(): + axi.set_ylabel('Frequency / Hz') + add_cbar(axes[0, 1], im, pad=0.02, width=0.03, xlabel='dB') + + +figsize = (4, 4) +figsize_subplots = (8, 10) +azideg = np.rad2deg(azi) +azilim = azideg[0], azideg[-1] +phiticks = np.arange(np.rad2deg(azimin), np.rad2deg(azimax)+90, 90) +flim = fmin, fmax +im_kw = {'cmap': 'Blues', 'vmin': -60, 'vmax': 20} +gridspec_kw = {'wspace': 0.1} + +# Beamformer output +fig, ax = plt.subplots(figsize=figsize) +im = ax.pcolormesh(azideg, f, db(Y.T), **im_kw) +ax.plot(np.zeros(N), f_xo, 'wx', alpha=0.5) +decorate_singleplot(ax, xticks=phiticks, xlabel='Azimuth / deg', + ylabel='Frequnecy / Hz') +add_cbar(ax, im, xlabel='dB') +plt.savefig('spatial-responses-fd.png', bbox_inches='tight') + +# each filterbank +fig, ax = plt.subplots(figsize=figsize_subplots, ncols=2, nrows=3, + sharex=True, sharey=True, gridspec_kw=gridspec_kw) +for i, (axi, Yi) in enumerate(zip(ax.flat, Y_band)): + im = axi.pcolormesh(azideg, f, db(Yi.T), **im_kw) + axi.set_title('Subband #{}'.format(i)) +decorate_subplots(ax, im, xticks=phiticks, xlim=azilim, ylim=flim) +plt.savefig('spatial-responses-subband-fd.png', bbox_inches='tight') + +# each order +fig, ax = plt.subplots(figsize=figsize_subplots, ncols=2, nrows=3, + sharex=True, sharey=True, gridspec_kw=gridspec_kw) +for i, (axi, Yi) in enumerate(zip(ax.flat, Y_order)): + im = axi.pcolormesh(azideg, f, db(Yi.T), **im_kw) + axi.set_title('Order #{}'.format(i)) +decorate_subplots(ax, im, xticks=phiticks, xlim=azilim, ylim=flim) +plt.savefig('spatial-responses-order-fd.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-modal-beamforming-td.py b/examples/linph-filterbank-modal-beamforming-td.py new file mode 100644 index 0000000..a7057de --- /dev/null +++ b/examples/linph-filterbank-modal-beamforming-td.py @@ -0,0 +1,111 @@ +"""Angular-temporal responses. + +- No spatial sampling + + 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, point_spread, modal_norm, db +from micarray.modal.radial import crossover_frequencies, tf_linph_filterbank + +c = 343 +fs = 44100 +N = 5 +R = 0.049 + +Nfft = 2048 +fmin, fmax, numf = 0, fs/2, int(Nfft/2)+1 +f = np.linspace(fmin, fmax, num=numf, endpoint=True) + +max_boost = 30 +f_xo = crossover_frequencies(N, R, max_boost, modal_weight=maxre_sph) +H_fbank = tf_linph_filterbank(f_xo, f) + +# Look directions +azimin, azimax, numazi = -np.pi, np.pi, 360 +azi = np.linspace(azimin, azimax, num=numazi) + +# Beamformer output +Y_band = np.zeros((N+1, numazi, numf), dtype='complex') +Y_order = np.zeros((N+1, numazi, numf), dtype='complex') +for i, Hi in enumerate(H_fbank): + ps = point_spread(i, azi, modal_weight=maxre_sph, equalization='diffuse') + Yi = ps[:, :, np.newaxis] * Hi + Y_order[:i+1, :] += Yi + Y_band[i, :] = np.sum(Yi, axis=0) +normN = modal_norm(maxre_sph(N)) +Y_band *= normN +Y_order *= normN +Y = np.sum(Y_order, axis=0) + +y_order = np.fft.irfft(Y_order, n=Nfft, axis=-1) +y_order = np.roll(y_order, int(Nfft/2), axis=-1) +y_band = np.fft.irfft(Y_band, n=Nfft, axis=-1) +y_band = np.roll(y_band, int(Nfft/2), axis=-1) +y = np.fft.irfft(Y, n=Nfft, axis=-1) +y = np.roll(y, int(Nfft/2), axis=-1) + + +# 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) + + +def decorate_singleplot(ax, **kwarg): + ax.axis('tight') + ax.set(**kwarg) + + +def decorate_subplots(axes, im, **kwarg): + for axi in axes.flat: + decorate_singleplot(axi, **kwarg) + if axi.is_last_row(): + axi.set_xlabel('Angle / deg') + if axi.is_first_col(): + axi.set_ylabel('Time / ms') + add_cbar(axes[0, 1], im, pad=0.02, width=0.03, xlabel='dB') + + +figsize = (4, 4) +figsize_subplots = (8, 10) +azideg = np.rad2deg(azi) +azilim = azideg[0], azideg[-1] +aziticks = np.arange(np.rad2deg(azimin), np.rad2deg(azimax)+90, 90) +tau = Nfft / 2 / fs * 1000 +tlim = [-3, 3] +im_kw = {'cmap': 'Blues', 'vmin': -60, 'vmax': 20, + 'extent': [azilim[0], azilim[-1], -tau, tau]} +gridspec_kw = {'wspace': 0.1} + +# Beamformer impulse responses +fig, ax = plt.subplots(figsize=figsize) +im = ax.imshow(db(y.T), **im_kw) +decorate_singleplot(ax, xticks=aziticks, xlim=azilim, ylim=tlim, + xlabel='Azimuth / deg', ylabel='Time / ms') +add_cbar(ax, im, xlabel='dB') +plt.savefig('spatial-responses-td.png', bbox_inches='tight') + +# each filterbank +fig, ax = plt.subplots(figsize=figsize_subplots, ncols=2, nrows=3, + sharex=True, sharey=True, gridspec_kw=gridspec_kw) +for i, (axi, yi) in enumerate(zip(ax.flat, y_band)): + im = axi.imshow(db(yi.T), **im_kw) + axi.set_title('Subband #{}'.format(i)) +decorate_subplots(ax, im, xticks=aziticks, xlim=azilim, ylim=tlim) +plt.savefig('spatial-responses-subband-td.png', bbox_inches='tight') + +# each order +fig, ax = plt.subplots(figsize=figsize_subplots, ncols=2, nrows=3, + sharex=True, sharey=True, gridspec_kw=gridspec_kw) +for i, (axi, yi) in enumerate(zip(ax.flat, y_order)): + im = axi.imshow(db(yi.T), **im_kw) + axi.set_title('Order #{}'.format(i)) +decorate_subplots(ax, im, xticks=aziticks, xlim=azilim, ylim=tlim) +plt.savefig('spatial-responses-order-td.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-modal-beamforming-zylia.py b/examples/linph-filterbank-modal-beamforming-zylia.py new file mode 100644 index 0000000..7a8f98b --- /dev/null +++ b/examples/linph-filterbank-modal-beamforming-zylia.py @@ -0,0 +1,119 @@ +"""Third-order Ambisonics (Zylia). + + 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 +from scipy.signal import unit_impulse, sosfilt, fftconvolve as conv,\ + kaiser, freqz +from scipy.special import sph_harm +from micarray.util import db +from micarray.modal.radial import crossover_frequencies, sos_radial_filter,\ + tf_equalized_radial_filters + +c = 343 +fs = 44100 +Nfft = 2048 + +array_order = 3 +azi_m, colat_m, R = np.loadtxt('zylia.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) +h_radial = np.stack([np.fft.irfft(Hi, n=Nfft, axis=-1) for Hi in H_radial]) +h_radial = np.roll(h_radial, int(Nfft/2), axis=-1) +h_radial *= kaiser(Nfft, beta=8.6) +pre_delay = -Nfft / 2 / fs + +# 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 + (Nfft + 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 / deg') +ax.set_ylabel('Time / ms') +add_cbar(ax, im, xlabel='dB') +plt.savefig('./zylina-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}'.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 / deg') +ax.set_ylabel('Frequency / Hz') +add_cbar(ax, im, xlabel='dB') +plt.savefig('./zylina-fd.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-noise-amplification.py b/examples/linph-filterbank-noise-amplification.py new file mode 100644 index 0000000..1311dc3 --- /dev/null +++ b/examples/linph-filterbank-noise-amplification.py @@ -0,0 +1,54 @@ +"""Noise amplification. + + Reference + --------- + Franz Zotter, "A linear-phase filter-bank approach to process + rigid spherical microphone array recordings", in Proc. IcETRAN, + Palic, Serbia, 2018. (see Fig. 11) +""" +import numpy as np +import matplotlib.pyplot as plt +from micarray.util import maxre_sph, modal_norm, db +from micarray.modal.radial import spherical_hn2, tf_linph_filterbank,\ + crossover_frequencies + +c = 343 +R = 0.042 +N = 4 + +fmin, fmax, numf = 10, 20000, 500 +f = np.logspace(np.log10(fmin), np.log10(fmax), num=numf) +kr = 2 * np.pi * f / c * R + +Max_boost = 0, 5, 10, 15, 20 +Noise_amp = np.zeros((numf, len(Max_boost))) +Freq_xo = np.zeros((N, len(Max_boost))) + +# Prototype radial filters +H_proto = np.stack([1j**(-n-1) * (kr)**2 + * spherical_hn2(n, kr, derivative=True) + for n in range(N+1)]) + +for k, max_boost in enumerate(Max_boost): + f_xo = crossover_frequencies(N, R, max_boost) + Freq_xo[:, k] = f_xo + H_fbank = tf_linph_filterbank(f_xo, f) + + 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] + Noise_amp[:, k] = (modal_norm(np.abs(H_radial.T)) / np.abs(H_proto[0]))**2 + +# Plot +fig, ax = plt.subplots() +ax.semilogx(f, db(Noise_amp, power=True)) +ax.set_xlim(fmin, fmax) +ax.set_ylim(-3, 23) +ax.set_xlabel('Frequency / Hz') +ax.set_ylabel('Magnitude / dB') +ax.grid(True) +ax.legend(Max_boost, title='max. boost / dB') +plt.savefig('./noise-amplification.png') diff --git a/examples/linph-filterbank-sidelobe-suppression.py b/examples/linph-filterbank-sidelobe-suppression.py new file mode 100644 index 0000000..2121245 --- /dev/null +++ b/examples/linph-filterbank-sidelobe-suppression.py @@ -0,0 +1,47 @@ +"""Sidelobe suppression with different equalization schemes + +- The target magnitude responses fo the filter-bank is designed + by using the zero-phase Butterworth responses. + (not to confused with typical Butterworth filters) +- modal weight: 3D max-rE +- equalization + (1) omnidirectional + (2) diffuse-field + (3) free-field + + Reference + --------- + Franz Zotter, "A linear-phase filter-bank approach to process + rigid spherical microphone array recordings", in Proc. IcETRAN, + Palic, Serbia, 2018. (see Fig. 8) +""" +import numpy as np +import matplotlib.pyplot as plt +from micarray.util import db, point_spread + +N = 4 +azi = np.linspace(0, np.pi, num=360) +equalization = ['omni', 'diffuse', 'free'] + +# Plot +rticks = np.arange(-36, 12, 12) +rlim = -36, 3 +rticklabels = ['-36', '-24', '-12', '0 dB'] +nn = np.arange(N+1) + +fig, ax = plt.subplots(figsize=(10, 6), ncols=3, subplot_kw={'polar': True}, + gridspec_kw={'wspace': 0.3}) +for axi, eq in zip(ax, equalization): + ps = np.stack([np.sum(point_spread(n, azi, equalization=eq), axis=0) + for n in range(N+1)]) + ps *= 1 / np.max(ps) + axi.plot(azi[:, np.newaxis] * (-1)**nn, db(ps.T), lw=3, alpha=0.5) + axi.set_title('{}'.format(eq), y=1.2) + axi.set_rlabel_position(135) + axi.set_rticks(rticks) + axi.set_rlim(*rlim) + axi.set_yticklabels(rticklabels) +for axi in ax[1:]: + axi.set_xticklabels([]) +ax[-1].legend(nn, title='subband', bbox_to_anchor=(1.1, 1)) +plt.savefig('./sidelobe-suppression.png', bbox_inches='tight') diff --git a/examples/zylia.txt b/examples/zylia.txt new file mode 100644 index 0000000..4a734ff --- /dev/null +++ b/examples/zylia.txt @@ -0,0 +1,19 @@ +0.000000000000000000e+00 0.000000000000000000e+00 4.900000000000000189e-02 +3.058094442459276599e-03 7.305422893435145060e-01 4.900000000000000189e-02 +2.096009867533636051e+00 7.306700739627713936e-01 4.900000000000000189e-02 +-2.093360581925928443e+00 7.299094216727589624e-01 4.900000000000000189e-02 +-1.434099592396968825e+00 1.231829149238461429e+00 4.900000000000000189e-02 +-6.564873917134568249e-01 1.231643393484136872e+00 4.900000000000000189e-02 +6.612328142115841967e-01 1.231937671113323418e+00 4.900000000000000189e-02 +1.436243081415389033e+00 1.231737410884538697e+00 4.900000000000000189e-02 +2.755459326219778848e+00 1.231628696190499195e+00 4.900000000000000189e-02 +-2.750632294631810915e+00 1.231514727261006081e+00 4.900000000000000189e-02 +-2.480359839378209141e+00 1.909654982476469920e+00 4.900000000000000189e-02 +-1.705349572174404083e+00 1.909855242705254419e+00 4.900000000000000189e-02 +-3.861333273700143232e-01 1.909963957399293921e+00 4.900000000000000189e-02 +3.909603589579822014e-01 1.910077926328787035e+00 4.900000000000000189e-02 +1.707493061192824513e+00 1.909763504351331909e+00 4.900000000000000189e-02 +2.485105261876336513e+00 1.909949260105656466e+00 4.900000000000000189e-02 +-3.138534559147334146e+00 2.411050364246278832e+00 4.900000000000000189e-02 +-1.045582786056157065e+00 2.410922579627021722e+00 4.900000000000000189e-02 +1.048232071663864895e+00 2.411683231917034487e+00 4.900000000000000189e-02 From 0aedc37073c37cb6ccd72d4c998ac13295cbd418 Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Fri, 13 Dec 2019 13:10:33 +0100 Subject: [PATCH 02/10] Update radial.py and util.py --- micarray/modal/radial.py | 228 +++++++++++++++++++++++++++++++++++++++ micarray/util.py | 62 ++++++++++- 2 files changed, 288 insertions(+), 2 deletions(-) diff --git a/micarray/modal/radial.py b/micarray/modal/radial.py index a2ad998..e3f5596 100644 --- a/micarray/modal/radial.py +++ b/micarray/modal/radial.py @@ -1,6 +1,7 @@ from __future__ import division import numpy as np from scipy import special +from scipy.signal import bilinear_zpk, zpk2sos from .. import util from functools import wraps from warnings import warn @@ -463,3 +464,230 @@ def circ_diagonal_mode_mat(bk): for k in range(K): Bk[k, :, :] = np.diag(bk[k, :]) return np.squeeze(Bk) + + +def spherical_hn2(n, z, derivative=False): + """Spherical Hankel function of the sedond kind. + + n : int, array_like + Order of the spherical Hankel function (n >= 0). + z : comiplex or float, array_like + Argument of the spherical Hankel function. + derivative : bool, optional + If True, the value of the derivative (rather than the function + itself) is returned. + + """ + return special.spherical_jn(n, z, derivative)\ + - 1j * special.spherical_yn(n, z, derivative) + + +def sos_radial_filter(N, r, setup, c=343, fs=44100, pzmap='mz'): + """Radial filter design for a plane wave. + + Parameters + ---------- + N : int + Maximum order. + r : float + Radius of microphone array + setup : {'rigid'} + Array configuration (e.g. rigid) + pzmap : {'mz', 'bt'} + Pole-zero mapping method (matched-z, bilinear transform) + + Returns + ------- + delay : float + Overall delay + sos : list of (L, 6) arrays + Second-order section filters + """ + sos = [] + if setup is 'rigid': + for n in range(N + 1): + s0 = c / r * np.zeros(n) + sinf = c / r * np.roots(derivative_bessel_poly(n)[::-1]) + if pzmap is 'mz': + z0 = np.exp(s0 / fs) + zinf = np.exp(sinf / fs) + elif pzmap is 'bt': + z0, zinf, _ = bilinear_zpk(s0, pre_warp(sinf, fs), 1, fs=fs) + z0 = np.delete(z0, -1) + fc = c * n / 2 / np.pi / r + k = normalize_digital_filter_gain( + s0, sinf, z0, zinf, fc, fs) * c / r + sos.append(zpk2sos(z0, zinf, k, pairing='nearest')) + return -r / c, sos + + +def tf_butter(order, w, wc, btype='low'): + """Butterworth responses. + + Parameters + ---------- + order : int + Butterworth order. + w : array_like + Evaluation frequencies in Hertz. + wc : float + Cut-off frequency in Hertz. + btype : {'low', 'high'}, optional + Response type. + + """ + x = w / wc + if btype == 'low': + return 1 / (1 + x**order) + elif btype == 'high': + return x**order / (1 + x**order) + else: + raise ValueError("'btype' must be either: 'low' or 'high'") + + +def pre_warp(s, fs): + """Pre-warping frequency axis for bilinear transform.""" + return np.real(s) + 1j * 2 * fs * np.tan(np.imag(s) / 2 / fs) + + +def normalize_digital_filter_gain(s0, sinf, z0, zinf, fc, fs): + """Match the digital filter gain at a control frequency. + + Parameters + ---------- + s0 : (N,) array_like + zeros in the Laplace domain + sinf : (N,) array_like + polse in the Laplace domain + z0 : (N,) array_like + zeros in the z-domain + zinf : (N,) array_like + zeros in the z-domain + fc : float + Control frequency in Hz + fs : int + Sampling frequency in Hz + + """ + k = 1 + s_c = 1j * 2 * np.pi * fc + z_c = 1 / np.exp(1j * 2 * np.pi * fc / fs) + k *= np.prod(s_c - s0) / np.prod(s_c - sinf) + k *= np.prod(1 - zinf * z_c) / np.prod(1 - z0 * z_c) + return np.abs(k) + + +def bessel_poly(n): + """Bessel polynomial coefficients""" + beta = np.zeros(n + 1) # n-th order polynomial has (n+1) coeffcieints + beta[n] = 1 # This is always 1! + for k in range(n - 1, -1, -1): # Recurrence relation + beta[k] = beta[k + 1] * (2 * n - k) * (k + 1) / (n - k) / 2 + return beta + + +def derivative_bessel_poly(n): + """Derivative of the Bessel polynomial""" + gamma = bessel_poly(n + 1) + gamma[:-1] -= n * decrease_bessel_order_by_one(gamma) + return gamma + + +def decrease_bessel_order_by_one(beta): + n = len(beta)-1 + alpha = np.zeros(n) + for k in range(n - 1): + alpha[k] = beta[k + 1] * (k + 1) / (2 * n - k - 1) + alpha[-1] = 1 # This is always one + return alpha + + +def crossover_frequencies(N, r_array, max_boost, modal_weight=util.maxre_sph, + c=343): + """Crossover frequencies for filter-bank design. + + The crossover frequencies are detemined in such a way + that the maximum boost caused by each radial filter below 'max_boost'. + The small argument approximation of the spheircal Hankel function + normalized by the gain for each band is used for the computation. + The returned array has '(N-1)' frequencies. + + Parameters + ---------- + N : int + Maximum spherical harmonic order (Ambisonic order). + r_array : float + Radius of spherical microphone array in meter. + max_boost : float + Maximum allowable boost by radial filters in decibel. + modal_weight : callable, optional + Gain for individual spherical harmonic order (n). + c : float, optional + Speed of sound in meter per second. + + """ + g = 10**(max_boost / 20) + band_gain = [1 / util.modal_norm(modal_weight(n)) for n in range(N+1)] + kr = [np.power(util.double_factorial(2*n-1) * (n+1) / g / np.sqrt(2) + * band_gain[n], 1/n) + for n in range(1, N+1)] + return c / 2 / np.pi / r_array * np.array(kr) + + +def tf_linph_filterbank(f_xo, f, type='butter'): + """Linear-phase filterbank transfer functions. + + f_xo : array_like + Crossover frequencies in Hertz. + f : array_like + Frequencies in Hertz for which the transfer functions are evaluated. + type : {'butter'} + Type of filter responses. + + """ + N = len(f_xo) + omega = 2 * np.pi * f + omega_xo = 2 * np.pi * f_xo + if type == 'butter': + H_lpf = np.array([tf_butter(n+2, omega, omega_xo[n], btype='low') + for n in range(N)]) + H_hpf = np.array([tf_butter(n+2, omega, omega_xo[n], btype='high') + for n in range(N)]) + H_bpf = np.vstack([H_lpf[0], H_lpf[1:] * H_hpf[:-1], H_hpf[-1]]) + H_bpf /= np.sum(H_bpf, axis=0) + return H_bpf + else: + raise ValueError("Only 'type' = 'butter' is available.") + + +def tf_equalized_radial_filters(N, R, f, max_boost, + modal_weight=util.maxre_sph, c=343): + """Tranfer functions of equalized radial filters. + + N : int + Highest spherical harmonic order (Ambisonic order). + R : float + Radius of spherical microphone array in meter. + f : array_like + Frequencies in Hertz. + max_boost : float + Maximum allowable boost by radial filters in decibel. + modal_weight : callable, optional + Modal weighting function. + c : float, optional + Speed of sound in m/s. + + """ + kr = 2 * np.pi * f / c * R + f_xo = crossover_frequencies(N, R, max_boost, modal_weight) + H_fbank = tf_linph_filterbank(f_xo, f) + 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 = util.maxre_sph(i) + ai *= 1 / util.modal_norm(ai) + for n, Hn in enumerate(H_proto[:i+1]): + H_radial[n] += Hn * Hi * ai[n] + return 1 / 4 / np.pi * util.modal_norm(util.maxre_sph(N)) * H_radial diff --git a/micarray/util.py b/micarray/util.py index 08eb60b..8c45d6a 100644 --- a/micarray/util.py +++ b/micarray/util.py @@ -1,5 +1,6 @@ import numpy as np from scipy import linalg +from scipy.special import eval_legendre as legendre def norm_of_columns(A, p=2): @@ -93,7 +94,7 @@ def matdiagmul(A, b): return C -def db(x, power=False): +def db(x, *, power=False): """Convert *x* to decibel. Parameters @@ -106,4 +107,61 @@ def db(x, power=False): """ with np.errstate(divide='ignore'): - return 10 if power else 20 * np.log10(np.abs(x)) + return (10 if power else 20) * np.log10(np.abs(x)) + + +def double_factorial(n): + """Double factorial""" + if n == 0: + return 1 + elif n == 1: + return 1 + else: + return n * double_factorial(n - 2) + + +def maxre_sph(N): + """max-reE modal weight for spherical harmonics expansion. + + Parameter + --------- + N : int + Highest spherical harmonic order (Ambisonic order). + + """ + theta = np.deg2rad(137.9 / (N + 1.52)) + return legendre(np.arange(N + 1), np.cos(theta)) + + +def point_spread(N, phi, modal_weight=maxre_sph, equalization='omni'): + """Directional response of a given modal weight function and + equalization scheme. + + Parameters + ---------- + N : int + Highest spherical harmonic order (Ambisonic order). + phi : array_like + Angular distance from the main axis in radian. + modal_weight : callable, optional + Modal weighting function. + equalization : {'omni', 'diffuse', 'free'}, optional + Equalization scheme. + + """ + a = modal_weight(N) + if equalization == 'omni': + pass + elif equalization == 'diffuse': + a *= 1 / modal_norm(a, ord=2) + elif equalization == 'free': + a *= 1 / modal_norm(a, ord=1) + return np.stack([(2*n+1) * a[n] * legendre(n, np.cos(phi)) + for n in range(N+1)]) + + +def modal_norm(a, ord=2): + """Norm of the coefficients in the spherical harmonics domain. + """ + num_degree = 2 * np.arange(a.shape[-1]) + 1 + return np.sum(num_degree * a**ord, axis=-1)**(1/ord) From cdad7c1949a56a550fd1dec8ea2335eb3273af12 Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Sun, 15 Dec 2019 14:21:50 +0100 Subject: [PATCH 03/10] Update docstring --- micarray/modal/radial.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/micarray/modal/radial.py b/micarray/modal/radial.py index e3f5596..0e655fb 100644 --- a/micarray/modal/radial.py +++ b/micarray/modal/radial.py @@ -607,7 +607,7 @@ def crossover_frequencies(N, r_array, max_boost, modal_weight=util.maxre_sph, """Crossover frequencies for filter-bank design. The crossover frequencies are detemined in such a way - that the maximum boost caused by each radial filter below 'max_boost'. + that the maximum boost caused by each radial filter is limited. The small argument approximation of the spheircal Hankel function normalized by the gain for each band is used for the computation. The returned array has '(N-1)' frequencies. From d3b1a1cee9f1cdfc2d28e6677eff7133b673b6b2 Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Sun, 15 Dec 2019 14:56:54 +0100 Subject: [PATCH 04/10] Fix typo in line label --- examples/linph-filterbank-butterworth-fd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/linph-filterbank-butterworth-fd.py b/examples/linph-filterbank-butterworth-fd.py index 9dfa116..35679f7 100644 --- a/examples/linph-filterbank-butterworth-fd.py +++ b/examples/linph-filterbank-butterworth-fd.py @@ -48,7 +48,7 @@ # Filter-bank fig, ax = plt.subplots() for i, Hi in enumerate(H_fbank): - ax.semilogx(f, db(Hi), lw=3, label='${}$'.format(n), alpha=0.5) + 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}'.format(fx), rotation=30, From f8f8b5aef987adf75c51eb7f09951b28c2840db6 Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Sun, 15 Dec 2019 15:03:57 +0100 Subject: [PATCH 05/10] Add new option for filterbank design: constant filter order --- micarray/modal/radial.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/micarray/modal/radial.py b/micarray/modal/radial.py index 0e655fb..3789969 100644 --- a/micarray/modal/radial.py +++ b/micarray/modal/radial.py @@ -655,13 +655,21 @@ def tf_linph_filterbank(f_xo, f, type='butter'): for n in range(N)]) H_bpf = np.vstack([H_lpf[0], H_lpf[1:] * H_hpf[:-1], H_hpf[-1]]) H_bpf /= np.sum(H_bpf, axis=0) - return H_bpf + elif type == 'butter_maxorder': + H_lpf = np.array([tf_butter(N+2, omega, omega_xo[n], btype='low') + for n in range(N)]) + H_hpf = np.array([tf_butter(N+2, omega, omega_xo[n], btype='high') + for n in range(N)]) + H_bpf = np.vstack([H_lpf[0], H_lpf[1:] * H_hpf[:-1], H_hpf[-1]]) + H_bpf /= np.sum(H_bpf, axis=0) else: raise ValueError("Only 'type' = 'butter' is available.") + return H_bpf def tf_equalized_radial_filters(N, R, f, max_boost, - modal_weight=util.maxre_sph, c=343): + modal_weight=util.maxre_sph, c=343, + type='butter'): """Tranfer functions of equalized radial filters. N : int @@ -680,7 +688,7 @@ def tf_equalized_radial_filters(N, R, f, max_boost, """ kr = 2 * np.pi * f / c * R f_xo = crossover_frequencies(N, R, max_boost, modal_weight) - H_fbank = tf_linph_filterbank(f_xo, f) + H_fbank = tf_linph_filterbank(f_xo, f, type) H_proto = np.stack([1j**(-n-1) * (kr)**2 * spherical_hn2(n, kr, derivative=True) for n in range(N+1)]) From 4c85846c92a8569a4062eba1d56fa5efc52eec1f Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Sun, 15 Dec 2019 17:23:27 +0100 Subject: [PATCH 06/10] Add filterbank option --- examples/linph-filterbank-butterworth-fd.py | 2 +- examples/linph-filterbank-butterworth-td.py | 4 ++-- examples/linph-filterbank-modal-beamforming-em32.py | 3 ++- examples/linph-filterbank-modal-beamforming-fd.py | 2 +- examples/linph-filterbank-modal-beamforming-td.py | 2 +- examples/linph-filterbank-modal-beamforming-zylia.py | 3 ++- examples/linph-filterbank-noise-amplification.py | 2 +- 7 files changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/linph-filterbank-butterworth-fd.py b/examples/linph-filterbank-butterworth-fd.py index 35679f7..be91dba 100644 --- a/examples/linph-filterbank-butterworth-fd.py +++ b/examples/linph-filterbank-butterworth-fd.py @@ -26,7 +26,7 @@ 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) +H_fbank = tf_linph_filterbank(f_xo, f, type='butter') H_tot = np.sum(H_fbank, axis=0) # Prototpye radial filters diff --git a/examples/linph-filterbank-butterworth-td.py b/examples/linph-filterbank-butterworth-td.py index 78fa10a..845880d 100644 --- a/examples/linph-filterbank-butterworth-td.py +++ b/examples/linph-filterbank-butterworth-td.py @@ -17,7 +17,7 @@ fs = 44100 Nfft = 2048 -R = 0.0429 +R = 0.049 N = 5 max_boost = 30 f_xo = crossover_frequencies(N, R, max_boost, modal_weight=maxre_sph) @@ -26,7 +26,7 @@ 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) +H_fbank = tf_linph_filterbank(f_xo, f, type='butter') # Prototype radial filters H_proto = np.stack([1j**(-n-1) * (kr)**2 diff --git a/examples/linph-filterbank-modal-beamforming-em32.py b/examples/linph-filterbank-modal-beamforming-em32.py index 192d4dd..e568991 100644 --- a/examples/linph-filterbank-modal-beamforming-em32.py +++ b/examples/linph-filterbank-modal-beamforming-em32.py @@ -45,7 +45,8 @@ 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) +H_radial = tf_equalized_radial_filters(array_order, R, f_dft, max_boost, + type='butter') h_radial = np.stack([np.fft.irfft(Hi, n=Nfft, axis=-1) for Hi in H_radial]) h_radial = np.roll(h_radial, int(Nfft/2), axis=-1) h_radial *= kaiser(Nfft, beta=8.6) diff --git a/examples/linph-filterbank-modal-beamforming-fd.py b/examples/linph-filterbank-modal-beamforming-fd.py index c01b06d..c50bd69 100644 --- a/examples/linph-filterbank-modal-beamforming-fd.py +++ b/examples/linph-filterbank-modal-beamforming-fd.py @@ -21,7 +21,7 @@ fmin, fmax, numf = 10, 20000, 500 f = np.logspace(np.log10(fmin), np.log10(fmax), num=numf) f_xo = crossover_frequencies(N, R, max_boost, modal_weight=maxre_sph) -H_fbank = tf_linph_filterbank(f_xo, f) +H_fbank = tf_linph_filterbank(f_xo, f, type='butter') # Look directions azimin, azimax, numazi = -np.pi, np.pi, 361 diff --git a/examples/linph-filterbank-modal-beamforming-td.py b/examples/linph-filterbank-modal-beamforming-td.py index a7057de..510681b 100644 --- a/examples/linph-filterbank-modal-beamforming-td.py +++ b/examples/linph-filterbank-modal-beamforming-td.py @@ -24,7 +24,7 @@ max_boost = 30 f_xo = crossover_frequencies(N, R, max_boost, modal_weight=maxre_sph) -H_fbank = tf_linph_filterbank(f_xo, f) +H_fbank = tf_linph_filterbank(f_xo, f, type='butter') # Look directions azimin, azimax, numazi = -np.pi, np.pi, 360 diff --git a/examples/linph-filterbank-modal-beamforming-zylia.py b/examples/linph-filterbank-modal-beamforming-zylia.py index 7a8f98b..4d3e031 100644 --- a/examples/linph-filterbank-modal-beamforming-zylia.py +++ b/examples/linph-filterbank-modal-beamforming-zylia.py @@ -45,7 +45,8 @@ 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) +H_radial = tf_equalized_radial_filters(array_order, R, f_dft, max_boost, + type='butter') h_radial = np.stack([np.fft.irfft(Hi, n=Nfft, axis=-1) for Hi in H_radial]) h_radial = np.roll(h_radial, int(Nfft/2), axis=-1) h_radial *= kaiser(Nfft, beta=8.6) diff --git a/examples/linph-filterbank-noise-amplification.py b/examples/linph-filterbank-noise-amplification.py index 1311dc3..0efc63a 100644 --- a/examples/linph-filterbank-noise-amplification.py +++ b/examples/linph-filterbank-noise-amplification.py @@ -32,7 +32,7 @@ for k, max_boost in enumerate(Max_boost): f_xo = crossover_frequencies(N, R, max_boost) Freq_xo[:, k] = f_xo - H_fbank = tf_linph_filterbank(f_xo, f) + H_fbank = tf_linph_filterbank(f_xo, f, type='butter') H_radial = np.zeros_like(H_proto) for i, Hi in enumerate(H_fbank): From 6eb1145d51d5bed0da2437d4f98b2447e5ad3aa3 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 18 Dec 2019 21:44:30 +0100 Subject: [PATCH 07/10] axis labeling, typos I'd vote for 'f / Hz' or 'frequency in Hz', but not for 'frequency / Hz' for axis labeling I changed all to 'frequency / Hz' style --- examples/linph-filterbank-butterworth-fd.py | 14 +++++++------- examples/linph-filterbank-butterworth-td.py | 4 ++-- examples/linph-filterbank-crossover-frequencies.py | 6 +++--- .../linph-filterbank-modal-beamforming-em32.py | 10 +++++----- examples/linph-filterbank-modal-beamforming-fd.py | 8 ++++---- examples/linph-filterbank-modal-beamforming-td.py | 6 +++--- .../linph-filterbank-modal-beamforming-zylia.py | 14 +++++++------- examples/linph-filterbank-noise-amplification.py | 6 +++--- examples/linph-filterbank-sidelobe-suppression.py | 2 +- 9 files changed, 35 insertions(+), 35 deletions(-) diff --git a/examples/linph-filterbank-butterworth-fd.py b/examples/linph-filterbank-butterworth-fd.py index be91dba..3a8768d 100644 --- a/examples/linph-filterbank-butterworth-fd.py +++ b/examples/linph-filterbank-butterworth-fd.py @@ -2,7 +2,7 @@ - The target magnitude responses fo the filter-bank is designed by using the zero-phase Butterworth responses. - (not to confused with typical Butterworth filters) + (not to confused with typical (minphase) Butterworth filters) Reference --------- @@ -51,15 +51,15 @@ 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}'.format(fx), rotation=30, + ax.text(fx, 0, '{:0.1f} Hz'.format(fx), rotation=30, horizontalalignment='left', verticalalignment='bottom') -ax.semilogx(f, db(H_tot), 'k:', label='Tot.') +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 / Hz') -ax.set_ylabel('Magnitude / dB') +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') @@ -74,7 +74,7 @@ ax.set_ylim(-23, 33) ax.set_xscale('log') ax.grid(True) -ax.set_xlabel('Frequency / Hz') -ax.set_ylabel('Magnitude / dB') +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') diff --git a/examples/linph-filterbank-butterworth-td.py b/examples/linph-filterbank-butterworth-td.py index 845880d..b82e8e5 100644 --- a/examples/linph-filterbank-butterworth-td.py +++ b/examples/linph-filterbank-butterworth-td.py @@ -55,9 +55,9 @@ def decorate_subplots(axes, **kwargs): for ax in axes.flat: if ax.is_first_col(): - ax.set_ylabel('Amplitude / dB') + ax.set_ylabel('Amplitude in dB') if ax.is_last_row(): - ax.set_xlabel('Time / ms') + ax.set_xlabel('Time in ms') ax.grid(True) ax.set(**kwargs) diff --git a/examples/linph-filterbank-crossover-frequencies.py b/examples/linph-filterbank-crossover-frequencies.py index 562629d..b030c7a 100644 --- a/examples/linph-filterbank-crossover-frequencies.py +++ b/examples/linph-filterbank-crossover-frequencies.py @@ -46,14 +46,14 @@ 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}'.format(fn), + 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 / Hz') -ax.set_ylabel('Magnitude / dB') +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') diff --git a/examples/linph-filterbank-modal-beamforming-em32.py b/examples/linph-filterbank-modal-beamforming-em32.py index e568991..d8a7417 100644 --- a/examples/linph-filterbank-modal-beamforming-em32.py +++ b/examples/linph-filterbank-modal-beamforming-em32.py @@ -96,8 +96,8 @@ def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg): ax.axis('tight') ax.set_ylim(tlim) ax.set_xticks(phiticks) -ax.set_xlabel('Azimuth / deg') -ax.set_ylabel('Time / ms') +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') @@ -106,7 +106,7 @@ def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg): 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}'.format(fi), va='bottom', ha='left', rotation=30) +[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='--') @@ -114,7 +114,7 @@ def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg): ax.set_xticks(phiticks) ax.set_ylim(flim) ax.set_yscale('log') -ax.set_xlabel('Azimuth / deg') -ax.set_ylabel('Frequency / Hz') +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') diff --git a/examples/linph-filterbank-modal-beamforming-fd.py b/examples/linph-filterbank-modal-beamforming-fd.py index c50bd69..387965f 100644 --- a/examples/linph-filterbank-modal-beamforming-fd.py +++ b/examples/linph-filterbank-modal-beamforming-fd.py @@ -58,9 +58,9 @@ def decorate_subplots(axes, im, **kwarg): for axi in axes.flat: decorate_singleplot(axi, **kwarg) if axi.is_last_row(): - axi.set_xlabel('Angle / deg') + axi.set_xlabel('angle in deg') if axi.is_first_col(): - axi.set_ylabel('Frequency / Hz') + axi.set_ylabel('frequency in Hz') add_cbar(axes[0, 1], im, pad=0.02, width=0.03, xlabel='dB') @@ -77,8 +77,8 @@ def decorate_subplots(axes, im, **kwarg): fig, ax = plt.subplots(figsize=figsize) im = ax.pcolormesh(azideg, f, db(Y.T), **im_kw) ax.plot(np.zeros(N), f_xo, 'wx', alpha=0.5) -decorate_singleplot(ax, xticks=phiticks, xlabel='Azimuth / deg', - ylabel='Frequnecy / Hz') +decorate_singleplot(ax, xticks=phiticks, xlabel='azimuth in deg', + ylabel='frequency in Hz') add_cbar(ax, im, xlabel='dB') plt.savefig('spatial-responses-fd.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-modal-beamforming-td.py b/examples/linph-filterbank-modal-beamforming-td.py index 510681b..98d51d6 100644 --- a/examples/linph-filterbank-modal-beamforming-td.py +++ b/examples/linph-filterbank-modal-beamforming-td.py @@ -67,9 +67,9 @@ def decorate_subplots(axes, im, **kwarg): for axi in axes.flat: decorate_singleplot(axi, **kwarg) if axi.is_last_row(): - axi.set_xlabel('Angle / deg') + axi.set_xlabel('angle in deg') if axi.is_first_col(): - axi.set_ylabel('Time / ms') + axi.set_ylabel('time in ms') add_cbar(axes[0, 1], im, pad=0.02, width=0.03, xlabel='dB') @@ -88,7 +88,7 @@ def decorate_subplots(axes, im, **kwarg): fig, ax = plt.subplots(figsize=figsize) im = ax.imshow(db(y.T), **im_kw) decorate_singleplot(ax, xticks=aziticks, xlim=azilim, ylim=tlim, - xlabel='Azimuth / deg', ylabel='Time / ms') + xlabel='azimuth in deg', ylabel='time in ms') add_cbar(ax, im, xlabel='dB') plt.savefig('spatial-responses-td.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-modal-beamforming-zylia.py b/examples/linph-filterbank-modal-beamforming-zylia.py index 4d3e031..22f1e5f 100644 --- a/examples/linph-filterbank-modal-beamforming-zylia.py +++ b/examples/linph-filterbank-modal-beamforming-zylia.py @@ -96,17 +96,17 @@ def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg): ax.axis('tight') ax.set_ylim(tlim) ax.set_xticks(phiticks) -ax.set_xlabel('Azimuth / deg') -ax.set_ylabel('Time / ms') +ax.set_xlabel('azimuth in deg') +ax.set_ylabel('time in ms') add_cbar(ax, im, xlabel='dB') -plt.savefig('./zylina-td.png', bbox_inches='tight') +plt.savefig('./zylia-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}'.format(fi), va='bottom', ha='left', rotation=30) +[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='--') @@ -114,7 +114,7 @@ def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg): ax.set_xticks(phiticks) ax.set_ylim(flim) ax.set_yscale('log') -ax.set_xlabel('Azimuth / deg') -ax.set_ylabel('Frequency / Hz') +ax.set_xlabel('azimuth in deg') +ax.set_ylabel('frequency in Hz') add_cbar(ax, im, xlabel='dB') -plt.savefig('./zylina-fd.png', bbox_inches='tight') +plt.savefig('./zylia-fd.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-noise-amplification.py b/examples/linph-filterbank-noise-amplification.py index 0efc63a..e79a4cb 100644 --- a/examples/linph-filterbank-noise-amplification.py +++ b/examples/linph-filterbank-noise-amplification.py @@ -47,8 +47,8 @@ ax.semilogx(f, db(Noise_amp, power=True)) ax.set_xlim(fmin, fmax) ax.set_ylim(-3, 23) -ax.set_xlabel('Frequency / Hz') -ax.set_ylabel('Magnitude / dB') +ax.set_xlabel('frequency in Hz') +ax.set_ylabel('magnitude in dB') ax.grid(True) -ax.legend(Max_boost, title='max. boost / dB') +ax.legend(Max_boost, title='max. boost in dB') plt.savefig('./noise-amplification.png') diff --git a/examples/linph-filterbank-sidelobe-suppression.py b/examples/linph-filterbank-sidelobe-suppression.py index 2121245..6ab7d8f 100644 --- a/examples/linph-filterbank-sidelobe-suppression.py +++ b/examples/linph-filterbank-sidelobe-suppression.py @@ -1,4 +1,4 @@ -"""Sidelobe suppression with different equalization schemes +"""Sidelobe suppression with different equalization schemes. - The target magnitude responses fo the filter-bank is designed by using the zero-phase Butterworth responses. From 5d374ba4bdb7747b9320eb2fbd252824ef054195 Mon Sep 17 00:00:00 2001 From: Frank Schultz Date: Wed, 18 Dec 2019 23:25:32 +0100 Subject: [PATCH 08/10] PEP style mods, improved comments --- examples/linph-filterbank-butterworth-fd.py | 4 +- ...linph-filterbank-modal-beamforming-em32.py | 4 +- ...inph-filterbank-modal-beamforming-zylia.py | 8 ++-- .../linph-filterbank-sidelobe-suppression.py | 2 +- micarray/modal/radial.py | 40 +++++++++++++------ micarray/util.py | 17 +++++--- 6 files changed, 49 insertions(+), 26 deletions(-) diff --git a/examples/linph-filterbank-butterworth-fd.py b/examples/linph-filterbank-butterworth-fd.py index 3a8768d..13999aa 100644 --- a/examples/linph-filterbank-butterworth-fd.py +++ b/examples/linph-filterbank-butterworth-fd.py @@ -1,7 +1,7 @@ """Linear-phase filterbank. -- The target magnitude responses fo the filter-bank is designed - by using the zero-phase Butterworth responses. +- 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 diff --git a/examples/linph-filterbank-modal-beamforming-em32.py b/examples/linph-filterbank-modal-beamforming-em32.py index d8a7417..8880bd0 100644 --- a/examples/linph-filterbank-modal-beamforming-em32.py +++ b/examples/linph-filterbank-modal-beamforming-em32.py @@ -1,4 +1,6 @@ -"""Fourth-order Ambisonics (Eigenmike em32). +"""Fourth-order Ambisonics micropphone Eigenmike em32. + +- filter bank design for Ambisonic encoding of em32 signals Reference --------- diff --git a/examples/linph-filterbank-modal-beamforming-zylia.py b/examples/linph-filterbank-modal-beamforming-zylia.py index 22f1e5f..6fdb43e 100644 --- a/examples/linph-filterbank-modal-beamforming-zylia.py +++ b/examples/linph-filterbank-modal-beamforming-zylia.py @@ -1,4 +1,6 @@ -"""Third-order Ambisonics (Zylia). +"""Third-order Ambisonics microphone Zylia ZM-1. + +- filter bank design for Ambisonic encoding of ZM-1 signals Reference --------- @@ -99,7 +101,7 @@ def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg): ax.set_xlabel('azimuth in deg') ax.set_ylabel('time in ms') add_cbar(ax, im, xlabel='dB') -plt.savefig('./zylia-td.png', bbox_inches='tight') +plt.savefig('./zylia-ZM1-td.png', bbox_inches='tight') # Transfer functions fig, ax = plt.subplots(figsize=(4, 4)) @@ -117,4 +119,4 @@ def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg): ax.set_xlabel('azimuth in deg') ax.set_ylabel('frequency in Hz') add_cbar(ax, im, xlabel='dB') -plt.savefig('./zylia-fd.png', bbox_inches='tight') +plt.savefig('./zylia-ZM1-fd.png', bbox_inches='tight') diff --git a/examples/linph-filterbank-sidelobe-suppression.py b/examples/linph-filterbank-sidelobe-suppression.py index 6ab7d8f..1a3cfe2 100644 --- a/examples/linph-filterbank-sidelobe-suppression.py +++ b/examples/linph-filterbank-sidelobe-suppression.py @@ -1,6 +1,6 @@ """Sidelobe suppression with different equalization schemes. -- The target magnitude responses fo the filter-bank is designed +- The target magnitude responses for the filter-bank is designed by using the zero-phase Butterworth responses. (not to confused with typical Butterworth filters) - modal weight: 3D max-rE diff --git a/micarray/modal/radial.py b/micarray/modal/radial.py index 3789969..2eee1d3 100644 --- a/micarray/modal/radial.py +++ b/micarray/modal/radial.py @@ -60,6 +60,7 @@ def spherical_pw(N, k, r, setup): ------- bn : (M, N+1) numpy.ndarray Radial weights for all orders up to N and the given wavenumbers. + """ kr = util.asarray_1d(k*r) n = np.arange(N+1) @@ -96,6 +97,7 @@ def spherical_ps(N, k, r, rs, setup): ------- bn : (M, N+1) numpy.ndarray Radial weights for all orders up to N and the given wavenumbers. + """ k = util.asarray_1d(k) krs = k*rs @@ -122,7 +124,8 @@ def weights(N, kr, setup): .. math:: - b_n(kr) = j_n(kr) - \frac{j_n^\prime(kr)}{h_n^{(2)\prime}(kr)}h_n^{(2)}(kr) + b_n(kr) = + j_n(kr) - \frac{j_n^\prime(kr)}{h_n^{(2)\prime}(kr)}h_n^{(2)}(kr) Parameters ---------- @@ -237,7 +240,6 @@ def regularize(dn, a0, method): hn : array_like """ - idx = np.abs(dn) > a0 if method == 'none': @@ -281,6 +283,7 @@ def diagonal_mode_mat(bk): Bk : (M, (N+1)**2, (N+1)**2) numpy.ndarray Multidimensional array containing diagnonal matrices with input vector on main diagonal. + """ bk = repeat_n_m(bk) if len(bk.shape) == 1: @@ -310,6 +313,7 @@ def repeat_n_m(v): ------- : (,(N+1)**2) numpy.ndarray Vector or stack of vectors containing repetated values. + """ krlist = [np.tile(v, (2*i+1, 1)).T for i, v in enumerate(v.T.tolist())] return np.squeeze(np.concatenate(krlist, axis=-1)) @@ -340,7 +344,8 @@ def circular_pw(N, k, r, setup): Returns ------- bn : (M, 2*N+1) numpy.ndarray - Radial weights for all orders up to N and the given wavenumbers. + Radial weights for all orders up to N and the given wavenumbers + """ kr = util.asarray_1d(k*r) n = np.roll(np.arange(-N, N+1), -N) @@ -376,7 +381,8 @@ def circular_ls(N, k, r, rs, setup): Returns ------- bn : (M, 2*N+1) numpy.ndarray - Radial weights for all orders up to N and the given wavenumbers. + Radial weights for all orders up to N and the given wavenumbers + """ k = util.asarray_1d(k) krs = k*rs @@ -400,7 +406,8 @@ def circ_radial_weights(N, kr, setup): .. math:: - b_n(kr) = J_n(kr) - \frac{J_n^\prime(kr)}{H_n^{(2)\prime}(kr)}H_n^{(2)}(kr) + b_n(kr) = + J_n(kr) - \frac{J_n^\prime(kr)}{H_n^{(2)\prime}(kr)}H_n^{(2)}(kr) Parameters ---------- @@ -456,6 +463,7 @@ def circ_diagonal_mode_mat(bk): Bk : (M, 2*N+1, 2*N+1) numpy.ndarray Multidimensional array containing diagnonal matrices with input vector on main diagonal. + """ if len(bk.shape) == 1: bk = bk[np.newaxis, :] @@ -502,6 +510,7 @@ def sos_radial_filter(N, r, setup, c=343, fs=44100, pzmap='mz'): Overall delay sos : list of (L, 6) arrays Second-order section filters + """ sos = [] if setup is 'rigid': @@ -578,7 +587,7 @@ def normalize_digital_filter_gain(s0, sinf, z0, zinf, fc, fs): def bessel_poly(n): - """Bessel polynomial coefficients""" + """Bessel polynomial coefficients.""" beta = np.zeros(n + 1) # n-th order polynomial has (n+1) coeffcieints beta[n] = 1 # This is always 1! for k in range(n - 1, -1, -1): # Recurrence relation @@ -587,13 +596,14 @@ def bessel_poly(n): def derivative_bessel_poly(n): - """Derivative of the Bessel polynomial""" + """Bessel polynomial derivative.""" gamma = bessel_poly(n + 1) gamma[:-1] -= n * decrease_bessel_order_by_one(gamma) return gamma def decrease_bessel_order_by_one(beta): + """Bessel polynomial decrease order.""" n = len(beta)-1 alpha = np.zeros(n) for k in range(n - 1): @@ -606,9 +616,9 @@ def crossover_frequencies(N, r_array, max_boost, modal_weight=util.maxre_sph, c=343): """Crossover frequencies for filter-bank design. - The crossover frequencies are detemined in such a way + The crossover frequencies are determined in such a way that the maximum boost caused by each radial filter is limited. - The small argument approximation of the spheircal Hankel function + The small argument approximation of the spherical Hankel function normalized by the gain for each band is used for the computation. The returned array has '(N-1)' frequencies. @@ -640,8 +650,8 @@ def tf_linph_filterbank(f_xo, f, type='butter'): f_xo : array_like Crossover frequencies in Hertz. f : array_like - Frequencies in Hertz for which the transfer functions are evaluated. - type : {'butter'} + Frequencies in Hertz at which the transfer functions are evaluated. + type : {'butter', 'butter_equal_slopes'} Type of filter responses. """ @@ -655,7 +665,9 @@ def tf_linph_filterbank(f_xo, f, type='butter'): for n in range(N)]) H_bpf = np.vstack([H_lpf[0], H_lpf[1:] * H_hpf[:-1], H_hpf[-1]]) H_bpf /= np.sum(H_bpf, axis=0) - elif type == 'butter_maxorder': + elif type == 'butter_equal_slopes': + # special case for development purpose + # all bands exhibit equal filter slopes of order N+2 H_lpf = np.array([tf_butter(N+2, omega, omega_xo[n], btype='low') for n in range(N)]) H_hpf = np.array([tf_butter(N+2, omega, omega_xo[n], btype='high') @@ -670,7 +682,7 @@ def tf_linph_filterbank(f_xo, f, type='butter'): def tf_equalized_radial_filters(N, R, f, max_boost, modal_weight=util.maxre_sph, c=343, type='butter'): - """Tranfer functions of equalized radial filters. + """Transfer functions of equalized radial filters. N : int Highest spherical harmonic order (Ambisonic order). @@ -684,6 +696,8 @@ def tf_equalized_radial_filters(N, R, f, max_boost, Modal weighting function. c : float, optional Speed of sound in m/s. + type : {'butter', 'butter_maxorder'} + Type of filter responses. """ kr = 2 * np.pi * f / c * R diff --git a/micarray/util.py b/micarray/util.py index 8c45d6a..7b28af3 100644 --- a/micarray/util.py +++ b/micarray/util.py @@ -17,6 +17,7 @@ def norm_of_columns(A, p=2): ------- array_like p-norm of each column of A. + """ _, N = A.shape return np.asarray([linalg.norm(A[:, j], ord=p) for j in range(N)]) @@ -35,7 +36,8 @@ def coherence_of_columns(A): Returns ------- array_like - Mutual coherence of columns of A. + Mutual coherence of columns of A + """ A = np.asmatrix(A) _, N = A.shape @@ -82,6 +84,7 @@ def matdiagmul(A, b): ------- array_like Result of matrix multiplication. + """ if len(b.shape) == 1: b = b[np.newaxis, :] @@ -111,7 +114,7 @@ def db(x, *, power=False): def double_factorial(n): - """Double factorial""" + """Double factorial.""" if n == 0: return 1 elif n == 1: @@ -128,14 +131,17 @@ def maxre_sph(N): N : int Highest spherical harmonic order (Ambisonic order). + see + F. Zotter, M. Frank: "All-Round Ambisonic Panning and Decoding." + J. Audio Eng. Soc. 60(10):207-820, October 2012 + """ theta = np.deg2rad(137.9 / (N + 1.52)) return legendre(np.arange(N + 1), np.cos(theta)) def point_spread(N, phi, modal_weight=maxre_sph, equalization='omni'): - """Directional response of a given modal weight function and - equalization scheme. + """Directional response for modal weight type and equalization scheme. Parameters ---------- @@ -161,7 +167,6 @@ def point_spread(N, phi, modal_weight=maxre_sph, equalization='omni'): def modal_norm(a, ord=2): - """Norm of the coefficients in the spherical harmonics domain. - """ + """Norm of the coefficients in the spherical harmonics domain.""" num_degree = 2 * np.arange(a.shape[-1]) + 1 return np.sum(num_degree * a**ord, axis=-1)**(1/ord) From 64a9d5822c5c578f337e975a1ee596ce6b77a68b Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Thu, 26 Dec 2019 18:43:22 +0100 Subject: [PATCH 09/10] fix typo --- examples/linph-filterbank-modal-beamforming-em32.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/linph-filterbank-modal-beamforming-em32.py b/examples/linph-filterbank-modal-beamforming-em32.py index 8880bd0..a66b0d2 100644 --- a/examples/linph-filterbank-modal-beamforming-em32.py +++ b/examples/linph-filterbank-modal-beamforming-em32.py @@ -1,4 +1,4 @@ -"""Fourth-order Ambisonics micropphone Eigenmike em32. +"""Fourth-order Ambisonics microphone Eigenmike em32. - filter bank design for Ambisonic encoding of em32 signals From 81efa0b9018fe27a10f23dc6f27d660dc0927032 Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Tue, 7 Jan 2020 17:15:00 +0100 Subject: [PATCH 10/10] 20kHz low-pass filtering; temporal windowing --- ...linph-filterbank-modal-beamforming-em32.py | 26 ++++++++++++++----- ...inph-filterbank-modal-beamforming-zylia.py | 26 ++++++++++++++----- micarray/util.py | 8 ++++++ 3 files changed, 46 insertions(+), 14 deletions(-) diff --git a/examples/linph-filterbank-modal-beamforming-em32.py b/examples/linph-filterbank-modal-beamforming-em32.py index a66b0d2..e6730cd 100644 --- a/examples/linph-filterbank-modal-beamforming-em32.py +++ b/examples/linph-filterbank-modal-beamforming-em32.py @@ -10,16 +10,18 @@ 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,\ - kaiser, freqz -from micarray.util import db + 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 = 2048 +Nfft = 4096 +Nfir = 1024 array_order = 4 azi_m, colat_m, R = np.loadtxt('em32.txt').T @@ -49,10 +51,20 @@ 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(Nfft/2), axis=-1) -h_radial *= kaiser(Nfft, beta=8.6) -pre_delay = -Nfft / 2 / fs +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 @@ -87,7 +99,7 @@ def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg): 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 + (Nfft + Nimp - 1)/fs * 1000 +tmax = tmin + (Nfir + Nimp - 1)/fs * 1000 tlim = -1.5, 1.5 flim = fmin, fmax diff --git a/examples/linph-filterbank-modal-beamforming-zylia.py b/examples/linph-filterbank-modal-beamforming-zylia.py index 6fdb43e..6402f3b 100644 --- a/examples/linph-filterbank-modal-beamforming-zylia.py +++ b/examples/linph-filterbank-modal-beamforming-zylia.py @@ -10,16 +10,18 @@ import numpy as np import matplotlib.pyplot as plt import micarray +import soundfile as sf from scipy.signal import unit_impulse, sosfilt, fftconvolve as conv,\ - kaiser, freqz + freqz, butter from scipy.special import sph_harm -from micarray.util import db +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 = 2048 +Nfft = 4096 +Nfir = 1024 array_order = 3 azi_m, colat_m, R = np.loadtxt('zylia.txt').T @@ -49,10 +51,20 @@ 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(Nfft/2), axis=-1) -h_radial *= kaiser(Nfft, beta=8.6) -pre_delay = -Nfft / 2 / fs +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('zm1-radial-filters.wav', h_radial.T, samplerate=fs, subtype='FLOAT') # beamforming bf_order = array_order @@ -87,7 +99,7 @@ def add_cbar(ax, im, pad=0.05, width=0.05, **kwarg): 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 + (Nfft + Nimp - 1)/fs * 1000 +tmax = tmin + (Nfir + Nimp - 1)/fs * 1000 tlim = -1.5, 1.5 flim = fmin, fmax diff --git a/micarray/util.py b/micarray/util.py index 7b28af3..1ca4a84 100644 --- a/micarray/util.py +++ b/micarray/util.py @@ -170,3 +170,11 @@ def modal_norm(a, ord=2): """Norm of the coefficients in the spherical harmonics domain.""" num_degree = 2 * np.arange(a.shape[-1]) + 1 return np.sum(num_degree * a**ord, axis=-1)**(1/ord) + + +def tapering_window(Nfull, Nedge): + w = np.ones(Nfull) + alpha = np.linspace(0, np.pi, Nedge, endpoint=False) + w[:Nedge] = 0.5 + 0.5 * np.cos(alpha[::-1]) + w[-Nedge:] = 0.5 + 0.5 * np.cos(alpha) + return w