LCAV / pyroomacoustics

Pyroomacoustics is a package for audio signal processing for indoor applications. It was developed as a fast prototyping platform for beamforming algorithms in indoor scenarios.
https://pyroomacoustics.readthedocs.io
MIT License
1.36k stars 419 forks source link

Help me understand energy scaling #352

Open Kanafani80 opened 1 month ago

Kanafani80 commented 1 month ago

Hi there! Thank you for this great repo.

I am having trouble understanding the following energy scaling that I see in a couple of places, for example in room.cpp when scattered energy is computed:

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) ;

This seems different from the lambert cosine formula (which is indeed applied earlier). In the code above, r_sq will generally be much larger than the mic radius squared, so it seems p_hit will be tiny, which will make energy blow up.

What is the physical explanation of this? I can’t find the equivalent formula in either the wayverb or Schroeder dissertation.

thanks!

fakufaku commented 4 weeks ago

Hi @Kanafani80 , you are doing a pretty deep dive into the code! 🤿

This is an excellent question. As is, p_hit is hard to understand and it would deserve to have some comment in the code explaining it. I had to dig myself to recall why it is like this... 🙈

This term is necessary to match the scales of the image source model and the ray tracing, as explained in Dirk Schroeder's thesis, Section 5.4 (in particular, page 75, equation 5.54). It is not called p_hit, but the value is the same, i.e. r ** 2 * (1 - cos(gamma)). The cos(gamma) can be computed as sqrt(1 - mic_dist **2 / r ** 2) (definition of cosine + pythagoras, see fig 5.20 for definition of gamma and r).

fakufaku commented 4 weeks ago

Added some comments in the code here and here.

Kanafani80 commented 3 weeks ago

Thanks for the reply! I had not read that far into the dissertation ;)

The scale matching idea makes sense to me. However, it seems that your calculation is using the cumulative distance travelled by the ray up to this point: *double r_sq = double(travel_dist_at_mic) travel_dist_at_mic;** If we go by geometry and Pythagoras, shouldn't the distance between the secondary source and the receiver be used instead (so, hop_dist instead of travel_dist_at_mic)?

Thanks again!

fakufaku commented 3 weeks ago

In my understanding, the quantity relates the total energy of the source, to the fraction that impacts the receiver. Thus, it should be computed according to the total distance traveled, not just from the last reflection.