Merge branch 'develop'

README.md

@@ -49,10 +49,10 @@ pip3 install .
 Usage
 --------
 
-See `examples` folder (they require `libsndfile` module installed). You need to `cd examples` and run `python3 mesh_sim.py` (we recommend starting with this one). This script demonstrates two equivalent ways to define the environment for sound propagation, and save the impulse response as an audio file. You can use a `.obj` file with an optional `.mtl` file with the same name to define the room geometry and materials. In this case, the `.mtl` file has two extra rows compared with conventional `.mtl` file used for visual rendering:
+See `examples` folder (extra modules may be required). You need to `cd examples` and run `python3 mesh_sim.py` (we recommend starting with this one). This script demonstrates two equivalent ways to define the environment for sound propagation, and save the impulse response as an audio file. You can use a `.obj` file with an optional `.mtl` file with the same name to define the room geometry and materials. In this case, the `.mtl` file has two extra rows compared with conventional `.mtl` file used for visual rendering:
 ```
 sound_a 0.5 0.6 0.6 0.7 0.75 0.8 0.9 0.9  # sound absorption coefficients, for 8 octave bands [62.5, 125, 250, 500, 1000, 2000, 4000, 8000]Hz
-sound_s 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 # sound scattering coefficients, if you don't know the details of diffuse/specular reflections, keep it 0.5
+sound_s 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 # sound scattering coefficients, if you don't know the details of diffuse/specular reflections, keep it low
 ```
 or directly create a shoebox shaped room using our API:
 ```

examples/compare_image.py

@@ -1,13 +1,14 @@
 import numpy as np
 import pygsound as ps
 import matplotlib.pyplot as plt
-import rirgenerator as rg
+import rir_generator as rir
+# using https://github.com/audiolabs/rir-generator/
 
 
 def main():
     l = 10
     w = 6
-    h = 2
+    h = 4
     absorb = 0.3
     reflec = np.sqrt(1.0 - absorb)
 
@@ -17,7 +18,7 @@ def main():
     ctx.channel_type = ps.ChannelLayoutType.mono
     scene = ps.Scene()
 
-    mesh = ps.createbox(l, w, h, absorb, 0.5)
+    mesh = ps.createbox(l, w, h, absorb, 0.1)
     scene.setMesh(mesh)
 
     src_coord = [1, 1, 0.5]
@@ -29,39 +30,34 @@ def main():
     lis = ps.Listener(lis_coord)
     lis.radius = 0.01
 
-    rir_gs = scene.computeIR(src, lis, ctx)
-    rir_img = rg.rir_generator(343, 16000, lis_coord, src_coord, [l, w, h], beta=[reflec]*6)
-    rir_img = rir_img / max(abs(rir_img[0]))
-
-    max_cnt = min(len(rir_gs['samples']), len(rir_img[0])) * 2
-
-    fig, axs = plt.subplots(4, 1, figsize=(10, 20))
+    rir_gs = scene.computeIR([src], [lis], ctx)
+    rir_img = rir.generate(
+        c=343,  # Sound velocity (m/s)
+        fs=ctx.sample_rate,  # Sample frequency (samples/s)
+        r=[lis_coord],
+        s=src_coord,  # Source position [x y z] (m)
+        L=[l, w, h],  # Room dimensions [x y z] (m)
+        beta=[reflec] * 6,  # Reflection coefficients
+    )
+    rir_img = rir_img[:, 0] / max(abs(rir_img))
+    rir_gs = rir_gs['samples'][0][0][0]
+    max_cnt = min(len(rir_gs), len(rir_img))
+
+    fig, axs = plt.subplots(2, 1, figsize=(16, 10))
     axs[0].set_title('Image Method (waveform)')
-    axs[0].plot(rir_img[0], linewidth=0.5)
+    axs[0].plot(rir_img, linewidth=0.5)
     axs[0].set_xlabel('Sample')
     axs[0].set_xlim(0, max_cnt)
     axs[0].set_ylim(-1, 1)
     axs[0].set_ylabel('Amplitude')
 
     axs[1].set_title('Geometric Sound Propagation (waveform)')
-    axs[1].plot(rir_gs['samples'], linewidth=0.5)
+    axs[1].plot(rir_gs, linewidth=0.5)
     axs[1].set_xlabel('Sample')
     axs[1].set_xlim(0, max_cnt)
     axs[1].set_ylim(-1, 1)
     axs[1].set_ylabel('Amplitude')
 
-    axs[2].set_title('Image Method (spectrogram)')
-    axs[2].specgram(rir_img[0], mode='magnitude', NFFT=1024, Fs=16000, noverlap=512)
-    axs[2].set_xlim(0, max_cnt / 16000)
-    axs[2].set_xlabel('Times (s)')
-    axs[2].set_ylabel('Frequency (Hz)')
-
-    axs[3].set_title('Geometric Sound Propagation (spectrogram)')
-    axs[3].specgram(rir_gs['samples'], mode='magnitude', NFFT=1024, Fs=16000, noverlap=512)
-    axs[3].set_xlim(0, max_cnt / 16000)
-    axs[3].set_xlabel('Times (s)')
-    axs[3].set_ylabel('Frequency (Hz)')
-
     fig.tight_layout()
     plt.savefig('img_comparison.png')
     plt.show()

examples/cube.mtl

@@ -11,4 +11,4 @@ Ni 1.000000
 d 1.000000
 illum 2
 sound_a 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
-sound_s 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
+sound_s 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

examples/custom_array.py

@@ -1,11 +1,10 @@
-import os
 import numpy as np
 import pygsound as ps
 from wavefile import WaveWriter, Format
 
 
 def compute_array(meshpath, src_coord, lis_coord, r, s, micarray):
-    mesh = ps.loadobj(meshpath, os.path.join(os.path.dirname(meshpath), ''), r, s)
+    mesh = ps.loadobj(meshpath, r, s)
     ctx = ps.Context()
     ctx.diffuse_count = 20000
     ctx.specular_count = 2000
@@ -19,14 +18,13 @@ def compute_array(meshpath, src_coord, lis_coord, r, s, micarray):
     res = {}
     res_buffer = []
     rate = 0
-    abandon_flag = False
     for offset in micarray:
         lis = ps.Listener((offset + lis_coord).tolist())
         lis.radius = 0.01
 
-        res_ch = scene.computeIR(src, lis, ctx)
+        res_ch = scene.computeIR([src], [lis], ctx)
         rate = res_ch['rate']
-        sa = res_ch['samples']
+        sa = res_ch['samples'][0][0][0]
         res['rate'] = rate
         res_buffer.append(sa)
     res['samples'] = np.zeros((len(res_buffer), np.max([len(ps) for ps in res_buffer])))
@@ -43,7 +41,7 @@ def main():
     src_coord = [1, 1, 1]
     lis_coord = [0.1, 0.1, 0.1]
 
-    res = compute_array("cube.obj", src_coord, lis_coord, 0.5, 0.5, custom_array)
+    res = compute_array("cube.obj", src_coord, lis_coord, 0.5, 0.1, custom_array)
 
     with WaveWriter('custom_array.wav', channels=np.shape(res['samples'])[0], samplerate=int(res['rate'])) as w:
         w.write(np.array(res['samples']))

examples/mesh_sim.py

@@ -10,7 +10,7 @@ def main():
     ctx.specular_count = 2000
     ctx.channel_type = ps.ChannelLayoutType.stereo
 
-    mesh1 = ps.loadobj("cube.obj", "")    # if the second argument is empty, the code will infer the .mtl name using .obj name
+    mesh1 = ps.loadobj("cube.obj")
     scene = ps.Scene()
     scene.setMesh(mesh1)
 
@@ -19,29 +19,26 @@ def main():
 
     src = ps.Source(src_coord)
     src.radius = 0.01
+    src.power = 1.0
 
     lis = ps.Listener(lis_coord)
     lis.radius = 0.01
 
-    res = scene.computeMultichannelIR(src, lis, ctx)
-
-    with WaveWriter('test1.wav', channels=np.shape(res['samples'])[0], samplerate=int(res['rate'])) as w1:
-        w1.write(np.array(res['samples']))
+    res = scene.computeIR([src], [lis], ctx)    # you may pass lists of sources and listeners to get N_src x N_lis IRs
+    audio_data = np.array(res['samples'][0][0])     # the IRs are indexed by [i_src, i_lis, i_channel]
+    with WaveWriter('test1.wav', channels=audio_data.shape[0], samplerate=int(res['rate'])) as w1:
+        w1.write(audio_data)
         print("IR using .obj input written to test1.wav.")
 
     # Simulation using a shoebox definition
-    ctx = ps.Context()
-    ctx.diffuse_count = 20000
-    ctx.specular_count = 2000
-    ctx.channel_type = ps.ChannelLayoutType.stereo
-
-    mesh2 = ps.createbox(10, 6, 2, 0.5, 0.5)
+    mesh2 = ps.createbox(10, 6, 2, 0.5, 0.1)
     scene = ps.Scene()
     scene.setMesh(mesh2)
 
-    res = scene.computeMultichannelIR(src, lis, ctx)
-    with WaveWriter('test2.wav', channels=np.shape(res['samples'])[0], samplerate=int(res['rate'])) as w:
-        w2.write(np.array(res['samples']))
+    res = scene.computeIR([src_coord], [lis_coord], ctx)    # use default source and listener settings if you only pass coordinates
+    audio_data = np.array(res['samples'][0][0])
+    with WaveWriter('test2.wav', channels=audio_data.shape[0], samplerate=int(res['rate'])) as w2:
+        w2.write(audio_data)
         print("IR using shoebox input written to test2.wav.")
 
 

setup.py

@@ -72,7 +72,7 @@ def build_extension(self, ext):
 
 setup(
     name='pygsound',
-    version='0.2',
+    version='0.3',
     author='Zhenyu Tang, Carl Schissler, Dinesh Manocha',
     author_email='[email protected]',
     description='A room impulse response simulator using for geometric sound propagation',
@@ -92,6 +92,5 @@ def build_extension(self, ext):
     install_requires=[
         'Cython',
         'numpy',
-        'wavefile',
     ],
 )

src/GSound/gsound/gsMeshRequest.cpp

@@ -74,7 +74,7 @@ MeshRequest:: MeshRequest()
 		weldTolerance( 0.00001f ),
 		simplifyTolerance( 0.0001f ),
 		minDiffractionEdgeAngle( 0.0f ),
-		minDiffractionEdgeLength( 0.0f ),
+		minDiffractionEdgeLength( 0.5f ),
 		diffuseResolution( 0.5f ),
 		edgeResolution( 0.5f ),
 		minRaysPerEdge( 1 ),

src/GSound/gsound/gsMeshRequest.h

@@ -125,7 +125,7 @@ class MeshRequest
 			  * If the angle between normals of two neighboring triangles
 			  * is less than this value, then the edge that they share is not diffracting.
 			  * Thus, a larger diffraction angle will result in less diffracting edges while
-			  * lower thresold will result in more edges.
+			  * lower threshold will result in more edges.
 			  */
 			Real minDiffractionEdgeAngle;
 

src/pygsound/src/Context.cpp

@@ -12,7 +12,7 @@ Context::Context()
     prop_request.flags.set(gs::PropagationFlags::DIFFUSE, true);
     prop_request.flags.set(gs::PropagationFlags::DIFFRACTION, true);
     prop_request.flags.set(gs::PropagationFlags::SOURCE_DIRECTIVITY, false);
-    prop_request.flags.set(gs::PropagationFlags::DOPPLER_SORTING, true);
+    prop_request.flags.set(gs::PropagationFlags::DOPPLER_SORTING, false);
     prop_request.flags.set(gs::PropagationFlags::ADAPTIVE_QUALITY, false);
     prop_request.flags.set(gs::PropagationFlags::AIR_ABSORPTION, true);
     prop_request.flags.set(gs::PropagationFlags::ADAPTIVE_IR_LENGTH, true);

src/pygsound/src/Listener.hpp

@@ -36,7 +36,7 @@ class Listener
 	friend std::ostream &operator<<( std::ostream &_strm,
 	                                 const Listener &_list );
 
-	gs::SoundListener   m_listener;
+	gs::SoundListener m_listener;
 };
 
 #endif // INC_LISTENER_HPP

src/pygsound/src/Scene.cpp

@@ -21,65 +21,109 @@ Scene::setMesh( SoundMesh &_mesh )
 }
 
 py::dict
-Scene::computeIR( SoundSource &_source, Listener &_listener, Context &_context )
+Scene::computeIR( std::vector<SoundSource> &_sources, std::vector<Listener> &_listeners, Context &_context )
 {
-	m_scene.addSource( &_source.m_source);
-	m_scene.addListener( &_listener.m_listener );
-	if (m_scene.getObjectCount() == 0){
-        std::cout << "object count is zero, cannot propagate sound!" << std::endl;
-	}
+    int n_src = _sources.size();
+    int n_lis = _listeners.size();
+
+    for (SoundSource& p : _sources){
+        m_scene.addSource(&p.m_source);
+    }
+    for (Listener& p : _listeners){
+        m_scene.addListener(&p.m_listener);
+    }
+
+    if (m_scene.getObjectCount() == 0){
+        std::cerr << "object count is zero, cannot propagate sound!" << std::endl;
+    }
+
     propagator.propagateSound(m_scene, _context.internalPropReq(), sceneIR);
-	const gs::SoundSourceIR& sourceIR = sceneIR.getListenerIR(0).getSourceIR(0);
-	gs::ImpulseResponse result;
-	result.setIR(sourceIR, _listener.m_listener, _context.internalIRReq());
 
-	auto rate = result.getSampleRate();
+    py::list IRPairs(n_src);
+    auto rate = _context.getSampleRate();
+    for (int i_src = 0; i_src < n_src; ++i_src){
+        py::list srcSamples(n_src);
+        for (int i_lis = 0; i_lis < n_lis; ++i_lis){
+            const gs::SoundSourceIR& sourceIR = sceneIR.getListenerIR(i_lis).getSourceIR(i_src);
+            gs::ImpulseResponse result;
+            result.setIR(sourceIR, *m_scene.getListener(i_lis), _context.internalIRReq());
+            auto numOfChannels = int(result.getChannelCount());
 
-    auto *sample = result.getChannel( 0 );
-	std::vector<float> samples(sample, sample+result.getLengthInSamples());
+            py::list samples;
+            for (int ch = 0; ch < numOfChannels; ch++)
+            {
+                auto *sample_ch = result.getChannel(ch);
+                std::vector<float> samples_ch(sample_ch, sample_ch+result.getLengthInSamples());
+                samples.append(samples_ch);
+            }
+            srcSamples[i_lis] = samples;
+        }
+        IRPairs[i_src] = srcSamples;
+    }
 
-	m_scene.clearSources();
-	m_scene.clearListeners();
+    m_scene.clearSources();
+    m_scene.clearListeners();
 
-	py::dict ret;
-	ret["rate"] = rate;
-	ret["samples"] = samples;
-	return ret;
+    py::dict ret;
+    ret["rate"] = rate;
+    ret["samples"] = IRPairs;   // index by [i_src, i_lis, i_channel]
 
-//    return result;
+    return ret;
 }
 
 py::dict
-Scene::computeMultichannelIR( SoundSource &_source, Listener &_listener, Context &_context )
+Scene::computeIR( std::vector<std::vector<float>> &_sources, std::vector<std::vector<float>> &_listeners, Context &_context,
+                        float src_radius, float src_power, float lis_radius)
 {
-	m_scene.addSource( &_source.m_source);
-	m_scene.addListener( &_listener.m_listener );
-	if (m_scene.getObjectCount() == 0){
-		std::cout << "object count is zero, cannot propagate sound!" << std::endl;
-	}
-	propagator.propagateSound(m_scene, _context.internalPropReq(), sceneIR);
-	const gs::SoundSourceIR& sourceIR = sceneIR.getListenerIR(0).getSourceIR(0);
-//	gs::ImpulseResponse result;
-//	result.setIR(sourceIR, _listener.m_listener, _context.internalIRReq());
-	auto result = std::make_shared<gs::ImpulseResponse>();
-	result->setIR(sourceIR, _listener.m_listener, _context.internalIRReq());
-	auto rate = result->getSampleRate();
-
-    auto numOfChannels = int(result->getChannelCount());
-    assert(numOfChannels > 0);
-
-    py::list samples;
-    for (int ch = 0; ch < numOfChannels; ch++)
-    {
-        auto *sample_ch = result->getChannel(ch);
-		std::vector<float> samples_ch(sample_ch, sample_ch+result->getLengthInSamples());
-        samples.append(samples_ch);
+    int n_src = _sources.size();
+    int n_lis = _listeners.size();
+
+    // listener propagation is most expensive, so swap them for computation if there are more listeners
+    bool swapBuffer = false;
+    auto src_pos = std::ref(_sources);
+    auto lis_pos = std::ref(_listeners);
+    if (n_src < n_lis){
+        swapBuffer = true;
+        src_pos = std::ref(_listeners);
+        lis_pos = std::ref(_sources);
+        std::swap(n_src, n_lis);
     }
+
+    std::vector<SoundSource> sources;
+    std::vector<Listener> listeners;
+
+    for (auto p : src_pos.get()){
+        auto source = new SoundSource(p);
+        source->setRadius(src_radius);
+        source->setPower(src_power);
+        sources.push_back( *source );
+    }
+    for (auto p : lis_pos.get()){
+        auto listener = new Listener(p);
+        listener->setRadius(lis_radius);
+        listeners.push_back( *listener );
+    }
+
+    py::dict _ret = computeIR(sources, listeners, _context);
+    py::list IRPairs = _ret["samples"];
+
     py::dict ret;
-    ret["rate"] = rate;
-    ret["samples"] = samples;
-    m_scene.clearSources();
-    m_scene.clearListeners();
+    ret["rate"] = _ret["rate"];
 
-    return ret;
+    // swap the IR buffer if sources and listeners have been swapped
+    if (swapBuffer){
+        py::list swapIRPairs;
+        for (int i_lis = 0; i_lis < n_lis; ++i_lis) {
+            py::list srcSamples;
+            for (int i_src = 0; i_src < n_src; ++i_src){
+                srcSamples.append(IRPairs[i_src].cast<py::list>()[i_lis]);
+            }
+            swapIRPairs.append(srcSamples);
+        }
+        ret["samples"] = swapIRPairs;
+    }else{
+        ret["samples"] = IRPairs;   // index by [i_src, i_lis, i_channel]
+    }
+
+	return ret;
 }