add minimal documentation on scattering

libtascar/include/acousticmodel.h

@@ -256,6 +256,7 @@ namespace TASCAR {
       bool muteonstop = false;
       // scatter filter parameters:
       uint32_t scatterreflections = 0;
+      float scatterspread = TASCAR_PIf * 0.125;
       float scatterstructuresize = 1;
       float scatterdamping = 0.0;
       // proxy parameters:

libtascar/include/fdn.h

@@ -215,6 +215,7 @@ namespace TASCAR {
     };
     void setpar_t60(float az, float daz, float t, float dt, float t60,
                     float damping, bool fixcirculantmat, bool truncate_forward);
+    void set_scatterpar(float daz, float t, float dt, float t60, float damping);
     void set_logdelays(bool ld) { logdelays_ = ld; };
     void set_zero()
     {

libtascar/src/acousticmodel.cc

@@ -592,7 +592,7 @@ receiver_t::receiver_t(tsccfg::node_t xmlsrc, const std::string& name,
                 "Number of reflections created by scattering filter");
   GET_ATTRIBUTE(scatterstructuresize, "m", "size of scatter structure");
   GET_ATTRIBUTE(scatterdamping, "", "damping of scatter reflection filter");
-
+  GET_ATTRIBUTE_DEG(scatterspread,"Spatial spread of scattering");
   // end proxy
   if(avgdist <= 0)
     avgdist = 0.5f * powf(volumetric.boxvolumef(), 0.33333f);
@@ -620,18 +620,17 @@ void receiver_t::configure()
   if(scatterreflections > 0) {
     scatterfilter = new TASCAR::fdn_t(scatterreflections, (uint32_t)f_sample,
                                       true, TASCAR::fdn_t::mean, false);
-    scatterfilter->setpar_t60(
-        0.0f, TASCAR_2PIf * 0.125f,
-        f_sample * (0.1f * scatterstructuresize / 340.0f),
+    scatterfilter->set_scatterpar(
+        scatterspread, f_sample * (0.1f * scatterstructuresize / 340.0f),
         f_sample * (scatterstructuresize / 340.0f), f_sample,
-        std::max(0.0f, std::min(0.999f, scatterdamping)), true, false);
+        std::max(0.0f, std::min(0.999f, scatterdamping)));
     scatterallpass_w.resize(scatterreflections);
     scatterallpass_x.resize(scatterreflections);
     scatterallpass_y.resize(scatterreflections);
     scatterallpass_z.resize(scatterreflections);
     size_t k = 1;
     for(auto& flt : scatterallpass_x) {
-      flt.set_allpass(0.9f, TASCAR_2PI * 0.25 * k / scatterreflections);
+      flt.set_allpass(0.89f, TASCAR_2PI * 0.25 * k / scatterreflections);
       ++k;
     }
     k = 1;
@@ -641,12 +640,12 @@ void receiver_t::configure()
     }
     k = 1;
     for(auto& flt : scatterallpass_z) {
-      flt.set_allpass(0.9f, TASCAR_2PI * 0.25 * k / scatterreflections);
+      flt.set_allpass(0.905f, TASCAR_2PI * 0.25 * k / scatterreflections);
       ++k;
     }
     k = 1;
     for(auto& flt : scatterallpass_w) {
-      flt.set_allpass(0.9f, TASCAR_2PI * 0.25 * k / scatterreflections);
+      flt.set_allpass(0.91f, TASCAR_2PI * 0.25 * k / scatterreflections);
       ++k;
     }
   }

libtascar/src/fdn.cc

@@ -157,6 +157,87 @@ void fdn_t::setpar_t60(float az, float daz, float t_min, float t_max, float t60,
   }
 }
 
+void fdn_t::set_scatterpar(float daz, float t_min, float t_max, float t60,
+                           float damping)
+{
+  // set delays:
+  set_zero();
+  float t_mean(0);
+  for(uint32_t tap = 0; tap < fdnorder_; ++tap) {
+    float t_(t_min);
+    if(logdelays_) {
+      // logarithmic distribution:
+      if(fdnorder_ > 1)
+        t_ =
+            t_min * powf(t_max / t_min, (float)tap / ((float)fdnorder_ - 1.0f));
+      ;
+    } else {
+      // squareroot distribution:
+      if(fdnorder_ > 1)
+        t_ = t_min + (t_max - t_min) *
+                         powf((float)tap / ((float)fdnorder_ - 1.0f), 0.5f);
+    }
+    uint32_t d((uint32_t)std::max(0.0f, t_));
+    fdnpath[tap].delay = std::max(2u, std::min(maxdelay_ - 1u, d));
+    fdnpath[tap].reflection.set_eta(0.87f * (float)tap /
+                                    ((float)fdnorder_ - 1.0f));
+    // eta[k] = 0.87f * (float)k / ((float)d1 - 1.0f);
+    t_mean += (float)(fdnpath[tap].delay);
+  }
+  // if feed forward model, then truncate delays:
+  if(!feedback) {
+    for(auto& path : fdnpath)
+      path.delay++;
+  }
+  t_mean /= (float)std::max(1u, fdnorder_);
+  float g(0.0f);
+  switch(gainmethod) {
+  case fdn_t::original:
+    g = expf(-4.2f * t_min / t60);
+    break;
+  case fdn_t::mean:
+    g = expf(-4.2f * t_mean / t60);
+    break;
+  case fdn_t::schroeder:
+    g = powf(10.0f, -3.0f * t_mean / t60);
+    break;
+  }
+  prefilt0.set_lp(g, damping);
+  prefilt1.set_lp(g, damping);
+  // set rotation:
+  for(uint32_t tap = 0; tap < fdnorder_; ++tap) {
+    // set reflection filters:
+    fdnpath[tap].reflection.set_lp(g, damping);
+    float laz = 0.0f;
+    if(fdnorder_ > 1)
+      laz = -daz + 2.0f * daz * (float)tap / (float)(fdnorder_ - 1);
+    fdnpath[tap].rotation.set_rotation(laz, TASCAR::posf_t(0, 0, 1));
+    TASCAR::quaternion_t q;
+    q.set_rotation(0.5f * daz * (float)(tap & 1) - 0.5f * daz,
+                   TASCAR::posf_t(0, 1, 0));
+    fdnpath[tap].rotation.rmul(q);
+    q.set_rotation(0.125f * daz * (float)(tap % 3) - 0.25f * daz,
+                   TASCAR::posf_t(1, 0, 0));
+    fdnpath[tap].rotation.rmul(q);
+  }
+  // set feedback matrix:
+  if(fdnorder_ > 1) {
+    TASCAR::fft_t fft(fdnorder_);
+    TASCAR::spec_t eigenv(fdnorder_ / 2 + 1);
+    for(uint32_t k = 0; k < eigenv.n_; ++k)
+      eigenv[k] = std::exp(i_f * TASCAR_2PIf *
+                           powf((float)k / (0.5f * (float)fdnorder_), 2.0f));
+    ;
+    fft.execute(eigenv);
+    for(uint32_t itap = 0; itap < fdnorder_; ++itap)
+      for(uint32_t otap = 0; otap < fdnorder_; ++otap)
+        feedbackmat[fdnorder_ * itap + otap] =
+            fft.w[(otap + fdnorder_ - itap) % fdnorder_];
+  } else {
+    feedbackmat[0] = 1.0;
+  }
+}
+
 reflectionfilter_t::reflectionfilter_t()
 {
   sy.set_zero();

manual/manual.tex

@@ -1470,6 +1470,9 @@ \subsection{Reflectors: {\tt <face .../>} and {\tt <facegroup .../>} elements}\l
 Polygons meshes are flattened by a projection on a plane which is
 orthogonal to the polygon normal vector.
 
+\subsubsection{Scattering}
+
+\tascar{} contains a very simple scattering model. For each reflector the amount of scattering can be controlled using the \attr{scattering} attribute. This is added to the diffuse sound field model path. By default, the spatial dispersion of the scattering is reproduced by the receiver's decorrelation stage, if enabled. Optionally, the scattering can be rendered explicitly using additional virtual sound sources added to the diffuse sound field model. This can be enabled in each receiver by setting the attribute \attr{scatterreflections} to a number greater than zero. Reasonable values with a reasonable trade-off between computational effort and spatial dispersion are 4 to 8.
 
 \subsection{Obstacles: {\tt <obstacle .../>} element}\index{obstacle}\label{sec:obstacle}
 

scripts/sdm/gen_ir.tsc

@@ -1,12 +1,13 @@
 <?xml version="1.0"?>
 <session>
   <scene>
-    <receiver type="amb1h1v" scatterreflections="5" channelorder="FuMa"/>
+    <receiver type="amb1h1v" scatterreflections="4" scatterspread="20" channelorder="FuMa"/>
     <source>
       <position>0 2 0 0</position>
       <sound/>
     </source>
     <facegroup name="walls" shoebox="7 6 3" scattering="0.1" mute="true"/>
-    <face name="w1" height="3" width="3" dorientation="0 0 0" dlocation="-2 -1 -1"/>
+    <face name="w1" height="3" width="3" dorientation="0 0 0" dlocation="-2 -1 -1" scattering="0.5"/>
+    <reverb type="simplefdn" volumetric="7 6 3" forwardstages="1"/>
   </scene>
 </session>