diff --git a/data/opt_lev_comb_drive.npy b/data/opt_lev_comb_drive.npy new file mode 100644 index 0000000..ae4a902 Binary files /dev/null and b/data/opt_lev_comb_drive.npy differ diff --git a/data/opt_lev_comb_response.npy b/data/opt_lev_comb_response.npy new file mode 100644 index 0000000..7ac5255 Binary files /dev/null and b/data/opt_lev_comb_response.npy differ diff --git a/data/opt_lev_drive.npy b/data/opt_lev_drive.npy new file mode 100644 index 0000000..93aacdb Binary files /dev/null and b/data/opt_lev_drive.npy differ diff --git a/data/opt_lev_response.npy b/data/opt_lev_response.npy new file mode 100644 index 0000000..518ce03 Binary files /dev/null and b/data/opt_lev_response.npy differ diff --git a/nb/Project_Optical_Tweezers_Calibration.ipynb b/nb/Project_Optical_Tweezers_Calibration.ipynb index 29f7dac..279d05d 100644 --- a/nb/Project_Optical_Tweezers_Calibration.ipynb +++ b/nb/Project_Optical_Tweezers_Calibration.ipynb @@ -2,12 +2,13 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", + "import scipy.signal as signal\n", "%matplotlib inline" ] }, @@ -17,23 +18,284 @@ "source": [ "# Project: Calibrating an Optical Tweezer\n", "\n", - "An [optical tweezer](https://en.wikipedia.org/wiki/Optical_tweezers) is essentially a tightly focused laser-beam, where the intensity gradient in three-dimensions tends to confine dielectric particles of order $100~\\rm{nm}$ to $10~\\mu\\rm{m}$ at the location of highest intensity, i.e. the central spot of the focus. When a dielectric particle is trapped within the optical tweezer, it acts like opto-mechanical system that, for certain conditions, is well-approximated by a mass-spring system. If some external force pushes or pulls on the dielectric particle, it gets displaced from its equilibrium position, which deflects some of the optical tweezer's light and results in an optical restoring force. The restoring force is proportional to the driving force for sufficiently small displacements.\n", + "An [optical tweezer](https://en.wikipedia.org/wiki/Optical_tweezers) is essentially a tightly focused laser-beam, where the intensity gradient in three-dimensions tends to confine dielectric particles (usually with radii of order $100~\\rm{nm}$ to $30~\\mu\\rm{m}$) at the location of highest intensity, i.e. the central spot of the focus. When a dielectric particle is trapped within the optical tweezer, it acts like opto-mechanical system that, for certain conditions, is well-approximated by a mass-spring system. If some external force pushes or pulls on the dielectric particle, it gets displaced from its equilibrium position, which deflects some of the optical tweezer's light and results in an optical restoring force. The restoring force is proportional to the driving force for sufficiently small displacements of the particle, i.e. $F_{\\rm opt} = -k_{\\rm trap} x$, where $x$ the is position of the particle along a single axis, relative to the equilibrium position. The system can thus be modeled as a **driven damped harmonic oscillator**, like the RLC circuit we've studied in Physics 81.\n", "\n", - "\n", + "\n", "\n", - "Making use of this property, it is possible to construct a force sensor with sensitivity to very miniscule forces, and perform precision physics experiments. Like all good appratuses\n", + "Making use of this property, it is possible to construct a force sensor with sensitivity to very miniscule forces, and perform precision physics experiments. Like all good appratuses, this one needs to be calibrated. This can be accomplished by pushing on the trapped particle with a known force at some fixed frequency, and then observing the response at that same frequency. We use an oscillating driving force because measurements are often difficult at DC. If we expect the system to have a different response at different frequencies, then we can drive it at multiple frequencies and observe the response at each frequency independently. This procedure measures the **transfer function** of the mechanical system, which is related to it's **mechanical susceptibility**.\n", + "\n", + "To properly analyze such a calibration measurement, we'll need to apply the principles of Fourier Analysis, and compare two different signals: the driving signal and the response signal.\n", + "\n", + "### The apparatus and the measurement\n", + "\n", + "The experiment and the data involved in this project are based on an earlier version of the apparatus discussed in [this research paper](https://pubs.aip.org/aip/rsi/article/91/8/083201/989394/High-sensitivity-levitated-microsphere-apparatus). Essentially, an infrared laser is brought to a very tight focus inside a vacuum chamber, and silica microspheres of radii $5~\\mu\\rm{m}$ to $15~\\mu\\rm{m}$ are trapped at this focus. The microsphere and optical trap are surrounded by 6 pyramidal electrodes that allow the experimenter to drive an electric field with any direction and a wide range of magnitudes. The microsphere's motion is monitored by an optical imaging system that produces some voltage that is proportional to the displacement of the microsphere from equilibrium. Before calibration, this signal is in arbitrary units as a result of signficant digital signal processing, but we'd like to know the actual displacement/force, so we need a calibration factor from (Arbitrary units) -> (Newtons).\n", + "\n", + "As part of the trapping process, the microspheres end up electrically charged. This charge can be monitored and controlled with single electron precision, and brought to a state of *exact* electrically neutrality (equal numbers of protons and electrons out of a total of $\\sim 10^{13}$). This measurement is detailed in [this other research paper](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.113.251801). To calibrate the apparatus, an electric field with known magnitude is applied to a microsphere with a single electron charge, yielding a known force via the expression $F = q E$, where $q = -e$ is the charge of a single electron. The detector response (in some arbitrary units) can then be calibrated to Newtons by comparing the driving force to the measured response.\n", + "\n", + "Once calibrated, the optically levitated force sensor can be used to test the neutrality of matter, or test the universal law of gravitation at short-range, where various Beyond the Standard Model (BSM) theories may manifest as short-range forces that couple to mass. We won't approach that complex of a problem in this project, rather, we'll work with some calibration data to understand the response of the force sensor.\n", + "\n", + "We've provided you with a few different time-series datafiles, all with a sampling rate of $f_{\\rm samp} = 5000~\\rm{Hz}$.\n", + "1. `opt_lev_drive.npy`: a 10-s long data file with a single sinusoidal drive, with data in units of Volts/meter (electric field)\n", + "2. `opt_lev_response.npy`: a 10-s long data file with the response of the microsphere, with data in arbitrary detector units\n", + "3. `opt_lev_comb_drive.npy`: a 50-s long data file with a frequency comb drive (i.e. many superposed sinusoids), again in V/m\n", + "4. `opt_lev_comb_response.npy`: a 50s long data file with the resposne of the microsphere to the frequency comb, again in arbitrary detector units\n", + "\n", + "This data is all taken for a single cartesian axis of drive and response, although the microsphere can be driven and measured fully in three dimensions.\n", + "\n", + "The response data files include the microsphere's response to stochastic driving forces from the environment (mainly intensity and pointing fluctuations in the trapping laser) so we will have to be aware of this in our analysis. To see the single-frequency drive and response by themselves, we'll have to filter our data in the time-domain to isolate the signal we're interested in. Fourier analysis will allow us to calibrate the data *without* explicitly filtering it.\n", "\n", "### Potential goals for this project\n", "\n", - "We expect you to complete all of the following goals for full credit. If you prefer to replace on of these goals with your own goal, please consult the instructor/TA.\n", + "We expect you to complete the first two of the following goals for full credit. If you prefer to replace on of these goals with your own goal, please consult the instructor/TA.\n", + "\n", + "1. Extract the amplitude and phase of the single-frequency response, relative to the drive. If we model our drive as a single tone, $f_{\\rm drive} = A_{\\rm drive} \\cos( \\omega t )$, then a general solution for a driven damped harmonic oscillator is given by $f_{\\rm resp} = A_{\\rm resp} \\cos( \\omega t + \\phi)$.\n", + " * If we take the Fourier transform of a signal, like with `np.fft.rfft`, the value of the transform at a particular frequency is generally a complex number, i.e. $\\tilde{f}(\\omega) = F[f(t)]$ with $F[\\dots]$ the Fourier transform and with $\\tilde{f}$ being complex\n", + " * This is because most implementations of an fast fourier transform (FFT) make use of complex phasors of the form $\\tilde{A}(t) = \\tilde{A} e^{-i \\omega t} = A e^{-i (\\omega t + \\phi)}$\n", + " * The amplitude and phase of the response, relative to the drive, can be evaluated in the frequency domain by taking the absolute value of the ratio of fourier transforms and the arg of the same ratio (where 'arg' is complex-plane equivalent to the arctangent function and corresponds to the phasor rotation angle between drive and response)\n", + " $$ A(\\omega) = \\frac{A_{\\rm resp}}{A_{\\rm drive}} = \\left| \\frac{\\tilde{f}_{\\rm resp}(\\omega)}{\\tilde{f}_{\\rm drive}(\\omega)} \\right| \\qquad \\text{and} \\qquad \\phi(\\omega) = \\phi_{\\rm resp} - \\phi_{\\rm drive} = \\text{Arg} \\left( \\frac{\\tilde{f}_{\\rm resp}(\\omega)}{\\tilde{f}_{\\rm drive}(\\omega)} \\right) $$\n", + "\n", + "2. Extend this analysis to a frequency-comb drive/response and thus qualitatively determine the functions $A(\\omega)$ and $\\phi(\\omega)$\n", + " * The frequency comb calibration allows us to measure many different drive frequencies simultaneously.\n", + " * Combining the amplitude ratios and phase differences for every drive frequency, we can find $A(\\omega)$ and $\\phi(\\omega)$\n", + "\n", + "3. (optional if interested) Use multiple calibration measurements and combine them to compute a mean value of the calibration factor and some standard deviation, corresponding to our certainty of the calibration factor, given the fluctuations present in the system." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Loading the data\n", + "\n", + "We can load the data into memory using NumPy. Note that the datafiles were stored as `.npy` files, which is essentially a binary file that takes advantage of compression to decrease file size (as opposed to a `.txt` file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Sampling frequency, common to all datafiles\n", + "fsamp = 5000.0\n", + "\n", + "drive = np.load('../data/opt_lev_drive.npy')\n", + "response = np.load('../data/opt_lev_response.npy')\n", + "\n", + "# Make an associated time vector for the drive/response\n", + "nsamp = len(drive)\n", + "tvec = np.arange(nsamp) / fsamp\n", + "\n", + "comb_drive = np.load('../data/opt_lev_comb_drive.npy')\n", + "comb_response = np.load('../data/opt_lev_comb_response.npy')\n", + "nsamp_comb = len(comb_drive)\n", + "\n", + "# Make an associated time vector for the frequency comb drive/response\n", + "tvec_comb = np.arange(nsamp_comb) / fsamp" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plotting the time-series\n", + "\n", + "It's always good practice to actually look at your raw data before doing any complicated analysis. In truth, the data we've provided to you has had some pre-processing done on it to make it more palatable, but it's essentially the raw drive signal and raw microsphere response.\n", + "\n", + "The drive frequency is sufficiently high, that it will be hard to see individual cycles of the driving sinusoid if we plot all 10-seconds of data at the same time, so we'll plot a subset of data: the first half-second." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "## Create a mask to only look at a subset of the data. The default value is \n", + "## the first 0.5 seconds of data, but you should change this and see what happens.\n", + "mask = tvec < 0.5\n", + "\n", + "fig, ax = plt.subplots(2, 1, figsize=(10, 6))\n", + "ax[0].plot(tvec[mask], drive[mask])\n", + "ax[1].plot(tvec[mask], response[mask])\n", + "\n", + "ax[0].set_ylabel('E-field Drive [V/m]')\n", + "ax[1].set_ylabel('Detector Response [arb]')\n", + "ax[1].set_xlabel('Time [s]')\n", + "\n", + "fig.tight_layout()\n", + "plt.show" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Things to think about:\n", + "\n", + "##### 1.1 Does the drive signal look like a single frequency sinusoid?\n", + "\n", + "##### 1.2 What about the microsphere response? Can you discern anything looking at just the time-domain?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Filtering the data\n", + "\n", + "At first glance it's hard to see what's going on in the time-domain. Let's take a fourier transform of our drive signal, robustly identify the drive frequency, and then filter our microsphere response around that frequency. In this way, we can verify our assumption that the response has the form $f_{\\rm resp} (t) = A_{\\rm resp} \\cos(\\omega t + \\phi)$, given a drive $f_{\\rm drive} (t) = A_{\\rm drive} \\cos(\\omega t)$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate the fourier transform of the drive signal\n", + "# Note that the frequencies we calculate are in units of Hz, whereas\n", + "# 'omega' in the definitions above are in units of rad/s, with the \n", + "# conversion factor omega = 2 * pi * freq\n", + "freqs = np.fft.rfftfreq(nsamp, 1/fsamp)\n", + "drive_fft = np.fft.rfft(drive)\n", + "\n", + "# Plot the spectrum of the drive signal\n", + "plt.figure(figsize=(6, 4))\n", + "plt.loglog(freqs, np.abs(drive_fft))\n", + "plt.xlabel('Frequency [Hz]')\n", + "plt.ylabel('Drive Spectrum [arb]')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Indeed, the Fourier transform of the drive signal looks like a single frequency together with some noise at a much lower amplitude. The single frequency looks like a delta-function at one frequency, while the rest of the spectrum is, more or less, flat across all frequencies.\n", + "\n", + "By rough estimate, it looks like our drive signal is $\\sim 40~\\text{Hz}$. \n", + "\n", + "Can you come up with a way to find this value more precisely than just reading it off the plot?\n", + "\n", + "Our estimate should be good enough to build a filter and filter the microsphere response." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Build a bandpass filter to isolate the drive frequency of ~40 Hz. We'll use a\n", + "# 4th order Butterworth filter for this, which has a pretty sharp rolloff.\n", + "sos = signal.butter(4, [30, 50], btype='band', fs=fsamp, output='sos')\n", + "\n", + "# Apply the filter to the response signal, using the signal.sosfiltfilt function\n", + "# which applies the filter twice (once forward in time, once backward in time)\n", + "# to eliminate phase distortion in the output.\n", + "response_filt = signal.sosfiltfilt(sos, response)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What does our filtered signal look like in the time-domain?\n", + "\n", + "Here, we'll compare it to the drive signal to see if we can see the sinusoid of the response." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(10, 4))\n", + "plt.plot(tvec[mask], drive[mask], label='Drive Signal')\n", + "plt.plot(tvec[mask], response_filt[mask]*4, label='Filtered Response (scaled for clarity)')\n", + "plt.xlabel('Time [s]')\n", + "plt.ylabel('Detector Response [arb]')\n", + "plt.legend(loc='upper right', framealpha=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Things to think about:\n", + "\n", + "##### 2.1 Does this look like what you expected?\n", + "\n", + "It turns out that this measurement has a lot of noise and some transient features, so even after filtering the data, noise at frequencies near to the drive signal make it through our bandpass filter and appear in superposition with our response signal.\n", + "\n", + "The presence of contaminating noise is somewhat expected: in this measurment, we're detecting the quasi-electrostatic response of a laboratory scale electric field (i.e. not that big) acting on a **single** electron, the smallest quanta of free-charge detectable, out of $\\sim 10^{13}$ total charge carriers (both positive and negative) present on the microsphere. It's actually quite remarkable we can see a response in the first place!\n", + "\n", + "##### 2.2 What does the sum of two sinusoids closely spaced in frequency look like in the time-domain? Do you see a 'beat-frequency'?\n", + "\n", + "This isn't particularly important to our calibration, but helps us to make sense of what we're seeing" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Examining the Fourier Transform of the Response\n", + "\n", + "The above plot begs the question: what does the Fourier transform of our signal actually look like?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Compute the FFT of the unfiltered response signal.\n", + "response_fft = np.fft.rfft(response)\n", + "\n", + "plt.figure(figsize=(6, 4))\n", + "plt.loglog(freqs, np.abs(response_fft))\n", + "plt.xlabel('Frequency [Hz]')\n", + "plt.ylabel('Drive Spectrum [arb]')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Clearly there's a lot more going on in the response than just the single frequency!\n", + "\n", + "There's multiple other 'tones' (single frequencies that look kind of like delta functions) in the response, as well as a lot of low-frequency noise, and some broad-spectrum features. For the calibration, we're only interested in the response at $\\sim 40~\\text{Hz}$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Suggestions moving forward\n", + "\n", + "1. Plot the filtered spectrum of the response\n", + " * Change the filter and see how this affects our response.\n", + "\n", + "2. Develop a way to programatically (i.e. without just guessing based on the plot) determine the drive frequency, and then select the Fourier component of drive and response at this determined drive frequency.\n", + " * Each element in the array `response_fft` is one such component, complex-valued, and corresponds to $\\tilde{A} = A e^{- i \\phi}$\n", + " * The same index in the array `freqs` gives you $f = \\omega / 2 \\pi$, the frequency of that paticular Fourier component\n", "\n", - "1. Come up with a way to evaluate the consistency of all of the data before subdividing into \"Early Time\" and \"Late Time\" measurements. Make sure you take the errors into account, quote a p-value (probability estimate), and explain why your results indicate the data should be subdivided. \n", + "3. Compute the amplitude ratio and phase difference using the formulas provided before\n", + " * The `np.angle()` function will be useful here.\n", "\n", - "2. Make a nice plot (or two) that summarizes the difference between \"Early Time\" results and \"Late Time\" results.\n", + "4. Make some similar plots of the frequency comb drive and response in both time and frequency domains\n", "\n", - "3. Come up with a way to evaluate the consistency of the data within each subset (\"Early Time\" and \"Late Time\"). Make sure you take the errors into account, quote a p-value (probability estimate), and interpret your result.\n", + "5. Extend your analysis to many freqencies simultaneously to analyze the entire frequency comb.\n", + " * Some clever array-masking will be useful, like we did to select times less than 0.5 seconds for plotting\n", "\n", - "4. Come up with a way to evaluate the consistency of the two different subsets (\"Early Time\" and \"Late Time\") with each other. Make sure you take the errors into account, quote a p-value (probability estimate), and interpret your result." + "6. Make plots of $A(\\omega)$ and $\\phi(\\omega)$, with appropriate axes labels" ] }, {