CoCalc -- chap10soln.ipynb

Path: examples/think-dsp/code / chap10soln.ipynb

Views: ⁷³⁴⁴

Kernel: Python 3

ThinkDSP

This notebook contains solutions to exercises in Chapter 10: Signals and Systems

License: Creative Commons Attribution 4.0 International

In [ ]:

# Get thinkdsp.py

import os

if not os.path.exists('thinkdsp.py'):
    !wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/thinkdsp.py

In [2]:

import numpy as np
import matplotlib.pyplot as plt

from thinkdsp import decorate

Exercise 1

In this chapter I describe convolution as the sum of shifted, scaled copies of a signal. Strictly speaking, this operation is linear convolution, which does not assume that the signal is periodic.

But when we multiply the DFT of the signal by the transfer function, that operation corresponds to circular convolution, which assumes that the signal is periodic. As a result, you might notice that the output contains an extra note at the beginning, which wraps around from the end.

Fortunately, there is a standard solution to this problem. If you add enough zeros to the end of the signal before computing the DFT, you can avoid wrap-around and compute a linear convolution.

Modify the example in chap10soln.ipynb and confirm that zero-padding eliminates the extra note at the beginning of the output.

Solution: I'll truncate both signals to $2^{16}$ elements, then zero-pad them to $2^{17}$ . Using powers of two makes the FFT algorithm most efficient.

Here's the impulse response:

In [3]:

if not os.path.exists('180960__kleeb__gunshot.wav'):
    !wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/180960__kleeb__gunshot.wav

In [4]:

from thinkdsp import read_wave

response = read_wave('180960__kleeb__gunshot.wav')

start = 0.12
response = response.segment(start=start)
response.shift(-start)

response.truncate(2**16)
response.zero_pad(2**17)

response.normalize()
response.plot()
decorate(xlabel='Time (s)')

And its spectrum:

In [5]:

transfer = response.make_spectrum()
transfer.plot()
decorate(xlabel='Frequency (Hz)', ylabel='Amplitude')

Here's the signal:

In [6]:

if not os.path.exists('92002__jcveliz__violin-origional.wav'):
    !wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/92002__jcveliz__violin-origional.wav

In [7]:

violin = read_wave('92002__jcveliz__violin-origional.wav')

start = 0.11
violin = violin.segment(start=start)
violin.shift(-start)

violin.truncate(2**16)
violin.zero_pad(2**17)

violin.normalize()
violin.plot()
decorate(xlabel='Time (s)')

And its spectrum:

In [8]:

spectrum = violin.make_spectrum()

Now we can multiply the DFT of the signal by the transfer function, and convert back to a wave:

In [9]:

output = (spectrum * transfer).make_wave()
output.normalize()

The result doesn't look like it wraps around:

In [10]:

output.plot()

And we don't hear the extra note at the beginning:

In [11]:

output.make_audio()

We should get the same results from np.convolve and scipy.signal.fftconvolve.

First I'll get rid of the zero padding:

In [12]:

response.truncate(2**16)
response.plot()

In [13]:

violin.truncate(2**16)
violin.plot()

Now we can compare to np.convolve:

In [14]:

output2 = violin.convolve(response)

The results are similar:

In [15]:

output2.plot()

And sound the same:

In [16]:

output2.make_audio()

But the results are not exactly the same length:

In [17]:

len(output), len(output2)

(131072, 131071)

scipy.signal.fftconvolve does the same thing, but as the name suggests, it uses the FFT, so it is substantially faster:

In [18]:

from thinkdsp import Wave

import scipy.signal
ys = scipy.signal.fftconvolve(violin.ys, response.ys)
output3 = Wave(ys, framerate=violin.framerate)

The results look the same.

In [19]:

output3.plot()

And sound the same:

In [20]:

output3.make_audio()

And within floating point error, they are the same:

In [21]:

output2.max_diff(output3)

1.7763568394002505e-14

Exercise 2

The Open AIR library provides a "centralized... on-line resource for anyone interested in auralization and acoustical impulse response data" (http://www.openairlib.net). Browse their collection of impulse response data and download one that sounds interesting. Find a short recording that has the same sample rate as the impulse response you downloaded.

Simulate the sound of your recording in the space where the impulse response was measured, computed two way: by convolving the recording with the impulse response and by computing the filter that corresponds to the impulse response and multiplying by the DFT of the recording.

Solution: I downloaded the impulse response of the Lady Chapel at St Albans Cathedral http://www.openairlib.net/auralizationdb/content/lady-chapel-st-albans-cathedral

Thanks to Audiolab, University of York: Marcin Gorzel, Gavin Kearney, Aglaia Foteinou, Sorrel Hoare, Simon Shelley.

In [22]:

if not os.path.exists('stalbans_a_mono.wav'):
    !wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/stalbans_a_mono.wav

In [23]:

response = read_wave('stalbans_a_mono.wav')

start = 0
duration = 5
response = response.segment(duration=duration)
response.shift(-start)

response.normalize()
response.plot()
decorate(xlabel='Time (s)')

Here's what it sounds like:

In [24]:

response.make_audio()

The DFT of the impulse response is the transfer function:

In [25]:

transfer = response.make_spectrum()
transfer.plot()
decorate(xlabel='Frequency (Hz)', ylabel='Amplitude')

Here's the transfer function on a log-log scale:

In [26]:

transfer.plot()
decorate(xlabel='Frequency (Hz)', ylabel='Amplitude',
         xscale='log', yscale='log')

Now we can simulate what a recording would sound like if it were played in the same room and recorded in the same way. Here's the trumpet recording we have used before:

In [27]:

if not os.path.exists('170255__dublie__trumpet.wav'):
    !wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/170255__dublie__trumpet.wav

In [28]:

wave = read_wave('170255__dublie__trumpet.wav')

start = 0.0
wave = wave.segment(start=start)
wave.shift(-start)

wave.truncate(len(response))
wave.normalize()
wave.plot()
decorate(xlabel='Time (s)')

Here's what it sounds like before transformation:

In [29]:

wave.make_audio()

Now we compute the DFT of the violin recording.

In [30]:

spectrum = wave.make_spectrum()

I trimmed the violin recording to the same length as the impulse response:

In [31]:

len(spectrum.hs), len(transfer.hs)

(110251, 110251)

In [32]:

spectrum.fs

array([0.00000e+00, 2.00000e-01, 4.00000e-01, ..., 2.20496e+04,
       2.20498e+04, 2.20500e+04])

In [33]:

transfer.fs

array([0.00000e+00, 2.00000e-01, 4.00000e-01, ..., 2.20496e+04,
       2.20498e+04, 2.20500e+04])

Now we can multiply in the frequency domain and the transform back to the time domain.

In [34]:

output = (spectrum * transfer).make_wave()
output.normalize()

Here's a comparison of the original and transformed recordings:

In [35]:

wave.plot()

In [36]:

output.plot()

And here's what it sounds like:

In [37]:

output.make_audio()

Now that we recognize this operation as convolution, we can compute it using the convolve method:

In [38]:

convolved2 = wave.convolve(response)
convolved2.normalize()
convolved2.make_audio()