{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Discrete Fourier Transform\n", "\n", "[return to main page](index.ipynb)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Five Cosines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Generate five [cosine tones](http://docs.scipy.org/doc/numpy/reference/generated/numpy.cos.html) with the following parameters:\n", "\n", "| | Amplitude | Frequency/Hz | Phase/Degree |\n", "|-------|-----------|-------------:|-------------:|\n", "| $x_1$ | 1.0 | 200 | 0 |\n", "| $x_2$ | 0.75 | 400 | 0 |\n", "| $x_3$ | 0.5 | 600 | 90 |\n", "| $x_4$ | 0.25 | 800 | 90 |\n", "| $x_5$ | 0.125 | 1000 | -90 |\n", "\n", "All five signals should have a duration of 1 second and a sampling rate of 44.1 kHz. *Hint:* Try to exploit [broadcasting](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Generate the sum signal $x_6 = x_1 + x_2 + x_3 + x_4 + x_5$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Plot all 6 signals at once." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Listen to the individual tones and to the sum signal. Use the function `tools.fade()` from [tools.py](tools.py) to avoid clicks at the beginning and/or the end. Watch for the volume of the sum signal, reduce the level if necessary. What happens if you don’t reduce the sound level?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that limiting the range from -1 to 1 is only relevant once you play back a signal over the sound card or save it as a sound file.\n", "For internal calculations you can of course use the full double precision (a.k.a. 'float64') range." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Naively Implementing the DFT Equation\n", "\n", "*Exercise:* Write a function named `naive_fourier_transform()`, that calculates the Discrete Fourier Transform $X[k]$ of a one-dimensional array $x[n]$ using the equation\n", "\n", "$$X[k] = \\sum_{n=0}^{N-1} x[n] \\text{e}^{-\\text{j}2\\pi kn/N},$$\n", "\n", "where $N$ is the number of samples of $x$ and $k$ is the frequency index with $0 \\le k \\lt N$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Call your function with $x_1$ as argument." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The calculation will take some time; in the meantime, estimate how many multiplications and additions are calculated and how often the exponential function is called." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If your calculation from above still isn't finished, stop it with the menu item \"Kernel\" $\\to$ \"Interrupt\" (or by clicking the \"stop\" icon).\n", "\n", "Call the function again, but this time only using the first 1000 values of $x_1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementing the DFT with Matrix Multiplication\n", "\n", "*Exercise:* Think about how to implement the same function as a matrix multiplication.\n", "Would that need less processing power?\n", "What about the memory requirements?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using the Fast Fourier Transform\n", "\n", "*Exercise:* Now, calculate the discrete Fourier transform of $x_1$ using the function [numpy.fft.fft()](http://docs.scipy.org/doc/numpy/reference/routines.fft.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That should be faster, right?\n", "\n", "*Exercise:* What algorithmic complexity does the FFT algorithm have?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you only computed the first 1000 values above, do the same here.\n", "\n", "*Exercise:* Calculate the difference of the results of `naive_fourier_transform()` and `numpy.fft.fft()`.\n", "What is the maximum absolute error?\n", "To assess the error, consider the maximum value of the result." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the maximum error (in this concrete example) is in the order of magnitude of\n", "around $10^{-10}$ to $10^{-8}$, everything is OK (what you are seeing are rounding errors).\n", "If it is significantly larger, something is wrong.\n", "If it is significantly smaller, as well.\n", "\n", "*Exercise:* How large is the error in decibels relative to the absolute maximum?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* In comparison, how large is the theoretically maximum signal-to-noise ratio (in decibel) of a linearly quantized 16-bit signal?\n", "Does this depend on the sampling frequency?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that you have checked (admittedly not very thoroughly) the correctness of the\n", "function `numpy.fft.fft()` and experienced its efficiency, you should also use it.\n", "\n", "*Exercise:* Apply the function to the signals $x_1$ to $x_6$ and plot the results. Plot $x_1$ to $x_5$ in a common plot and $x_6$ separately. Don't forget the frequency axis!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should have noticed that the result is complex - it is called the *complex spectrum*.\n", "The DFT of a real signal is in general complex, by the way.\n", "In order to recognize something in the plots, you should plot magnitude\n", "and phase of the complex spectrum separately.\n", "\n", "*Exercise:* Calculate the inverse FFT ([numpy.fft.ifft()](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.ifft.html)) of the FFT of $x_1$.\n", "Is there a difference between the result and $x_1$?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Write a function named `plot_fft()`, that calculates the DFT of a signal and plots a figure with two subplots with magnitude and phase (in degrees).\n", "Pass the sampling frequency to the function, in order to be able to label the x-axis in Hertz." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## DFT of a Real Signal\n", "\n", "One property of the DFT is, that the transform of a real signal is symmetrical.\n", "\n", "*Exercise:* What does this mean with regards to the plots?\n", "Especially to the magnitude plot?\n", "How can we \"simplify\" our plots with that?\n", "Use this simplification in a new version of your function `plot_fft()`.\n", "If you don't remember the symmetries, you can look them up in any book about signal processing (or on Wikipedia)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since half of the data is redundant, we wouldn't even have to calculate it.\n", "NumPy provides the functions [numpy.fft.rfft()](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfft.html) and [numpy.fft.irfft()](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.irfft.html) to provide only the non-redundant result.\n", "\n", "*Bonus Exercise:* Make yet another variant of your function `plot_fft()` that uses the \"real\" FFT functions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Plot all signals from the table at the top of the page in the time domain (with `matplotlib.pyplot.plot()`) and in the frequency domain (with your function `plot_fft()`).\n", "Which of the given parameters can you recognize in which plots? Which are not visible?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "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](dft-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": { "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.4" } }, "nbformat": 4, "nbformat_minor": 1 }