Convolution Part 2

Chris Tralie

In [2]:
import numpy as np
import matplotlib.pyplot as plt
import IPython.display as ipd
import librosa
from scipy.signal import fftconvolve

Let's look at how to define convolution mathematically. We'll write it in code first, and then we'll translate it into math

In [ ]:
# We have an underlying signal x that holds
# our original sound
# We want to hear what it sounds like in an
# envrionment whose impulse response is stored
# in the array h

# Assume I have an ouput array y that's big
# enough to hold all of the echoes

# Fill in samples in y one at a time
for i in range(len(y)):
    for k in range(len(h)):
        if i - k > 0 and i - k < len(x):
            y[i] += h[k]*x[i-k]

Here's how to write this in mathematical notation

$ y[n] = \sum_{k = 0}^N h[k]x[n-k] $

Amazingly, if we switch the roles of $x$ and $h$, we get the same answer. In other words, convolution is "commutative"

$y[n] = \sum_{k=0}^{M} x[k]h[n-k]$

Designing Impulse Responses

The easiest way to design an impulse response is to measure it! If we record ourselves clapping or popping a balloon in some environment, then we will get one little impulse/clap for every echo that happens in that environment, and that is our impulse response! So, for example, we can hear what "Jessie's girl" sounds like from the car in front of us when we're stuck in the JFK tunnel.

Here's "Jessie's girl"

In [3]:
x, sr = librosa.load("jessiesgirl.mp3", sr=44100)
ipd.Audio(x, rate=sr)
/home/ctralie/anaconda3/lib/python3.7/site-packages/librosa/core/audio.py:162: UserWarning: PySoundFile failed. Trying audioread instead.
  warnings.warn("PySoundFile failed. Trying audioread instead.")
Out[3]:

Here's the impulse response of the JFK tunnel:

In [4]:
h, sr = librosa.load("imp_JFKTunnel.mp3", sr=44100)
ipd.Audio(h, rate=sr)
/home/ctralie/anaconda3/lib/python3.7/site-packages/librosa/core/audio.py:162: UserWarning: PySoundFile failed. Trying audioread instead.
  warnings.warn("PySoundFile failed. Trying audioread instead.")
Out[4]:

And here's the result of the convolution

In [5]:
y = fftconvolve(x, h)
ipd.Audio(y, rate=sr)
Out[5]: