Open an interactive online version by clicking the badge Binder badge or download the notebook.

Assignment 4: Early Reflections and Wave Domain Modeling#

Sebastian J. Schlecht (1), Nils Meyer-Kahlen (2)
Notebook: Cristóbal Andrade (1)

(1) Friedrich-Alexander-Universität Erlangen-Nürnberg (2) Aalto University

Contact: sebastian.schlecht@fau.de, cristobal.andrade@fau.de

This assignment is divided into two parts. In the first part, you will develop a simple yet effective model of early reflections in a shoebox-shaped room. In the second part, you will model the low-frequency modes of a room with rigid boundaries.

Duration: 12 Hours

Requirements: Basics of digital signal processing, room acoustics and wave equation.

References
Allen, J. B., & Berkley, D. A. (1979). Image method for efficiently simulating small‐room acoustics. The Journal of the Acoustical Society of America, 65(4), 943–950.
Bilbao, S. (2009). Numerical Sound Synthesis: Finite Difference Schemes and Simulation in Musical Acoustics. Wiley.

Dependencies pip install numpy scipy matplotlib soundfile pyfar sofar jupyter ipython pillow watermark

[ ]:

import soundfile as sf
import pyfar as pf
import numpy as np
import pooch
import matplotlib.pyplot as plt
import os
import time
from IPython.display import Audio, display
from scipy.fft import fft, fftfreq

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_py = 'https://github.com/pyfar/open-educational-resources/raw/main/courses/Virtual_Acoustics_Lab_FAU/Assignment4'
url_large_files = 'https://github.com/pyfar/files/raw/main/education/VAL_FAU'
url_hrtfs = 'https://github.com/pyfar/files/raw/main/education/VAR_TUB'

# Get current working directory (where the notebook was started)
notebook_dir = os.getcwd()

# Create a Pooch object using that directory
my_pooch_py = pooch.create(
    path=notebook_dir,
    base_url=url_py,
    registry={
        "Scene.py" : None,
        "Binaural_DSP.py" : None,
        "Early_DSP.py" : None,
        "Direct_DSP.py" : None,
        "DSP.py" : None,
        "WaveDomain_DSP.py" : None,
        "Config.py" : None,
        "BlockConvolver_DSP.py" : None,
        "hrirsDiffuseFieldEQ.py" : None,
        "air_absorption_iso.py" : None,
        "cart2sph_vec.py" : None,
        "VariableDelay_DSP.py" : None,
        "VariableSOS_DSP.py" : None,
        "design_one_pole_filter.py" : None,
        "Propagation_DSP.py" : None,
    }
)

my_pooch_hrtfs = pooch.create(
    path=notebook_dir,
    base_url=url_hrtfs,
    registry={
        "FABIAN_HRIR_measured_HATO_0.sofa" : None
    }
)


# Download all files
for fname in my_pooch_py.registry:
    fpath = my_pooch_py.fetch(fname)
    print(f"Downloaded: {fpath}")

for fname in my_pooch_hrtfs.registry:
    fpath = my_pooch_hrtfs.fetch(fname)
    print(f"Downloaded: {fpath}")

1. Early Reflections Modeling#

The main steps involved in early reflections modeling in the assigment are:

  1. Computing the image source indices

  2. Computing the image sources depending on the source position

  3. Computing the relative position to the listener position

  4. Rendering the propagation path ofeach image source including air absorptions.

[ ]:

from Config import Config
from Scene import Scene
from Early_DSP import Early_DSP
from Direct_DSP import Direct_DSP
from Binaural_DSP import Binaural_DSP
from WaveDomain_DSP import WaveDomain_DSP
from DSP import mag2db

# Set up matplotlib for interactive plotting and animations
%matplotlib inline
plt.style.use('default')
plt.rcParams['figure.figsize'] = (12, 8)

Configuration and Impulse signal generation#

[ ]:

# Define Config = Fixed during runtime
config = Config()
config.Early_maxImageSourceOrder = 3
config.roomSize = [10, 7]

scene = Scene()

# Create impulse signal
numBlocks = 30
signalLength = config.blockSize * numBlocks
signal = np.zeros((signalLength, 1))
signal[100, 0] = 1
time = np.arange(signalLength) / config.fs

1.1 TASK: Prepare Image Sources#

Implement Early_DSP.prepareImageSources(). This method is called once when the Early_DSP object is initialized and is responsible for generating a list of image source (IS) indices. Specifically, it should create the image source indices up to maxImageSourceOrder.

Compute imageSourceList a numberOfImageSources x 3 matrix with image room indices. For example: [2 -1 0] means:

  • two image rooms further in positive X direction

  • one image room further in negative Y direction

  • no displacement in Z direction

[ ]:

def prepareImageSources(self, maxImageSourceOrder):
    """
    Create the image source indices up to maxImageSourceOrder.

    Args:
        maxImageSourceOrder (int): Maximum image source order

    Returns:
        np.ndarray: Image source list matrix
    """

    # YOUR CODE HERE
    raise NotImplementedError()
    return imageSourceList

Early_DSP.prepareImageSources = prepareImageSources

early = Early_DSP(config)
print("Image source list:")
print(early.imageSourceList)

1.2 TASK: Calculate Image Source Positions#

Implement Early_DSP.setPosition(). The final IS positions should be calculated in setPosition(), whenever a new sourcePosition or listenerPosition is set during runtime.

Briefly describe the 3D image source plot.

[ ]:

def setPosition(self, sourcePosition, listenerPosition):
    """
    Compute image source positions.

    Args:
        sourcePosition (np.ndarray): Source position in Cartesian coordinates
        listenerPosition (np.ndarray): Listener position in Cartesian coordinates
    """
    # Compute image source locations
    # YOUR CODE HERE
    raise NotImplementedError()
Early_DSP.setPosition = setPosition

early.setPosition(scene.sourcePosition, scene.listenerPosition)

%matplotlib widget
early.plotImageSources()

1.3 TASK: Implement and Test Early_DSP.process Routine#

Process the signal thorugh early reflections. HINT: See propagation_DSP

[ ]:

def process(self, input_signal):
    """
    Process input signal through early reflections.

    Args:
        input_signal (np.ndarray): Input signal

    Returns:
        np.ndarray: Output signal with one channel per image source
    """
    # YOUR CODE HERE
    raise NotImplementedError()
    return output

Early_DSP.process = process
# Initialize the output signal buffer
output = np.zeros((signalLength, early.numberOfOutputs))

# Perform block by block processing of the signal
blockSize = config.blockSize
for itBlock in range(numBlocks):
    blockIndex = slice(itBlock * blockSize, (itBlock + 1) * blockSize)
    block = signal[blockIndex, :]
    output[blockIndex, :] = early.process(block)

outputMono = np.sum(output, axis=1)

# Save audio file
sf.write('assignment_Early_DSP.wav', outputMono, config.fs)

1.4 TASK: Analyze Time-Domain Response and Animation#

Describe the time-domain response and the 2D animation. Also compare the 2D animation with the WaveDomain_DSP animation.

[ ]:

%matplotlib inline
# Also compare the 2D animation with the WaveDomain_DSP animation
plt.figure(figsize=(12, 6))
plt.plot(time, output)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.title('Early Reflections Response')
plt.grid(True)
plt.show()

early.animationImageSources(use_animation=True, save_gif=True) #Saved as image_sources_animation.gif

plt.show()

1.5 TASK: Combine direct sound, early reflections and binaural rendering#

  1. Write a full pipeline processing loop the input signal with block processing

  2. Use Direct_DSP, Binaural_DSP, Early_DSP

  3. The scene is static so you can initialize all DOAs before the loop

  4. Implement and use Early_DSP.getDoa()

[ ]:

# Load input signal
inputSignal = pf.signals.files.speech(voice='male', sampling_rate=config.fs).time.transpose()

# try:
#     inputSignal, fs = sf.read('december_tour.wav')
#     if len(inputSignal.shape) == 1:
#         inputSignal = inputSignal.reshape(-1, 1)
# except FileNotFoundError:
#     print("Warning: december_tour.wav not found, using impulse signal instead")
#     inputSignal = signal


direct = Direct_DSP(config)
binaural = Binaural_DSP(config, early.numberOfOutputs + 1)
direct.setScene(scene)

doaDirect = direct.getDoa()
doaReflections = early.getDoa()
doaAll = np.vstack([doaDirect, doaReflections])
binaural.set_doa(doaAll[:, 0], doaAll[:, 1])

outputBinaural = np.zeros((len(inputSignal), 2))
numBlocks = len(inputSignal) // blockSize

# Implement for loop for block processing
# Combine direct sound, early reflections and binaural rendering
# YOUR CODE HERE
raise NotImplementedError()

sf.write('assignment_Early_DSP_binaural.wav', outputBinaural, config.fs)
# You can listen to your ouput signal the following widget. Please check your sound set-up volume before playing.
display(Audio(outputBinaural.transpose(), rate=config.fs))

1.6 TASK: Analyze Sonic Results with Different Room Sizes#

Describe the sonic result. Does it sound like a room?

  1. Try at least two different room sizes

  2. Save sound examples with telling names

[ ]:

# Dummy room sizes
room_sizes = [[3, 2]]
# YOUR CODE HERE
raise NotImplementedError()
direct = Direct_DSP(config)
binaural = Binaural_DSP(config, early.numberOfOutputs + 1)
direct.setScene(scene)

# Implement for loop for block processing as in the previous task, now also with different room sizes
# YOUR CODE HERE
raise NotImplementedError()

[ ]:

# Write your description here
# YOUR CODE HERE
raise NotImplementedError()

2. Wave Domain Modeling#

In this assignment, you will create a basic modeling of the low-frequency modes in a room with rigid boundaries. To simplify the implementation only 2D simulation are considered. The main steps involved are:

  1. Creating an input signal

  2. Implementing a FDTD scheme

  3. Manipulating the room shape

To keep the processing reasonable, the maximum valid frequency is set at 1000 Hz and sampling rate is ~8400 Hz. Please start with watching a small excerpt of this tutorial video (5:30 - 12:00): www.youtube.com/watch?v=xgJJwmrX568

In the following you implement the same time step function.

Define Config and Setup FDTD scheme#

[ ]:

# Define Config = Fixed during runtime
config = Config()
config.roomSize = [10, 7, 4]
scene = Scene()

2.1 TASK: Implement Wave Domain Processing#

  1. Add the FDTD update to WaveDomain_DSP process. Hint: Use the following equation:

  • FDTD update including boundary nodes

    \[u_{l,m}^{n+1}=\left(2-\frac{{}_{1}}{2}K_{l,m}\right)u_{l,m}^{n}+\textstyle{\frac12}\bigl(u_{l+1,m}^{n}+u_{l-1,m}^{n}+u_{l,m+1}^{n}+u_{l,m-1}^{n}\bigr)-u_{l,m}^{n-1}\]

  • \(K_{l,m}\) number of interior neighbours ( \(0<K_{l,m} \leq 4 )\))

  1. Inject source (input) and extract receiver (output)

  2. Check 2D animation (for development animation can be disabled with property draw = false)

[ ]:

def process(self, input_signal):
    """
    Process input signal through FDTD simulation

    Args:
        input_signal: Input signal array [t, 1]

    Returns:
        output: Output signal array [t, 1]
    """
    assert input_signal.shape[1] == self.numberOfInputs, 'Number of Inputs incorrect.'

    # Initialize state arrays
    Nxy = np.ceil(np.array(self.roomSize) / self.dx).astype(int) + 2  # number of points in x-dir and y-dir
    u0 = np.zeros(Nxy, dtype=np.float64)  # values at next time position
    u1 = np.zeros(Nxy, dtype=np.float64)  # values at current time position
    u2 = np.zeros(Nxy, dtype=np.float64)  # values at previous time position

    xv = (np.arange(Nxy[0]) * self.dx) - 0.5 * self.dx  # x-sampling points
    yv = (np.arange(Nxy[1]) * self.dx) - 0.5 * self.dx  # y-sampling points

    # Create interior mask
    inMask = self.createInteriorMask(xv, yv)
    KMap = self.computeNumberOfNeighbors(inMask)

    # Set grid forcing grid point indices and ensure that it is inside the domain
    inxy = np.round(self.xyIn / self.dx + 1.5).astype(int)
    assert inMask[inxy[0], inxy[1]], "Source position is outside domain"

    # Set grid extraction grid point indices and ensure that it is inside the domain
    outxy = np.round(self.xyOut / self.dx + 1.5).astype(int)
    assert inMask[outxy[0], outxy[1]], "Receiver position is outside domain"

    # FDTD Simulation
    #start_time = time.time()
    output = np.zeros_like(input_signal[:, 0])  # match MATLAB: output = 0*input
    Nt = len(input_signal)

    for nt in range(Nt):
        iX = slice(1, Nxy[0]-1)
        iY = slice(1, Nxy[1]-1)

        # Apply FDTD update including boundary nodes
        # YOUR CODE HERE
        raise NotImplementedError()

        # Enforce the boundary condition by setting the values outside the exterior grid points to zero
        # YOUR CODE HERE
        raise NotImplementedError()

        # Add input and record output at the grid indices inxy and outxy;
        # YOUR CODE HERE
        raise NotImplementedError()

        # Step forward in time
        u2 = u1.copy()
        u1 = u0.copy()

        # Plotting and Recording
        if (self.draw or self.record) and nt % 5 == 0:  # Only every 10th frame
            self.drawRoom(u2, inMask, xv, yv)
            if self.record:
                self.recordVideo(self.videoFileName)

        #elapsed = time.time() - start_time
        #print(f'\rProgress: nt={nt+1} out of {Nt} time-steps, {elapsed:.2f} s elapsed', end='')

    if self.record:
        self.recordVideo(None)  # close video writing

    return output.reshape(-1, 1)

WaveDomain_DSP.process = process

# Setup FDTD scheme
waveDomain = WaveDomain_DSP(config)
waveDomain.setPosition(scene.sourcePosition, scene.listenerPosition)

# Create simple excitation signal
duration = 0.2  # seconds
fs = 1 / waveDomain.dt
Nt = int(np.ceil(duration * fs))  # number of time-steps to compute
time = np.arange(Nt) / fs
input_signal = np.zeros((Nt, 1))
input_signal[0] = 1

waveDomain.draw = False
waveDomain.videoFileName = 'ftdtAnimation.gif' # Animation is saved as gif
waveDomain.record = True  # record animation
output = waveDomain.process(input_signal)

2.2 Describe Time-domain Response#

  1. Plot and describe problems with impulse signal (if not bandlimited)

  2. Describe the numerical dispersion in the time-domain animation

[ ]:

# Plot Output
# YOUR CODE HERE
raise NotImplementedError()

# Write your description here
# YOUR CODE HERE
raise NotImplementedError()

2.3 TASK: Improve the excitation signal#

  1. Create a bandlimited one period sine excitation

  2. Describe the improved result

[ ]:

# Setup FDTD scheme
waveDomain = WaveDomain_DSP(config)
waveDomain.setPosition(scene.sourcePosition, scene.listenerPosition)
# YOUR CODE HERE
raise NotImplementedError()

waveDomain.draw = False  # show animation
waveDomain.videoFileName = 'ftdtAnimationBandlimited.gif' # Animation is saved as gif
waveDomain.record = True  # record animation
output = waveDomain.process(input_signal)

plt.figure(figsize=(12, 6))
plt.plot(time, output)
plt.xlabel('Time in samples')
plt.ylabel('Amplitude')
plt.title('Bandlimited Excitation Response - Time Domain')
plt.grid(True)
plt.show()

2.4 Compare time-domain response and animation with image source model#

[ ]:

# Write your description here
# YOUR CODE HERE
raise NotImplementedError()

2.5 TASK: Modal frequencies#

  1. Determine lowest ~5 room mode frequency (e.g., https://amcoustics.com/tools/amroc)

  2. Plot the magnitude response of the impulse response (e.g. up to 50 Hz)

  3. Plot the room mode frequencies (e.g., use function xline)

  4. Describe the plot

[ ]:

# YOUR CODE HERE
raise NotImplementedError()

2.6 TASK: Modal excitation#

  1. Determine a room mode frequency (e.g., https://amcoustics.com/tools/amroc)

  2. Create an input signal at this frequency

  3. Describe the output

[ ]:

waveDomain.draw = False  # show animation
waveDomain.record = False  # record animation - set to true when ready

# YOUR CODE HERE
raise NotImplementedError()

2.7 TASK: Manipulate the room shape#

Change the inMask definition in createInteriorMask to create a new room shape. Observe again the resulting wave propagation animation. Submit a video of your new room shape and give a description here.

[ ]:

waveDomain.resetVideoRecording()
# 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,matplotlib,soundfile,pyfar,sofar,jupyter,ipython,pillow,watermark

Python implementation: CPython
Python version       : 3.12.9
IPython version      : 8.12.3

numpy     : 1.26.4
scipy     : 1.15.2
matplotlib: 3.10.1
soundfile : 0.13.1
pyfar     : 0.7.1
sofar     : 1.2.1
jupyter   : 1.1.1
ipython   : not installed
pillow    : not installed
watermark : 2.5.0

Compiler    : Clang 13.0.0 (clang-1300.0.29.30)
OS          : Darwin
Release     : 24.5.0
Machine     : arm64
Processor   : arm
CPU cores   : 12
Architecture: 64bit