From 783aae3cc12b4da039434dcf6229976611661e97 Mon Sep 17 00:00:00 2001 From: jaseg Date: Mon, 15 Mar 2021 13:44:59 +0100 Subject: fix up interval-based frequency estimator --- .../Accelerometer Data Analysis.ipynb | 302 ++++++++++++--------- 1 file changed, 167 insertions(+), 135 deletions(-) (limited to 'prototype/sensor-analysis') diff --git a/prototype/sensor-analysis/Accelerometer Data Analysis.ipynb b/prototype/sensor-analysis/Accelerometer Data Analysis.ipynb index 75f0751..b7cf70d 100644 --- a/prototype/sensor-analysis/Accelerometer Data Analysis.ipynb +++ b/prototype/sensor-analysis/Accelerometer Data Analysis.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 703, "metadata": {}, "outputs": [], "source": [ @@ -24,7 +24,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 704, "metadata": {}, "outputs": [], "source": [ @@ -33,19 +33,30 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 705, + "metadata": {}, + "outputs": [], + "source": [ + "#last_run, = db.execute('SELECT MAX(run_id) FROM packets').fetchone()\n", + "# Override\n", + "#last_run = 40\n", + "last_run=33" + ] + }, + { + "cell_type": "code", + "execution_count": 706, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Last run was ID #42 with 1745 packets total, 256 distinct over 217.80828976631165s\n" + "Last run was ID #33 with 3178 packets total, 125 distinct over 99.94346904754639s\n" ] } ], "source": [ - "last_run, = db.execute('SELECT MAX(run_id) FROM packets').fetchone()\n", "num_packets, = db.execute('SELECT COUNT(*) FROM packets WHERE run_id=?', (last_run,)).fetchone()\n", "num_packets_distinct, = db.execute('SELECT COUNT(*) FROM (SELECT DISTINCT data FROM packets WHERE run_id=?)', (last_run,)).fetchone()\n", "timespan_start, timespan_end = db.execute('SELECT MIN(timestamp_us)/1e6, MAX(timestamp_us)/1e6 FROM packets WHERE run_id=?', (last_run,)).fetchone()\n", @@ -55,17 +66,7 @@ }, { "cell_type": "code", - "execution_count": 39, - "metadata": {}, - "outputs": [], - "source": [ - "# Override\n", - "last_run = 40" - ] - }, - { - "cell_type": "code", - "execution_count": 40, + "execution_count": 707, "metadata": {}, "outputs": [], "source": [ @@ -76,7 +77,7 @@ }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 708, "metadata": {}, "outputs": [ { @@ -1040,7 +1041,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -1052,10 +1053,10 @@ { "data": { "text/plain": [ - "[]" + "[]" ] }, - "execution_count": 41, + "execution_count": 708, "metadata": {}, "output_type": "execute_result" } @@ -1067,8 +1068,10 @@ }, { "cell_type": "code", - "execution_count": 42, - "metadata": {}, + "execution_count": 709, + "metadata": { + "scrolled": true + }, "outputs": [ { "name": "stdout", @@ -1087,7 +1090,7 @@ }, { "cell_type": "code", - "execution_count": 43, + "execution_count": 710, "metadata": {}, "outputs": [ { @@ -1106,7 +1109,7 @@ }, { "cell_type": "code", - "execution_count": 44, + "execution_count": 711, "metadata": {}, "outputs": [], "source": [ @@ -1119,7 +1122,7 @@ }, { "cell_type": "code", - "execution_count": 45, + "execution_count": 712, "metadata": { "scrolled": false }, @@ -1140,14 +1143,14 @@ }, { "cell_type": "code", - "execution_count": 46, + "execution_count": 713, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Sequence number range: 739 ... 1457\n" + "Sequence number range: 35 ... 161\n" ] } ], @@ -1158,7 +1161,7 @@ }, { "cell_type": "code", - "execution_count": 47, + "execution_count": 714, "metadata": {}, "outputs": [], "source": [ @@ -1170,7 +1173,7 @@ }, { "cell_type": "code", - "execution_count": 97, + "execution_count": 715, "metadata": { "scrolled": false }, @@ -2136,7 +2139,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -2149,63 +2152,73 @@ "name": "stdout", "output_type": "stream", "text": [ - "Found sensor offset: 0.57 g / 5.55 m/s^2\n", + "Found sensor offset: 1.00 g / 9.81 m/s^2\n", "\n", "Centrifugal acceleration at 6.49 Hz:\n", " Theory: 2.16 g / 21.14 m/s^2\n", - " Measurement: 1.90 g / 18.586194313627637 m/s^2\n", - " Rel. Error: 13.72 %\n", - " Abs. Error: 0.26 g / 2.55 m/s^2\n", + " Measurement: 55.57 g / 544.9912192549019 m/s^2\n", + " Rel. Error: -96.12 %\n", + " Abs. Error: -53.42 g / -523.85 m/s^2\n", "\n", "Centrifugal acceleration at 6.49 Hz:\n", " Theory: 6.82 g / 66.88 m/s^2\n", - " Measurement: 6.60 g / 64.67817413725541 m/s^2\n", - " Rel. Error: 3.41 %\n", - " Abs. Error: 0.22 g / 2.20 m/s^2\n", + " Measurement: -0.58 g / -5.674407202881151 m/s^2\n", + " Rel. Error: -1278.66 %\n", + " Abs. Error: 7.40 g / 72.56 m/s^2\n", "\n", "Centrifugal acceleration at 6.49 Hz:\n", " Theory: 14.89 g / 146.00 m/s^2\n", - " Measurement: 14.82 g / 145.35641444809548 m/s^2\n", - " Rel. Error: 0.44 %\n", - " Abs. Error: 0.07 g / 0.64 m/s^2\n", + " Measurement: nan g / nan m/s^2\n", + " Rel. Error: nan %\n", + " Abs. Error: nan g / nan m/s^2\n", "\n", "Centrifugal acceleration at 6.49 Hz:\n", " Theory: 23.04 g / 225.90 m/s^2\n", - " Measurement: 23.11 g / 226.6766272456559 m/s^2\n", - " Rel. Error: -0.34 %\n", - " Abs. Error: -0.08 g / -0.77 m/s^2\n", + " Measurement: nan g / nan m/s^2\n", + " Rel. Error: nan %\n", + " Abs. Error: nan g / nan m/s^2\n", "\n", "Centrifugal acceleration at 6.49 Hz:\n", " Theory: 34.60 g / 339.27 m/s^2\n", - " Measurement: 34.80 g / 341.26606334638393 m/s^2\n", - " Rel. Error: -0.59 %\n", - " Abs. Error: -0.20 g / -2.00 m/s^2\n", + " Measurement: nan g / nan m/s^2\n", + " Rel. Error: nan %\n", + " Abs. Error: nan g / nan m/s^2\n", "\n", "Centrifugal acceleration at 6.49 Hz:\n", " Theory: 53.88 g / 528.41 m/s^2\n", - " Measurement: 51.75 g / 507.51233891199473 m/s^2\n", - " Rel. Error: 4.12 %\n", - " Abs. Error: 2.13 g / 20.90 m/s^2\n", + " Measurement: nan g / nan m/s^2\n", + " Rel. Error: nan %\n", + " Abs. Error: nan g / nan m/s^2\n", "\n", "Centrifugal acceleration at 6.49 Hz:\n", " Theory: 81.62 g / 800.43 m/s^2\n", - " Measurement: 82.09 g / 805.0345119987545 m/s^2\n", - " Rel. Error: -0.57 %\n", - " Abs. Error: -0.47 g / -4.60 m/s^2\n", + " Measurement: nan g / nan m/s^2\n", + " Rel. Error: nan %\n", + " Abs. Error: nan g / nan m/s^2\n", "\n", "Centrifugal acceleration at 6.49 Hz:\n", " Theory: 29.79 g / 292.17 m/s^2\n", - " Measurement: 29.71 g / 291.36536349393447 m/s^2\n", - " Rel. Error: 0.28 %\n", - " Abs. Error: 0.08 g / 0.81 m/s^2\n", + " Measurement: nan g / nan m/s^2\n", + " Rel. Error: nan %\n", + " Abs. Error: nan g / nan m/s^2\n", "\n", "Centrifugal acceleration at 6.49 Hz:\n", " Theory: 9.33 g / 91.46 m/s^2\n", - " Measurement: 9.21 g / 90.28568973827844 m/s^2\n", - " Rel. Error: 1.30 %\n", - " Abs. Error: 0.12 g / 1.17 m/s^2\n", + " Measurement: nan g / nan m/s^2\n", + " Rel. Error: nan %\n", + " Abs. Error: nan g / nan m/s^2\n", "\n" ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + ":46: RuntimeWarning: Mean of empty slice.\n", + " ivl_avg = (reassembled_values / mems_lsb_per_g)[idx].mean()\n", + "/usr/lib/python3.9/site-packages/numpy/core/_methods.py:188: RuntimeWarning: invalid value encountered in double_scalars\n", + " ret = ret.dtype.type(ret / rcount)\n" + ] } ], "source": [ @@ -2277,7 +2290,7 @@ "print(f'Found sensor offset: {sensor_offx:.2f} g / {sensor_offx*g:.2f} m/s^2')\n", "print()\n", "\n", - "for theory, meas in zip(acc_theory, acc_meas):\n", + "for theory, meas, interval in zip(acc_theory, acc_meas, interval_speeds):\n", " ax.axhline(theory - sensor_offx, color='orange', alpha=1, zorder=1)\n", " meas += sensor_offx\n", " \n", @@ -2291,34 +2304,8 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 716, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Average speed of rotation: 19.12 Hz / 1147 rpm\n" - ] - } - ], - "source": [ - "speed_ivl_min, speed_ivl_max = ivl_start, ivl_end\n", - "\n", - "def fun(x, args):\n", - " deltas, = args # poor api\n", - " return np.sqrt(np.mean([ ((val + 0.5*x[0]) % x[0] - 0.5*x[0])**2 for val in deltas ]))\n", - "res = scipy.optimize.minimize(fun, 0.05, args=[sorted(deltas[speed_ivl_min:speed_ivl_max])[:-2]])\n", - "interval = np.abs(res.x[0])\n", - "print(f'Average speed of rotation: {1/interval:.2f} Hz / {60 / interval:.0f} rpm')" - ] - }, - { - "cell_type": "code", - "execution_count": 15, - "metadata": { - "scrolled": false - }, "outputs": [ { "data": { @@ -3281,7 +3268,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -3290,41 +3277,84 @@ "metadata": {}, "output_type": "display_data" }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Largest peak at 4.31 Hz / 259 rpm\n", + "Mixing product 1 at 5.69 Hz / 341 rpm\n", + "Mixing product 2 at 14.31 Hz / 859 rpm\n", + "Mixing product 3 at 15.69 Hz / 941 rpm\n", + "Mixing product 4 at 24.31 Hz / 1459 rpm\n", + "Mixing product 5 at 25.69 Hz / 1541 rpm\n", + "Mixing product 6 at 34.31 Hz / 2059 rpm\n" + ] + }, { "data": { "text/plain": [ - "[]" + "" ] }, - "execution_count": 15, + "execution_count": 716, "metadata": {}, "output_type": "execute_result" } ], "source": [ + "tsa = np.array(timestamps)\n", + "#s_min, s_max = 70, 120\n", + "s_min, s_max = 40, 90\n", + "speed_idx = (tsa > s_min) & (tsa < s_max)\n", + "\n", + "ts = np.arange(0, len(reassembled_values)) / sampling_rate\n", + "fft_idx = (ts > s_min) & (ts < s_max)\n", + "\n", + "N = fft_idx.sum()\n", + "T = 1/sampling_rate\n", + "x = np.linspace(0.0, N*T, N)\n", + "y = reassembled_values[fft_idx] / mems_lsb_per_g # cut out beginning and that time we tapped the thing\n", + "y *= scipy.signal.windows.blackmanharris(len(y))\n", + "yf = scipy.fftpack.fft(y)\n", + "xf = np.linspace(0.0, 1/(2*T), N//2)\n", + "mag = 2/N * np.abs(yf[:N//2])\n", + "\n", "fig, ax = plt.subplots()\n", "ax.grid()\n", - "for i in range(int(max(deltas[speed_ivl_min:speed_ivl_max])//interval)):\n", - " ax.axhline(i*interval, color='orange')\n", - "ax.plot(sorted(deltas[speed_ivl_min:speed_ivl_max])[:-2])" + "ax.plot(xf, mag)\n", + "ax.set_ylabel('Magnitude [g]')\n", + "ax.set_xlabel('Frequency [Hz]')\n", + "\n", + "peaks, _ = scipy.signal.find_peaks(mag, height=.1, distance=1/T)\n", + "assert peaks.any()\n", + "\n", + "peak_data = sorted([ (-mag[idx], xf[idx]) for idx in peaks ])\n", + "largest_peak_f = peak_data[0][1]\n", + "print(f'Largest peak at {largest_peak_f:.2f} Hz / {largest_peak_f * 60:.0f} rpm')\n", + "for i in range(1,4):\n", + " mix1 = i*sampling_rate - largest_peak_f\n", + " mix2 = i*sampling_rate + largest_peak_f\n", + " print(f'Mixing product {2*i-1} at {mix1:.2f} Hz / {mix1 * 60:.0f} rpm')\n", + " print(f'Mixing product {2*i} at {mix2:.2f} Hz / {mix2 * 60:.0f} rpm')\n", + "\n", + "ax.axvline(largest_peak_f, color='orange')" ] }, { "cell_type": "code", - "execution_count": 77, - "metadata": {}, + "execution_count": 717, + "metadata": { + "scrolled": false + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Largest peak at 1.68 Hz / 101 rpm\n", - "Mixing product 1 at 8.32 Hz / 499 rpm\n", - "Mixing product 2 at 11.68 Hz / 701 rpm\n", - "Mixing product 3 at 18.32 Hz / 1099 rpm\n", - "Mixing product 4 at 21.68 Hz / 1301 rpm\n", - "Mixing product 5 at 28.32 Hz / 1699 rpm\n", - "Mixing product 6 at 31.68 Hz / 1901 rpm\n" + "interval: -0.06357678079289661\n", + "scores: [0.0005478345240435783, 0.003917745560629984, 0.00382343194410793, 0.006241784699377939, 0.0074863442856081385, 0.008964752044997151, 0.008813512329706083, 0.008780553883453538, 0.008530974440876703]\n", + "argmin: 1\n", + "Average speed of rotation: 15.73 Hz / 944 rpm\n" ] }, { @@ -4288,7 +4318,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -4300,47 +4330,49 @@ { "data": { "text/plain": [ - "Text(0.5, 0, 'Frequency [Hz]')" + "[]" ] }, - "execution_count": 77, + "execution_count": 717, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "N = ivl_end - ivl_start\n", - "T = 1/sampling_rate\n", - "x = np.linspace(0.0, N*T, N)\n", - "y = reassembled_values[ivl_start:ivl_end] / mems_lsb_per_g # cut out beginning and that time we tapped the thing\n", - "y *= scipy.signal.windows.blackmanharris(len(y))\n", - "yf = scipy.fftpack.fft(y)\n", - "xf = np.linspace(0.0, 1/(2*T), N//2)\n", - "mag = 2/N * np.abs(yf[:N//2])\n", - "\n", - "peaks, _ = scipy.signal.find_peaks(mag, height=.1, distance=1/T)\n", - "assert peaks.any()\n", + "def calc_rspeed_deltas(target_deltas):\n", + " target_deltas = np.array(target_deltas)\n", + " target_deltas = target_deltas[0:-int(len(target_deltas)*0.1)]\n", + " target_deltas = target_deltas[target_deltas > 1/30]\n", + " def fun(x):\n", + " return np.sqrt(np.mean([ ((val + 0.5*x[0]) % x[0] - 0.5*x[0])**2 for val in target_deltas ]))\n", "\n", - "peak_data = sorted([ (-mag[idx], xf[idx]) for idx in peaks ])\n", - "largest_peak_f = peak_data[0][1]\n", - "print(f'Largest peak at {largest_peak_f:.2f} Hz / {largest_peak_f * 60:.0f} rpm')\n", - "for i in range(1,4):\n", - " mix1 = i*sampling_rate - largest_peak_f\n", - " mix2 = i*sampling_rate + largest_peak_f\n", - " print(f'Mixing product {2*i-1} at {mix1:.2f} Hz / {mix1 * 60:.0f} rpm')\n", - " print(f'Mixing product {2*i} at {mix2:.2f} Hz / {mix2 * 60:.0f} rpm')\n", + " #def accept(x_old, x_new, **kwargs):\n", + " # return 1/30 < x_new[0] and x_new[0] < 1\n", + " res = scipy.optimize.minimize(fun, 0.1)\n", + " #res = scipy.optimize.basinhopping(fun, max(1/largest_peak_f, 0.5), accept_test=accept, stepsize=0.2)\n", + " \n", + " interval = res.x[0]\n", + " scores = [ fun([i*interval]) * (1/(i+10)) for i in range(1, 10) ]\n", + " print('interval:', interval)\n", + " print('scores:', scores)\n", + " argmin = np.argmin(scores)+1\n", + " print('argmin:', argmin)\n", + " interval = np.abs(interval) * argmin\n", + " print(f'Average speed of rotation: {1/interval:.2f} Hz / {60 / interval:.0f} rpm')\n", + " return interval\n", "\n", + "target_deltas = sorted(np.array(deltas)[speed_idx[:-1]])\n", + "interval = calc_rspeed_deltas(target_deltas)\n", "fig, ax = plt.subplots()\n", "ax.grid()\n", - "ax.axvline(xf[peaks], color='orange')\n", - "ax.plot(xf, mag)\n", - "ax.set_ylabel('Magnitude [g]')\n", - "ax.set_xlabel('Frequency [Hz]')" + "for i in range(int(max(target_deltas)//interval)):\n", + " ax.axhline(i*interval, color='orange')\n", + "ax.plot(target_deltas)" ] }, { "cell_type": "code", - "execution_count": 78, + "execution_count": 718, "metadata": { "scrolled": false }, @@ -5306,7 +5338,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -5324,15 +5356,15 @@ }, { "cell_type": "code", - "execution_count": 80, + "execution_count": 719, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Centrifugal acceleration at 1.68 Hz: 6.11 m/s^2 / 0.62 g\n", - "Centrifugal acceleration at 8.35 Hz: 151.53 m/s^2 / 15.45 g\n" + "Centrifugal acceleration at 4.31 Hz: 40.42 m/s^2 / 4.12 g\n", + "Centrifugal acceleration at 15.73 Hz: 537.19 m/s^2 / 54.78 g\n" ] } ], @@ -5352,7 +5384,7 @@ }, { "cell_type": "code", - "execution_count": 81, + "execution_count": 720, "metadata": { "scrolled": false }, @@ -5402,7 +5434,7 @@ }, { "cell_type": "code", - "execution_count": 83, + "execution_count": 721, "metadata": {}, "outputs": [ { @@ -6366,7 +6398,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" -- cgit