How to use the reflection module? #249

Activity

1. FRidh commented on Jun 29, 2021

   

   FRidh

   on Jun 29, 2021

   Member

   More actions

   The reflection factor returned consists of the positive frequencies. You could indeed do a convolution if you create an impulse response out of the spectrum. Using np.fft.irfft should give you a usable impulse response. Alternatively, when stacking positive and negative frequencies you can use np.fft.ifft https://github.com/FRidh/auraliser/blob/master/auraliser/propagation.py#L38.
2. dizcza commented on Jun 30, 2021

   

   dizcza

   on Jun 30, 2021

   Author

   More actions

   Thanks for the feedback.

   I have several comments.

   1. Isn't ir_reflection function you mentioned equivalent to impulse_response for any spectrum? I mean, since you assume that the frequencies are mirrored, they should give identical results.

   2. Bounday should not accept a zero frequency or? Therefore, I cannot use np.fft.rfftfreq(ntaps, d=1. / Fs) as the input, because it starts with zero. How to alleviate this? Currently, I set a very small number, say 1e-10, as the first frequency component.

   3. reflection_factor_spherical_wave is numerically unstable for large frequencies.

      b = Boundary(np.arange(1, 20_000), flow_resistivity=2e5, angle=np.deg2rad(60), reflection_model='spherical', distance=10)
      b.plot_reflection_factor()

      The issue is in this line
      

      python-acoustics/acoustics/reflection.py

      Line 338 in fbc8745

      F = 1.0 - 1j * np.sqrt(np.pi) * w * np.exp(-w**2.0) * erfc(1j * w)

      I'll come up with a PR that substitutes np.exp(-w**2.0) * erfc(1j * w) by its numerical stable equivalent erfcx(1j * w).

   4. What to do with the phase information? I guess, my current script accounts for the amplitude change only.

   import matplotlib.pyplot as plt
   import numpy as np
   
   from acoustics.reflection import Boundary
   from acoustics.signal import impulse_response_real_even
   
   
   def friedlanderMW(t_interval, ta, amplitude, tau=0.05):
       """
       Friedlander model to calculate a muzzle blast wave.
   
       Parameters
       ----------
       ta -- Arrival time, s (float)
       tau -- Positive phase duration, s (float)
       """
       x = (t_interval - ta) / tau
       Pmb = amplitude * (1 - x) * np.exp(-x)
       Pmb[t_interval <= ta] = 0
       return Pmb
   
   
   Fs = 96000
   duration = 0.5  # s
   distance = 1000  # m
   t_interval = np.linspace(0, duration, int(Fs * duration))
   pmb = friedlanderMW(t_interval, ta=0.1, amplitude=300, tau=0.05)
   plt.plot(t_interval, pmb, label='orig')
   
   freq = np.fft.rfftfreq(len(pmb), d=1. / Fs)
   freq[0] = 1e-10  # avoid zero
   b = Boundary(freq, flow_resistivity=2e5, angle=np.deg2rad(60), reflection_model='spherical', distance=1)
   rf = b.reflection_factor.squeeze()  # why is it 2d?
   ir = impulse_response_real_even(rf, ntaps=len(pmb))
   
   pmb = np.convolve(pmb, ir, mode='same')
   plt.plot(t_interval, pmb, label='reflected')
   plt.legend()
   plt.show()