revision of signal analysis notebook

notebooks/2_tgv_signal.ipynb

@@ -11,7 +11,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -27,16 +27,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Rocking motion of a TGV"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Our data are stored in a compressed file `TGV_data.csv.bz2`. Each row of the uncompressed file contains entries separated by commas and the first row contains labels explaining the content of the respective column.\n",
-    "\n",
-    "NumPy provides rather universal tools to import data from files into NumPy arrays. We will use `genfromtxt` which allows to deal with `bz2` compressed data files and also handles the labels in the first row of the data."
+    "### Importing the data"
    ]
   },
   {
@@ -45,23 +36,24 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "np.info(np.genfromtxt)"
+    "data = np.genfromtxt('data/TGV_data.csv.bz2', delimiter=',', names=True)\n",
+    "time = data['Time_s']\n",
+    "omega_y = data['Gyroscope_y_rads']\n",
+    "_ = plt.plot(time, omega_y)"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "data = np.genfromtxt('data/TGV_data.csv.bz2', delimiter=',', names=True)"
+    "### Time spacing"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "There are five columns identified by names:"
+    "Before analyzing the data any further, let us take a look at the time elapsed between subsequent measurements."
    ]
   },
   {
@@ -70,7 +62,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "data.dtype.names"
+    "time_intervals = time[1:]-time[:-1]"
    ]
   },
   {
@@ -79,16 +71,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "time = data['Time_s']\n",
-    "omega_x = data['Gyroscope_x_rads']\n",
-    "omega_y = data['Gyroscope_y_rads']"
+    "_ = plt.hist(time_intervals, bins=100)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Let us first get an idea of the data."
+    "We will ignore these small differences and assume the data to be equally spaced in time with the mean time difference between subsequent data points."
    ]
   },
   {
@@ -97,22 +87,25 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(time, omega_x)\n",
-    "plt.plot(time, omega_y)"
+    "delta_t = np.mean(time_intervals)\n",
+    "delta_t"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Spectral properties"
+    "### Spectral analysis"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Before analyzing the data any further, let us take a look at the time elapsed between subsequent measurements."
+    "Fourier transformation of a discrete real signal with `rfft`:\n",
+    "$$y_j = \\sum_{k=0}^{n-1} x_k \\exp\\left(-\\text{i}\\frac{2\\pi jk}{n}\\right)$$\n",
+    "With exception of $y_0$ and, for even $n$, $y_{n/2}$ all $y_j$ are complex.\n",
+    "`rfft` returns an array $y_0, \\text{Re}(y_1), \\text{Im}(y_1),\\ldots$"
    ]
   },
   {
@@ -121,7 +114,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "time_intervals = time[1:]-time[:-1]"
+    "fft_y = fftpack.rfft(omega_y)"
    ]
   },
   {
@@ -130,14 +123,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "_ = plt.hist(time_intervals, bins=100)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We will ignore these small differences and assume the data to be equally spaced in time with the mean time difference between subsequent data points."
+    "_ = plt.plot(fft_y[::2])"
    ]
   },
   {
@@ -146,15 +132,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "delta_t = np.mean(time_intervals)\n",
-    "delta_t"
+    "_ = plt.plot(fft_y[1::2])"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
+   "outputs": [],
    "source": [
-    "We will do the spectral analysis with the data for $\\omega_y$."
+    "freqs = fftpack.rfftfreq(omega_y.shape[0], delta_t)"
    ]
   },
   {
@@ -163,7 +150,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(*signal.periodogram(omega_y, 1/delta_t))"
+    "_ = plt.plot(freqs[::2], fft_y[::2])"
    ]
   },
   {
@@ -172,25 +159,23 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "z = fftpack.rfft(omega_y)"
+    "_ = plt.plot(freqs[1::2], fft_y[1::2])"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "plt.plot(z[::2])"
+    "<div style=\"color: DarkBlue;background-color: LightYellow\">\n",
+    "<i>Exercise:</i> Given the Fourier transform, how does one get back to the original data?\n",
+    "</div>"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "plt.plot(z[1::2])"
+    "### Power spectral density"
    ]
   },
   {
@@ -199,7 +184,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "freqs = fftpack.rfftfreq(omega_y.shape[0], delta_t)"
+    "fft_y = np.insert(fft_y, 0, 0)\n",
+    "fft_square = fft_y*fft_y\n",
+    "power = fft_square[::2]+fft_square[1::2]\n",
+    "power = 2*delta_t/fft_square.shape[0]*power\n",
+    "_ = plt.plot(freqs[::2], power)"
    ]
   },
   {
@@ -208,23 +197,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(freqs[::2], z[::2])"
+    "_ = plt.plot(*signal.periodogram(omega_y, 1/delta_t))"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "plt.plot(freqs[1::2], z[1::2])"
+    "### Filter out the rocking signal"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We are interested in the signal around 1.4 Hz. Filter out frequencies beyond 2 Hz."
+    "We are interested in the signal around 1.4 Hz. Filter out frequencies beyond 3 Hz."
    ]
   },
   {
@@ -242,7 +229,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "freqs, response = signal.freqz(filter_coeffs, fs=1/delta_t)"
+    "_ = plt.plot(filter_coeffs, '.')"
    ]
   },
   {
@@ -251,7 +238,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(freqs, 20*np.log10(np.abs(response)))"
+    "freqs, response = signal.freqz(filter_coeffs, fs=1/delta_t)"
    ]
   },
   {
@@ -260,7 +247,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.polar(np.angle(response), np.abs(response))"
+    "_ = plt.plot(freqs, 20*np.log10(np.abs(response)))"
    ]
   },
   {
@@ -269,7 +256,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(time, omega_y)"
+    "_ = plt.polar(np.angle(response), np.abs(response))"
    ]
   },
   {
@@ -296,7 +283,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(time[150:-150], omega_y[150:-150])"
+    "_ = plt.plot(time[150:-150], omega_y[150:-150])"
    ]
   },
   {
@@ -305,7 +292,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(time[150:-150], omega_y_filtered)"
+    "_ = plt.plot(time[150:-150], omega_y_filtered)"
    ]
   },
   {
@@ -339,7 +326,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(freqs, 20*np.log10(abs(response)))"
+    "_ = plt.plot(freqs, 20*np.log10(abs(response)))"
    ]
   },
   {
@@ -357,8 +344,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(time[150:-150], omega_y[150:-150])\n",
-    "plt.plot(time[150:-150], omega_y_filtered)"
+    "_ = plt.plot(time[150:-150], omega_y[150:-150])\n",
+    "_ = plt.plot(time[150:-150], omega_y_filtered)"
    ]
   },
   {
@@ -367,8 +354,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(fftpack.rfftfreq(omega_y_filtered.shape[0], delta_t)[::2],\n",
-    "         fftpack.rfft(omega_y_filtered)[::2])"
+    "_ = plt.plot(fftpack.rfftfreq(omega_y_filtered.shape[0], delta_t)[::2],\n",
+    "             fftpack.rfft(omega_y_filtered)[::2])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Spectrogram"
    ]
   },
   {
@@ -377,8 +371,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "freq, sp_time, Sxx = signal.spectrogram(omega_y, fs=1/delta_t, nperseg=512)\n",
-    "plt.pcolormesh(freq, sp_time, Sxx)"
+    "freq, sp_time, Sxx = signal.spectrogram(omega_y, fs=1/delta_t)\n",
+    "_ = plt.pcolormesh(sp_time, freq, Sxx)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<div style=\"color: DarkBlue;background-color: LightYellow\">\n",
+    "<i>Exercise:</i> Try to obtain more structure by plotting the logarithm of the spectrogram.\n",
+    "</div>"
    ]
   },
   {
@@ -389,7 +392,7 @@
    "source": [
     "fig = plt.figure()\n",
     "ax = Axes3D(fig)\n",
-    "ax.plot_surface(freq[:, np.newaxis], time, Sxx, cmap=cm.viridis)"
+    "_ = ax.plot_surface(sp_time, freq[:, np.newaxis], Sxx, cmap=cm.viridis)"
    ]
   }
  ],