{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Unit 8: Signals and Sampling\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "
\n", "

## Overview and Learning Objectives

\n", "\n", "In technical fields such as engineering or computer science, a signal may be defined as a function that conveys information about the state or behavior of a physical system. For example, a signal may describe the time-varying sound pressure at some place, the motion of a particle through some space, the distribution of light on a screen representing an image, or the sequence of images as in the case of a video signal. In this unit, we consider the case of sound or audio signals, which may be represented graphically by a plot that shows the relative air pressure (with respect to a reference air pressure) over time. In some way, a sinusoid of a given frequency may be thought of as a prototype of an audio signal. We start this unit by formally defining a signal as a mathematical function. Interleaving theory with practice, we then provide a Python function to generate a sinusoid with different parameters (duration, amplitude, frequency, phase, sampling rate), which you can look at via a plot and listen to via audio playback. To digitally compute with signals, one typically needs to convert a signal given on a continuous-time axis into one given on a discrete-time axis—a process commonly referred to as sampling. While providing a Python function that performs some equidistant sampling, we also discuss the kind of information that may be lost in the sampling process. This leads us to topics related to aliasing and the sampling theorem. Finally, using sinusoidal signals as an example, we cover the phenomena of constructive and destructive interference. In the exercises, we further deepen these aspects by studying the phenomenon known as beating and by looking at further examples illustrating aliasing. We assume that you are familiar with these fundamental concepts that are typically taught in an introductory signal processing course. Besides reviewing these concepts, another main learning objective of this unit is to show you how to generate and use multimedia objects (e.g., audio, image, and video objects) within the Jupyter notebook framework. In particular, in fields such as multimedia engineering, the interaction with concrete multimedia objects (e.g., sound examples) will help you move from recalling and reciting theoretical concepts towards comprehension and application. Following similar educational considerations, you may find additional material in the FMP Notebooks for teaching and learning fundamentals of music processing. In particular, the FMP Notebook on Sampling and the FMP Notebook on Interference and Beating have served as basis for this unit. \n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "## Signals\n", "\n", "A **sound** is generated by a vibrating object such as the vocal cords of a singer, the string and soundboard of a violin, or the diaphragm of a kettledrum. In signal processing, such a sound is typically represented by a function or **signal** $f\\colon\\mathbb{R}\\to\\mathbb{R}$, which encodes the sound's air pressure changes over time. The signal is called **periodic** if its values repeat at regular intervals. Intuitively speaking, the **period** of the signal is defined as the time required to complete a cycle. The **frequency**, measured in **Hertz** (Hz), is the reciprocal of the period. The prototype of such a periodic signal is a **sinusoid**, which is specified by its **frequency**, its **amplitude** (the peak deviation of the sinusoid from its mean), and its **phase** (determining where in its cycle the sinusoid is at time zero). In the following code cell, we provide a function for generating a sinusoid. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "\n", "def generate_sinusoid(dur=1, amp=1, freq=1, phase=0, Fs=100):\n", " \"\"\"Generation of sinusoid\n", "\n", " Notebook: PCP_08_signal.ipynb\n", "\n", " Args:\n", " dur: Duration (in seconds) of sinusoid (Default value = 1)\n", " amp: Amplitude of sinusoid (Default value = 1)\n", " freq: Frequency (in Hertz) of sinusoid (Default value = 1)\n", " phase: Phase (relative to interval [0,1)) of sinusoid (Default value = 0)\n", " Fs: Sampling rate (in samples per second) (Default value = 100)\n", "\n", " Returns:\n", " x: Signal\n", " t: Time axis (in seconds)\n", " \"\"\"\n", " num_samples = int(Fs * dur)\n", " t = np.arange(num_samples) / Fs\n", " x = amp * np.sin(2 * np.pi * (freq * t - phase))\n", " return x, t\n", "\n", "amp = 1.2\n", "freq = 2\n", "phase = 0.2\n", "x, t = generate_sinusoid(dur=5, amp=amp, freq=freq, phase=phase, Fs=1000)\n", "\n", "plt.figure(figsize=(8, 2.5))\n", "plt.plot(t, x, color='blue', linewidth=2.0, linestyle='-')\n", "plt.xlim([0, t[-1]])\n", "plt.ylim([-amp*1.2, amp*1.2])\n", "plt.xlabel('Time (seconds)')\n", "plt.ylabel('Amplitude')\n", "plt.title('Sinusoid with freq = %.1f, amp = %.1f, and phase = %.1f' % (freq, amp, phase))\n", "plt.grid()\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Interpreted as audio signal, one may listen to the sound of a sinusoid as long as its frequency lies in the human hearing range (roughly between $20$ and $20,000~\\mathrm{Hz}$). To integrate an audio signal into a Jupyter notebook, we can use the module [IPython.display](https://ipython.readthedocs.io/en/stable/api/generated/IPython.display.html), which is an application programming interface (API) for displaying various tools in IPython. As for audio, one can use the class [IPython.display.Audio](https://ipython.org/ipython-doc/3/api/generated/IPython.display.html), which creates an in-browser audio player. In the following code cell, this functionality is demonstrated using a sinusoid of $440~\\mathrm{Hz}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import IPython.display as ipd\n", "\n", "Fs = 4000\n", "dur = 1\n", "x, t = generate_sinusoid(dur=dur, amp=1, freq=440, Fs=Fs)\n", "\n", "ipd.display(ipd.Audio(x, rate=Fs))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "## Sampling\n", "\n", "In digital signal processing, one often reduces a **continuous-time** (CT) signal $f\\colon\\mathbb{R}\\to\\mathbb{R}$ into a **discrete-time** (DT) signal $x\\colon\\mathbb{Z}\\to\\mathbb{R}$ by a procedure known as **equidistant sampling**. Fixing a positive real number $T>0$, the DT-signal $x$ is obtained by setting \n", "\n", "\\begin{equation}\n", " x(n):= f(n \\cdot T)\n", "\\end{equation}\n", "\n", "for $n\\in\\mathbb{Z}$. The value $x(n)$ is called the **sample** taken at time $t=n\\cdot T$ of the original analog signal $f$. In short, this procedure is also called **$T$-sampling**. The number $T$ is referred to as the **sampling period** and the inverse $F_\\mathrm{s}:=1/T$ as the **sampling rate**. The sampling rate specifies the number of samples per second and is measured in Hertz (Hz). \n", "\n", "In the following code cell, we start with a CT-signal $f$ that is defined via linear interpolation of a DT-signal sampled at a high sampling rate (generated by the generate_example_function). In the plot, this CT-signal is visualized as black curve. Applying equidistant sampling (sampling_equidistant), we obtain the DT-signal visualized as the red stem plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def generate_example_signal(dur=1, Fs=100):\n", " \"\"\"Generate example signal\n", "\n", " Notebook: PCP_08_signal.ipynb\n", "\n", " Args:\n", " dur: Duration (in seconds) of signal (Default value = 1)\n", " Fs: Sampling rate (in samples per second) (Default value = 100)\n", "\n", " Returns:\n", " x: Signal\n", " t: Time axis (in seconds)\n", " \"\"\"\n", " N = int(Fs * dur)\n", " t = np.arange(N) / Fs\n", " x = 1 * np.sin(2 * np.pi * (1.9 * t - 0.3))\n", " x += 0.5 * np.sin(2 * np.pi * (6.1 * t - 0.1))\n", " x += 0.1 * np.sin(2 * np.pi * (20 * t - 0.2))\n", " return x, t\n", "\n", "def sampling_equidistant(x_1, t_1, Fs_2, dur=None):\n", " \"\"\"Equidistant sampling of interpolated signal\n", "\n", " Notebook: PCP_08_signal.ipynb\n", "\n", " Args:\n", " x_1: Signal to be interpolated and sampled\n", " t_1: Time axis (in seconds) of x_1\n", " Fs_2: Sampling rate used for equidistant sampling\n", " dur: Duration (in seconds) of sampled signal (Default value = None)\n", "\n", " Returns:\n", " x_2: Sampled signal\n", " t_2: time axis (in seconds) of sampled signal\n", " \"\"\"\n", " if dur is None:\n", " dur = len(t_1) * t_1\n", " N = int(Fs_2 * dur)\n", " t_2 = np.arange(N) / Fs_2\n", " x_2 = np.interp(t_2, t_1, x_1)\n", " return x_2, t_2\n", " \n", "Fs_1 = 100\n", "x_1, t_1 = generate_example_signal(Fs=Fs_1, dur=2)\n", "\n", "Fs_2 = 20\n", "x_2, t_2 = sampling_equidistant(x_1, t_1, Fs_2)\n", " \n", "plt.figure(figsize=(8, 2.2))\n", "plt.plot(t_1, x_1, 'k')\n", "plt.title('Original CT-signal')\n", "plt.xlabel('Time (seconds)')\n", "plt.ylim([-1.8, 1.8])\n", "plt.xlim([t_1, t_1[-1]])\n", "plt.tight_layout()\n", "\n", "plt.figure(figsize=(8, 2.2))\n", "plt.stem(t_2, x_2, linefmt='r', markerfmt='ro', basefmt='None', use_line_collection=True)\n", "plt.plot(t_1, x_1, 'k', linewidth=1, linestyle='dotted')\n", "plt.title(r'DT-Signal obtained by sampling with $F_\\mathrm{s} = %.0f$' % Fs_2)\n", "plt.xlabel('Time (seconds)')\n", "plt.ylim([-1.8, 1.8])\n", "plt.xlim([t_1, t_1[-1]])\n", "plt.tight_layout()\n", "\n", "plt.figure(figsize=(8, 2.2))\n", "plt.stem(x_2, linefmt='r', markerfmt='ro', basefmt='None', use_line_collection=True)\n", "plt.title(r'DT-Signal')\n", "plt.xlabel('Time (samples)')\n", "plt.ylim([-1.8, 1.8])\n", "plt.xlim([0, len(t_2)])\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "## Aliasing\n", "\n", "In general, sampling is a **lossy** operation in the sense that information is lost in this process and that the original CT-signal cannot be recovered from its sampled version. Only if the CT-signal has additional properties in terms of its frequency spectrum (it needs to be **bandlimited**), a perfect reconstruction is possible. This is the assertion of the famous **sampling theorem**. The sampling theorem also shows how the original CT-signal can be reconstructed by superimposing suitably shifted $\\mathrm{sinc}$-functions weighted by the samples of the DT-signal.\n", "\n", "
\n", "\n", "
\n", "Recap: Sampling Theorem
\n", "The sampling theorem, which is often associated with the names Harry Nyquist and Claude Shannon, states that a continuous-time (CT) signal that is bandlimited can be reconstructed perfectly under certain conditions. More precisely, a CT-signal $f\\in L^2(\\mathbb{R})$ is called $\\Omega$-bandlimited if the Fourier transform $\\hat{f}$ vanishes for $|\\omega|>\\Omega$ (i.e., $\\hat{f}(\\omega) = 0$ for $|\\omega|>\\Omega$). Let $f\\in L^2(\\mathbb{R})$ be an $\\Omega$-bandlimited function and let $x$ be the $T$-sampled version of $f$ with $T:=1/(2\\Omega)$ (i.e., $x(n)=f(nT)$, $n\\in\\mathbb{Z}$). Then $f$ can be reconstructed from $x$ by\n", "\n", "$$\n", " f(t)=\\sum_{n\\in\\mathbb{Z}}x(n)\\mathrm{sinc}\\left(\\frac{t-nT}{T}\\right)\n", " =\\sum_{n\\in\\mathbb{Z}}f\\left(\\frac{n}{2\\Omega}\\right) \\mathrm{sinc}\\left(2\\Omega t-n\\right),\n", "$$\n", "\n", "where the $\\mathrm{sinc}$-function is defined as\n", "\n", "\\begin{equation}\n", " \\mathrm{sinc}(t):=\\left\\{\\begin{array}{ll}\n", " \\frac{\\sin \\pi t}{\\pi t},&\\mbox{ if $t\\not= 0$,}\\\\\n", " 1,&\\mbox{ if $t= 0$.}\n", "\\end{array}\\right.\n", "\\end{equation}\n", "\n", "In other words, the CT-signal $f$ can be perfectly reconstructed from the DT-signal obtained by equidistant sampling if the bandlimit is no greater than half the sampling rate. This reconstruction based on the $\\mathrm{sinc}$-function has been used in the function reconstruction_sinc.\n", "
\n", "\n", "Without additional properties, sampling may cause an effect known as **aliasing** where certain frequency components of the signal become indistinguishable. This effect is illustrated by the following figure. Using a high sampling rate, the original CT-signal can be reconstructed with high accuracy. However, when decreasing the sampling rate, the higher-frequency components are not captured well and only a coarse approximation of the original signal remains. This phenomenon is demonstrated by the following example." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def reconstruction_sinc(x, t, t_sinc):\n", " \"\"\"Reconstruction from sampled signal using sinc-functions\n", "\n", " Notebook: PCP_08_signal.ipynb\n", "\n", " Args:\n", " x: Sampled signal\n", " t: Equidistant discrete time axis (in seconds) of x\n", " t_sinc: Equidistant discrete time axis (in seconds) of signal to be reconstructed\n", "\n", " Returns:\n", " x_sinc: Reconstructed signal having time axis t_sinc\n", " \"\"\"\n", " Fs = 1 / t\n", " x_sinc = np.zeros(len(t_sinc))\n", " for n in range(0, len(t)):\n", " x_sinc += x[n] * np.sinc(Fs * t_sinc - n)\n", " return x_sinc\n", "\n", "def plot_signal_reconstructed(t_1, x_1, t_2, x_2, t_sinc, x_sinc, figsize=(8, 2.2)):\n", " \"\"\"Plotting three signals\n", "\n", " Notebook: PCP_08_signal.ipynb\n", "\n", " Args:\n", " t_1: Time axis of original signal\n", " x_1: Original signal\n", " t_2: Time axis for sampled signal\n", " x_2: Sampled signal\n", " t_sinc: Time axis for reconstructed signal\n", " x_sinc: Reconstructed signal\n", " figsize: Figure size (Default value = (8, 2.2))\n", " \"\"\"\n", " plt.figure(figsize=figsize)\n", " plt.plot(t_1, x_1, 'k', linewidth=1, linestyle='dotted', label='Orignal signal')\n", " plt.stem(t_2, x_2, linefmt='r:', markerfmt='r.', basefmt='None', label='Samples', use_line_collection=True)\n", " plt.plot(t_sinc, x_sinc, 'b', label='Reconstructed signal')\n", " plt.title(r'Sampling rate $F_\\mathrm{s} = %.0f$' % (1/t_2))\n", " plt.xlabel('Time (seconds)')\n", " plt.ylim([-1.8, 1.8])\n", " plt.xlim([t_1, t_1[-1]])\n", " plt.legend(loc='upper right', framealpha=1)\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "Fs_2 = 40\n", "x_2, t_2 = sampling_equidistant(x_1, t_1, Fs_2)\n", "t_sinc = t_1\n", "x_sinc = reconstruction_sinc(x_2, t_2, t_sinc)\n", "plot_signal_reconstructed(t_1, x_1, t_2, x_2, t_sinc, x_sinc)\n", "\n", "Fs_2 = 20\n", "x_2, t_2 = sampling_equidistant(x_1, t_1, Fs_2)\n", "t_sinc = t_1\n", "x_sinc = reconstruction_sinc(x_2, t_2, t_sinc)\n", "plot_signal_reconstructed(t_1, x_1, t_2, x_2, t_sinc, x_sinc)\n", "\n", "Fs_2 = 10\n", "x_2, t_2 = sampling_equidistant(x_1, t_1, Fs_2)\n", "t_sinc = t_1\n", "x_sinc = reconstruction_sinc(x_2, t_2, t_sinc)\n", "plot_signal_reconstructed(t_1, x_1, t_2, x_2, t_sinc, x_sinc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "## Interference\n", "\n", "In signal processing, **interference** occurs when a wave is superimposed with another wave that has the same or nearly the same frequency. When a crest of one wave meets a crest of the other wave at some point, then the individual magnitudes add up for a certain period of time, which is known as **constructive interference**. Vice versa, when a crest of one wave meets a trough of the other wave, then the magnitudes cancel out for a certain period of time, which is known as **destructive interference**. We illustrate these phenomena in the following code cell. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_interference(t, x1, x2, figsize=(8, 2), xlim=None, ylim=None, title=''):\n", " \"\"\"Plotting two signals and its superposition\n", "\n", " Notebook: PCP_08_signal.ipynb\n", "\n", " Args:\n", " t: Time axis\n", " x1: Signal 1\n", " x2: Signal 2\n", " figsize: Figure size (Default value = (8, 2))\n", " xlim: x-Axis limits (Default value = None)\n", " ylim: y-Axis limits (Default value = None)\n", " title: Figure title (Default value = '')\n", " \"\"\" \n", " plt.figure(figsize=(8, 2))\n", " plt.plot(t, x1, color='gray', linewidth=1.0, linestyle='-', label='x1')\n", " plt.plot(t, x2, color='cyan', linewidth=1.0, linestyle='-', label='x2')\n", " plt.plot(t, x1+x2, color='red', linewidth=2.0, linestyle='-', label='x1+x2')\n", " if xlim is None:\n", " plt.xlim([0, t[-1]])\n", " else:\n", " plt.xlim(xlim)\n", " if ylim is not None:\n", " plt.ylim(ylim)\n", " plt.xlabel('Time (seconds)')\n", " plt.ylabel('Amplitude')\n", " plt.title(title)\n", " plt.legend(loc='upper right')\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "dur = 5\n", "x1, t = generate_sinusoid(dur=dur, Fs=100, amp=1, freq=1.05, phase=0.0)\n", "x2, t = generate_sinusoid(dur=dur, Fs=100, amp=1, freq=0.95, phase=0.8)\n", "plot_interference(t, x1, x2, xlim=[0, dur], ylim=[-2.2,2.2], title='Constructive Interference')\n", "\n", "dur = 5\n", "x1, t = generate_sinusoid(dur=dur, Fs=100, amp=1, freq=1.05, phase=0.0)\n", "x2, t = generate_sinusoid(dur=dur, Fs=100, amp=1, freq=1.00, phase=0.4)\n", "plot_interference(t, x1, x2, xlim=[0, dur], ylim=[-2.2,2.2], title='Destructive Interference')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises and Results" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import libpcp.signal\n", "show_result = True" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "
\n", "Exercise 1: Beating
\n", "Let $f_1(t)=\\sin(2\\pi \\omega_1 t)$ and $f_2(t)=\\sin(2\\pi \\omega_2 t)$ be two sinusoids with distinct but nearby frequencies $\\omega_1\\approx\\omega_2$. The superposition $f_1+f_2$ of these two sinusoids results in a function that looks like a single sine wave with a slowly varying amplitude, a phenomenon also known as beating. Mathematically, this phenomenon results from a trigonometric identity yielding \n", "\n", "$$\n", "\\sin(2\\pi \\omega_1 t)+\\sin(2\\pi \\omega_2 t)=\n", "2\\cos\\left(2\\pi\\frac{\\omega_1-\\omega_2}{2}t\\right)\\sin\\left(2\\pi\\frac{\\omega_1+\\omega_2}{2}t\\right).\n", "$$\n", "\n", "If the difference $\\omega_1-\\omega_2$ is small, the cosine term has a low frequency compared to the sine term. As a result the signal $f_1+f_2$ can be seen as a sine wave of frequency $(\\omega_1+\\omega_2)/2$ with a slowly varying amplitude envelope of frequency $|\\omega_1-\\omega_2|$. Note that this rate is twice the frequency $(\\omega_1-\\omega_2)/2$ of the cosine term. \n", "
\n", "
• Use the functions generate_sinusoid and plot_interference to illustrate the beating effect for various parameters $\\omega_1$ and $\\omega_2$.
• \n", "
• In particular, consider $\\omega_1=200$ and $\\omega_2=203$ to generate a superimposed signal x at a sampling rate of Fs=4000. Listen to x using ipd.display(ipd.Audio(x, rate=Fs)).
• \n", "
\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#\n", "# Your Solution\n", "#" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "libpcp.signal.exercise_beating(show_result=show_result) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "
\n", "Exercise 2: Aliasing with Sinsuoids
\n", "First, using x, t = generate_sinusoid(dur=2, Fs=128, freq=10), generate a sinusoid of frequency $\\omega=10~\\mathrm{Hz}$ at Fs=128 to mimic a CT-signal. Then, use the function sampling_equidistant to successively decrease the sampling rate using [64, 32, 20, 16, 12, 8, 4]. Visualize the result via plot_signal_reconstructed and discuss your observations.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#\n", "# Your Solution\n", "#" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "libpcp.signal.exercise_aliasing_sinus(show_result=show_result) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "
\n", "Exercise 3: Visual Aliasing
\n", "Aliasing effects occur also in other areas such as the visual domain. One famous example is the wagon-wheel effect, where a wheel appears to rotate at a slower speed or even backwards when the speed is actually increasing. Identify some YouTube videos that illustrate this effect. Furthermore, integrate the video into your notebooks. To this end, do the following:\n", "
\n", "
• Import the API for display tools in IPython: import IPython.display as ipd.
• \n", "
• Use the following commands: ipd.display, ipd.YouTubeVideo
• \n", "
\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#\n", "# Your Solution\n", "#" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "libpcp.signal.exercise_aliasing_visual(show_result=show_result) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "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.9.17" } }, "nbformat": 4, "nbformat_minor": 1 }