diff --git a/examples/directivities/simulation_with_measured_directivity.py b/examples/directivities/simulation_with_measured_directivity.py index f4291666..55d4d416 100644 --- a/examples/directivities/simulation_with_measured_directivity.py +++ b/examples/directivities/simulation_with_measured_directivity.py @@ -50,7 +50,6 @@ """ import matplotlib.pyplot as plt - import pyroomacoustics as pra from pyroomacoustics.directivities import ( Cardioid, diff --git a/pyroomacoustics/beamforming.py b/pyroomacoustics/beamforming.py index 42cf2b19..3595f0c2 100644 --- a/pyroomacoustics/beamforming.py +++ b/pyroomacoustics/beamforming.py @@ -418,6 +418,7 @@ def record(self, signals, fs): raise NameError("The signals should be a 2D array.") if fs != self.fs: + self.signals = u.resample(signals, fs, self.fs) try: import samplerate diff --git a/pyroomacoustics/parameters.py b/pyroomacoustics/parameters.py index 978de47f..6e5c1b25 100644 --- a/pyroomacoustics/parameters.py +++ b/pyroomacoustics/parameters.py @@ -82,6 +82,7 @@ def get_num_threads(): "octave_bands_n_fft": 512, "octave_bands_base_freq": 125.0, "octave_bands_keep_dc": False, + "resample_backend": "soxr", } diff --git a/pyroomacoustics/tests/test_resample.py b/pyroomacoustics/tests/test_resample.py new file mode 100644 index 00000000..fd13bead --- /dev/null +++ b/pyroomacoustics/tests/test_resample.py @@ -0,0 +1,49 @@ +import matplotlib.pyplot as plt +import numpy as np +import pyroomacoustics as pra +import pytest + + +@pytest.mark.parametrize( + "fs_in, fs_out, backend", + [ + (240, 160, "soxr"), + (240, 160, "samplerate"), + (240, 160, "scipy"), + ], +) +def test_downsample(fs_in, fs_out, backend): + """Idea use a sine above Nyquist of fs_out. It should disappear.""" + assert fs_in > fs_out + f_sine = fs_out / 2.0 + (fs_in - fs_out) / 2.0 * 0.95 + time = np.arange(fs_in * 10) / fs_in + signal_in = np.sin(2.0 * np.pi * time * f_sine) + signal_out = pra.resample(signal_in, fs_in, fs_out, backend=backend) + + assert abs(signal_out).mean() < 1e-3 + + +if __name__ == "__main__": + + # Reads the file containing the Eigenmike's directivity measurements + eigenmike = pra.MeasuredDirectivityFile("EM32_Directivity") + + fs_tgt = 16000 + fs_file = eigenmike.fs + + print(f"{fs_file=} {fs_tgt=}") + + rir_original = eigenmike.impulse_responses[10, 20] + time_file = np.arange(rir_original.shape[0]) / fs_file + + rirs = {} + for backend in ["soxr", "samplerate", "scipy"]: + rirs[backend] = pra.resample(rir_original, fs_file, fs_tgt, backend=backend) + + fig, ax = plt.subplots(1, 1) + ax.plot(time_file, rir_original, label="Original") + for idx, (backend, rir) in enumerate(rirs.items()): + time_rir = np.arange(rir.shape[0]) / fs_tgt + ax.plot(time_rir, rir, label=backend, linewidth=(3 - idx)) + ax.legend() + plt.show() diff --git a/pyroomacoustics/utilities.py b/pyroomacoustics/utilities.py index cd5cc3d9..3f6ab51f 100644 --- a/pyroomacoustics/utilities.py +++ b/pyroomacoustics/utilities.py @@ -23,11 +23,12 @@ # not, see . from __future__ import division +import fractions import functools import itertools +import warnings import numpy as np -import soxr from scipy import signal from scipy.io import wavfile @@ -35,6 +36,20 @@ from .parameters import constants, eps from .sync import correlate +try: + import soxr + + _has_soxr = True +except ImportError: + _has_soxr = False + +try: + import samplerate + + _has_samplerate = True +except ImportError: + _has_samplerate = False + def requires_matplotlib(func): @functools.wraps(func) # preserves name, docstrings, and signature of function @@ -821,25 +836,68 @@ def angle_function(s1, v2): return np.vstack((az, co)) -def resample(data, old_fs, new_fs): +def resample(data, old_fs, new_fs, backend=None, *args, **kwargs): """ Resample an ndarray from ``old_fs`` to ``new_fs`` along the last axis. Parameters ---------- data : numpy array - Input data to be resampled. + Input data to be resampled expected in shape (..., num_samples). old_fs : int Original sampling rate. new_fs : int New sampling rate. + backend: str + The resampling backend to use. Options are as follows. + All extra arguments are passed to the backend. + + - `soxr`: The default backend. It is the fastest and most + accurate. It is not installed by default, but can be installed + via `pip install python-soxr`. + - `samplerate`: It is the first fallback backend. It is slower, + but as accurate as `soxr`. It is not installed by default, but can + be installed by `pip install samplerate`. + - `scipy`: It is the fallback when none of the other libraries + are installed. This uses `scipy.signal.resample_poly` and is not as + good as the other backend. This will generate a warning unless + specified explicitely. + + The backend used package-wide is set via the constants, + e.g., `pra.constants.set("resample_backend", "soxr")`. Returns ------- - numpy array + The resampled signal. """ + + if backend is None: + # get the package-wide default backend + backend = constants.get("resample_backend") + + if backend not in ("soxr", "samplerate", "scipy"): + raise ValueError( + "Possible choices for the resampling backend are " + "soxr | samplerate | scippy." + ) + + # select the backend + if backend == "soxr" and not _has_soxr: + backend = "samplerate" + + if backend == "samplerate" and not _has_samplerate: + backend = "scipy" + warnings.warn( + "Neither of the resampling backends `soxr` or `samplerate` are installed. " + "Falling back to scipy.signal.resample_poly. To silence this warning, " + "specify `backend=scipy` explicitely." + ) + + # format the data ndim = data.ndim + # for samplerate and soxr the data needs to be in format + # (num_samples, num_channels) if ndim == 1: data = data[:, None] elif ndim == 2: @@ -848,8 +906,25 @@ def resample(data, old_fs, new_fs): shape = data.shape data = data.reshape(-1, data.shape[-1]).T - resampled_data = soxr.resample(data, old_fs, new_fs) + if backend == "soxr": + resampled_data = soxr.resample(data, old_fs, new_fs, *args, **kwargs) + elif backend == "samplerate": + resampled_data = samplerate.resample( + data, new_fs / old_fs, "sinc_best", *args, **kwargs + ) + else: + # first, simplify the fraction + rate_frac = fractions.Fraction(int(new_fs), int(old_fs)) + resampled_data = signal.resample_poly( + data, + up=rate_frac.numerator, + down=rate_frac.denominator, + axis=0, + *args, + **kwargs + ) + # restore the original shape of the data if ndim == 1: resampled_data = resampled_data[:, 0] elif ndim == 2: