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

Assignment 5: Artificial Reverberation#

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

In this assignment, your task is to implement a late reverberation algorithm based on feedback delay networks (FDNs).

Duration: 12 Hours

Requirements: DSP

Dependencies pip install numpy scipy matplotlib soundfile pyfar jupyter pooch watermark

[ ]:

import numpy as np
import matplotlib.pyplot as plt
import pooch
import os
from scipy.signal import freqz
from scipy.io import wavfile
from scipy.linalg import hadamard

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/courses/Virtual_Acoustics_Lab_FAU/Assignment5'

# 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={
        "SOSarray_DSP.py" : None,
        "Config.py" : None,
        "DSP.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.

[ ]:

from Config import Config
from DSP import m2smp
from SOSarray_DSP import SOSarray_DSP


%matplotlib inline
config = Config()

Create impulse signal#

[ ]:

numBlocks = 300
signalLength = config.blockSize * numBlocks
signal = np.zeros((signalLength, 1))
signal[0, 0] = 1
time = np.arange(1, signalLength + 1) / config.fs

1 TASK: Compute delay-proportional absorption filters#

  1. Compute filter for 100 samples delay. For this implement onePoleAbsorption and designOnePoleFilter.
    Hint: Recall the frequency response for a one-pole filter:
    \[H(e^{j\omega T}) = \frac{b_0}{1 + a_1 e^{-j\omega T}}\]

    Think about a way of solving this for \(a_1\) and \(b_0\) given the magnitude values at DC and Nyquist.

  2. Plot the filter magnitude response
  3. Convert filter magnitude response to attenuation per sample
  4. Plot reverberation time corresponding to this attenuation per sample
  5. Repeat a)-d) for 1000 samples delay

  6. Explain the plots

[ ]:

def designOnePoleFilter(HDc, HNyq):
    """
    Design one-pole filter

    Inputs:
        HDc  - Linear magnitude at DC of size [1, number of filters] or [number of filters]
        HNyq - Linear magnitude at Nyquist of size [1, number of filters] or [number of filters]

    Outputs:
        sos - sos filters of size [6 x number of filters]

    Example:
        sos = designOnePoleFilter(np.random.rand(3), np.random.rand(3))
    """

    # Ensure inputs are numpy arrays
    HDc = np.atleast_1d(HDc)
    HNyq = np.atleast_1d(HNyq)

    numFilters = HDc.size
    sos = np.zeros((6, numFilters))
    sos[3, :] = 1  # a0 = 1

    # TASK: Implement one-pole filter
    # a) See: https://ccrma.stanford.edu/~jos/fp/One_Pole.html
    # b) Solve for a1 and b0 (e.g. by plugging in HDc and HNyq and rearranging)
    # c) Implement the filter coefficients into second-order sections (SOSs)
    # [b0; b1; b2; a0; a1; a2], where a0 = 1.
    # d) Vectorize for multiple filters.

    # YOUR CODE HERE
    raise NotImplementedError()

    return sos

def onePoleAbsorption(RT_DC, RT_NY, delays, fs):
    """
    Design one-pole absorption filters according to specified T60

    Syntax:  sos = onePoleAbsorption(RT_DC, RT_NY, delays, fs)

    Inputs:
        RT_DC  - Reverberation time in [seconds] at DC frequency
        RT_NY  - Reverberation time in [seconds] at Nyquist frequency
        delays - Array of delay times in [samples]
        fs     - Sampling frequency

    Outputs:
        sos - Filter coefficients of size [6, 1, N] where N is number of filters

    Example:
        sos = onePoleAbsorption(2, 1, np.array([100, 130]), 48000)
    """

    # Ensure delays is a numpy array
    delays = np.atleast_1d(delays)

    # Dummy
    N = delays.size
    sos = np.zeros((6, 1, N))
    sos[[0, 3], :, :] = 1

    # Calculate magnitude response at DC and Nyquist
    # Hint: RT60 to slope conversion: -60 dB / (RT60 * fs) gives dB per sample
    # YOUR CODE HERE
    raise NotImplementedError()

    # Design one-pole filters
    sos = designOnePoleFilter(HDc, HNyq)

    # conform to shape [6, 1, N]
    sos = sos[:, np.newaxis, :]

    return sos

delayLengths = np.array([100])

# Compute filter for 100 samples delay
# YOUR CODE HERE
raise NotImplementedError()

# Convert filter magnitude response to attenuation per sample
# YOUR CODE HERE
raise NotImplementedError()

delayLengths = np.array([1000])

# Repeat the same steps for 1000 samples delay
# YOUR CODE HERE
raise NotImplementedError()

# Plot the filter magnitude response
# YOUR CODE HERE
raise NotImplementedError()

# Plot the reverberation time
# YOUR CODE HERE
raise NotImplementedError()

2 TASK: Implement a lossless feedback matrix#

  1. Generate a lossless matrix of size numberOfDelays

  2. Prove numerically the losslessness

  3. Repeat a) and b) for another type of matrix

  4. Explain differences between those feedback matrices, e.g., density, computational cost.

[ ]:

numberOfDelays = config.Late.numberOfDelays

# YOUR CODE HERE
raise NotImplementedError()

# Set feedback matrix for a later task
config.Late.feedbackMatrix = A

# Explain diffrences
# YOUR CODE HERE
raise NotImplementedError()

3 TASK: Initialize the FDN parameters#

  1. Compute delay line lengths from config.roomSize

  2. Initialize and set a multichannel delay with correct lengths. Use dsp.Delay(delayLengths), where delayLengths is [1xnCh] vector in samples.

  3. Compute absorption filters

  4. Initialize a SOSarray_DSP with filters

  5. Generate input and output gains (can be ones or random)

  6. Set the lossless feedback matrix

[ ]:

blockSize = config.blockSize
RT60 = config.RT60
fs = config.fs

# Compute delay line lengths
# YOUR CODE HERE
raise NotImplementedError()

# Init delay buffers and indices arrays
# YOUR CODE HERE
raise NotImplementedError()

# Compute absorption filters and set them in the SOSarray_DSP
# YOUR CODE HERE
raise NotImplementedError()

# Set input and output gains
# YOUR CODE HERE
raise NotImplementedError()

# Set feedback matrix
# YOUR CODE HERE
raise NotImplementedError()

4 TASK: Implement the late reverberation processing with a Feedback Delay Network.#

Before implementing a full Feedback Delay Network (FDN), we’ll start with a simpler building block: a one-pole feedback comb filter. This filter has the transfer function:

\[H(z) = \frac{1}{1 - g z^{-D}}\]

where \(g\) is the feedback gain and \(D\) is the delay length.

The processing structure we use here (circular buffers, block processing, feedback paths) is very similar to what you’ll implement for the FDN.

[ ]:


def comb_filter(signal: np.ndarray,
                delayLength: int,
                feedbackGain: float,
                blockSize: int = 256) -> np.ndarray:
    """
    Process a 1D feedback comb filter using block processing and a circular buffer.

    Args:
        signal: Input signal array of shape [num_samples, 1] or [num_samples].
        delayLength: Delay line length in samples (integer, >= 1).
        feedbackGain: Feedback gain scalar (typically |g| < 1 for stability).
        blockSize: Block size for processing.

    Returns:
        Output signal array of shape [num_samples, 1].
    """
    if signal.ndim == 1:
        signal = signal.reshape(-1, 1)

    numSamples = signal.shape[0]

    # Delay state
    delayBuffer = np.zeros(delayLength, dtype=np.float64)
    delayWriteIndex = 0

    # Processing buffers
    lastBlock = np.zeros(blockSize, dtype=np.float64)
    output = np.zeros((numSamples, 1), dtype=np.float64)

    numBlocks = int(np.ceil(numSamples / blockSize))

    for it_block in range(numBlocks):
        block_index = slice(it_block * blockSize, (it_block + 1) * blockSize)

        # Skip if not full block
        if block_index.stop > numSamples:
            break

        block = signal[block_index, :]

        # Output mixing: direct reference since we don't modify lastBlock
        out = lastBlock

        # Feedback and input injection
        feedback_paths = lastBlock * feedbackGain  # Direct scalar multiplication
        delay_input = block[:, 0] + feedback_paths  # Direct 1D operations

        # Circular-buffer delay processing
        delay_output = np.zeros(blockSize, dtype=np.float64)

        # Single delay line processing
        for i in range(blockSize):
            read_idx = (delayWriteIndex - delayLength) % delayLength
            delay_output[i] = delayBuffer[read_idx]
            delayBuffer[delayWriteIndex] = delay_input[i]
            delayWriteIndex = (delayWriteIndex + 1) % delayLength

        # Update state for next block
        lastBlock = delay_output

        # Write output
        output[block_index, 0] = out

    return output

# Example usage
fs = 48000
duration = 1.0
t = np.arange(int(fs * duration)) / fs
# Impulse input
x = np.zeros_like(t)
x[0] = 1.0

y = comb_filter(x, delayLength=int(0.03 * fs), feedbackGain=0.6, blockSize=256)


# Plot input and output (impulse response)

plt.figure()
plt.plot(t, y, label='Output', alpha=0.9)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Feedback Comb Filter - Impulse Response')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Using a similar processing structure as the comb filter above

  1. Implement a FDN processing loop that handles multiple delay lines simultaneously

  2. Use lastBlock to retain signal between the for loop iterations

FDN

[ ]:

lastBlock = np.zeros((blockSize, numberOfDelays))
output = np.zeros((signalLength, 1))

for it_block in range(numBlocks):
    block_index = slice(it_block * blockSize, (it_block + 1) * blockSize)
    block = signal[block_index, :]

    # Compute output using lastBlock and outputGains
    # YOUR CODE HERE
    raise NotImplementedError()

    # Compute feedback paths and delay input
    # YOUR CODE HERE
    raise NotImplementedError()

    # Implement delay lines using circular buffer
    delay_output = np.zeros((blockSize, numberOfDelays))
    # YOUR CODE HERE
    raise NotImplementedError()

    # Process delay lines through absorption filters
    # YOUR CODE HERE
    raise NotImplementedError()

    lastBlock = delay_output

    output[block_index, :] = out

5 TASK: describe the result#

  1. Highlight sonic features such as smoothness, density, ringing.

  2. How does the time-domain plot conform to your expectation of the reverb design. Highlight features such as early reflections and late reverberation.

  3. How does the spectrogram conform to your expectation of the reverb design. Highlight features such as reverberation time.

[ ]:

# Save audio file
output_normalized = output / np.max(np.abs(output))
wavfile.write('assignment_Late_DSP.wav', config.fs, output_normalized.astype(np.float32))

# Plot and describe your results
# YOUR CODE HERE
raise NotImplementedError()

6 TASK: make one more FDN design#

  1. Replace the input signal with something more interesting

  2. Change at least one of these parameters: feedback matrix, room size, reverberation time, number of delays

  3. Explain your change and the result

[ ]:

# 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,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
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