Impulse Response Project Log 1: Raytracing around a Ring

Goal

• Simulate acoustic impulse responses (IRs) for convolutional reverb in different geometries.

• Start with a simple ring geometry, then move on to more exotic spaces like a Möbius strip.

What I Did

• Set up a ring geometry with a listener/microphone on one side of the ring and a source on the other

• Used raytracing to simulate the propagation of sound waves around the ring, allowing for reflection of the walls

• Tracked cumulative distance travelled by each ray before reaching the listener to add dissipation to the final signal (the rays which travelled further before reaching the listener have a weaker amplitude)

• Turned it into a soundfile with a small amount of processing

import numpy as np
import scipy.signal as signal
import soundfile as sf

# Constants
c = 344  # speed of sound (m/s)
dt = 0.001  # simulation step size in seconds
t_max = 5.0  # total simulation time
R = 35.0  # source radius
fs = 44100  # sampling rate for IR output

def ring(r, phi):
    x = r * np.cos(phi)
    y = r * np.sin(phi)
    return np.array([x, y])

# Main ray tracing function with cumulative distance tracking
def raytrace():
    rays = []
    for j in np.arange(0, 360, 0.005):
        t = 0
        pos = ring(R, 0)
        theta = np.deg2rad(j)
        direction = np.array([np.cos(theta), np.sin(theta)])
        ray_path = [pos.copy()]
        cumulative_distance = [0.0]
        total_dist = 0.0

        while t < t_max:
            step = direction * c * dt
            pos += step
            total_dist += np.linalg.norm(step)
            t += dt

            r = np.linalg.norm(pos)
            if r <= 31 or r >= 39:
                radial = pos / r
                scatter_angle = np.random.normal(0, np.deg2rad(2))
                cos_a = np.cos(scatter_angle)
                sin_a = np.sin(scatter_angle)
                R_mat = np.array([[cos_a, -sin_a], [sin_a, cos_a]])
                direction = R_mat @ (direction - 2 * np.dot(direction, radial) * radial)

            ray_path.append(pos.copy())
            cumulative_distance.append(total_dist)

        rays.append((ray_path, cumulative_distance))
    return rays

# Build the impulse response using distance-aware hits
def impulseresponse(rays, listener_pos=np.array([0.0, 35.0]), listener_radius=2,
                    fs=44100, t_max=5.0):
    n_samples = int(t_max * fs)
    ir = np.zeros(n_samples)
    impulse_data = []

    for ray_path, distances in rays:
        for i, pos in enumerate(ray_path):
            if np.linalg.norm(pos - listener_pos) < listener_radius:
                distance = distances[i]
                arrival_time = distance / c
                impulse_data.append((arrival_time, distance))
                # allow multiple hits

    print(f"Number of impulses detected: {len(impulse_data)}")

    for arrival_time, distance in impulse_data:
        i = int(arrival_time * fs)
        if i < n_samples:
            T60 = 2.5
            amp = (1.0 / (distance + 1e-6)) * np.exp(-6.91 * arrival_time / T60)
            kernel = amp * signal.windows.gaussian(11, std=1.5)
            start = max(i - 5, 0)
            end = min(i + 6, n_samples)
            ir[start:end] += kernel[:end - start]

    return ir

# Final cleanup and export
def prepare_ir_for_convolution(ir, fs, t_max=5.0, t60=2.5):
    threshold = 1e-6
    nonzero_indices = np.where(np.abs(ir) > threshold)[0]
    if len(nonzero_indices) == 0:
        print("No impulses detected in IR.")
        return ir
    first = nonzero_indices[0]
    ir = np.roll(ir, -first)
    ir[-first:] = 0

    t = np.linspace(0, t_max, len(ir))
    decay_env = np.exp(-6.91 * t / t60)
    ir *= decay_env

    max_amp = np.max(np.abs(ir))
    if max_amp > 0:
        ir /= max_amp

    fade_len = int(0.1 * fs)
    ir[-fade_len:] *= np.linspace(1, 0, fade_len)

    return ir

# Run full pipeline
rays = raytrace()
ir = impulseresponse(rays, t_max=t_max)
ir = prepare_ir_for_convolution(ir, fs)
sf.write("impulse_response_ready3.wav", ir, fs, subtype='PCM_16')

Notes

• Code was ineffective due to hardware limitation (my laptop is a decade old).

• Saving the cumulative distance tracking, which is important for the dissipation of the sound waves, and the path of each ray uses too much RAM (around 8GB) which caused my code to run for several minutes before throwing me a memory error.

• Need dissipation for it to sound realistic and need a high number of rays to generate a decent impulse response.

• I will need to implement a way to have dissipation but reduce RAM. One way I will considering is only keeping rays in memory which intersect with the listener or are near the listener.

• Another core issue for the project is how to get an impulse response out of a mobius strip using raytracing since the half twist in the Mobius strip causes phase to invert as waves move around the strip, but rays don’t have phase…

What Next

• Incorporate a dissipation model which doesn’t eat all my RAM

• Figure out how to give rays a phase or explore other numerical wave simulation methods, possibly Finite Difference Time Domain (FDTD) (computationally expensive though)