|
| 1 | +""" |
| 2 | +# Octave Band Filter Banks |
| 3 | +
|
| 4 | +Multi-frequency simulation relies on octave-band filterbanks |
| 5 | +to run the simulation in different perceptually relevant frequency bands |
| 6 | +before merging the results into a single room impulse response. |
| 7 | +
|
| 8 | +This scripts demonstrate some of the octave band filters available |
| 9 | +in pyroomacoustics. |
| 10 | +
|
| 11 | +Two ocatave band filter banks are implemented in pyroomacoustics. |
| 12 | +
|
| 13 | +## Cosine Filterbank |
| 14 | +
|
| 15 | +This filterbank uses a number of overlapping cosine filters to |
| 16 | +cover the octaves. |
| 17 | +It guarantees perfect reconstruction, but does not conserve the energy |
| 18 | +in the bands. |
| 19 | +
|
| 20 | +## Antoni's Orthogonal-like Fractional Octave Bands |
| 21 | +
|
| 22 | +This class implements a type of fractional octave filter bank with |
| 23 | +both perfect reconstruction and energy conservation. |
| 24 | +
|
| 25 | +J. Antoni, Orthogonal-like fractional-octave-band filters, J. Acoust. Soc. |
| 26 | +Am., 127, 2, February 2010 |
| 27 | +""" |
| 28 | + |
| 29 | +import matplotlib.pyplot as plt |
| 30 | +import numpy as np |
| 31 | +from scipy.signal import chirp |
| 32 | + |
| 33 | +import pyroomacoustics as pra |
| 34 | + |
| 35 | +if __name__ == "__main__": |
| 36 | + |
| 37 | + # Test unit energy. |
| 38 | + fs = 16000 |
| 39 | + n_fft = 2**10 # 1024 |
| 40 | + base_freq = 125.0 # Hertz, default frequency in pyroomacoustics. |
| 41 | + |
| 42 | + # The cosine filter bank |
| 43 | + octave_bands = pra.OctaveBandsFactory( |
| 44 | + fs=fs, |
| 45 | + n_fft=n_fft, |
| 46 | + keep_dc=True, |
| 47 | + base_frequency=base_freq, |
| 48 | + ) |
| 49 | + |
| 50 | + # The orthogonal-like filterbankd with perfect reconstruction and energy |
| 51 | + # conservation. |
| 52 | + # The `band_overlap_ratio` and `slope` parameters control the transition |
| 53 | + # between adjacent bands. |
| 54 | + antoni_octave_bands = pra.AntoniOctaveFilterBank( |
| 55 | + fs=fs, |
| 56 | + base_frequency=base_freq, |
| 57 | + band_overlap_ratio=0.5, |
| 58 | + slope=0, |
| 59 | + n_fft=n_fft, |
| 60 | + third=False, |
| 61 | + ) |
| 62 | + |
| 63 | + fig, axes = plt.subplots(2, 2, figsize=(10, 5)) |
| 64 | + |
| 65 | + filters_orth = np.zeros(n_fft) |
| 66 | + filters_orth[n_fft // 2] = 1.0 |
| 67 | + filters_orth = antoni_octave_bands.analysis(filters_orth) |
| 68 | + |
| 69 | + for i, (lbl, ob) in enumerate( |
| 70 | + {"Cosine": octave_bands.filters, "Antoni": filters_orth}.items() |
| 71 | + ): |
| 72 | + filters_spectrum = np.fft.rfft(ob, axis=0) |
| 73 | + n_freq = filters_spectrum.shape[0] |
| 74 | + freq = np.arange(n_freq) / n_freq * (fs / 2) |
| 75 | + time = np.arange(ob.shape[0]) / fs |
| 76 | + |
| 77 | + sum_reconstruction = np.sum(abs(filters_spectrum), axis=1) |
| 78 | + |
| 79 | + for b in range(ob.shape[-1]): |
| 80 | + c = antoni_octave_bands.centers[b] |
| 81 | + line_label = ( |
| 82 | + (f"{c if c < 1000 else c / 1000:.0f}{'k' if c >= 1000 else ''}Hz") |
| 83 | + if lbl == "Antoni" |
| 84 | + else None |
| 85 | + ) |
| 86 | + axes[0, i].plot(freq, abs(filters_spectrum[:, b]) ** 2, label=line_label) |
| 87 | + axes[1, i].plot(time, ob[:, b]) |
| 88 | + axes[0, i].plot( |
| 89 | + freq, sum_reconstruction, label="sum magnitude" if lbl == "Antoni" else None |
| 90 | + ) |
| 91 | + axes[0, i].set_title(f"{lbl} - energy response") |
| 92 | + axes[1, i].set_title(f"{lbl} - impulse response") |
| 93 | + axes[1, i].set_xlim(0.025, 0.04) |
| 94 | + |
| 95 | + fig.legend(loc="lower right") |
| 96 | + fig.tight_layout() |
| 97 | + |
| 98 | + # Test octave bands interpolation |
| 99 | + coeffs = np.arange(antoni_octave_bands.n_bands)[::-1] + 1 |
| 100 | + |
| 101 | + mat_interp = { |
| 102 | + "Cosine": octave_bands.synthesis(coeffs, min_phase=True), |
| 103 | + "Cosine, min. phase": octave_bands.synthesis(coeffs, min_phase=False), |
| 104 | + "Antoni": antoni_octave_bands.synthesis(coeffs, min_phase=False), |
| 105 | + "Antoni, min. phase": antoni_octave_bands.synthesis(coeffs, min_phase=True), |
| 106 | + } |
| 107 | + |
| 108 | + # Compare the energy of the original coefficients to those after filtering. |
| 109 | + energy = (octave_bands.get_bw() / octave_bands.fs * 2.0) * coeffs**2 |
| 110 | + bar_labels = [ |
| 111 | + f"{c if c < 1000 else c / 1000:.0f}{'k' if c >= 1000 else ''}" |
| 112 | + for c in antoni_octave_bands.centers |
| 113 | + ] + ["Total"] |
| 114 | + bar_energy_original = energy.tolist() + [energy.sum()] |
| 115 | + |
| 116 | + bar_width = 0.9 / (len(mat_interp) + 2) |
| 117 | + |
| 118 | + def bar_x(idx): |
| 119 | + bar_space = 0.1 / (len(mat_interp) + 2) |
| 120 | + x = np.arange(len(bar_energy_original)) |
| 121 | + return bar_width / 2.0 + idx * (bar_width + bar_space) + x |
| 122 | + |
| 123 | + fig, axes = plt.subplots(1, 3, figsize=(15, 5)) |
| 124 | + axes[0].set_title("Impulse response") |
| 125 | + axes[1].set_title("Magnitude response") |
| 126 | + axes[2].set_title("Per-band energy") |
| 127 | + for idx, (lbl, values) in enumerate(mat_interp.items()): |
| 128 | + axes[0].plot(np.arange(values.shape[-1]) / fs, values, label=lbl) |
| 129 | + axes[0].set_xlabel("Time (s)") |
| 130 | + |
| 131 | + H = abs(np.fft.rfft(values)) |
| 132 | + axes[1].plot(np.arange(H.shape[-1]) * fs / values.shape[-1], H, label=lbl) |
| 133 | + axes[1].set_xlabel("Frequency (Hz)") |
| 134 | + |
| 135 | + energy_bands = antoni_octave_bands.energy(values, oversampling=4).tolist() |
| 136 | + energy_bands += [sum(energy_bands)] |
| 137 | + axes[2].bar(bar_x(idx), energy_bands, label=lbl, width=bar_width) |
| 138 | + |
| 139 | + axes[2].bar( |
| 140 | + bar_x(len(mat_interp)), bar_energy_original, label="True", width=bar_width |
| 141 | + ) |
| 142 | + axes[2].set_xticks(np.arange(len(bar_energy_original)) + 0.5, bar_labels) |
| 143 | + axes[2].set_xlabel("Band centers (Hz)") |
| 144 | + axes[2].legend() |
| 145 | + fig.tight_layout() |
| 146 | + |
| 147 | + # Test reconstruction |
| 148 | + time = np.arange(fs * 5) / fs |
| 149 | + f0 = 1.0 |
| 150 | + x = np.zeros_like(time) |
| 151 | + x[fs : 3 * fs] = np.sin(2 * np.pi * f0 * time[fs : 3 * fs]) |
| 152 | + x = chirp(time, 100, time[-1], 7500, method="linear") |
| 153 | + |
| 154 | + x_bands = octave_bands.analysis(x) |
| 155 | + x_rec = x_bands.sum(axis=-1) |
| 156 | + |
| 157 | + band_energy = antoni_octave_bands.energy(x, oversampling=4) |
| 158 | + print( |
| 159 | + f"Energy of input signal: {np.square(x).sum()}, " |
| 160 | + f"sum of band energies: {band_energy.sum()}" |
| 161 | + ) |
| 162 | + print(f"Reconstruction error: {abs(x - x_rec).max()}") |
| 163 | + |
| 164 | + low = None |
| 165 | + high = None |
| 166 | + |
| 167 | + low = 0 if low is None else low |
| 168 | + high = x.shape[-1] if high is None else high |
| 169 | + |
| 170 | + num_plots = x_bands.shape[-1] + 1 |
| 171 | + freq = np.arange(x_bands.shape[-2] // 2 + 1) / x_bands.shape[-2] * fs |
| 172 | + time = np.arange(x.shape[-1]) / fs |
| 173 | + fig, axes = plt.subplots(num_plots, 2, figsize=(10, 6)) |
| 174 | + axes[0, 0].plot(time[low:high], x_rec[low:high], label="reconstructed") |
| 175 | + axes[0, 0].plot(time[low:high], x[low:high], label="original") |
| 176 | + axes[0, 0].plot(time[low:high], x[low:high] - x_rec[low:high], label="error") |
| 177 | + L = axes[0, 1].magnitude_spectrum(x, label="reconstructed", Fs=fs) |
| 178 | + L = axes[0, 1].magnitude_spectrum(x_rec, label="reconstructed", Fs=fs) |
| 179 | + axes[0, 0].legend(fontsize="xx-small") |
| 180 | + |
| 181 | + ylim_time = x.min(), x.max() |
| 182 | + ylim = -0.01 * max(L[0]), max(L[0]) |
| 183 | + axes[0, 0].set_title("Filtered signal") |
| 184 | + axes[0, 1].set_title("Magnitude response") |
| 185 | + axes[0, 1].set_xlabel("") |
| 186 | + axes[0, 1].set_ylabel("") |
| 187 | + axes[0, 1].set_ylim(ylim) |
| 188 | + axes[0, 1].set_xticks([]) |
| 189 | + axes[0, 0].set_ylim(ylim_time) |
| 190 | + axes[0, 0].set_xticks([]) |
| 191 | + axes[0, 0].set_ylabel("Sum") |
| 192 | + for b in range(1, num_plots): |
| 193 | + axes[b, 0].plot(time[low:high], x_bands[low:high, b - 1]) |
| 194 | + axes[b, 0].set_ylim(ylim_time) |
| 195 | + axes[b, 1].magnitude_spectrum(x_bands[low:high, b - 1], Fs=fs) |
| 196 | + axes[b, 1].set_ylim(ylim) |
| 197 | + axes[b, 1].set_xlabel("") |
| 198 | + axes[b, 1].set_ylabel("") |
| 199 | + if b < num_plots - 1: |
| 200 | + axes[b, 0].set_xticks([]) |
| 201 | + axes[b, 1].set_xticks([]) |
| 202 | + else: |
| 203 | + axes[b, 0].set_xlabel("Time (s)") |
| 204 | + ticks = np.arange(0, fs / 2 + 1, 1000) |
| 205 | + ticklabels = [ |
| 206 | + f"{c if c < 1000 else c / 1000:.0f}{'k' if c > 1000 else ''} Hz" |
| 207 | + for c in ticks |
| 208 | + ] |
| 209 | + axes[b, 1].set_xticks(ticks, ticklabels) |
| 210 | + if b > 0: |
| 211 | + c = antoni_octave_bands.centers[b - 1] |
| 212 | + axes[b, 0].set_ylabel( |
| 213 | + f"{c if c < 1000 else c / 1000:.0f}{'k' if c > 1000 else ''} Hz" |
| 214 | + ) |
| 215 | + axes[b, 0].set_xlabel("Frequency (Hz)") |
| 216 | + fig.tight_layout() |
| 217 | + |
| 218 | + plt.show() |
0 commit comments