{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-02-09T09:53:16.947062400Z", "start_time": "2024-02-09T09:53:16.683055600Z" }, "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "import numpy as np\n", "from qutip import *\n", "import numpy as np\n", "import matplotlib.transforms as transforms\n", "from matplotlib.animation import FuncAnimation\n", "from IPython.display import HTML\n", "\n", "from scipy import signal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Lecture 1: Noise in measurements\n", "\n", "```{admonition} Expected prior knowledge\n", ":class: tip\n", "Before the start of this lecture, you should be able to:\n", "\n", "- recall the definition of a Fourier transform \n", "```\n", "\n", "```{admonition} Learning goals\n", ":class: important\n", "After this lecture you will be able to:\n", "\n", "- estimate average quantities and their statistical error from noisy measurements\n", "- characterize noise via its power spectral density\n", "- identify common kinds of noise spectra\n", "``` \n", "\n", "Consider the following problem. We have a small voltage that we want\n", "to measure. We assume the voltage to be constant in time. To measure it, we connect the small voltage to an amplifier and then\n", "connect our amplified signal to an oscilloscope that records the voltage at its\n", "input as a function of time. \n", "\n", "[
](scope.png)\n", "\n", "On our oscilloscope screen, we see the following trace:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-02-09T09:53:22.645987100Z", "start_time": "2024-02-09T09:53:22.462972800Z" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "mu, sigma = 0, 1\n", "ts = np.linspace(0,100,100)\n", "vs = np.asarray([np.random.normal(mu,sigma) for t in ts])\n", "\n", "plt.plot(ts, vs)\n", "plt.title(\"Signal\")\n", "plt.xlabel(\"Time (s)\")\n", "plt.ylim(-3,3)\n", "plt.ylabel(\"Voltage (mV)\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What is approximately the value of the voltage?\n", "- the average of the trace seems to be about 0 mV\n", "- How can we quantify how close it is to 0 mV?\n", "\n", "A good first step: create a histogram of the data, which might look like" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-09-27T19:15:29.463059Z", "start_time": "2021-09-27T19:15:29.301727Z" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "count, bins, ignored = plt.hist(vs, 20, density=True)\n", "plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (bins - mu)**2 / (2 * sigma**2) ),\n", " linewidth=2, color='r')\n", "plt.title(\"Histogram of measured voltages and fit to a Gaussian probability distribution\")\n", "plt.xlabel(\"Voltage [mV]\")\n", "plt.ylabel(\"Probability of voltage value\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(In this histogram, we have normalized the bin counts to match the Gaussian fit.)\n", "\n", "Fitting to a Gaussian, we find a full-width half-maximum (FWHM) of 2 mV. We thus estimate that, on average, the statistical error in each measurement point is about $\\pm 1$ mV.\n", "\n", "However, because we assumed that the voltage we are trying to measure is constant in time, we can reduce the error in our estimate of the \n", "unknown small voltage by averaging points together. We know how to calculate the average, but what is the uncertainty in this average value? \n", "\n", "```{admonition} Standard error of the mean\n", ":class: note\n", "The error $\\sigma_N$ of our average of $N$ points is related to the error $\\sigma_1$ of a single point by \n", "\n", "$$\n", "\\sigma_N^2 = \\frac{\\sigma_1^2}{N}\n", "$$\n", "\n", "The quantity $\\sigma_N$ is known as the standard error of the mean. How is this quantity related to our histogram above? Adding\n", "more points will not change the *shape* of the distribution, it will only give \n", "more counts in each bin. \n", "```\n", "\n", "The reduction of $\\sigma_N$ with respect to $\\sigma_1$, however, has nothing to do with the width of the histogram. Rather, it determines how accurately we can determine the center of the peak (the average): The more points we add, the smaller the (relative) statistical fluctuations in the height of each bin will become, and the more accurately we will be able to find the center of the peak.\n", "\n", "```{admonition} Caveat\n", ":class: warning\n", "The above equation relating the standard error of the mean $\\sigma_N$ to the standard error of a single measurement $\\sigma_1$ is only true if the noise in your data is *uncorrelated*, which means that each data point in your trace represents an independent sampling of the random noise of your amplifier. \n", "```\n", "\n", "```{admonition} Autocorrelation function\n", ":class: note\n", "How can one know if the noise is uncorrelated? By calculating its autocorrelation function. For a function $v(t)$ it is defined as\n", "\n", "$$\n", "R_{vv}(\\tau) = \\langle v(t+\\tau) v(t) \\rangle\n", "$$\n", "\n", "where the angled brackets denote the expectation value/ensemble average.\n", "```\n", "\n", "A common example of correlated noise arises when your signal has some \"memory\" that is slower than your sampling rate such that consecutive samples will have similar values. For example, if you put a 1 kHz low-pass\n", "filter after your amplifier, but configure your oscilloscope to record the signal\n", "with a 10 kHz sampling rate, then $R_{vv}$ will look something like this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-09-27T19:15:29.836673Z", "start_time": "2021-09-27T19:15:29.465757Z" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "ts = np.linspace(0,3,31)\n", "ys = np.exp(-ts)\n", "plt.plot(ts, ys)\n", "plt.title(\"Autocorrelation function\")\n", "plt.ylabel(\"$R_{vv}$\")\n", "plt.xlabel(\"Time difference between data points (ms)\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we can see that each measured value displayed on our oscilloscope will be strongly correlated with several values in its direct neighborhood. In this case, the values are not independently drawn from the random distribution, and the error in our estimate of the average voltage will *not* scale with $\\frac{1}{N}$ if we average $N$ such points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Power Spectral Density\n", "\n", "To summarize so far:\n", "\n", "- The error in estimating the average of a noisy measurement depends on how many (uncorrelated) points we measure\n", "- If we record a dataset of uncorrelated points for a time $T_m$, the error in our estimate of the mean value will go down as $\\sigma_T^2 \\propto \\frac{1}{T_m}$\n", "\n", "However, if we measure for longer periods of time, then our measurements are\n", "slower. This becomes an issue if we want to measure a signal that is not constant in time, because increasing $T_m$ reduces\n", "the speed at which we can record a changing signal. Thus, there is a tradeoff\n", "between the noise level $\\sigma_v$ and the *bandwidth* of your measurement.\n", "\n", "```{admonition} Question\n", ":class: question\n", "if you want buy an amplifier, does it make sense to look at its noise level without looking at its bandwidth? Why (not)?\n", "```\n", "\n", "So then, how do we characterize noisy, time-varying signals? We can compute\n", "the power spectral density $S_{xx}(\\omega)$, which tells us how the \"noise\n", "power\" is distributed across different frequencies.\n", "\n", "```{admonition} Definitions\n", ":class: discussion\n", "1. For a signal $V(t)$, we first define a (finite-time, truncated) Fourier \n", "transform:\n", "\n", "$$\n", "\\hat{V}(\\omega) = \\frac{1}{\\sqrt{T}} \\int_0^T V(t)e^{-i \\omega t} dt\n", "$$\n", "\n", "where the prefactor $\\frac{1}{\\sqrt{T}}$ is for normalization.\n", "\n", "2. The power spectral density is then defined as the \"ensemble average\" of the magnitude\n", "of $\\hat{V}(\\omega)$:\n", "\n", "$$\n", "S_{VV}(\\omega) = \\big \\langle | \\hat{V}(\\omega) |^2 \\big \\rangle\n", "$$\n", "\n", "3. What do we mean by \"ensemble average\"?\n", "\n", "- If the noise in our measurement is random, its Fourier transform is also \n", "random\n", "- To reduce the \"noise\" in our estimate of the noise, we average over the \n", "$|\\hat{V}(\\omega)|^2$ of many different measurement traces\n", "\n", "The power spectral density describes how the power in a signal is distributed over different frequencies. It is also called the power spectrum. \n", "```\n", "\n", "How could we measure a power spectral density in practice?\n", "\n", "[
](filter.png)\n", "\n", "Here, the bandpass filter transmits the part of the signal within a frequency band $\\Delta \\omega$ centered at $\\omega_0$. The rectifier converts the transmitted signal into a DC voltage $V_\\text{out}$. As a result, $V_\\text{out}$ is proportional to the amount of power present in $v(t)$ in the frequency band $\\omega_0\\pm\\Delta\\omega/2$, given by $ V_{out}^2 = S_\\text{vv}(\\omega_0) {\\Delta \\omega} $\n", "\n", "## Types of noise specta\n", "\n", "### White noise\n", "Uncorrelated noise is called 'white noise'. Its power spectral density is flat (independent of frequency)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-09-27T19:15:30.241593Z", "start_time": "2021-09-27T19:15:29.839261Z" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "fs = 1e3\n", "N = 1e5\n", "amp = 2*np.sqrt(2)\n", "freq = 1234.0\n", "noise_power = 0.001 * fs / 2\n", "time = np.arange(N) / fs\n", "x = 0\n", "# x = amp*np.sin(2*np.pi*freq*time)\n", "x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)\n", "\n", "f, Pxx_den = signal.welch(x, fs, nperseg=len(x)//4)\n", "plt.semilogy(f, Pxx_den)\n", "plt.ylim([1e-7, 1e2])\n", "plt.title(\"White Noise\")\n", "plt.xlabel('frequency [Hz]')\n", "plt.ylabel('PSD [V**2/Hz]')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1/f noise\n", "Common in many experiments, particularly for f < 100kHz" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-09-27T19:15:30.739820Z", "start_time": "2021-09-27T19:15:30.246532Z" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "plt.semilogy(f[1:], Pxx_den[1:]/f[1:])\n", "#plt.ylim([1e-7, 1e2])\n", "plt.xlim([0,200])\n", "plt.title(\"1/f Noise (also called 'pink' noise)\")\n", "plt.xlabel('frequency [Hz]')\n", "plt.ylabel('PSD [V**2/Hz]')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "plt.semilogy(f[1:], Pxx_den[1:]/f[1:])\n", "#plt.ylim([1e-7, 1e2])\n", "plt.xlim([0,200])\n", "plt.title(\"1/f Noise (also called 'pink' noise)\")\n", "plt.xlabel('frequency [Hz]')\n", "plt.ylabel('PSD [V**2/Hz]')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 0 }, "source": [ "### \"Peaked\" noise" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2021-09-27T19:15:31.237172Z", "start_time": "2021-09-27T19:15:30.742009Z" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "plt.semilogy(f, Pxx_den/(1+(f-100)**2))\n", "plt.ylim([1e-8, 1e-2])\n", "plt.xlim([0,200])\n", "plt.title(\"Peaked noise\")\n", "plt.xlabel('frequency [Hz]')\n", "plt.ylabel('PSD [V**2/Hz]')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{admonition} Calculating the noise in a measurement from a known power spectral density\n", ":class: example\n", "Suppose the spec sheet of our amplifier specifies $S_\\text{vv}=1 \\, \\text{nV}^2/\\text{Hz}$ within the amplifier bandwidth (0-1 MHz). We can then calculate the expected fluctuations $\\sigma_v$ of the data points in a timetrace $V(t)$ that we would measure with an oscilloscope, using\n", "\n", "$$\n", "\\sigma_v^2 = \\int_{\\omega_1}^{\\omega_2} S_\\text{vv}(\\omega) d\\omega\n", "$$\n", "\n", "What determines $\\omega_1$ and $\\omega_2$? If we measure for a total time $T_m$, we are insensitive to frequencies below $\\omega_1 = 1/T_m$. Furthermore, $\\omega_2$ is given by the upper limit of the bandwidth of our amplifier (1 MHz), assuming that our oscilloscope does not limit the measurement bandwidth. If $T_m= 1$ s, then \n", "\n", "$$\n", "\\begin{align}\n", "\\sigma_\\text{v}^{2} &= \\int_{1}^{10^6} \\frac{1 \\, \\text{nV}^2}{\\text{Hz}} df\\\\\n", "&= 10^{-18} (10^6 - 1) \\approx 10^{-12} \\, \\text{V}^2\\\\\n", "\\sigma_\\text{v} &= 1 \\, \\mu\\text{V}\n", "\\end{align}\n", "$$\n", "\n", "so the RMS noise on the measurement will be 1$\\mu $V\n", "\n", "What if we put a 1 kHZ low-pass filter on the output of our amplifier?\n", "\n", "$$ \n", "\\begin{align}\n", "\\sigma_\\text{v}^{2} &= \\int_{1}^{10^3} \\frac{1 \\, \\text{nV}^2}{\\text{Hz}} df \\\\\n", "= & 10^{-15} \\, \\text{V}^2\\\\\n", "\\sigma_\\text{v} &= 31 \\, \\text{nV}\n", "\\end{align}\n", "$$\n", "```\n", "\n", "```{admonition} Conclusions\n", ":class: important\n", "1. Measuring longer decreases the uncertainty in the average of the measured quantity according to $\\sigma_N^2= \\sigma_1^2/N $, provided the N data points are uncorrelated.\n", "2. The power spectral density $S_{xx}(\\omega)$ is the Fourier transform of the autocorrelation function $R_{xx}(t)$:\n", "\n", "$$\n", "S_{xx}(\\omega) = \\int_{-\\infty}^{+\\infty} R_{xx}(t) e^{-i\\omega t} dt\n", "$$\n", "\n", "3. For uncorrelated noise, $R_{xx}(\\tau) = \\delta(\\tau)$, such that $S_{xx}(\\omega)$ is a constant (white noise)\n", "4. The expected noise in a data trace can be calculated by integrating over the power spectral density\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Demonstration: Noise after averaging \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# Increase the animation embed limit\n", "plt.rcParams['animation.embed_limit'] = 60\n", "\n", "# Start with Gaussian noise\n", "N = 10000\n", "x = np.random.normal(size=N)\n", "\n", "Navg = np.geomspace(1, N / 20, 100).astype(int) # 100 points between 1 and N/20\n", "Navg = np.unique(Navg) # remove duplicates\n", "\n", "def average_points(x, n):\n", " n_out = len(x) // n\n", " # A trick for using reshape\n", " x = x[0:n * n_out]\n", " x = np.reshape(x, [n, n_out])\n", " return np.average(x, axis=0)\n", "\n", "averages = []\n", "sigma = np.zeros(len(Navg))\n", "# calculates the averages for different values of n and stores them in the averages list.\n", "# also calculates the standard deviation of these averages and stores them in the sigma array\n", "for i in range(len(Navg)):\n", " averages.append(average_points(x, Navg[i]))\n", " sigma[i] = np.std(averages[i])\n", "\n", "i = 0\n", "N = 2000\n", "c = 'orange'\n", "\n", "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 12), dpi=200)\n", "ax1.set_xscale('log')\n", "ax1.set_yscale('log')\n", "ax1.set_title(\"N_avg = %d\" % Navg[i])\n", "l1, = ax1.plot(Navg, sigma)\n", "cir1 = ax1.scatter(Navg, sigma)\n", "cir2 = ax1.scatter(Navg[i:i+1], sigma[i:i+1], color=c)\n", "\n", "l2, = ax2.plot(range(N), x[0:N])\n", "cir3 = ax2.scatter(range(N), x[0:N])\n", "xax = np.linspace(0, N, N // Navg[i])\n", "yax = averages[i][0:N // Navg[i]]\n", "l3, = ax2.plot(xax, yax, color=c)\n", "cir4 = ax2.scatter(xax, yax, color=c, facecolors='none')\n", "\n", "def update_plot(i):\n", " n = Navg[i] # Navg[i] is the number of points to average\n", " ax1.set_title(\"N_avg = %d\" % n)\n", " cir2.set_offsets(np.column_stack((Navg[i:i+1], sigma[i:i+1])))\n", " xax = np.linspace(0, N, N // n)\n", " yax = averages[i][0:N // n]\n", " l3.set_data(xax, yax)\n", " cir4.set_offsets(np.column_stack((xax, yax)))\n", "\n", "frames = range(len(Navg)) # include all frames\n", "ani = FuncAnimation(fig, update_plot, frames=frames)\n", "\n", "plt.close(fig)\n", "# show the animation\n", "HTML(ani.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Demonstration: \"Ensemble averaging\" PSD (reducing the noise of the noise...)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# Increase the animation embed limit\n", "plt.rcParams['animation.embed_limit'] = 60\n", "\n", "# Number of points in an individual data trace\n", "N = 1000\n", "\n", "# Number of \"ensemble\" members to make\n", "M = 500\n", "\n", "v = np.random.normal(size=(N, M))\n", "\n", "# Take the Fourier transform\n", "vt = np.fft.fft(v, axis=0)\n", "f = np.fft.fftfreq(len(v))\n", "\n", "# Make our noise spectrum \"peaky\" so it looks a bit interesting\n", "f0 = max(f) / 3\n", "Q = 10\n", "delta = f - f0\n", "filt = f0 ** 2 / (f0 ** 2 - f ** 2 + 1j * f * f0 / Q)\n", "vt2 = np.zeros(vt.shape) + 0j\n", "for i in range(vt2.shape[1]):\n", " vt2[:, i] = vt[:, i] * filt\n", "\n", "# Create the figure and axes\n", "fig, ax = plt.subplots(figsize=(8, 6) , dpi=150)\n", "line1, = ax.plot(f[:N // 3], np.abs(vt2[:N // 3, 0]), color='blue', lw=1)\n", "line2, = ax.plot(f[:N // 3], np.abs(vt2[:N // 3, 0]), 'ro', markersize=1.2)\n", "\n", "\n", "# Set the axis labels and title\n", "ax.set_xlabel('Frequency')\n", "ax.set_ylabel('Amplitude')\n", "ax.set_title('Ensemble Averaging')\n", "\n", "# Update function for the animation\n", "def update(frame):\n", " avg = np.average(np.abs(vt2[:N // 3, :frame + 1]), axis=1)\n", " line1.set_ydata(avg)\n", " line2.set_ydata(np.abs(vt2[:N // 3, -1]))\n", " ax.set_title('Ensemble Averaging, N_avg = %d' % (frame + 1))\n", " return line1, line2\n", "\n", "# Create the animation\n", "animation = FuncAnimation(fig, update, frames=M, interval=200, blit=True)\n", "\n", "# Display the animation\n", "plt.close(fig)\n", "HTML(animation.to_jshtml())\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.12.5" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }