Correcting the sparse rir problem, putting the distance attenuation d…

…irectly in c++

pyroomacoustics/libroom_src/room.cpp

@@ -629,13 +629,13 @@ std::vector<entry> Room::get_rir_entries(size_t nb_phis,
 
 	for (size_t i(0); i<nb_phis; ++i)  {
 		//std::cout << "\n===============\ni="<< i << std::endl;
-		float phi = 2*M_PI/nb_phis;
+		float phi = i*2*M_PI/nb_phis;
 
 		for (size_t j(0); j<nb_thetas; ++j){
 
 			//std::cout << "j=" << j << std::endl;
 
-			float theta = M_PI/nb_thetas;
+			float theta = j*M_PI/nb_thetas;
 
 			// For 2D, this parameter means nothing, but we set it to
 			// PI/2 to be consistent
@@ -656,6 +656,11 @@ std::vector<entry> Room::get_rir_entries(size_t nb_phis,
 
 	}
 
+	// Now we must apply the distance attenuation for every ray
+	for (size_t k(0); k < output.size(); k++){
+		output[k][1] /= (output[k][0] * sound_speed);
+	}
+
 	return output;
 
 

pyroomacoustics/libroom_src/tests/test_ray_tracing.py

@@ -93,15 +93,17 @@ def test_room_construct():
 def compute_rir(log, time_thres, fs, plot=True):
 
 
+    TIME = 0
+    ENERGY = 1
+
+    # ======= WITH FRACTIONAL PART =======
+
     # the python utilities to compute the rir
     fdl = pra.constants.get('frac_delay_length')
     fdl2 = (fdl - 1) // 2  # Integer division
 
     ir = np.zeros(int(time_thres * fs) + fdl)
 
-    TIME = 0
-    ENERGY = 1
-
     for entry in log:
         time_ip = int(np.floor(entry[TIME] * fs))
 
@@ -111,7 +113,7 @@ def compute_rir(log, time_thres, fs, plot=True):
         time_fp = (entry[TIME] * fs) - time_ip
 
         # Distance attenuation
-        ir[time_ip - fdl2:time_ip + fdl2 + 1] += (entry[ENERGY] * fractional_delay(time_fp)) / (sound_speed * entry[TIME])
+        ir[time_ip - fdl2:time_ip + fdl2 + 1] += (entry[ENERGY] * fractional_delay(time_fp))
 
 
     if plot :
@@ -124,7 +126,7 @@ def compute_rir(log, time_thres, fs, plot=True):
 
     return ir
 
-def apply_rir(rir, wav_file_name, fs, result_name="aaa.wav"):
+def apply_rir(rir, wav_file_name, fs, cutoff=200, result_name="aaa.wav"):
 
 
     fs0, audio_anechoic = wavfile.read(wav_file_name)
@@ -138,20 +140,27 @@ def apply_rir(rir, wav_file_name, fs, result_name="aaa.wav"):
     # Compute the convolution and set all coefficients between -1 and 1 (range for float32 .wav files)
     result = scipy.signal.fftconvolve(rir, audio_anechoic)
 
+    if cutoff > 0:
+        result = highpass(result, fs, cutoff)
+
     result /= np.abs(result).max()
     result -= np.mean(result)
     wavfile.write(result_name, rate=fs, data=result.astype('float32'))
 
-
+def highpass(audio, fs, cutoff=200, butter_order=5):
+    nyq = 0.5 * fs
+    fc_norm = cutoff / nyq
+    b, a = signal.butter(butter_order, fc_norm, btype="high", analog=False)
+    return signal.lfilter(b, a, audio)
 
 
 if __name__ == '__main__':
     room = test_room_construct()
 
     # parameters
-    nb_phis = 20
-    nb_thetas = 20
-    source_pos = [9,9,7]
+    nb_phis = 50
+    nb_thetas = 50
+    source_pos = [9,7,9]
     mic_radius = 0.05
     scatter_coef = 0.1
     time_thres = 0.5 #s
@@ -168,7 +177,10 @@ def apply_rir(rir, wav_file_name, fs, result_name="aaa.wav"):
 
     rir = compute_rir(log, time_thres, fs, plot=True)
 
-    apply_rir(rir, "0riginal.wav", fs)
+    if len(rir) == 0:
+        raise ValueError("The room impulse response is empty !")
+
+    apply_rir(rir, "0riginal.wav", fs, cutoff=200)