Merge with master

.github/workflows/pythonpackage.yml

@@ -10,7 +10,13 @@ jobs:
       fail-fast: false
       max-parallel: 12
       matrix:
-        os: [ubuntu-latest, macos-13, macos-latest, windows-latest]
+        # switch temporarily to windows-2019 due to a breaking
+        # change making mutex constexpr
+        # a crash occurs if python and the extension are built
+        # with different versions
+        # over time, this should hopefully resolves itself
+        # ref: https://github.com/microsoft/STL/issues/4875
+        os: [ubuntu-latest, macos-13, macos-latest, windows-2019]
         python-version: [3.8, 3.9, "3.10", "3.11", "3.12"]
         exclude:
           - os: macos-latest
@@ -29,7 +35,7 @@ jobs:
       uses: actions/setup-python@v4
       with:
         python-version: ${{ matrix.python-version }}
-    - name: Install dependencies
+    - name: Upgrade pip
       run: |
         python -m pip install --upgrade pip
     - name: Workaround for windows and python 3.7

CHANGELOG.rst

@@ -43,6 +43,35 @@ Changed
 `Unreleased`_
 -------------
 
+Nothing yet
+
+`0.7.7`_ - 2024-09-09
+---------------------
+
+Bugfix
+~~~~~~
+
+- Fixes a bug when plotting an anechoic room
+- Fixes a buggy assert from ``rir_build.cpp``
+- Makes the build system consistent for all source files in ``pyroomacoustics/libroom_src``
+- Fixes a call to a deprecated function of scipy.integrate in ``pyroomacoustics/denoise/iterative_wiener`` (issue #362)
+
+`0.7.6`_ - 2024-08-05
+---------------------
+
+Bugfix
+~~~~~~
+
+- Fixes the default value of ``energy_thres`` in ``experimental.rt60`` to match the
+  previous behavior of the function (issue #358)
+- Further fixes issue with cast reflections delays to float32 in room.py (#353)
+  which was not fully fixed by the previous update
+- Fixes calls to deprecated API of ``scipy.linalg.eigh``
+- Fixes use of deprecated feature of numpy (conversion of singleton array to scalar)
+
+`0.7.5`_ - 2024-06-18
+---------------------
+
 Changed
 ~~~~~~~
 
@@ -58,6 +87,8 @@ Bugfix
   ticks when called with ``kind="tf"``
 - Fixes an issue where the ``visibility`` attribute of the room is not
   set if there are no visible source or image source in the room. (#313)
+- Fixes issue with cast reflections delays to float32 in room.py (#353)
+- Fixes calls to ``numpy.linalg.solve`` with Numpy 2.0 API
 
 `0.7.4`_ - 2024-04-25
 ---------------------
@@ -603,7 +634,10 @@ Changed
    ``pyroomacoustics.datasets.timit``
 
 
-.. _Unreleased: https://github.com/LCAV/pyroomacoustics/compare/v0.7.4...master
+.. _Unreleased: https://github.com/LCAV/pyroomacoustics/compare/v0.7.7...master
+.. _0.7.7: https://github.com/LCAV/pyroomacoustics/compare/v0.7.6...v0.7.7
+.. _0.7.6: https://github.com/LCAV/pyroomacoustics/compare/v0.7.5...v0.7.6
+.. _0.7.5: https://github.com/LCAV/pyroomacoustics/compare/v0.7.4...v0.7.5
 .. _0.7.4: https://github.com/LCAV/pyroomacoustics/compare/v0.7.3...v0.7.4
 .. _0.7.3: https://github.com/LCAV/pyroomacoustics/compare/v0.7.2...v0.7.3
 .. _0.7.2: https://github.com/LCAV/pyroomacoustics/compare/v0.7.1...v0.7.2

pyroomacoustics/beamforming.py

@@ -1209,7 +1209,7 @@ def rake_max_udr_filters(
         SINR, v = la.eigh(
             A,
             b=B,
-            eigvals=(self.M * self.Lg - 1, self.M * self.Lg - 1),
+            subset_by_index=(self.M * self.Lg - 1, self.M * self.Lg - 1),
             overwrite_a=True,
             overwrite_b=True,
             check_finite=False,
@@ -1309,7 +1309,7 @@ def rake_max_sinr_filters(self, source, interferer, R_n, epsilon=5e-3, delay=0.0
         SINR, v = la.eigh(
             K_s,
             b=K_nq,
-            eigvals=(self.M * self.Lg - 1, self.M * self.Lg - 1),
+            subset_by_index=(self.M * self.Lg - 1, self.M * self.Lg - 1),
             overwrite_a=True,
             overwrite_b=True,
             check_finite=False,

pyroomacoustics/bss/auxiva.py

@@ -227,7 +227,7 @@ def demix(Y, X, W):
             )
 
             WV = np.matmul(W_hat, V)
-            W[:, s, :] = np.conj(np.linalg.solve(WV, eyes[:, :, s]))
+            W[:, s, :] = np.conj(np.linalg.solve(WV, eyes[:, :, [s]]))[..., 0]
 
             # normalize
             denom = np.matmul(

pyroomacoustics/bss/fastmnmf.py

@@ -175,8 +175,8 @@ def separate():
 
                 try:
                     tmp_FM = np.linalg.solve(
-                        np.matmul(Q_FMM, V_FMM), np.eye(n_chan)[None, m]
-                    )
+                        np.matmul(Q_FMM, V_FMM), np.eye(n_chan)[None, :, [m]]
+                    )[..., 0]
                 except np.linalg.LinAlgError:
                     # If Gaussian elimination fails due to a singlular matrix, we
                     # switch to the pseudo-inverse solution

pyroomacoustics/bss/fastmnmf2.py

@@ -165,8 +165,8 @@ def separate():
                     np.einsum("ftij, ft -> fij", XX_FTMM, 1 / Y_FTM[..., m]) / n_frames
                 )
                 tmp_FM = np.linalg.solve(
-                    np.matmul(Q_FMM, V_FMM), np.eye(n_chan)[None, m]
-                )
+                    np.matmul(Q_FMM, V_FMM), np.eye(n_chan)[None, :, [m]]
+                )[..., 0]
                 Q_FMM[:, m] = (
                     tmp_FM
                     / np.sqrt(

pyroomacoustics/bss/ilrma.py

@@ -159,7 +159,7 @@ def demix(Y, X, W):
             C = np.matmul((X * iR[s, :, None, :]), np.conj(X.swapaxes(1, 2))) / n_frames
 
             WV = np.matmul(W, C)
-            W[:, s, :] = np.conj(np.linalg.solve(WV, eyes[:, :, s]))
+            W[:, s, :] = np.conj(np.linalg.solve(WV, eyes[:, :, [s]]))[..., 0]
 
             # normalize
             denom = np.matmul(

pyroomacoustics/bss/sparseauxiva.py

@@ -148,7 +148,7 @@ def demixsparse(Y, X, S, W):
             W_H = np.conj(np.swapaxes(W, 1, 2))
             WV = np.matmul(W_H, V[:, s, :, :])
             rhs = I[None, :, s][[0] * WV.shape[0], :]
-            W[:, :, s] = np.linalg.solve(WV, rhs)
+            W[:, :, s] = np.linalg.solve(WV, rhs[..., None])[..., 0]
 
             # normalize
             P1 = np.conj(W[:, :, s])

pyroomacoustics/denoise/iterative_wiener.py

@@ -337,7 +337,14 @@ def _lpc_all_pole(omega):
     d_omega = 2 * np.pi / 1000
     omega_vals = np.arange(-np.pi, np.pi, d_omega)
     vec_integrand = np.vectorize(_lpc_all_pole)
-    integral = integrate.trapz(vec_integrand(omega_vals), omega_vals)
+
+    try:
+        integral = integrate.trapezoid(vec_integrand(omega_vals), omega_vals)
+    except AttributeError:
+        # older versions of scipy do not have the function 'trapezoid'
+        # fall back to the legacy function name
+        integral = integrate.trapz(vec_integrand(omega_vals), omega_vals)
+
     return rhs * 2 * np.pi / N / integral
 
 

pyroomacoustics/denoise/tests/test_iterative_wiener.py

@@ -0,0 +1,38 @@
+import numpy as np
+
+import pyroomacoustics as pra
+
+
+def test_iterative_wiener():
+    """
+    A simple functional test that the call does not produce any errors.
+    """
+    # parameters
+    num_blocks = 20
+    nfft = 512
+    hop = nfft // 2
+
+    # create a dummy signal
+    blocks = np.random.randn(num_blocks, hop)
+
+    # initialize STFT and IterativeWiener objects
+    stft = pra.transform.STFT(nfft, hop=hop, analysis_window=pra.hann(nfft))
+    scnr = pra.denoise.IterativeWiener(
+        frame_len=nfft, lpc_order=20, iterations=2, alpha=0.8, thresh=0.01
+    )
+
+    # apply block-by-block
+    processed_blocks = []
+    for n in range(num_blocks):
+
+        # go to frequency domain, 50% overlap
+        stft.analysis(blocks[n])
+
+        # compute wiener output
+        X = scnr.compute_filtered_output(
+            current_frame=stft.fft_in_buffer, frame_dft=stft.X
+        )
+
+        # back to time domain
+        mono_denoised = stft.synthesis(X)
+        processed_blocks.append(mono_denoised)

pyroomacoustics/doa/music.py

@@ -150,7 +150,7 @@ def _compute_spatial_spectrum(self, cross, k):
             Dc = np.array(self.mode_vec[k, :, n], ndmin=2).T
             Dc_H = np.conjugate(np.array(self.mode_vec[k, :, n], ndmin=2))
             denom = np.linalg.multi_dot([Dc_H, cross, Dc])
-            P[n] = 1 / abs(denom)
+            P[n] = 1 / abs(denom[0, 0])
 
         return P
 

pyroomacoustics/experimental/localization.py

@@ -208,7 +208,7 @@ def edm_line_search(R, tdoa, bounds, steps):
     for i in range(d.shape[0]):
         D[-1, :-1] = D[:-1, -1] = (dif + d[i]) ** 2
         w = np.sort(np.abs(la.eigh(D, eigvals_only=True)))
-        # w = la.eigh(D, eigvals_only=True, eigvals=(D.shape[0]-6,D.shape[0]-6))
+        # w = la.eigh(D, eigvals_only=True, subset_by_index=(D.shape[0]-6,D.shape[0]-6))
         cost[i] = np.sum(w[: D.shape[0] - 5])
 
     return cost, d

pyroomacoustics/experimental/rt60.py

@@ -33,7 +33,7 @@
 import numpy as np
 
 
-def measure_rt60(h, fs=1, decay_db=60, energy_thres=0.95, plot=False, rt60_tgt=None):
+def measure_rt60(h, fs=1, decay_db=60, energy_thres=1.0, plot=False, rt60_tgt=None):
     """
     Analyze the RT60 of an impulse response. Optionaly plots some useful information.
 

pyroomacoustics/libroom_src/rir_builder.cpp

@@ -73,8 +73,6 @@ void threaded_rir_builder_impl(
   int max_sample = int(std::ceil(fs * t_max)) + fdl2;
   int min_sample = int(std::floor(fs * t_min)) - fdl2;
 
-  assert(min_samples >= 0);
-
   if (min_sample < 0)
     throw std::runtime_error("minimum time recorded is less than 0");
   if (max_sample >= int(rir_len))

pyroomacoustics/libroom_src/rir_builder.hpp

@@ -16,4 +16,6 @@ void delay_sum(const py::buffer irs, const py::buffer delays, py::buffer output,
 void fractional_delay(py::buffer out, const py::buffer time, size_t lut_gran,
                       size_t num_threads);
 
+#include "rir_builder.cpp"
+
 #endif  // __RIR_BUILDER_HPP__

pyroomacoustics/libroom_src/room.cpp

@@ -832,9 +832,10 @@ bool Room<D>::scat_ray(
       // of scat_per_slot
       if (travel_dist_at_mic < distance_thres && scat_trans.maxCoeff() > energy_thres)
       {
-
-        //output[k].push_back(Hit(travel_dist_at_mic, scat_trans));        
-        //microphones[k].log_histogram(output[k].back(), hit_point);
+        // The (r_sq * p_hit) normalization factor is necessary to equalize the energy
+        // of the IR computed with ray tracing to that of the image source method.
+        // Ref: D. Schroeder, "Physically based real-time auralization of interactive virtual environments",
+        // section 5.4, eq. 5.54.
         double r_sq = double(travel_dist_at_mic) * travel_dist_at_mic;
         auto p_hit = (1 - sqrt(1 - mic_radius_sq / std::max(mic_radius_sq, r_sq)));
         Eigen::ArrayXf energy = scat_trans / (r_sq * p_hit) ;
@@ -944,6 +945,10 @@ void Room<D>::simul_ray(
           //   because the ray will continue its way          
           float travel_dist_at_mic = travel_dist + distance;
 
+          // The (r_sq * p_hit) normalization factor is necessary to equalize the energy
+          // of the IR computed with ray tracing to that of the image source method.
+          // Ref: D. Schroeder, "Physically based real-time auralization of interactive virtual environments",
+          // section 5.4, eq. 5.54.
           double r_sq = double(travel_dist_at_mic) * travel_dist_at_mic;
           auto p_hit = (1 - sqrt(1 - mic_radius_sq / std::max(mic_radius_sq, r_sq)));
           energy = transmitted / (r_sq * p_hit);