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..13999aa
--- /dev/null
+++ b/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')
diff --git a/examples/linph-filterbank-butterworth-td.py b/examples/linph-filterbank-butterworth-td.py
new file mode 100644
index 0000000..b82e8e5
--- /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.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')
diff --git a/examples/linph-filterbank-crossover-frequencies.py b/examples/linph-filterbank-crossover-frequencies.py
new file mode 100644
index 0000000..b030c7a
--- /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} 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')
diff --git a/examples/linph-filterbank-modal-beamforming-em32.py b/examples/linph-filterbank-modal-beamforming-em32.py
new file mode 100644
index 0000000..e6730cd
--- /dev/null
+++ b/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')
diff --git a/examples/linph-filterbank-modal-beamforming-fd.py b/examples/linph-filterbank-modal-beamforming-fd.py
new file mode 100644
index 0000000..387965f
--- /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, type='butter')
+
+# 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 in deg')
+        if axi.is_first_col():
+            axi.set_ylabel('frequency in 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 in deg',
+                    ylabel='frequency in 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..98d51d6
--- /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, type='butter')
+
+# 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 in deg')
+        if axi.is_first_col():
+            axi.set_ylabel('time in 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 in deg', ylabel='time in 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..6402f3b
--- /dev/null
+++ b/examples/linph-filterbank-modal-beamforming-zylia.py
@@ -0,0 +1,134 @@
+"""Third-order Ambisonics microphone Zylia ZM-1.
+
+- filter bank design for Ambisonic encoding of ZM-1 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.signal import unit_impulse, sosfilt, fftconvolve as conv,\
+                         freqz, butter
+from scipy.special import sph_harm
+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 = 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,
+                                       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('zm1-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('./zylia-ZM1-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('./zylia-ZM1-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..e79a4cb
--- /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, type='butter')
+
+    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 in Hz')
+ax.set_ylabel('magnitude in dB')
+ax.grid(True)
+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
new file mode 100644
index 0000000..1a3cfe2
--- /dev/null
+++ b/examples/linph-filterbank-sidelobe-suppression.py
@@ -0,0 +1,47 @@
+"""Sidelobe suppression with different equalization schemes.
+
+- 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
+- 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
diff --git a/micarray/modal/radial.py b/micarray/modal/radial.py
index a2ad998..2eee1d3 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
@@ -59,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)
@@ -95,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
@@ -121,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
     ----------
@@ -236,7 +240,6 @@ def regularize(dn, a0, method):
     hn : array_like
 
     """
-
     idx = np.abs(dn) > a0
 
     if method == 'none':
@@ -280,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:
@@ -309,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))
@@ -339,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)
@@ -375,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
@@ -399,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
     ----------
@@ -455,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, :]
@@ -463,3 +472,244 @@ 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):
+    """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):
+        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 determined in such a way
+    that the maximum boost caused by each radial filter is limited.
+    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.
+
+    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 at which the transfer functions are evaluated.
+    type : {'butter', 'butter_equal_slopes'}
+        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)
+    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')
+                          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,
+                                type='butter'):
+    """Transfer 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.
+    type : {'butter', 'butter_maxorder'}
+        Type of filter responses.
+
+    """
+    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, type)
+    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..1ca4a84 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):
@@ -16,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)])
@@ -34,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
@@ -81,6 +84,7 @@ def matdiagmul(A, b):
     -------
     array_like
         Result of matrix multiplication.
+
     """
     if len(b.shape) == 1:
         b = b[np.newaxis, :]
@@ -93,7 +97,7 @@ def matdiagmul(A, b):
     return C
 
 
-def db(x, power=False):
+def db(x, *, power=False):
     """Convert *x* to decibel.
 
     Parameters
@@ -106,4 +110,71 @@ 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).
+
+    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 for modal weight type 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)
+
+
+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