This question is a part of a more general question the answer of which I don't know - How to apply a filter in the freq domain and then convert the filtered signal back to the time domain? Well, I partially googled the answer that I need to
- convert the signal in FFT
- multiply by the filter
- convert back to the time domain
but I'm not entirely sure if I've applied this idea correctly to atmospheric absorption filtering (see below).
I'm also not sure whether I can do the filtering entirely in the time domain (as suggested here by convolution?) somehow alleviating the need to switch back and forth to the freq domain.
Input: sound wave (gunshot sound pressure subtracted atmospheric pressure) and atmospheric conditions. Output: the same sound wave attenuated with an atmospheric absorption filter.
I'm using python-acoustics to find the atmospheric attenuation coefficient as described in Engineering Acoustics/Outdoor Sound Propagation.
I've came up with the following code:
def atmosphericAttenuation(signal, distance, Fs, **kwargs):
"""
Apply atmospheric absorption to the `signal` for all its FFT frequencies.
It does not account for the geometrical attenuation.
Parameters
----------
signal - a pressure waveform (time domain)
distance - the travelled distance, m
Fs - sampling frequency of the `signal`, Hz
kwargs - passed to `Atmosphere` class
Returns
-------
signal_attenuated - attenuated signal in the original time domain
"""
# pip install acoustics
from acoustics.atmosphere import Atmosphere
atm = Atmosphere(**kwargs)
signal_rfft = np.fft.rfft(signal)
freq = np.fft.rfftfreq(n=len(signal), d=1. / Fs)
# compute the attenuation coefficient for each frequency
a_coef = atm.attenuation_coefficient(freq)
# (option 2) signal_rfft *= 10 ** (-a_coef * distance / 20)
signal_rfft *= np.exp(-a_coef * distance)
signal_attenuated = np.fft.irfft(signal_rfft)
return signal_attenuated
Am I doing it right? Which one is correct:
signal_rfft *= np.exp(-a_coef * distance)<- $P(r) = P(0) \exp (-\alpha r)$signal_rfft *= 10 ** (-a_coef * distance / 20)<- $A_a = -20 \log_{10} \frac{P(r)}{P(0)} = \alpha r$
If neither, please describe how it should be done. Thank you.
ir = atm.impulse_response(distance=distance, fs=Fs, ntaps=len(signal)); signal_attenuated = np.convolve(signal, ir, mode='same')? Indeed, this looks better. – dizcza Jun 28 '21 at 12:14