{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Impulse Responses, Convolution\n", "\n", "[return to main page](index.ipynb)\n", "\n", "In this notebook we will investigate an unknown digital system.\n", "The only information we have about this system is that it is LTI (linear and time invariant) and that it is defined by the function `tools.blackbox()` in the file [tools.py](tools.py).\n", "\n", "To be able to use it, you have to import it:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import tools" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After that, you should have a look at the documentation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "tools.blackbox?" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "In this notebook, we'll try to find out as many things as possible about this system, without looking at its source code.\n", "Later, we'll also have a quick look at a non-linear system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Listening to the System\n", "\n", "*Exercise:* Load the audio file [data/xmas.wav](data/xmas.wav) and apply the function `tools.blackbox()` to it." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Listen to the result. Compare the original signal and the resulting signal.\n", "Describe the differences.\n", "What does the system realized by `tools.blackbox()` sound like?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Obtaining the Impulse Response\n", "\n", "An LTI system can be completely described by its impulse response (which may be infinitely long, however).\n", "\n", "*Exercise:* Determine the impulse response (or at least an estimation of it) of `tools.blackbox()`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "To do that, use a unit impulse as input to the system.\n", "The resulting output is the impulse response.\n", "To get a meaningful response, append zeros (this is called *zero-padding*) to your unit impulse signal until it has a total length of $\\frac{1}{10}$ of a second.\n", "\n", "*Exercise:* Plot the impulse response (with the time in seconds along the x-axis).\n", "Note that the amplitude seems to be zero in the end, but it's not!\n", "Zoom into the plot until you can see the non-zero values." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To obtain more insight about the parts of the impulse response with very small amplitudes, we can try logarithmic scaling.\n", "\n", "*Exercise:* Plot the impulse response in decibels (i.e. with logarithmic amplitude scaling).\n", "Have a look in [tools.py](tools.py), you might find a useful function there ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Use the function [scipy.signal.freqz()](http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.freqz.html) to calculate the frequency response given the impulse response." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Plot the magnitude of the frequency response in decibels on a logarithmic frequency axis." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Try all combinations of logarithmic/linear scaling on the x-/y-axis (that's 4 combinations)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Naive Convolution\n", "\n", "*Exercise:* Write a function called `naive_convolution()` that computes the convolution of two one-dimensional arrays by means of two nested loops according to the equation\n", "\n", "$$y[n] = x[n] \\ast h[n] = \\sum_{m=-\\infty}^{\\infty} x[m] \\cdot h[n-m],$$\n", "\n", "where $x$ and $h$ are one-dimensional arrays of finite lengths.\n", "The infinite sum can be changed to a finite sum by assuming that all values before index 0 and all values after the last array element are equal to zero.\n", "\n", "Following this assumption, at which indices $n$ does $y[n]$ have its first and last non-zero value?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Use the function `naive_convolution()` to convolve the audio signal with the impulse response from the previous exercise.\n", "\n", "WARNING: The calculation may take a very long time!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Listen to the result (only if you had the patience to wait for the computation to finish)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* How long is the output signal compared to the input signals?\n", "Does this make sense according to the equation above?\n", "Compare the output signal with the result from the function\n", "`tools.blackbox()`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Matrix Convolution\n", "\n", "The implementation used in the previous exercise is very inefficient.\n", "\n", "*Exercise:* Write a function called `matrix_convolution()`, that does the same thing, but using matrix multiplication." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Call this function as well (using the same input signal and impulse response as\n", "before) and check if there is a difference in computation time.\n", "Warning: Depending on the order of arguments, this may need huge amounts of memory and bring your computer to a halt.\n", "You should try it first with a small part of the signal (say, the first 1000 samples or so) and then gradually increase the length of the input signal until you know which order is the \"right\" one." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Listen to the resulting output signal.\n", "The length of the output signal should be the same as in the previous exercise." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Switch the two input parameters and check if this has an influence on the computation time.\n", "Warning: See above, you should try this only with a small part of the input signal." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fast Convolution\n", "\n", "As you know, convolution in the time domain is equivalent to element-wise multiplication in the frequency domain.\n", "\n", "*Exercise:* Write a function named `fft_convolution()`, that transforms both inputs to the frequency domain using [numpy.fft.rfft()](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfft.html), does an element-wise multiplication of the (complex) spectra and transforms the result back with [numpy.fft.irfft()](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.irfft.html).\n", "Deliberate whether the result is supposed to be real or complex.\n", "Check the data types and sizes of the intermediate arrays." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* What would happen if you use `fft()` instead of `rfft()`?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "WARNING: The fast convolution using FFT implements *cyclic convolution*!\n", "\n", "Take care to use *zero-padding* to increase the length of the FFT to avoid components from the cyclic convolution in the output signal.\n", "By the way, the equation for cyclic convolution looks like this:\n", "\n", "$$y[n] = x[n] \\circledast_N h[n] = \\sum_{m=0}^{N-1} x[m] \\cdot h[(n-m) \\bmod N].$$\n", "\n", "If you want, you can also choose the next bigger power of two as FFT length - this is more efficient in many (but not all) cases.\n", "After the inverse transform you should trim the result to the appropriate length.\n", "\n", "Is the calculation with `fft_convolution()` faster than with the previously used functions?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Using Existing Functions\n", "\n", "Convolution is very important in many areas of signal processing, thus it is not surprising that a function for it is available in NumPy: [numpy.convolve()](http://docs.scipy.org/doc/numpy/reference/generated/numpy.convolve.html).\n", "This function does basically the same as our `naive_convolution()`.\n", "Just much more efficient.\n", "\n", "*Exercise:* Try it!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Have a look at the `mode` argument.\n", "Which \"mode\" did we implement before?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "But that's not everything yet!\n", "\n", "Let's import the signal processing module from SciPy and see if we find something useful there:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from scipy import signal" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "*Exercise:* Do the same convolution as before, but now using [scipy.signal.fftconvolve()](http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "TODO: longer input signal!\n", "\n", "TODO: Have a look at the documentation, especially about the sizes of the two input arrays.\n", "What happens if you switch the two input arrays (i.e. if you switch the input signal and the impulse response)?\n", "\n", "TODO: `fftfilt()`:\n", "\n", "* https://github.com/scipy/scipy/issues/1364\n", "\n", "* https://github.com/ColumbiaCMB/kid_readout/blob/master/kid_readout/analysis/timeseries/fftfilt.py\n", "\n", "* http://cnx.org/contents/9479779f-bd46-4d3a-a1d4-30580ae8aacc@10/Convolution_Algorithms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What if the System is Non-Linear?\n", "\n", "The function `tools.blackbox_nonlinear()` provides a non-linear system.\n", "\n", "*Exercise:* Listen how it sounds if you send the signal from `xmas.wav` (or some other\n", "audio signal) through this system." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Check if this system can be described by an impulse response.\n", "Use the unit impulse from before to obtain the impulse response of `tools.blackbox_nonlinear()`.\n", "\n", "Convolve our audio signal with this impulse response and listen to it.\n", "Do you hear any difference between the signal convolved with the impulse response and the signal sent directly through the non-linear system?\n", "\n", "Which leads us to the final question:\n", "Can a non-linear system be described completely with an impulse response?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dynamic Range Control\n", "Speaking of non-linear systems, it is worth noting that a non-linear processing can be quite useful. Dynamic signal processors, like limiter, compressor, expander and gate, are being frequently used to control the dynamic range of a signal. It allows to increase the loudness without causing amplitude clipping.\n", "\n", "*Exercise:* Load the audio file [data/singing.wav](data/singing.wav) and apply the function `tools.compressor()` to it. Try with `threshold=-30`, `ratio=3`, `makeup_gain=12`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Plot the original and the processed signal. Compare the signals in terms of maximum value and energy." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Listen to the signals, and compare the loudness." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Not bad, huh?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solutions\n", "\n", "If you had problems solving some of the exercises, don't despair!\n", "Have a look at the [example solutions](ir-solutions.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", " \n", " \"CC0\"\n", " \n", "
\n", " To the extent possible under law,\n", " the person who associated CC0\n", " with this work has waived all copyright and related or neighboring\n", " rights to this work.\n", "

" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [default]", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.2" } }, "nbformat": 4, "nbformat_minor": 1 }