{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# ForcePlates IV: Calibration of Sensor Depth\n", "\n", "Falk Mielke, FunMorph, Universiteit Antwerpen\n", "\n", "You can download this notebook [here](user://notebooks/ForcePlates4_Calibration2.ipynb?download&target=_blank).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is one part of a series of blog posts on Force Plates.\n", "- [ForcePlates I: Notes on Equipment](../fpequipment)\n", "- [ForcePlates II: Calculations](../fpcalculations)\n", "- [ForcePlates III: Probabilistic Calibration of Magnitude](../fpcalibration1)\n", "- [ForcePlates IV: Calibration of Sensor Depth](../fpcalibration2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Sensor Depth Calibration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall that my original goal was to use contact sockets to improve the probability of recording single limb forces of quadruped vertebrates that locomote. This looks as follows:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now because we add _padding_, i.e. the socked that increases the height of the contact surface above the internal sensor, the situation of the force recording changes. \n", "Whereas the forces stay the same (given that the socket itself is perfectly solid and inelastic), the moments are transferred via the additional stiff element. \n", "\n", "But what is the value for padding that has to be set in the kistler software? \n", "Should it be found on calibration sheets that were long lost?\n", "Does Kistler software know, from the serial number of the force plate, what the basic padding is and transform all data to \"surface level\" values? \n", "Last but not least: our force plate had a mounting adapter (metal plate with screw holes) attached on top. Is that included, or not?\n", "\n", "So many questions, one answer: calibration. **Let's find out definitely how deep below the surface the piezo sensors are located**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I have previously demonstrated how to get the right force units out of a raw voltage recording, which required the calibration of a uniform conversion factor. This unit conversion is no prerequisite for what follows.\n", "\n", "Here, I will demonstrate how to regress the contact point = centre of pressure (CoP). In principle, the calculation is simple, but it requires one of the most relevant tools of force measurement.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is a collection of `python` libraries used below, as well as some helper functions.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#######################################################################\n", "### Libraries ###\n", "#######################################################################\n", "# data management\n", "import numpy as NP\n", "import pandas as PD\n", "import scipy.stats as STATS\n", "\n", "PD.set_option('precision', 5) # for display\n", "\n", "# symbolic maths\n", "import sympy as SYM\n", "import sympy.physics.mechanics as MECH\n", "from sympy.matrices import rot_axis3 as ROT3\n", "\n", "# plots and graphics\n", "import matplotlib as MP\n", "import matplotlib.pyplot as MPP\n", "from IPython.core.display import SVG\n", "\n", "%matplotlib inline\n", "\n", "# formula printing\n", "SYM.init_printing(use_latex=False)\n", "\n", "# vector (de)composition helpers\n", "VectorComponents = lambda vec, rf: [vec.dot(rf.x), vec.dot(rf.y), vec.dot(rf.z)]\n", "MakeVector = lambda components, coords, n = 3: sum(NP.multiply(components[:n], coords[:n]))\n", "ChangeAFrame = lambda vec, rf_old, rf_new: MakeVector(VectorComponents(vec, rf_old), [rf_new.x, rf_new.y, rf_new.z])\n", "\n", "SmallAngleApprox = lambda ang: [(SYM.sin(ang), ang), (SYM.cos(ang), 1), (SYM.tan(ang), ang)]\n", "\n", "def PrintSolution(solution_dict):\n", " # print and substitute in known values\n", " for param, eqn in solution_dict.items():\n", " print ('\\n', '_'*20)\n", " print (param)\n", " print ('_'*20)\n", " SYM.pprint(eqn.subs(fp_dimensions))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What Has Three Legs and a Rectum on the Top Side?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... a drum stool. \n", "\n", "This old joke is based on the observation that drum seats are essentially **tripods**. \n", "For our purpose, let's take away the drummer (just because our kind can't sit still). \n", "Tripods are _very_ useful, and every student seriously aiming at force recordings should have one. \n", "\n", "Tripods are part of regular photo equipment, drumming equipment, etc.\n", "But most importantly: tripods are well determined systems due to their symmetry. The reasons why tripods are so useful are that they have a predictable horizontal force component in the static case, that the leg angle and length can be adjusted (and with it the horizontal force), and that mounting some extra weight on top is easy.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "For testing, I put one leg of a tripod of a given mass on the force plate. In fact, the mass does not matter much, only that its legs are fixed and remain fixed during the recording. \n", "\n", "\n", "I did this with the Kistler BioWare software, but it is independent of DAQ/Computer assembly. \n", "Kistler software allows the export of the desired parameters, but I was not sure whether I entered everything correctly. \n", "Since I did not know the padding height (\$az_0\$ in Kistler terminology; has to be entered in the software), I thought it's worth comparing theoretical calculations, Kistler output, and my own conversion.\n", "\n", "Fun fact: one can order support by a Kistler technician, who comes to help you with the software settings at a daily rate of approximately 1800.-EUR. Must be a clever guy. But enough Kistler bashing for now. Always know that you can have me for 999,-EUR (plus travel expenses).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below an overview of the test setup." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This was one of my first installations - today, I consider it more a kind of art than relevant measurement. \n", "I had carefully measured tripod angles, weighted the tripod, had meticulously piled books and papers to level out the feet, fixed a wooden socket on the force plate, and even made sure that all legs of the tripod get the same share of its weight. \n", "\n", "So what do you guess is irrelevant in this setting (excluding the measurement tools)?\n", "\n", "\$\$\\ldots\$\$\n", "\n", "\n", "Correct! All except the tripod. I learned a lot with these measurements.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What can benefit our purpose though is a point tip on one of the legs. I tend to tape screws to it - simple and efficient. The pointier, the better, and it should be tightly fixed, but careful not to hurt yourself or anyone." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The Tripod Trick" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You [are already able to](../fpcalculations) calculate moments and contact points from the raw Kistler output.\n", "Put the tripod with one pointy leg on the force plate. \n", "Imagine you knew precisely the correct \$z\$-offset = padding = \$az_0\$. \n", "The calculated contact point should be independent of the direction of the tripod. \n", "\n", "So put a tripod on an arbitrary point on the forceplate, take a calibration measurement (load, begin recording, unload after a time, end recording). Then turn the tripod, place it on exactly the same point, measure again." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The challenge is to hit the point repeatedly. But that can be overcome by statistics, because the placement error should be fairly normally distributed around the target. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculation Code - Rehearsal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we need our `sympy` code to calculate parameters from channels. However, this time we leave one force plate dimension unknown. I will just use formulas as before - refer to the [previous post](../fpcalculations) for details. \n", "\n", "You can also [skip this section](#skip_calculation)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "### reference frame of the force plate\n", "world = MECH.ReferenceFrame('N')\n", "origin = MECH.Point('O')\n", "origin.set_vel(world, 0)\n", "\n", "phi = SYM.symbols('phi_{x:z}') # force plate rotation in the world\n", "forceplate = world.orientnew('fp', 'Body', phi[::-1], 'ZYX')\n", "\n", "# center point = force plate center of mass\n", "forceplate_offset = SYM.symbols('c_{x:z}')\n", "center = origin.locatenew('C', MakeVector(forceplate_offset, [world.x, world.y, world.z]))\n", "\n", "\n", "# force plate coordinate system\n", "x = forceplate.x\n", "y = forceplate.y\n", "z = forceplate.z\n", "coordinates = [x, y, z]\n", "coord_strings = ['x', 'y', 'z']\n", "\n", "# notation: \n", "# i ... leg index\n", "# j ... coordinate index\n", "\n", "# impact point\n", "pj = SYM.symbols('p_{x:z}')\n", "impact = center.locatenew('P', MakeVector(pj, coordinates))\n", "\n", "# impact position vector\n", "rc = impact.pos_from(center)\n", "\n", "# leg points, relative to force plate center\n", "n_legs = 4\n", "leg_ids = range(n_legs)\n", "leg_positions = NP.array([(k,l) for k, l in [[+1,+1], [-1,+1], [-1,-1], [+1,-1]] ])\n", "\n", "# leg vectors\n", "qj = SYM.symbols('q_{x:y}') # plate dimensions\n", "rq = [leg_positions[i,0]*qj*x + leg_positions[i,1]*qj*y + 0*z for i in leg_ids] \n", "\n", "\n", "# leg points\n", "legs = [center.locatenew('Q_{%i}' % (i), rq[i]) for i in leg_ids]\n", "\n", "# vector from leg to impact\n", "sq = [leg.pos_from(impact) for leg in legs] # leg - impact" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "fp_dimensions = [(qj, 0.035), (qj, 0.075)]\n", "# padding (pj) has to be estimated from the data" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "### Balance of Forces\n", "\n", "fij = NP.array(SYM.symbols('f_{:4x:z}')).reshape(n_legs, -1)\n", "\n", "Fi = [MakeVector(fij[i,:], coordinates) for i in leg_ids]\n", "Fj = SYM.symbols('F_{x:z}')\n", "\n", "impact_components = Fj\n", "reactn_components = VectorComponents(sum(Fi), forceplate)\n", "\n", "# resulting equations\n", "force_balances = [SYM.Eq(impact_components[coord], reactn_components[coord]) for coord, _ in enumerate(coordinates) ]\n", "\n", "force_subs = [(imp, rcn) for imp, rcn in zip(impact_components, reactn_components)]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "### Balance of Moments\n", "# free moment\n", "Tz = SYM.symbols('T_{z}')\n", "\n", "\n", "# moments of the impact force on the whole object\n", "impact_moments = rc.cross(MakeVector(Fj, coordinates))\n", "impact_moments += Tz*z\n", "\n", "\n", "# moments of the reaction forces on the legs on the whole object (i.e. center)\n", "reactn_moments = [rq[i].cross(MakeVector(fij[i, :], coordinates )) \\\n", " for i in leg_ids]\n", "\n", "reactn_moments = [vc.factor(pj+qj) for vc in VectorComponents(sum(reactn_moments), forceplate)]\n", "\n", "# resulting equations\n", "moment_equations = [SYM.Eq(imp, rcn).subs(force_subs).simplify() \\\n", " for imp, rcn in zip(VectorComponents(impact_moments, forceplate), reactn_moments)]" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "### solution\n", "fx01, fx23, fy03, fy12 = SYM.symbols('fx_{01}, fx_{23}, fy_{03}, fy_{12}')\n", "fz = SYM.symbols('fz_{:4}')\n", "\n", "force_components = [fx01, fx23, fy03, fy12] + [f for f in fz]\n", "\n", "group_subs = [ (fij[0,0]+fij[1,0], fx01) \\\n", " , (fij[2,0]+fij[3,0], fx23) \\\n", " , (fij[0,1]+fij[3,1], fy03) \\\n", " , (fij[1,1]+fij[2,1], fy12) \\\n", " ] + [ \\\n", " (fij[i,2], fz[i]) for i in leg_ids \\\n", " ]\n", "\n", "moment_equations = [mmnt.factor(qj+pj).subs(group_subs) for mmnt in moment_equations]\n", "\n", "# (1) contact point\n", "solvents = [pj, pj]\n", "\n", "p_solutions = [ \\\n", " {solvents[s]: sln.factor(pj+qj) \\\n", " for s, sln in enumerate(solution)} \\\n", " for solution in SYM.nonlinsolve(moment_equations, solvents) \\\n", " ]\n", "\n", "tz_equation = moment_equations.subs([(pnt, sol) for pnt, sol in p_solutions.items()]).factor(qj+pj)\n", "tz_solution = [sol for sol in SYM.solveset(tz_equation, Tz)]\n", "\n", "cop_solutions = {Tz: tz_solution.simplify(), **p_solutions}\n", "\n", "# (2) force resubstitution\n", "force_solutions = {imp: rcn.subs(group_subs) for imp, rcn in force_subs}\n", "\n", "#(3) moment resubstitution\n", "Mj = SYM.symbols('M_{x:z}')\n", "moment_solutions = VectorComponents(impact_moments.subs([(param, sol) for param, sol in cop_solutions.items()]), forceplate)\n", "\n", "moment_solutions = {Mj[j]: moment_solutions[j].subs( \\\n", " [(force, sol) for force, sol in force_solutions.items()] \\\n", " ).simplify() \\\n", " for j, _ in enumerate(coordinates)}\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " ____________________\n", "F_{x}\n", "____________________\n", "fx_{01} + fx_{23}\n", "\n", " ____________________\n", "F_{y}\n", "____________________\n", "fy_{03} + fy_{12}\n", "\n", " ____________________\n", "F_{z}\n", "____________________\n", "fz_{0} + fz_{1} + fz_{2} + fz_{3}\n", "\n", " ____________________\n", "M_{x}\n", "____________________\n", "0.075⋅fz_{0} + 0.075⋅fz_{1} - 0.075⋅fz_{2} - 0.075⋅fz_{3}\n", "\n", " ____________________\n", "M_{y}\n", "____________________\n", "-0.035⋅fz_{0} + 0.035⋅fz_{1} + 0.035⋅fz_{2} - 0.035⋅fz_{3}\n", "\n", " ____________________\n", "M_{z}\n", "____________________\n", "-0.075⋅fx_{01} + 0.075⋅fx_{23} + 0.035⋅fy_{03} - 0.035⋅fy_{12}\n", "\n", " ____________________\n", "T_{z}\n", "____________________\n", "2⋅(-0.075⋅fx_{01}⋅fz_{2} - 0.075⋅fx_{01}⋅fz_{3} + 0.075⋅fx_{23}⋅fz_{0} + 0.075\n", "──────────────────────────────────────────────────────────────────────────────\n", " fz\n", "\n", "⋅fx_{23}⋅fz_{1} + 0.035⋅fy_{03}⋅fz_{1} + 0.035⋅fy_{03}⋅fz_{2} - 0.035⋅fy_{12}⋅\n", "──────────────────────────────────────────────────────────────────────────────\n", "_{0} + fz_{1} + fz_{2} + fz_{3} \n", "\n", "fz_{0} - 0.035⋅fy_{12}⋅fz_{3})\n", "──────────────────────────────\n", " \n", "\n", " ____________________\n", "p_{x}\n", "____________________\n", "0.035⋅fz_{0} - 0.035⋅fz_{1} - 0.035⋅fz_{2} + 0.035⋅fz_{3} + p_{z}⋅(fx_{01} + f\n", "──────────────────────────────────────────────────────────────────────────────\n", " fz_{0} + fz_{1} + fz_{2} + fz_{3} \n", "\n", "x_{23})\n", "───────\n", " \n", "\n", " ____________________\n", "p_{y}\n", "____________________\n", "0.075⋅fz_{0} + 0.075⋅fz_{1} - 0.075⋅fz_{2} - 0.075⋅fz_{3} + p_{z}⋅(fy_{03} + f\n", "──────────────────────────────────────────────────────────────────────────────\n", " fz_{0} + fz_{1} + fz_{2} + fz_{3} \n", "\n", "y_{12})\n", "───────\n", " \n" ] } ], "source": [ "### combined solution\n", "\n", "all_solutions = {**force_solutions, **moment_solutions, **cop_solutions}\n", "\n", "PrintSolution(all_solutions)\n", "\n", "parameters = list(all_solutions.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, here come the conversion functions, but this time including the padding as a free parameter." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# generate conversion functions\n", "ComponentsToParameters = { str(param): \\\n", " SYM.lambdify(force_components + [pj], eqn.subs(fp_dimensions), \"numpy\") \\\n", " for param, eqn in all_solutions.items()\\\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Empirical Regression " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's first [load test data](user://notebooks/cop_data.csv?target=_blank&download) that I acquired when calibrating a Kistler force plate.\n", "\n", "I took the four corner points, measured five repetitions for three directions." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
daterec_nrpt_nrgainFx12Fx34Fy14Fy23Fz1Fz2Fz3Fz4target_pxtarget_py
02019060412111001.15986-0.57157-1.56003-1.261003.205570.56842-0.461100.795750.0350.075
12019060412211000.79159-0.55819-1.04109-1.232263.197960.62729-0.485340.780920.0350.075
22019060412311000.89389-0.56254-0.83282-1.261613.254820.61714-0.479650.770700.0350.075
32019060412411001.06494-0.51332-1.43133-1.153573.170080.59010-0.477800.811280.0350.075
42019060412511001.04886-0.52417-1.34062-1.158753.184970.59545-0.484770.805680.0350.075
\n", "