or
download
the notebook.
Assignment 1: Sound Propagation Under Free Field Conditions#
(1) Friedrich-Alexander-Universität Erlangen-Nürnberg (2) Aalto University
Contact: sebastian.schlecht@fau.de, cristobal.andrade@fau.de
The rendered scene in this assignment is a emergency car driving by a static listener. You will learn to render theses components: stereo panning, propagation delay, the air absorption, distance attenuation. This assignment also shows you the basic workflow for rendering dynamic time-varying scene using block based processing and interpolation.
Duration: 12 Hours
Requirements: Basics of Digital Filtering
Dependencies pip install matplotlib==3.10.3 numpy==2.2.6 pyfar==0.7.1 scipy==1.15.3 watermark
[1]:
import pyfar as pf
import numpy as np
import numpy.matlib as npm
import matplotlib.pyplot as plt
from IPython.display import Audio, display
import pooch
import os
Let’s start by downloading all the necessary files for running this assignment.
[ ]:
# Leave this as it is: This is the URL from which the data will be downloaded
url = 'https://github.com/pyfar/open-educational-resources/tree/main/docs/oer/courses/Virtual_Acoustics_Lab_FAU/Assignment1'
# Get current working directory (where the notebook was started)
notebook_dir = os.getcwd()
# Create a Pooch object using that directory
my_pooch = pooch.create(
path=notebook_dir,
base_url=url, # Change this to your actual URL
registry={
"Scene.py" : None,
"VariableDelay_DSP.py" : None,
"VariableSOS_DSP.py" : None,
"DSP.py" : None,
"Config.py" : None,
"medium_attenuation.py" : None,
}
)
# Download all files
for fname in my_pooch.registry:
fpath = my_pooch.fetch(fname)
print(f"Downloaded: {fpath}")
Now let’s import all the necessary files.
[3]:
from Config import Config
from Scene import getDriveByScene
from medium_attenuation import air_attenuation
from VariableDelay_DSP import VariableDelay_DSP
from VariableSOS_DSP import VariableSOS_DSP
from DSP import m2smp, call112
%matplotlib inline
config = Config()
Create source signal#
First, we define the signal length and then create a source signal. The source signal is an approximate rendering sound of an emergency vehicle.
[ ]:
signalLength = 6 * config.fs; # samples corresponding to 6 seconds
signal = call112(signalLength, config.fs)
# You can listen with the following widget. Please check your sound set-up volume before playing.
display(Audio(signal.time, rate=config.fs))
1 TASK: Compute panning curve#
Define power-preserving panning gains gainLeft and gainRight according to the direction of arrival (doa). Recall how azimuth is defined in the appropriate coordinate.
Hint: \(g_l^2(\phi) + g_r^2(\phi) = const\)
Plot the panning gains
Plot the summed power
[ ]:
doa_azimuth = np.linspace(np.pi/2, -np.pi/2,100)
# YOUR CODE HERE
raise NotImplementedError()
2 TASK: Pan the source signal to 25 degrees on the left#
Apply the panning gains for the source signal
Write the result to a 2-channel output called signalPanned of size(signalLength, 2)
Listen to the signalPanned
Hint: use deg2rad for conversion
[ ]:
# YOUR CODE HERE
raise NotImplementedError()
# You can listen to your ouput signal the following widget. Please check your sound set-up volume before playing.
#display(Audio(singalPanned, rate=config.fs))
Create air absorption filters#
Next, we implement the subcomponents for the transmission through air. For this, we compute the absorption depending on atmospheric parameters at each frequency.
[7]:
f = np.linspace(0, config.fs/2, 2 ** 10)
T = config.temperature
hr = config.relativeHumidity
alpha_iso = air_attenuation(T, f, hr)
3 TASK: plot the air absorption#
Plot on a logarithmic frequency scale
Add the correct labels
[ ]:
# YOUR CODE HERE
raise NotImplementedError()
4 TASK: Decribe the main features of the plot here#
[9]:
# Write your description here
# YOUR CODE HERE
raise NotImplementedError()
5 TASK: compute the air absorption IIR filters#
Next, we have to convert the air absorption specification into a digital filter which can be applied efficiently to the source signal. A good representation of a digital filter is the second-order section (SOS). We use this framework here, although this first filter is only a first order filter.
5.1 Implement the function designOnePoleFilter#
Solve for a1 and b0 (e.g. by plugging in HDc and HNyq and rearranging)
Implement the filter coefficients into second-order sections (SOSs) [b0; b1; b2; a0; a1; a2], where a0 = 1.
Vectorize for multiple filters. For example: you want to compute the filter coefficients for two filters at the same time. Let’s say the input values are HDc = [1, 0.9] and HNyq = [0.7, 0.5]. The output shoud be a matrix of sos coefficient of size [6 x 2]. The first column of sos is the first filter corresponding to gains 1 and 0.7. The second column is the second filter corresponding to gains 0.9 and 0.5.
[10]:
def designOnePoleFilter(HDc, HNyq, fs):
#designOnePoleFilter - compute one pole filter
#
# Inputs:
# HDc - Linear magnitude at DC of size [1, number of filters]
# HNyq - Linear magnitude at Nyquist of size [1, number of filters]
#
# Outputs:
# sos - sos filters of size [6 x number of filters]
numFilters = np.size(HDc)
sos = pf.FilterSOS(np.zeros((numFilters, 6)), sampling_rate=fs)
sos.coefficients[0,:, 3] = 1
# YOUR CODE HERE
raise NotImplementedError()
return sos
5.2 Compute Filters#
Compute filter for distance of 1 meter
Plot filter magnitude response and ideal response
Compute filter for distance of 7 meter
Plot filter magnitude response and ideal response
Hint: You can use pyfar’s sos.process()
[ ]:
# YOUR CODE HERE
raise NotImplementedError()
6 TASK: describe and comment the plots#
WRITE YOUR ANSWER HERE
Note:#
The simple one-pole filter design is quite inaccurate for larger distances, and it can be highly improved with more dedicated filter designs for second or third order filters. For further reference, please see Kates, J. M. & Brandewie, E. J. Adding air absorption to simulated room acoustic models. J Acoust Soc Am 148, EL408-EL413 (2020).
7 TASK: Filter ‘signal’ according to 10m distance of air absorption.#
[ ]:
# YOUR CODE HERE
raise NotImplementedError()
8 TASK: Describe briefly what you see and hear.#
WRITE YOUR ANSWER HERE
Now we make the scene dynamic#
Define Scene = Source position is dynamic during runtime
The Scene struct contains parameters of the acoustic scene, such as source and listener positions. Typically such scene descriptions can vary over time, where the parameters can change at each block, e.g., every 256 samples. Here, we create a scene with 1000 blocks.
[13]:
scenes = getDriveByScene(1000)
9 TASK: Implement Processing#
The signal processing is performed over signal blocks, which is typical for real-time and time-varying audio. Each block provides 256 samples of audio and a few parameters from the scene description. Your task is to implement direct sound processing, including air absorption and the stereo panning for the right impression of the sound direction. You can find more detailed tasks below.
We initialize the processing blocks. Please read the DSP block implementation and familiarize yourself with the corresponding python documentation.
Example usage can be seen in
VariableSOS_DSP
VariableDelay_DSP
[14]:
variablePropagationDelay = VariableDelay_DSP(1, config)
variableAirAbsorption = VariableSOS_DSP(config)
Implement the processing of the direct sound
Compute start and end distance
Compute and process distance gain (hint: 1/r law)
Set and process propagation delay from distance (hint: VariableDelay_DSP)
Compute air absorption filter coefficients
Set and process air absorption filter (hint: VariableSOS_DSP)
Pan the source
[15]:
# Initialize some variables that might be helpful
lastRelativePosition = np.array([100, 0, 0])
outputSignal = np.zeros((signalLength, 2))
for itBlock in range(len(scenes)):
# Get sample block and current scene instance
# YOUR CODE HERE
raise NotImplementedError()
# Calculate distance gain, delay, air absorption, and panning
# YOUR CODE HERE
raise NotImplementedError()
# Apply panning gains to the processed block and store in output signal
# YOUR CODE HERE
raise NotImplementedError()
10 TASK: Describe the sonic result#
Describe features such as Doppler shift, air absorption, and direction. Please also include the .wav file in your submission (make sure to not clip the audio file).
Use Audio Widget for playback
Use audiowrite for writing the assignment_Direct.wav file
[ ]:
# YOUR CODE HERE
raise NotImplementedError()
# You can listen to your ouput signal the following widget. Please check your sound set-up volume before playing.
#display(Audio(outputSignal.transpose(), rate=config.fs))
11 TASK: Describe the time-domain response#
How does it conform to your expectation of direct sound processing? Describe features such as periodicity shift and envelope.
Use plot
Include correct axis labels
[ ]:
# YOUR CODE HERE
raise NotImplementedError()
12 TASK : Describe the spectrogram#
How does it conform to your expectation of direct sound processing? Describe features such as Doppler shift and air absorption.
use spectrogram
scale the frequency axis logarithmically, see XScale property
[ ]:
# YOUR CODE HERE
raise NotImplementedError()
License notice#
This notebook is licensed under CC BY 4.0
Watermark#
The following watermark might help others to install specific package versions that might be required to run the notebook.
[1]:
%load_ext watermark
%watermark -v -m -p numpy,scipy,pyfar,sofar,nbgrader,watermark
Python implementation: CPython
Python version : 3.12.9
IPython version : 8.12.3
numpy : 2.1.3
scipy : 1.15.2
pyfar : 0.7.1
sofar : 1.2.1
nbgrader : 0.9.5
watermark: 2.5.0
Compiler : Clang 13.0.0 (clang-1300.0.29.30)
OS : Darwin
Release : 24.4.0
Machine : arm64
Processor : arm
CPU cores : 12
Architecture: 64bit