dtk Package

bicycle Module

dtk.bicycle.basu_sig_figs()

Returns the number of significant figures reported in Table 1 of [BasuMandal2007].

Examples

>>> from pprint import pprint
>>> from dtk.bicycle import basu_sig_figs
>>> pprint(basu_sig_figs())
{'betaf': 0,
 'betafd': 14,
 'betafdd': 13,
 'betar': 0,
 'betard': 13,
 'betardd': 14,
 'phi': 14,
 'phid': 12,
 'phidd': 13,
 'psi': 13,
 'psid': 13,
 'psidd': 14,
 'psif': 13,
 'psifd': 13,
 'psifdd': 14,
 'theta': 0,
 'thetad': 13,
 'thetadd': 13,
 'x': 0,
 'xd': 14,
 'xdd': 13,
 'y': 0,
 'yd': 13,
 'ydd': 13,
 'z': 13,
 'zd': 13,
 'zdd': 13}
dtk.bicycle.basu_table_one_input()

Examples

>>> from pprint import pprint
>>> from dtk.bicycle import basu_table_one_input
>>> pprint(basu_table_one_input())
{'betaf': 0.0,
 'betafd': 8.0133620584155,
 'betar': 0.0,
 'betard': 8.912989661489,
 'phi': 3.1257073014894,
 'phid': -0.0119185528069,
 'psi': 0.9501292851472,
 'psid': 0.6068425835418,
 'psif': 0.2311385135743,
 'psifd': 0.4859824687093,
 'theta': 0.0,
 'thetad': 0.7830033527065,
 'x': 0.0,
 'xd': -2.8069345714545,
 'y': 0.0,
 'yd': -0.1480982396001,
 'z': 0.2440472102925,
 'zd': 0.1058778746261}
dtk.bicycle.basu_table_one_output()

Examples

>>> from pprint import pprint
>>> from dtk.bicycle import basu_table_one_output
>>> pprint(basu_table_one_output())
{'betafdd': 2.454807290455,
 'betardd': 1.8472554144217,
 'phidd': 0.1205543897884,
 'psidd': -7.8555281128244,
 'psifdd': -4.6198904039403,
 'thetadd': 0.8353281706379,
 'xdd': -0.5041626315047,
 'ydd': -0.3449706619454,
 'zdd': -1.460452833298}
dtk.bicycle.basu_to_moore_input(basu, rr, lam)

Returns the coordinates and speeds of the [Moore2012] derivation of the Whipple bicycle model as a function of the states and speeds of the [BasuMandal2007] coordinates and speeds.

Parameters:
basudictionary

A dictionary containing the states and speeds of the Basu-Mandal formulation. The states are represented with words corresponding to the greek letter and the speeds are the words with d appended, e.g. psi and psid.

rrfloat

Rear wheel radius.

lamfloat

Steer axis tilt.

Returns:
mooredictionary

A dictionary with the coordinates, q’s, and speeds, u’s, for the Moore formulation.

Examples

>>> import numpy as np
>>> from pprint import pprint
>>> from dtk.bicycle import basu_table_one_input, basu_to_moore_input
>>> vars = basu_table_one_input()
>>> pprint(basu_to_moore_input(vars, 0.3, np.pi/10))
{'q1': -0.0,
 'q2': -0.17447337661787718,
 'q3': -0.0,
 'q4': 0.6206670416476966,
 'q5': 0.3300446174593725,
 'q6': -0.0,
 'q7': -0.2311385135743,
 'q8': -0.0,
 'u1': 2.6703213326046784,
 'u2': -2.453592884421596e-14,
 'u3': -0.7830033527065,
 'u4': -0.6068425835418,
 'u5': 0.0119185528069,
 'u6': -8.912989661489,
 'u7': -0.4859824687093,
 'u8': -8.0133620584155}
dtk.bicycle.benchmark_matrices()

Returns the entries to the M, C1, K0, and K2 matrices for the benchmark parameter set printed in [Meijaard2007].

Returns:
Mndarray, shape(2,2)

The mass matrix.

C1ndarray, shape(2,2)

The speed proportional damping matrix.

K0ndarray, shape(2,2)

The gravity proportional stiffness matrix.

K2ndarray, shape(2,2)

The speed squared proportional stiffness matrix.

Notes

The equations of motion take this form:

M * q’’ + v * C1 * q’ + [g * K0 + v**2 * K2] * q’ = f

where q = [roll angle,

steer angle]

and f = [roll torque,

steer torque]

Examples

>>> from dtk.bicycle import benchmark_matrices
>>> M, C1, K0, K2 = benchmark_matrices()
>>> M
array([[80.81722   ,  2.31941332],
       [ 2.31941332,  0.29784188]])
>>> C1
array([[ 0.        , 33.86641391],
       [-0.85035641,  1.68540397]])
>>> K0
array([[-80.95      ,  -2.59951685],
       [ -2.59951685,  -0.80329488]])
>>> K2
array([[ 0.        , 76.5973459 ],
       [ 0.        ,  2.65431524]])
dtk.bicycle.benchmark_par_to_canonical(p)

Returns the canonical matrices of the Whipple bicycle model linearized about the upright constant velocity configuration. It uses the parameter definitions from [Meijaard2007].

Parameters:
pdictionary

A dictionary of the benchmark bicycle parameters. Make sure your units are correct, best to ue the benchmark paper’s units!

Returns:
Mndarray, shape(2,2)

The mass matrix.

C1ndarray, shape(2,2)

The damping like matrix that is proportional to the speed, v.

K0ndarray, shape(2,2)

The stiffness matrix proportional to gravity, g.

K2ndarray, shape(2,2)

The stiffness matrix proportional to the speed squared, v**2.

Examples

>>> from dtk.bicycle import benchmark_parameters, benchmark_par_to_canonical
>>> M, C1, K0, K2 = benchmark_par_to_canonical(benchmark_parameters())
>>> M
array([[80.81722   ,  2.31941332],
       [ 2.31941332,  0.29784188]])
>>> C1
array([[ 0.        , 33.86641391],
       [-0.85035641,  1.68540397]])
>>> K0
array([[-80.95      ,  -2.59951685],
       [ -2.59951685,  -0.80329488]])
>>> K2
array([[ 0.        , 76.5973459 ],
       [ 0.        ,  2.65431524]])
dtk.bicycle.benchmark_parameters()

Returns the benchmark bicycle parameters from [Meijaard2007].

Examples

>>> from pprint import pprint
>>> from dtk.bicycle import benchmark_parameters
>>> pprint(benchmark_parameters())
{'IBxx': 9.2,
 'IBxz': 2.4,
 'IByy': 11.0,
 'IBzz': 2.8,
 'IFxx': 0.1405,
 'IFyy': 0.28,
 'IHxx': 0.05892,
 'IHxz': -0.00756,
 'IHyy': 0.06,
 'IHzz': 0.00708,
 'IRxx': 0.0603,
 'IRyy': 0.12,
 'c': 0.08,
 'g': 9.81,
 'lam': 0.3141592653589793,
 'lambda': 0.3141592653589793,
 'mB': 85.0,
 'mF': 3.0,
 'mH': 4.0,
 'mR': 2.0,
 'rF': 0.35,
 'rR': 0.3,
 'w': 1.02,
 'xB': 0.3,
 'xH': 0.9,
 'zB': -0.9,
 'zH': -0.7}
dtk.bicycle.benchmark_state_space(M, C1, K0, K2, v, g)

Calculate the A and B matrices for the Whipple bicycle model linearized about the upright configuration.

Parameters:
Mndarray, shape(2,2)

The mass matrix.

C1ndarray, shape(2,2)

The damping like matrix that is proportional to the speed, v.

K0ndarray, shape(2,2)

The stiffness matrix proportional to gravity, g.

K2ndarray, shape(2,2)

The stiffness matrix proportional to the speed squared, v**2.

vfloat

Forward speed.

gfloat

Acceleration due to gravity.

Returns:
Andarray, shape(4,4)

System dynamic matrix.

Bndarray, shape(4,2)

Input matrix.

The states are [roll angle,

steer angle, roll rate, steer rate]

The inputs are [roll torque,

steer torque]

Examples

>>> from dtk.bicycle import benchmark_matrices, benchmark_state_space
>>> A, B = benchmark_state_space(*benchmark_matrices(), 5.2, 9.81)
>>> A
array([[  0.        ,   0.        ,   1.        ,   0.        ],
       [  0.        ,   0.        ,   0.        ,   1.        ],
       [  9.48977445, -24.66951001,  -0.54871674,  -1.71868007],
       [ 11.71947687, -22.40642251,  19.11938721, -16.04130074]])
>>> B
array([[ 0.        ,  0.        ],
       [ 0.        ,  0.        ],
       [ 0.01593498, -0.12409203],
       [-0.12409203,  4.32384018]])
dtk.bicycle.benchmark_state_space_vs_speed(M, C1, K0, K2, speeds=None, v0=0.0, vf=10.0, num=50, g=9.81)

Returns the state and input matrices for a set of speeds.

Parameters:
Marray_like, shape(2,2)

The mass matrix.

C1array_like, shape(2,2)

The speed proportional damping matrix.

K0array_like, shape(2,2)

The gravity proportional stiffness matrix.

K2array_like, shape(2,2)

The speed squared proportional stiffness matrix.

speedsarray_like, shape(n,), optional

An array of speeds in meters per second at which to compute the state and input matrices. If none, the v0, vf, and num parameters are used to generate a linearly spaced array.

v0float, optional, default: 0.0

The initial speed.

vffloat, optional, default: 10.0

The final speed.

numint, optional, default: 50

The number of speeds.

gfloat, optional, default: 9.81

Acceleration due to gravity in meters per second squared.

Returns:
speedsndarray, shape(n,)

An array of speeds in meters per second.

Asndarray, shape(n,4,4)

The state matrices evaluated at each speed in speeds.

Bsndarray, shape(n,4,2)

The input matrices

Notes

The second order equations of motion take this form:

M * q’’ + v * C1 * q’ + [g * K0 + v**2 * K2] * q’ = f

where q = [roll angle,

steer angle]

and f = [roll torque,

steer torque]

The first order equations of motion take this form:

x’ = A * x + B * u

where x = [roll angle,

steer angle, roll rate, steer rate]

and u = [roll torque,

steer torque]

Examples

>>> from dtk.bicycle import benchmark_matrices, benchmark_state_space_vs_speed
>>> M, C1, K0, K2 = benchmark_matrices()
>>> vs, As, Bs = benchmark_state_space_vs_speed(M, C1, K0, K2, num=3)
>>> vs
array([ 0.,  5., 10.])
>>> As
array([[[   0.        ,    0.        ,    1.        ,    0.        ],
        [   0.        ,    0.        ,    0.        ,    1.        ],
        [   9.48977445,   -0.57152317,   -0.        ,   -0.        ],
        [  11.71947687,   30.90875339,   -0.        ,   -0.        ]],

       [[   0.        ,    0.        ,    1.        ,    0.        ],
        [   0.        ,    0.        ,    0.        ,    1.        ],
        [   9.48977445,  -22.85146663,   -0.52761225,   -1.65257699],
        [  11.71947687,  -18.38412373,   18.38402617,  -15.42432764]],

       [[   0.        ,    0.        ,    1.        ,    0.        ],
        [   0.        ,    0.        ,    0.        ,    1.        ],
        [   9.48977445,  -89.69129698,   -1.0552245 ,   -3.30515399],
        [  11.71947687, -166.26275511,   36.76805233,  -30.84865527]]])
>>> Bs
array([[[ 0.        ,  0.        ],
        [ 0.        ,  0.        ],
        [ 0.01593498, -0.12409203],
        [-0.12409203,  4.32384018]],

       [[ 0.        ,  0.        ],
        [ 0.        ,  0.        ],
        [ 0.01593498, -0.12409203],
        [-0.12409203,  4.32384018]],

       [[ 0.        ,  0.        ],
        [ 0.        ,  0.        ],
        [ 0.01593498, -0.12409203],
        [-0.12409203,  4.32384018]]])
dtk.bicycle.benchmark_to_moore(benchmarkParameters, oldMassCenter=False)

Returns the parameters for the Whipple model as derived by Jason K. Moore.

Parameters:
benchmarkParametersdictionary

Contains the set of parameters for the Whipple bicycle model as presented in [Meijaard2007].

oldMassCenterboolean

If true it returns the fork mass center dimensions, l3 and l4, with respect to the rear offset intersection with the steer axis, otherwise the dimensions are with respect to the front wheel.

Returns:
mooreParametersdictionary

The parameter set for the Moore derivation of the whipple bicycle model as presented in [Moore2012].

Examples

>>> from pprint import pprint
>>> from dtk.bicycle import benchmark_parameters, benchmark_to_moore
>>> par = benchmark_parameters()
>>> pprint(benchmark_to_moore(par))
{'d1': 0.9534570696121849,
 'd2': 0.2676445084476887,
 'd3': 0.03207142672761929,
 'g': 9.81,
 'ic11': 7.178169776497895,
 'ic12': 0.0,
 'ic22': 11.0,
 'ic23': 0.0,
 'ic31': 3.8225535938357873,
 'ic33': 4.821830223502103,
 'id11': 0.0603,
 'id22': 0.12,
 'id33': 0.0603,
 'ie11': 0.05841337700152972,
 'ie12': 0.0,
 'ie22': 0.06,
 'ie23': 0.0,
 'ie31': 0.009119225261946298,
 'ie33': 0.007586622998470264,
 'if11': 0.1405,
 'if22': 0.28,
 'if33': 0.1405,
 'l1': 0.4707271515135145,
 'l2': -0.47792881146460797,
 'l3': -0.00597083392418685,
 'l4': -0.3699518200282974,
 'mc': 85.0,
 'md': 2.0,
 'me': 4.0,
 'mf': 3.0,
 'rf': 0.35,
 'rr': 0.3}
dtk.bicycle.front_contact(q1, q2, q3, q4, q7, d1, d2, d3, rr, rf, guess=None)

Returns the location in the ground plane of the front wheel contact point.

Parameters:
q1float

The location of the rear wheel contact point with respect to the inertial origin along the 1 axis (forward).

q2float

The location of the rear wheel contact point with respect to the inertial origin along the 2 axis (right).

q3float

The yaw angle.

q4float

The roll angle.

q7float

The steer angle.

d1float

The distance from the rear wheel center to the steer axis.

d2float

The distance between the front and rear wheel centers along the steer axis.

d3float

The distance from the front wheel center to the steer axis.

rrfloat

The radius of the rear wheel.

rffloat

The radius of the front wheel.

guessfloat, optional

A guess for the pitch angle. This may be only needed for extremely large steer and roll angles.

Returns:
q9float

The location of the front wheel contact point with respect to the inertial origin along the 1 axis.

q10float

The location of the front wheel contact point with respect to the inertial origin along the 2 axis.

Examples

>>> import numpy as np
>>> from dtk.bicycle import front_contact
>>> front_contact(0.0, 0.0,
...               np.deg2rad(5.0), np.deg2rad(5.0), np.deg2rad(5.0),
...               0.6, 0.3, 0.03, 0.3, 0.3)
(0.6987001194987257, 0.05266663513621053)
dtk.bicycle.lambda_from_abc(rF, rR, a, b, c)

Returns the steer axis tilt, lamba, for the parameter set based on the offsets from the steer axis.

Parameters:
rFfloat

Front wheel radius.

rRfloat

Rear wheel radius.

afloat

The rear wheel offset from the steer axis.

bfloat

The front wheel offset from the steer axis.

cfloat

The distance along the steer axis between the front wheel and rear wheel.

Returns:
lamfloat

The steer axis tilt as described in [Meijaard2007].

Examples

>>> from dtk.bicycle import lambda_from_abc
>>> import numpy as np
>>> float(np.rad2deg(lambda_from_abc(0.31, 0.29, 1.0, 0.1, 0.5)))
25.392364580504058
dtk.bicycle.meijaard_figure_four(time, rollRate, steerRate, speed)

Returns a figure that matches Figure #4 in [Meijaard2007].

import numpy as np
from scipy.signal import lti, lsim
from dtk.bicycle import (benchmark_matrices, benchmark_state_space,
                         meijaard_figure_four)
t = np.linspace(0.0, 5.0)
x0 = np.array([0.0, 0.0, 0.5, 0.0])
speed = 4.6  # m/s
A, B = benchmark_state_space(*benchmark_matrices(), speed, 9.81)
C, D = np.eye(4), np.zeros((4, 2))
system = lti(A, B, C, D)
t, y, _ = lsim(system, 0.0, t, X0=x0)
meijaard_figure_four(t, y[:, 2], y[:, 3], speed*np.ones_like(t))

(Source code, png, hires.png, pdf)

_images/dtk-1.png
dtk.bicycle.moore_to_basu(moore, rr, lam)

Returns the coordinates, speeds, and accelerations in [BasuMandal2007]’s convention.

Parameters:
mooredictionary

A dictionary containg values for the q’s, u’s and u dots.

rrfloat

Rear wheel radius.

lamfloat

Steer axis tilt.

Returns:
basudictionary

A dictionary containing the coordinates, speeds and accelerations.

Examples

>>> import numpy as np
>>> from pprint import pprint
>>> from dtk.bicycle import basu_table_one_input
>>> from dtk.bicycle import basu_to_moore_input, moore_to_basu
>>> rr, lam = 0.3, np.pi/10
>>> basu = basu_table_one_input()
>>> pprint(basu)
{'betaf': 0.0,
 'betafd': 8.0133620584155,
 'betar': 0.0,
 'betard': 8.912989661489,
 'phi': 3.1257073014894,
 'phid': -0.0119185528069,
 'psi': 0.9501292851472,
 'psid': 0.6068425835418,
 'psif': 0.2311385135743,
 'psifd': 0.4859824687093,
 'theta': 0.0,
 'thetad': 0.7830033527065,
 'x': 0.0,
 'xd': -2.8069345714545,
 'y': 0.0,
 'yd': -0.1480982396001,
 'z': 0.2440472102925,
 'zd': 0.1058778746261}
>>> moore = basu_to_moore_input(basu, rr, lam)
>>> pprint(moore)
{'q1': -0.0,
 'q2': -0.17447337661787718,
 'q3': -0.0,
 'q4': 0.6206670416476966,
 'q5': 0.3300446174593725,
 'q6': -0.0,
 'q7': -0.2311385135743,
 'q8': -0.0,
 'u1': 2.6703213326046784,
 'u2': -2.453592884421596e-14,
 'u3': -0.7830033527065,
 'u4': -0.6068425835418,
 'u5': 0.0119185528069,
 'u6': -8.912989661489,
 'u7': -0.4859824687093,
 'u8': -8.0133620584155}
>>> moore['u1p'] = 1.0
>>> moore['u2p'] = 2.0
>>> moore['u3p'] = 3.0
>>> moore['u4p'] = 4.0
>>> moore['u5p'] = 5.0
>>> moore['u6p'] = 6.0
>>> moore['u7p'] = 7.0
>>> moore['u8p'] = 8.0
>>> pprint(moore_to_basu(moore, rr, lam))
{'betaf': 0.0,
 'betafd': 8.0133620584155,
 'betafdd': -8.0,
 'betar': 0.0,
 'betard': 8.912989661489,
 'betardd': -6.0,
 'phi': 3.1257073014894,
 'phid': -0.0119185528069,
 'phidd': -5.0,
 'psi': 0.9501292851472,
 'psid': 0.6068425835418,
 'psidd': -4.0,
 'psif': 0.2311385135743,
 'psifd': 0.4859824687093,
 'psifdd': -7.0,
 'theta': 0.0,
 'thetad': 0.7830033527065,
 'thetadd': -3.0,
 'x': 0.0,
 'xd': -2.8069345714545,
 'xdd': -0.24465703387278925,
 'y': 0.0,
 'yd': -0.14809823960010002,
 'ydd': 2.804969014148545,
 'z': 0.2440472102925096,
 'zd': 0.10587787462605407,
 'zdd': -0.7877658248084111}
dtk.bicycle.pitch_from_roll_and_steer(q4, q7, rF, rR, d1, d2, d3, guess=None)

Returns the pitch angle of the bicycle frame for a given roll, steer and geometry.

Parameters:
q4float

Roll angle.

q5float

Steer angle.

rFfloat

Front wheel radius.

rRfloat

Rear wheel radius.

d1float

The rear wheel offset from the steer axis.

d2float

The distance along the steer axis between the intersection of the front and rear offset lines.

d3float

The front wheel offset from the steer axis.

guessfloat, optional

A good guess for the pitch angle. If not specified, the program will make a good guess for most roll and steer combinations.

Returns:
q5float

Pitch angle.

Notes

All of the geometry parameters should be expressed in the same units.

Examples

>>> import numpy as np
>>> from dtk.bicycle import pitch_from_roll_and_steer
>>> from dtk.bicycle import benchmark_parameters, benchmark_to_moore
>>> steer, roll = np.deg2rad(10.0), np.deg2rad(-5.0)
>>> p = benchmark_to_moore(benchmark_parameters())
>>> float(np.rad2deg(pitch_from_roll_and_steer(steer, roll,
...                                            p['rf'], p['rr'],
...                                            p['d1'], p['d2'], p['d3'])))
18.062710178550127
dtk.bicycle.sort_modes(evals, evecs)

Sort eigenvalues and eigenvectors into weave, capsize, caster modes.

Parameters:
evalsndarray, shape (n, 4)

eigenvalues

evecsndarray, shape (n, 4, 4)

eigenvectors

Returns:
weave[‘evals’]ndarray, shape (n, 2)

The eigen value pair associated with the weave mode.

weave[‘evecs’]ndarray, shape (n, 4, 2)

The associated eigenvectors of the weave mode.

capsize[‘evals’]ndarray, shape (n,)

The real eigenvalue associated with the capsize mode.

capsize[‘evecs’]ndarray, shape(n, 4, 1)

The associated eigenvectors of the capsize mode.

caster[‘evals’]ndarray, shape (n,)

The real eigenvalue associated with the caster mode.

caster[‘evecs’]ndarray, shape(n, 4, 1)

The associated eigenvectors of the caster mode.

This only works on the standard bicycle eigenvalues, not necessarily on any
general eigenvalues for the bike model (e.g. there isn’t always a distinct
weave, capsize and caster). Some type of check unsing the derivative of the
curves could make it more robust.

Examples

>>> from dtk.bicycle import (benchmark_matrices,
...                          benchmark_state_space_vs_speed, sort_modes)
>>> from dtk.control import eig_of_series
>>> M, C1, K0, K2 = benchmark_matrices()
>>> _, A, _ = benchmark_state_space_vs_speed(M, C1, K0, K2)
>>> weave, capsize, caster = sort_modes(*eig_of_series(A))
>>> weave['evals'][0:2]
array([[-5.53094372+0.j, -3.13164325+0.j],
       [-5.8702391 +0.j, -3.11751166+0.j]])
>>> capsize['evals'][0:2]
array([3.13164325+0.j, 3.16834073+0.j])
>>> caster['evals'][0:2]
array([5.53094372+0.j, 5.16831044+0.j])
dtk.bicycle.trail(rF, lam, fo)

Returns the trail and mechanical trail.

Parameters:
rF: float

The front wheel radius

lam: float

The steer axis tilt (pi/2 - headtube angle). The angle between the headtube and a vertical line.

fo: float

The fork offset

Returns:
c: float

Trail

cm: float

Mechanical Trail

Examples

>>> import numpy as np
>>> from dtk.bicycle import trail
>>> trail(0.3, np.deg2rad(10.0), 0.05)
(0.002126763618252235, 0.0020944533000790966)

control Module

class dtk.control.Bode(frequency, *args, **kwargs)

Bases: object

A class for creating Bode plots and the associated data.

Parameters:
frequencyndarray, shape(n,)

An array of frequencies at which to evaluate the system frequency reponse in radians per second. Use numpy.logspace to generate them.

*argssequence of dtk.control.StateSpace objects

One or more state space systems. If more than one system is provided, they must all have the same inputs and outputs.

Examples

import numpy as np
from dtk.bicycle import benchmark_matrices, benchmark_state_space
from dtk.control import StateSpace, Bode

speed = 4.6  # m/s
A, B = benchmark_state_space(*benchmark_matrices(), speed, 9.81)
C, D = np.eye(4), np.zeros((4, 2))

states = ['Roll Angle', 'Steer Angle', 'Roll Rate', 'Steer Rate']
inputs = ['Roll Torque', 'Steer Torque']

sys = StateSpace(A, B, C, D,
    name='Carvallo-Whipple Bicycle',
    stateNames=states,
    inputNames=inputs,
    outputNames=states,
)

freqs = np.logspace(0.0, 3.0, num=400)

bode = Bode(freqs, sys)

bode.plot()

(Source code)

_images/dtk-2_00.png

(png, hires.png, pdf)

_images/dtk-2_01.png

(png, hires.png, pdf)

_images/dtk-2_02.png

(png, hires.png, pdf)

_images/dtk-2_03.png

(png, hires.png, pdf)

_images/dtk-2_04.png

(png, hires.png, pdf)

_images/dtk-2_05.png

(png, hires.png, pdf)

_images/dtk-2_06.png

(png, hires.png, pdf)

_images/dtk-2_07.png

(png, hires.png, pdf)

mag_phase()

Computes the magnitude and phase for all the systems in the Bode object. This is called on instantiation.

mag_phase_system(system)

Returns the magnitude and phase for a single system.

Parameters:
systemdtk.control.StateSpace

A state space system.

Returns:
magnitudendarray, shape(n, m, p)

An array with the magnitude of the input-output transfer functions for each frequency.

phasendarray, shape(n, m, p)

An array with the phase of the in input-output transfer functions for each frequency in radians per second.

Notes

  • n : number of frequencies

  • m : number of outputs

  • p : number of inputs

Examples

>>> import numpy as np
>>> from dtk.bicycle import benchmark_matrices, benchmark_state_space
>>> from dtk.control import StateSpace, Bode
>>> speed = 4.6  # m/s
>>> A, B = benchmark_state_space(*benchmark_matrices(), speed, 9.81)
>>> C, D = np.eye(4), np.zeros((4, 2))
>>> states = ['Roll Angle', 'Steer Angle', 'Roll Rate', 'Steer Rate']
>>> inputs = ['Roll Torque', 'Steer Torque']
>>> sys = StateSpace(A, B, C, D,
...     name='Carvallo-Whipple Bicycle',
...     stateNames=states,
...     inputNames=inputs,
...     outputNames=states,
... )
>>> freqs = np.logspace(0.0, 3.0, num=400)
>>> bode = Bode(freqs, sys)
>>> mag, phase = bode.mag_phase_system(sys)
>>> mag[:3]
array([[[0.01169673, 0.38514231],
        [0.00676025, 0.21053334],
        [0.01169673, 0.38514231],
        [0.00676025, 0.21053334]],

       [[0.01158059, 0.38127272],
        [0.0067131 , 0.20907207],
        [0.01178282, 0.38793104],
        [0.00683034, 0.21272318]],

       [[0.01146535, 0.37743215],
        [0.00666676, 0.20763614],
        [0.01186929, 0.39072976],
        [0.00690164, 0.21495153]]])
>>> phase[:3]
array([[[-0.9832055 ,  2.09569516],
        [-1.00369635,  1.99821008],
        [ 0.58759082, -2.61669382],
        [ 0.56709997, -2.7141789 ]],

       [[-0.99052413,  2.08728955],
        [-1.01180201,  1.98810804],
        [ 0.5802722 , -2.62509943],
        [ 0.55899431, -2.72428094]],

       [[-0.99778294,  2.07892517],
        [-1.01988136,  1.97801793],
        [ 0.57301338, -2.63346381],
        [ 0.55091497, -2.73437105]]])
plot(**kwargs)

Plots the Bode plots for all systems in the Bode object.

Parameters:
**kwargsdictionary

Sets the color and linestyle attributes on this object and passes the rest through to plot_system.

plot_system(system, magnitude, phase, decibel=True, degree=True, **kwargs)

Plots the Bode plots of a single system. If a system for this object has already been plotted, it will add new lines to the existing plots.

Parameters:
systemdtk.control.StateSpace

The state space system.

magnitudendarray, shape(n, m, p)

An array with the magnitude of the input-output transfer functions for each frequency.

phasendarray, shape(n, m, p)

An array with the phase of the in input-output transfer functions for each frequency in radians per second.

Examples

import numpy as np
from dtk.bicycle import benchmark_matrices, benchmark_state_space
from dtk.control import StateSpace, Bode

speed = 4.6  # m/s
A, B = benchmark_state_space(*benchmark_matrices(), speed, 9.81)
C, D = np.array([1.0, 0.0, 0.0, 0.0]).reshape(1, 4), np.zeros((1, 1))

states = ['Roll Angle', 'Steer Angle', 'Roll Rate', 'Steer Rate']
inputs = ['Roll Torque', 'Steer Torque']
outputs = ['Roll Angle']

sys = StateSpace(A, B, C, D,
    name='Carvallo-Whipple Bicycle',
    stateNames=states,
    inputNames=inputs,
    outputNames=outputs,
)

freqs = np.logspace(0.0, 3.0, num=400)

bode = Bode(freqs, sys)

mag, phase = bode.mag_phase_system(sys)

bode.plot_system(sys, mag, phase, decibel=False, degree=False)

(Source code)

_images/dtk-3_00.png

(png, hires.png, pdf)

_images/dtk-3_01.png

(png, hires.png, pdf)

show()

Shows all figures stored in the object.

class dtk.control.StateSpace(A, B, C, D, **kwargs)

Bases: object

A linear time invariant system described by its state space.

Parameters:
Andarray, shape(n,n)

The state matrix.

Bndarray, shape(n,p)

The input matrix.

Cndarray, shape(m,n)

The output matrix.

Dndarray, shape(m,p)

The feedforward matrix.

namestring, optional

A name of the system.

stateNameslist, len(n), optional

A list of names of each state in order corresponding to A.

inputNameslist, len(p), optional

A list of names of each input in order corresponding to B.

outputNameslist, len(m), optional

A list of names of each output in order corresponding to C.

Examples

>>> import numpy as np
>>> from dtk.bicycle import benchmark_matrices, benchmark_state_space
>>> from dtk.control import StateSpace, Bode
>>> speed = 4.6  # m/s
>>> A, B = benchmark_state_space(*benchmark_matrices(), speed, 9.81)
>>> C, D = np.eye(4), np.zeros((4, 2))
>>> states = ['Roll Angle', 'Steer Angle', 'Roll Rate', 'Steer Rate']
>>> inputs = ['Roll Torque', 'Steer Torque']
>>> sys = StateSpace(A, B, C, D,
...     name='Carvallo-Whipple Bicycle',
...     stateNames=states,
...     inputNames=inputs,
...     outputNames=states,
... )
>>> print(sys)
A Carvallo-Whipple Bicycle system with 4 states, 2 inputs, and 4 outputs.
>>> sys.numStates, sys.numInputs, sys.numOutputs
(4, 2, 4)
>>> sys.A
array([[  0.        ,   0.        ,   1.        ,   0.        ],
       [  0.        ,   0.        ,   0.        ,   1.        ],
       [  9.48977445, -19.42926731,  -0.48540327,  -1.52037084],
       [ 11.71947687, -10.81273781,  16.91330407, -14.19038143]])
>>> sys.B
array([[ 0.        ,  0.        ],
       [ 0.        ,  0.        ],
       [ 0.01593498, -0.12409203],
       [-0.12409203,  4.32384018]])
>>> sys.C
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])
>>> sys.D
array([[0., 0.],
       [0., 0.],
       [0., 0.],
       [0., 0.]])
dtk.control.bode(system, frequency, fig=None, label=None, title=None, color=None)

Creates a Bode plot of the given system.

Parameters:
systemtuple

The system can be defined as a state space or the numerator and denominator of a transfer function. If defined in state space form it should include ndarrays for the state, input, output and feed-forward matrices, in that order. These should only be defined for a single input and single output. If in transfer function form the ndarrays for the numerator and denomonator coefficients must be provided.

frequencyndarray

An array of frequencies at which to evaluate the system frequency reponse in radians per second.

figmatplotlib Figure instance, optional
Returns:
magnitudendarray

The magnitude in dB of the frequency response.

phasendarray

The phase in degrees of the frequency response.

figmatplotlib Figure instance

The Bode plot.

Examples

import numpy as np
from dtk.bicycle import benchmark_matrices, benchmark_state_space
from dtk.control import bode

speed = 4.6  # m/s
A, B = benchmark_state_space(*benchmark_matrices(), speed, 9.81)
C, D = np.array([1.0, 0.0, 0.0, 0.0]), np.zeros(1)

freqs = np.logspace(0.0, 3.0, num=301)

bode((A, B[:, 0].reshape(4, 1), C, D), freqs)

(Source code, png, hires.png, pdf)

_images/dtk-4.png
bode((A, B[:, 0].reshape(4, 1), C, D), freqs,
    label='Nice Curve',
    title='My Bode Plot',
    color='black',
)

(Source code, png, hires.png, pdf)

_images/dtk-5.png
dtk.control.eig_of_series(matrices)

Returns the eigenvalues and eigenvectors for a series of matrices.

Parameters:
matricesarray_like, shape(n, m, m)

A series of square matrices.

Returns:
eigenvaluesndarray, shape(n, m)

The eigenvalues of the matrices.

eigenvectorsndarray, shape(n, m, m)

The eigenvectors of the matrices.

Examples

>>> import matplotlib.pyplot as plt
>>> from dtk.bicycle import (benchmark_matrices,
...                          benchmark_state_space_vs_speed)
>>> from dtk.control import eig_of_series, sort_modes, plot_phasor
>>> M, C1, K0, K2 = benchmark_matrices()
>>> v, A, B = benchmark_state_space_vs_speed(M, C1, K0, K2)
>>> evals, evecs = eig_of_series(A)
>>> evals[0:5]
array([[ 5.53094372+0.j        ,  3.13164325+0.j        ,
        -5.53094372+0.j        , -3.13164325+0.j        ],
       [ 5.16831044+0.j        ,  3.16834073+0.j        ,
        -5.8702391 +0.j        , -3.11751166+0.j        ],
       [ 4.76080633+0.j        ,  3.24904588+0.j        ,
        -6.19610903+0.j        , -3.11594235+0.j        ],
       [ 4.22644752+0.j        ,  3.45546901+0.j        ,
        -6.51427475+0.j        , -3.12094054+0.j        ],
       [ 3.67619382+0.52184908j,  3.67619382-0.52184908j,
        -3.12830521+0.j        , -6.82848078+0.j        ]])
dtk.control.plot_phasor(eigenvalues, eigenvectors, components=None, compNames=None, show=False)

Returns a phasor plot of the given eigenvalues and eigenvectors.

Parameters:
eigenvaluesarray_like, shape(n, )

The eigenvalues.

eigenvectorsarray_like, shape(n, n)

The eigenvectors where each column corresponds to the eigenvalues.

componentsarray_like, optional

The indices of the eigenvector components to plot.

showboolean, optional, default False

If true the plots will be displayed.

Returns:
figslist

A list of matplotlib figures.

Notes

Plots are not produced for zero eigenvalues.

Examples

import matplotlib.pyplot as plt
from dtk.bicycle import (benchmark_matrices,
                         benchmark_state_space_vs_speed)
from dtk.control import eig_of_series, sort_modes, plot_phasor

M, C1, K0, K2 = benchmark_matrices()
v, A, B = benchmark_state_space_vs_speed(M, C1, K0, K2)
evals, evecs = sort_modes(*eig_of_series(A))
plot_phasor(evals[25], evecs[25])

(Source code)

_images/dtk-6_00.png

(png, hires.png, pdf)

_images/dtk-6_01.png

(png, hires.png, pdf)

_images/dtk-6_02.png

(png, hires.png, pdf)

_images/dtk-6_03.png

(png, hires.png, pdf)

dtk.control.plot_root_locus(parvalues, eigenvalues, skipZeros=False, fig=None, parName=None, parUnits=None, **kwargs)

Returns a root locus plot of a series of eigenvalues with respect to a series of values.

Parameters:
parvaluesarray_like, shape(n,)

The parameter values corresponding to each eigenvalue.

eigenvaluesarray_like, shape(n,m)

The m eigenvalues for each parameter value.

skipZerosboolean, optional, default = False

If true any eigenvalues close to zero will not be plotted.

figmatplotlib.Figure, optional, default = None

Pass in a figure to plot on.

parNamestring, optional

Specify the name or abbreviation of the parameter name.

parUnitsstring, optional

Specify the units of the parameter.

**kwargsvaries

Any option keyword argments for a matplotlib scatter plot.

Returns:
figmatplotlib.Figure

Examples

from dtk.bicycle import (benchmark_matrices,
                         benchmark_state_space_vs_speed)
from dtk.control import eig_of_series, plot_root_locus

M, C1, K0, K2 = benchmark_matrices()
v, A, B = benchmark_state_space_vs_speed(M, C1, K0, K2)
evals, evecs = eig_of_series(A)
plot_root_locus(v, evals, parName='Speed', parUnits='[m/s]')

(Source code, png, hires.png, pdf)

_images/dtk-7.png
dtk.control.plot_root_locus_components(parvalues, eigenvalues, parts='both', parName=None, parUnits=None, skipZeros=True, ax=None, **kwargs)

Returns a root locus plot of a series of eigenvalues with respect to a series of values.

Parameters:
parvaluesarray_like, shape(n,)

The parameter values corresponding to each eigenvalue.

eigenvaluesarray_like, shape(n,m)

The m eigenvalues for each parameter value.

partsstring, optional, {'both'|'real'|'imaginary'}

Specify whether both the real and imaginary lines should be plotted or one or the other. Default is 'both'.

parNamestring, optional

Specify the name or abbreviation of the parameter name.

parUnitsstring, optional

Specify the units of the parameter.

skipZerosboolean, optional

If true (default) any eigenvalues close to zero will not be plotted.

**kwargsvaries

Any option keyword argments for the matplotlib plot function. This will be applied to all lines.

Returns:
figmatplotlib.Figure

Nothing is returned if an axis is provided.

Examples

import matplotlib.pyplot as plt
from dtk.bicycle import (benchmark_matrices,
                         benchmark_state_space_vs_speed)
from dtk.control import (eig_of_series, sort_modes,
                         plot_root_locus_components)

M, C1, K0, K2 = benchmark_matrices()
v, A, B = benchmark_state_space_vs_speed(M, C1, K0, K2)
evals, evecs = sort_modes(*eig_of_series(A))
fig, ax = plt.subplots(layout='constrained')
plot_root_locus_components(v, evals, parName='Speed',
                           parUnits='[m/s]', ax=ax)

(Source code, png, hires.png, pdf)

_images/dtk-8.png
dtk.control.sort_modes(evals, evecs)

Sort a series of eigenvalues and eigenvectors into modes.

Parameters:
evalsarray_like, shape (n, m)

eigenvalues

evecsarray_like, shape (n, m, m)

eigenvectors

Examples

>>> import matplotlib.pyplot as plt
>>> from dtk.bicycle import (benchmark_matrices,
...                          benchmark_state_space_vs_speed)
>>> from dtk.control import eig_of_series, sort_modes, plot_phasor
>>> M, C1, K0, K2 = benchmark_matrices()
>>> v, A, B = benchmark_state_space_vs_speed(M, C1, K0, K2)
>>> evals, evecs = eig_of_series(A)
>>> evals[0:5]
array([[ 5.53094372+0.j        ,  3.13164325+0.j        ,
        -5.53094372+0.j        , -3.13164325+0.j        ],
       [ 5.16831044+0.j        ,  3.16834073+0.j        ,
        -5.8702391 +0.j        , -3.11751166+0.j        ],
       [ 4.76080633+0.j        ,  3.24904588+0.j        ,
        -6.19610903+0.j        , -3.11594235+0.j        ],
       [ 4.22644752+0.j        ,  3.45546901+0.j        ,
        -6.51427475+0.j        , -3.12094054+0.j        ],
       [ 3.67619382+0.52184908j,  3.67619382-0.52184908j,
        -3.12830521+0.j        , -6.82848078+0.j        ]])
>>> evals, evecs = sort_modes(evals, evecs)
>>> evals[0:5]
array([[ 5.53094372+0.j        ,  3.13164325+0.j        ,
        -5.53094372+0.j        , -3.13164325+0.j        ],
       [ 5.16831044+0.j        ,  3.16834073+0.j        ,
        -5.8702391 +0.j        , -3.11751166+0.j        ],
       [ 4.76080633+0.j        ,  3.24904588+0.j        ,
        -6.19610903+0.j        , -3.11594235+0.j        ],
       [ 4.22644752+0.j        ,  3.45546901+0.j        ,
        -6.51427475+0.j        , -3.12094054+0.j        ],
       [ 3.67619382+0.52184908j,  3.67619382-0.52184908j,
        -6.82848078+0.j        , -3.12830521+0.j        ]])

inertia Module

dtk.inertia.compound_pendulum_inertia(m, g, l, T)

Returns the moment of inertia for an object hung as a compound pendulum.

Parameters:
mfloat

Mass of the pendulum.

gfloat

Acceration due to gravity.

lfloat

Length of the pendulum.

Tfloat

The period of oscillation.

Returns:
float

Moment of interia of the pendulum.

Examples

>>> from dtk.inertia import compound_pendulum_inertia
>>> compound_pendulum_inertia(3.0, 9.81, 0.2, 1.4)
0.1722244785902121
dtk.inertia.cylinder_inertia(l, m, ro, ri)

Calculate the moment of inertia for a hollow cylinder (or solid cylinder) where the x axis is aligned with the cylinder’s axis.

Parameters:
lfloat

The length of the cylinder.

mfloat

The mass of the cylinder.

rofloat

The outer radius of the cylinder.

rifloat

The inner radius of the cylinder. Set this to zero for a solid cylinder.

Returns:
Ixfloat

Moment of inertia about cylinder axis.

Iy, Izfloat

Moment of inertia about axis perpendicular to cylinder axis.

Examples

>>> from dtk.inertia import cylinder_inertia
>>> cylinder_inertia(1.0, 0.4, 0.02, 0.015)
(0.000125, 0.03339583333333333, 0.03339583333333333)
>>> cylinder_inertia(1.0, 0.4, 0.02, 0.0)
(8e-05, 0.03337333333333334, 0.03337333333333334)
dtk.inertia.euler_123(angles)

Returns the direction cosine matrix as a function of the Euler 123 angles.

Parameters:
anglesnumpy.array or list or tuple, shape(3,)

Three angles (in units of radians) that specify the orientation of a new reference frame with respect to a fixed reference frame. The first angle, phi, is a rotation about the fixed frame’s x-axis. The second angle, theta, is a rotation about the new y-axis (which is realized after the phi rotation). The third angle, psi, is a rotation about the new z-axis (which is realized after the theta rotation). Thus, all three angles are “relative” rotations with respect to the new frame. Note: if the rotations are viewed as occuring in the opposite direction (z, then y, then x), all three rotations are with respect to the initial fixed frame rather than “relative”.

Returns:
Rndarray, shape(3,3)

Three dimensional rotation matrix about three different orthogonal axes.

Examples

>>> import numpy as np
>>> from dtk.inertia import euler_123
>>> euler_123(np.deg2rad([12.0, 22.0, 45.0]))
array([[ 0.65561799, -0.65561799,  0.37460659],
       [ 0.74672788,  0.63658173, -0.19277236],
       [-0.11208268,  0.40611422,  0.90692266]])
dtk.inertia.euler_rotation(angles, order)

Returns a rotation matrix for a reference frame, B, in another reference frame, A, where the B frame is rotated relative to the A frame via body fixed rotations (Euler angles).

Parameters:
anglesarray_like

An array of three angles in radians that are in order of rotation.

ordertuple

A three tuple containing a combination of 1, 2, and 3 where 1 is about the x axis of the first reference frame, 2 is about the y axis of the this new frame and 3 is about the z axis. Note that (1, 1, 1) is a valid entry and will give you correct results, but combinations like this are not necessarily useful for describing a general configuration.

Returns:
Rndarray, shape(3,3)

A rotation matrix.

Notes

The rotation matrix is defined such that a R times a vector v equals the vector expressed in the rotated reference frame.

v’ = R * v

Where v is the vector expressed in the original reference frame and v’ is the same vector expressed in the rotated reference frame.

Examples

>>> import numpy as np
>>> from dtk.inertia import euler_rotation
>>> angles = [np.pi, np.pi / 2., -np.pi / 4.]
>>> rotMat = euler_rotation(angles, (3, 1, 3))
>>> rotMat
array([[-7.07106781e-01,  1.29893408e-16, -7.07106781e-01],
       [-7.07106781e-01,  4.32978028e-17,  7.07106781e-01],
       [ 1.22464680e-16,  1.00000000e+00,  6.12323400e-17]])
>>> v = np.array([[1.], [0.], [0.]])
>>> vp = rotMat @ v
>>> vp
array([[-7.07106781e-01],
       [-7.07106781e-01],
       [ 1.22464680e-16]])
dtk.inertia.inertia_components(jay, beta)

Returns the 2D orthogonal inertia tensor.

When at least three moments of inertia and their axes orientations are known relative to a common inertial frame of a planar object, the orthoganal moments of inertia relative the frame are computed.

Parameters:
jayndarray, shape(n,)

An array of at least three moments of inertia. (n >= 3)

betandarray, shape(n,)

An array of orientation angles corresponding to the moments of inertia in jay.

Returns:
ndarray, shape(3,)

Ixx, Ixz, Izz

Examples

>>> import numpy as np
>>> from dtk.inertia import inertia_components
>>> inertia_components([1.2, 0.5, 3.1], np.deg2rad([45.0, 90.0, 135.0]))
array([3.8 , 0.95, 0.5 ])
>>> inertia_components([1.2, 0.5, 0.51, 3.1],
...                    np.deg2rad([45.0, 90.0, 90.2, 135.0]))
array([3.79833626, 0.95000581, 0.50166378])
dtk.inertia.parallel_axis(Ic, m, d)

Returns the moment of inertia of a body about a different point.

Parameters:
Icndarray, shape(3,3)

The moment of inertia about the center of mass of the body with respect to an orthogonal coordinate system.

mfloat

The mass of the body.

dndarray, shape(3,)

The distances along the three ordinates that located the new point relative to the center of mass of the body.

Returns:
Indarray, shape(3,3)

The moment of inertia of a body about a point located by the distances in d.

Examples

>>> import numpy as np
>>> from dtk.bicycle import benchmark_parameters
>>> from dtk.inertia import parallel_axis
>>> p = benchmark_parameters()
>>> Ic = np.array([[p['IBxx'], 0.0, p['IBxz']],
...                [0.0, p['IByy'], 0.0],
...                [p['IBxz'], 0.0, p['IBzz']]])
>>> d = np.array([-p['xB'], 0.0, -p['zB']])  # about rear wheel contact
>>> parallel_axis(Ic, p['mB'], d)
array([[78.05,  0.  , 25.35],
       [ 0.  , 87.5 ,  0.  ],
       [25.35,  0.  , 10.45]])
dtk.inertia.principal_axes(I)

Returns the principal moments of inertia and the orientation.

Parameters:
Indarray, shape(3,3)

An inertia tensor.

Returns:
Ipndarray, shape(3,)

The principal moments of inertia. This is sorted smallest to largest.

Cndarray, shape(3,3)

The rotation matrix.

Examples

>>> import numpy as np
>>> from dtk.bicycle import benchmark_parameters
>>> from dtk.inertia import principal_axes
>>> p = benchmark_parameters()
>>> Ic = np.array([[p['IBxx'], 0.0, p['IBxz']],
...                [0.0, p['IByy'], 0.0],
...                [p['IBxz'], 0.0, p['IBzz']]])
>>> Ip, C = principal_axes(Ic)
>>> Ip
array([ 2., 10., 11.])
>>> C
array([[-0.31622777,  0.        ,  0.9486833 ],
       [ 0.9486833 ,  0.        ,  0.31622777],
       [ 0.        ,  1.        ,  0.        ]])
>>> C @ Ic @ C.T
array([[ 2.00000000e+00, -5.28515252e-17,  0.00000000e+00],
       [-3.40171594e-16,  1.00000000e+01,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  1.10000000e+01]])
dtk.inertia.rotate3(angles)

Produces a three-dimensional rotation matrix as rotations around the three cartesian axes.

Parameters:
anglesarray_like, shape(3,)

Three angles (in units of radians) that specify the orientation of a new reference frame with respect to a fixed reference frame. The first angle is a pure rotation about the x-axis, the second about the y-axis, and the third about the z-axis. All rotations are with respect to the initial fixed frame, and they occur in the order x, then y, then z.

Returns:
Rndarray, shape(3,3)

Three dimensional rotation matrix about three different orthogonal axes.

Examples

>>> import numpy as np
>>> from dtk.inertia import rotate3
>>> rotate3(np.deg2rad([12.0, 22.0, 45.0]))
array([[ 0.65561799, -0.63658173,  0.40611422],
       [ 0.65561799,  0.74672788,  0.11208268],
       [-0.37460659,  0.19277236,  0.90692266]])
dtk.inertia.rotate3_inertia(RotMat, relInertia)

Rotates an inertia tensor. A derivation of the formula in this function can be found in Crandall 1968, Dynamics of mechanical and electromechanical systems. This function only transforms an inertia tensor for rotations with respect to a fixed point. To translate an inertia tensor, one must use the parallel axis analogue for tensors. An inertia tensor contains both moments of inertia and products of inertia for a mass in a cartesian (xyz) frame.

Parameters:
RotMatarray_like, shape(3,3)

Three-dimensional rotation matrix specifying the coordinate frame that the input inertia tensor is in, with respect to a fixed coordinate system in which one desires to express the inertia tensor.

relInertiaarray_like, shape(3,3)

Three-dimensional cartesian inertia tensor describing the inertia of a mass in a rotated coordinate frame.

Returns:
Inertiandarray, shape(3,3)

Inertia tensor with respect to a fixed coordinate system (“unrotated”).

Examples

>>> import numpy as np
>>> from dtk.bicycle import benchmark_parameters
>>> from dtk.inertia import principal_axes, rotate3_inertia
>>> p = benchmark_parameters()
>>> Ic = np.array([[p['IBxx'], 0.0, p['IBxz']],
...                [0.0, p['IByy'], 0.0],
...                [p['IBxz'], 0.0, p['IBzz']]])
>>> Ip, C = principal_axes(Ic)
>>> Ip
array([ 2., 10., 11.])
>>> C
array([[-0.31622777,  0.        ,  0.9486833 ],
       [ 0.9486833 ,  0.        ,  0.31622777],
       [ 0.        ,  1.        ,  0.        ]])
>>> rotate3_inertia(C, Ic)
array([[ 2.00000000e+00, -5.28515252e-17,  0.00000000e+00],
       [-3.40171594e-16,  1.00000000e+01,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  1.10000000e+01]])
dtk.inertia.rotate_inertia_about_y(I, angle)

Returns inertia tensor rotated through angle about the Y axis.

Parameters:
Indarray, shape(3, 3)

An inertia tensor.

anglefloat

Angle in radians about the positive Y axis of which to rotate the inertia tensor.

Returns:
ndarray, shape(3, 3)

Rotated inerita tensor.

Examples

>>> import numpy as np
>>> from dtk.inertia import rotate_inertia_about_y
>>> rotate_inertia_about_y(np.diag([1.0, 2.0, 3.0]), np.deg2rad(45.0))
array([[ 2.,  0., -1.],
       [ 0.,  2.,  0.],
       [-1.,  0.,  2.]])
dtk.inertia.torsional_pendulum_inertia(k, T)

Calculate the moment of inertia for an ideal torsional pendulum.

Parameters:
kfloat

Torsional stiffness.

Tfloat

Period of oscillation.

Returns:
float

Moment of inertia.

Examples

>>> from dtk.inertia import torsional_pendulum_inertia
>>> torsional_pendulum_inertia(50.0, 1.0)
1.2665147955292222
dtk.inertia.total_com(coordinates, masses)

Returns the center of mass of a group of objects if the indivdual centers of mass and mass is provided.

coordinatesndarray, shape(3,n)

The rows are the x, y and z coordinates, respectively and the columns are for each object.

massesndarray, shape(3,)

An array of the masses of multiple objects, the order should correspond to the columns of coordinates.

Returns:
mTfloat

Total mass of the objects.

cTndarray, shape(3,)

The x, y, and z coordinates of the total center of mass.

Examples

>>> import numpy as np
>>> from dtk.bicycle import benchmark_parameters
>>> from dtk.inertia import total_com
>>> par = benchmark_parameters()
>>> coordinates = np.array([[par['xB'], par['xH']],
...                         [0.0, 0.0],
...                         [par['zB'], par['zH']]])
...
>>> coordinates
array([[ 0.3,  0.9],
       [ 0. ,  0. ],
       [-0.9, -0.7]])
>>> masses = np.array([par['mB'], par['mH']])
>>> masses
array([85.,  4.])
>>> total_com(coordinates, masses)
(89.0, array([ 0.32696629,  0.        , -0.89101124]))
dtk.inertia.tube_inertia(l, m, ro, ri)

Calculate the moment of inertia for a tube (or rod) where the x axis is aligned with the tube’s axis.

Parameters:
lfloat

The length of the tube.

mfloat

The mass of the tube.

rofloat

The outer radius of the tube.

rifloat

The inner radius of the tube. Set this to zero if it is a rod instead of a tube.

Returns:
Ixfloat

Moment of inertia about tube axis.

Iy, Izfloat

Moment of inertia about normal axis.

Examples

>>> from dtk.inertia import tube_inertia
>>> tube_inertia(1.0, 0.4, 0.02, 0.015)
(0.000125, 0.03339583333333333, 0.03339583333333333)
dtk.inertia.x_rot(angle)

Returns the rotation matrix for a reference frame rotated through an angle about the x axis.

Parameters:
anglefloat

The angle in radians.

Returns:
Rxndarray, shape(3,3)

The rotation matrix.

Notes

v’ = Rx * v where v is the vector expressed the reference in the original reference frame and v’ is the vector expressed in the new rotated reference frame.

Examples

>>> import numpy as np
>>> from dtk.inertia import x_rot
>>> x_rot(np.deg2rad(45.0))
array([[ 1.        ,  0.        ,  0.        ],
       [ 0.        ,  0.70710678,  0.70710678],
       [ 0.        , -0.70710678,  0.70710678]])
dtk.inertia.y_rot(angle)

Returns the rotation matrix for a reference frame rotated through an angle about the y axis.

Parameters:
anglefloat

The angle in radians.

Returns:
Rxndarray, shape(3,3)

The rotation matrix.

Notes

v’ = Rx * v where v is the vector expressed the reference in the original reference frame and v’ is the vector expressed in the new rotated reference frame.

Examples

>>> import numpy as np
>>> from dtk.inertia import y_rot
>>> y_rot(np.deg2rad(45.0))
array([[ 0.70710678,  0.        , -0.70710678],
       [ 0.        ,  1.        ,  0.        ],
       [ 0.70710678,  0.        ,  0.70710678]])
dtk.inertia.z_rot(angle)

Returns the rotation matrix for a reference frame rotated through an angle about the z axis.

Parameters:
anglefloat

The angle in radians.

Returns:
Rxndarray, shape(3,3)

The rotation matrix.

Notes

v’ = Rx * v where v is the vector expressed the reference in the original reference frame and v’ is the vector expressed in the new rotated reference frame.

Examples

>>> import numpy as np
>>> from dtk.inertia import z_rot
>>> z_rot(np.deg2rad(45.0))
array([[ 0.70710678,  0.70710678,  0.        ],
       [-0.70710678,  0.70710678,  0.        ],
       [ 0.        ,  0.        ,  1.        ]])

process Module

dtk.process.butterworth(data, cutoff, samplerate, order=2, axis=-1, btype='lowpass', **kwargs)

Returns the data filtered by a two-pass (forward/backward) Butterworth filter.

Parameters:
dataarray_like, shape(n, ) or shape(n, m)

The data to filter. Only handles 1D and 2D arrays.

cutofffloat

The filter cutoff frequency in Hertz.

sampleratefloat

The sample rate of the data in Hertz. Sample rate must be constant.

orderint

The order of the Butterworth filter.

axisint

The axis to filter along.

btype{‘lowpass’|’highpass’}

The type of filter. Default is ‘lowpass’.

kwargs

Any extra arguments to get passed to scipy.signal.sosfiltfilt.

Returns:
filtered_datandarray

The low pass filtered version of data.

Notes

The provided cutoff frequency is corrected by a multiplicative factor to ensure the double pass filter cutoff frequency matches that of a single pass filter, see [Winter2009].

References

[Winter2009]

David A. Winter (2009) Biomechanics and motor control of human movement. 4th edition. Hoboken: Wiley.

Examples

import numpy as np
import matplotlib.pyplot as plt
from dtk.process import butterworth, freq_spectrum

sample_rate = 1000  # Hz
duration = 10.0  # seconds
time = np.linspace(0.0, duration, num=int(sample_rate*duration) + 1)
white_noise = np.random.normal(0.0, 1.0, size=len(time))
cutoff = 200  # Hz
order = 4
filtered = butterworth(white_noise, cutoff, sample_rate, order=order)

freq, amp = freq_spectrum(white_noise, sample_rate)
freq_filt, amp_filt = freq_spectrum(filtered, sample_rate)

fig, ax = plt.subplots(layout='constrained')
ax.plot(freq, amp, label='Unfiltered')
ax.plot(freq_filt, amp_filt, alpha=0.75, label='Filtered')
ax.axvline(cutoff, color='black')
ax.set_ylabel('Amplitude of White Noise with STD=1')
ax.set_xlabel('Frequency [Hz]')
msg = 'Sample rate: {} Hz, Cutoff: {} Hz, Order: {}'
ax.set_title(msg.format(sample_rate, cutoff, order))
ax.legend()

(Source code, png, hires.png, pdf)

_images/dtk-9.png
dtk.process.coefficient_of_determination(measured, predicted)

Computes the coefficient of determination with respect to a measured and predicted array.

Parameters:
measuredarray_like, shape(n,)

The observed or measured values.

predictedarray_like, shape(n,)

The values predicted by a model.

Returns:
r_squaredfloat

The coefficient of determination.

Notes

The coefficient of determination [also referred to as R^2 and VAF (variance accounted for)] is computed either of these two ways:

      sum( [predicted - mean(measured)] ** 2 )
R^2 = ----------------------------------------
      sum( [measured - mean(measured)] ** 2 )

or:

          sum( [measured - predicted] ** 2 )
R^2 = 1 - ---------------------------------------
          sum( [measured - mean(measured)] ** 2 )

Examples

>>> import numpy as np
>>> from dtk.process import coefficient_of_determination
>>> np.random.seed(10)
>>> t = np.linspace(0.0, 10.0, num=1001)
>>> predicted = np.sin(t)
>>> measured = predicted + np.random.normal(0.01, 0.1, size=len(t))
>>> coefficient_of_determination(measured, predicted)
0.980225686442542
dtk.process.cumulative_power_spectrum(data, sample_rate, relative=True, remove_dc_component=False)

Return the cumulative power spectrum of a signal:

S(f) = \sum_{k=0}^f |X(k)|^2
Parameters:
datandarray, shape (m,) or shape(n,m)

The array of time signals where n is the number of variables and m is the number of time steps.

sample_rateint

The signal sampling rate in Hertz.

relativebool, optional

If True, the returned amplitued is expressed relative to the total power. The default is True.

remove_dc_componentbool, optional

If True, the DC component (f = 0) is not included in the returned spectrum ]0,f_N[. If False the returned spectrum covers [0, f_N[. The default is False.

Returns:
frequencyndarray, shape (p,)

The frequencies where p is a power of 2 close to m.

cumulative_powerndarray, shape (p,n)

The cumulative power up to each frequency.

Notes

  • cumulative_power_spectrum() performs zero-padding. Parseval’s theorem is satisfied for the padded input signal. Provide input signals with 2^p samples to prevent zero-padding.

  • The power contributions of positive and negative frequencies are combined in the positive half spectrum so that the results satisfy Parseval’s theoreom on the interval [0, f_N].

  • If the dc component is removed with remove_dc_component=True, the results do not satisfy Parseval’s theorem.

Examples

Create the cumulative power spectrum of a rect pulse and plot in time and frequency domain.

import numpy as np
import matplotlib.pyplot as plt
from dtk.process import cumulative_power_spectrum

# sampling parameters
N = 64  # signal period
f_s = 10  # sample rate
T = N/f_s

t = np.arange(0, T, 1/f_s)

# rect test signal
A = 3  # amplitude
tau = 0.2*T  # "on"-time

x = np.zeros_like(t)
x[0:int(tau*f_s)] = A

# power spectrum
freq, amp = cumulative_power_spectrum(x, f_s)

# plot
fig, ax = plt.subplots(2,1, layout="constrained")
ax[0].stem(t, x, )
ax[0].set_xlabel("$t$ in s")
ax[0].set_ylabel("$x(t)$")
ax[1].stem(freq,amp)
ax[1].set_xlabel("$f$ in Hz")
ax[1].set_ylabel("cumulative avg. power")
plt.suptitle((f"Sample rate: {f_s} Hz, Signal period: {T} s,"
             " relative=True"))

(Source code, png, hires.png, pdf)

_images/dtk-10.png
dtk.process.curve_area_stats(x, y)

Return the box plot stats of a curve based on area.

Parameters:
xndarray, shape (n,)

The x values

yndarray, shape (n,m)

The y values n are the time steps m are the various curves

Returns:
A dictionary containing:
medianndarray, shape (m,)

The x value corresponding to 0.5*area under the curve

lqndarray, shape (m,)

lower quartile

uqndarray, shape (m,)

upper quartile

98pndarray, shape (m,)

98th percentile

2pndarray, shape (m,)

2nd percentile

Examples

>>> from pprint import pprint
>>> import numpy as np
>>> from dtk.process import curve_area_stats
>>> x = np.linspace(0.0, 10.0, num=1001)
>>> y = np.vstack((np.exp(x), 0.5*x)).T
>>> pprint(curve_area_stats(x, y))
{'2p': array([6.09, 1.41]),
 '98p': array([9.97, 9.89]),
 'lq': array([8.61, 5.  ]),
 'median': array([9.3 , 7.07]),
 'uq': array([9.71, 8.66])}
dtk.process.derivative(x, y, method='forward', padding=None)

Returns the derivative of y with respect to x.

Parameters:
xndarray, shape(n,)

The monotonically increasing independent variable.

yndarray, shape(n,) or shape(n, m)

The dependent variable(s).

methodstring, optional
‘forward’

Use the forward difference method.

‘backward’

Use the backward difference method.

‘central’

Use the central difference method.

‘combination’

This is equivalent to method='central', padding='second order' and is in place for backwards compatibility. Selecting this method will ignore and user supplied padding settings.

paddingNone, float, ‘adjacent’ or ‘second order’, optional

The default, None, will result in the derivative vector being n-a in length where a=1 for forward and backward and a=2 for central. If you provide a float this value will be used to pad the result so that len(dydx) == n. If ‘adjacent’ is used, the nearest neighbor will be used for padding. If ‘second order’ is chosen second order foward and backward difference are used to pad the end points.

Returns:
dydxndarray, shape(n,) or shape(n-1,)

for combination else shape(n-1,)

Examples

import numpy as np
import matplotlib.pyplot as plt
from dtk.process import derivative

x = np.linspace(-10.0, 10.0, num=201)
y = x**2

fig, axes = plt.subplots(2, 1, layout='constrained')

axes[0].plot(x, y)
axes[1].plot(x, derivative(x, y, method='combination'))

(Source code, png, hires.png, pdf)

_images/dtk-11.png
dtk.process.find_timeshift(signal1, signal2, sample_rate, guess=None, plot=False)

Returns the timeshift, tau, of the second signal relative to the first signal.

Parameters:
signal1array_like, shape(n, )

The base signal.

signal2array_like, shape(n, )

A signal shifted relative to the first signal. The second signal should be leading the first signal.

sample_rateinteger or float

Sample rate of the signals. This should be the same for each signal.

guessfloat, optional, default=None

If you’ve got a good guess for the time shift then supply it here.

plotboolean, optional, defaul=False

If true, a plot of the error landscape will be shown.

Returns:
taufloat

The timeshift between the two signals.

Examples

import numpy as np
import matplotlib.pyplot as plt
from dtk.process import find_timeshift

t = np.linspace(0.0, 4.0, num=401)
sig1 = np.sin(2.0*t) + np.random.normal(0.0, 0.1, size=len(t))
sig2 = np.sin(2.0*t + 0.3) + np.random.normal(0.0, 0.1, size=len(t))

tau = find_timeshift(sig1, sig2, 100, guess=0.2)

fig, ax = plt.subplots(layout='constrained')
ax.plot(t, sig1, t, sig2)
ax.legend(['Signal 1', 'Signal 2'])
ax.set_title('Shift: {:1.3f} s'.format(tau))
ax.set_xlabel('Time [s]')

(Source code, png, hires.png, pdf)

_images/dtk-12.png
dtk.process.fit_goodness(ym, yp)

Calculate the goodness of fit.

Parameters:
ymndarray, shape(n,)

The vector of measured values.

ypndarry, shape(n,)

The vector of predicted values.

Returns:
rsqfloat

The r squared value of the fit.

SSEfloat

The error sum of squares.

SSTfloat

The total sum of squares.

SSRfloat

The regression sum of squares.

Notes

SST = SSR + SSE

Examples

>>> import numpy as np
>>> from dtk.process import fit_goodness
>>> np.random.seed(10)
>>> t = np.linspace(0.0, 10.0, num=1001)
>>> predicted = np.sin(t)
>>> measured = predicted + np.random.normal(0.01, 0.1, size=len(t))
>>> fit_goodness(measured, predicted)
(0.9862246380225299, 6.197685397793293, 449.9108922095648, 443.7132068117715)
dtk.process.freq_spectrum(data, sampleRate, norm='forward', remove_dc_component=True)

Return the frequency spectrum of a data set. Combines negative and positive frequencies in the positive frequency range. Returns frequencies up until the Nyquist frequency f_N.

Parameters:
dataarray_like, shape (m, ) or shape(n, m)

The array of time signals where n is the number of variables and m is the number of time steps.

sampleRateint

The signal sampling rate in Hertz.

normstr, optional

Normalization of the returned spectrum. See https://numpy.org/doc/stable/reference/routines.fft.html#normalization for explanation. The default is “forward”, which normalizes the frequency spectrum by 1/N.

remove_dc_componentbool, optional

If True, the DC component (f = 0) is not included in the returned spectrum ]0,f_N[. If False the returned spectrum covers [0, f_N[. The default is True.

Returns:
frequencyndarray, shape (p,)

Frequencies where p is a power of 2 close to m, in Hertz.

amplitudendarray, shape (p, n)

Amplitude at each frequency.

Examples

Create a sum of two sinusoids with the low frequency sinusoid having an amplitude of 2 and frequency of 5 Hz and the high frequency sinusoid having an amplitude of 1 and frequency of 50 Hz. Plot the frequency spectrum of the sum.

import numpy as np
import matplotlib.pyplot as plt
from dtk.process import freq_spectrum

sample_rate = 1000
duration = 1.0
time = np.linspace(0.0, duration, num=int(duration*sample_rate) + 1)

low_freq = 2.0*np.sin(5.0*2.0*np.pi*time)  # 5 Hz * 2 pi rad / cycle
high_freq = np.sin(50.0*2.0*np.pi*time)  # 50 Hz * 2 pi rad / cycle

fig, ax = plt.subplots(layout='constrained')
ax.plot(time, low_freq + high_freq)
ax.set_xlim((0.0, 0.4))
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')

(Source code, png, hires.png, pdf)

_images/dtk-13.png
freqs, amps = freq_spectrum(low_freq + high_freq, sample_rate)

fig, ax = plt.subplots(layout='constrained')
ax.plot(freqs, amps)
ax.set_xlim((0.0, 100.0))
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Amplitude')

(Source code, png, hires.png, pdf)

_images/dtk-14.png
dtk.process.least_squares_variance(A, sum_of_residuals)

Returns the variance in the ordinary least squares fit and the covariance matrix of the estimated parameters.

Parameters:
Andarray, shape(n,d)

The left hand side matrix in Ax=B.

sum_of_residualsfloat

The sum of the residuals (residual sum of squares).

Returns:
variancefloat

The variance of the fit.

covariancendarray, shape(d,d)

The covariance of x in Ax = b.

dtk.process.normalize(sig, hasNans=False)

Normalizes the vector with respect to the maximum value.

Parameters:
signdarray, shape(n,)
hasNansboolean, optional

If your data has nans use this flag if you want to ignore them.

Returns:
normSigndarray, shape(n,)

The signal normalized with respect to the maximum value.

Examples

>>> import numpy as np
>>> from dtk.process import normalize
>>> t = np.linspace(0.0, 2*np.pi, num=11)
>>> y = 5.0*np.sin(t)
>>> float(np.max(normalize(y)))
1.0
dtk.process.power_spectrum(data, sample_rate, remove_dc_component=False)

Return the power spectrum of a signal:

S(f) = |X(f)|^2
Parameters:
datandarray, shape (m,) or shape(n,m)

The array of time signals where n is the number of variables and m is the number of time steps.

sample_rateint

The signal sampling rate in Hertz.

remove_dc_componentbool, optional

If True, the DC component (f = 0) is not included in the returned spectrum ]0,f_N[. If False the returned spectrum covers [0, f_N[. The default is True.

Returns:
frequencyndarray, shape (p,)

The frequencies where p is a power of 2 close to m.

powerndarray, shape (p,n)

The power at each frequency.

Notes

  • power_spectrum() performs zero-padding. Parseval’s theorem is satisfied for the padded input signal. Provide input signals with 2^p samples to prevent zero-padding.

  • The power contributions of positive and negative frequencies are combined in the positive half spectrum so that the results satisfy Parseval’s theoreom on the interval [0, f_N].

  • If the dc component is removed with remove_dc_component=True, the results do not satisfy Parseval’s theorem.

Examples

Create the power spectrum of a rect pulse and plot in time and frequency domain. Note how the power of frequencies f>0 is larger then f=0 because the positive frequencies inlclude the contribution of negative frequencies. As a result, the mean power in the displayed half spectrum equals the the mean power of the input signal.

import numpy as np
import matplotlib.pyplot as plt
from dtk.process import power_spectrum

# sampling parameters
N = 64   # signal period
f_s = 10  # sample rate
T = N/f_s

t = np.arange(0, T, 1/f_s)

# rectangle test signal
A = 3  # amplitude
tau = 0.2*T # "on"-time

x = np.zeros_like(t)
x[0:int(tau*f_s)] = A

# power spectrum
freq, amp = power_spectrum(x, f_s)

# check Parseval's theorem
energy_time = np.mean(np.abs(x)**2)
energy_freq = np.sum(amp)

print(f"Mean power in time domain: {energy_time:.6f}")
print(f"Mean power in frequency domain: {energy_freq:.6f}")

# plot
fig, ax = plt.subplots(2, 1, layout="constrained")
ax[0].stem(t, x)
ax[0].set_xlabel("$t$ in s")
ax[0].set_ylabel("$x(t)$")
ax[1].stem(freq,amp)
ax[1].set_xlabel("$f$ in Hz")
ax[1].set_ylabel("$|X(f)|^2$")
plt.suptitle(f"Sample rate: {f_s} Hz, Signal period: {T} s")

(Source code, png, hires.png, pdf)

_images/dtk-15.png
dtk.process.spline_over_nan(x, y)

Returns a vector of which a cubic spline is used to fill in gaps in the data from nan values.

Parameters:
xarray_like, shape(n,)

This x values should not contain nans.

yarray_like, shape(n,)

The y values may contain nans.

Returns:
ndarray, shape(n,)

The splined y values. If y doesn’t contain any nans the y is returned unmodified.

Notes

The splined data is identical to the input data, except that the nan’s are replaced by new data from the spline fit.

Examples

import numpy as np
from dtk.process import spline_over_nan

x = np.linspace(0.0, 2.0, num=201)
y = np.sin(3*2*np.pi*x) + np.random.normal(0.0, 0.1, size=len(x))

y[78:89] = np.nan
y[95:102] = np.nan
y[189:192] = np.nan

y_splined = spline_over_nan(x, y)

fig, ax = plt.subplots(layout='constrained')
ax.fill_between(x, np.min(y_splined) - 0.5, np.max(y_splined) + 0.5,
                where=np.isnan(y), alpha=0.5, color='grey',
                transform=ax.get_xaxis_transform())
ax.plot(x, y_splined, linewidth=4, color='black', label='Filled')
ax.plot(x, y, linestyle='', marker='o', color='red', label='Missing')
ax.legend()

(Source code, png, hires.png, pdf)

_images/dtk-16.png
dtk.process.subtract_mean(sig, hasNans=False)

Subtracts the mean from a signal with nanmean.

Parameters:
signdarray, shape(n,)
hasNansboolean, optional

If your data has nans use this flag if you want to ignore them.

Returns:
ndarray, shape(n,)

sig minus the mean of sig

Examples

>>> import numpy as np
>>> from dtk.process import subtract_mean
>>> t = np.linspace(0.0, 2*np.pi, num=11)
>>> y1 = np.sin(t)
>>> y2 = np.sin(t) + 0.3
>>> print(np.allclose(y1, subtract_mean(y2)))
True
>>> y2[2:5] = np.nan
>>> subtract_mean(y2, hasNans=True)
array([ 0.31123729,  0.89902254,         nan,         nan,         nan,
        0.31123729, -0.27654797, -0.63981923, -0.63981923, -0.27654797,
        0.31123729])
dtk.process.sync_error(tau, signal1, signal2, time, plot=False)

Returns the error between two signal time histories given a time shift, tau.

Parameters:
taufloat

The time shift.

signal1ndarray, shape(n,)

The signal that will be interpolated. This signal is typically “cleaner” that signal2 and/or has a higher sample rate.

signal2ndarray, shape(n,)

The signal that will be shifted to syncronize with signal 1.

timendarray, shape(n,)

The time vector for the two signals

plotboolean, optional, default=False

If true a plot will be shown of the resulting signals.

Returns:
errorfloat

Error between the two signals for the given tau.

Examples

import numpy as np
import matplotlib.pyplot as plt
from dtk.process import find_timeshift, sync_error

t = np.linspace(0.0, 4.0, num=401)
sig1 = np.sin(2.0*t) + np.random.normal(0.0, 0.1, size=len(t))
sig2 = np.sin(2.0*t + 0.3) + np.random.normal(0.0, 0.1, size=len(t))

tau = find_timeshift(sig1, sig2, 100, guess=0.2)

sync_error(tau, sig1, sig2, t, plot=True)

(Source code, png, hires.png, pdf)

_images/dtk-17.png
dtk.process.time_vector(num_samples, sample_rate, start_time=0.0)

Returns a time vector starting at zero.

Parameters:
num_samplesint

Total number of samples.

sample_ratefloat

Sample rate of the signal in hertz.

start_timefloat, optional, default=0.0

The start time of the time series.

Returns:
timendarray, shape(numSamples,)

Time vector starting at zero.

Examples

>>> from dtk.process import time_vector
>>> time_vector(101, 100)
array([0.  , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1 ,
       0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2 , 0.21,
       0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3 , 0.31, 0.32,
       0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4 , 0.41, 0.42, 0.43,
       0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5 , 0.51, 0.52, 0.53, 0.54,
       0.55, 0.56, 0.57, 0.58, 0.59, 0.6 , 0.61, 0.62, 0.63, 0.64, 0.65,
       0.66, 0.67, 0.68, 0.69, 0.7 , 0.71, 0.72, 0.73, 0.74, 0.75, 0.76,
       0.77, 0.78, 0.79, 0.8 , 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87,
       0.88, 0.89, 0.9 , 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98,
       0.99, 1.  ])
dtk.process.truncate_data(tau, signal1, signal2, sample_rate)

Returns the truncated vectors with respect to the time shift tau. It assume you’ve found the time shift between two signals with find_time_shift or something similar.

Parameters:
taufloat

The time shift.

signal1array_like, shape(n, )

A time series.

signal2array_like, shape(n, )

A time series.

sample_rateinteger

The sample rate of the two signals.

Returns:
truncated1ndarray, shape(m, )

The truncated time series.

truncated2ndarray, shape(m, )

The truncated time series.

Examples

import numpy as np
import matplotlib.pyplot as plt
from dtk.process import find_timeshift, truncate_data

t = np.linspace(0.0, 4.0, num=401)
sig1 = np.sin(2.0*t) + np.random.normal(0.0, 0.1, size=len(t))
sig2 = np.sin(2.0*t + 0.3) + np.random.normal(0.0, 0.1, size=len(t))

tau = find_timeshift(sig1, sig2, 100, guess=0.2)

sigtr1, sigtr2 = truncate_data(tau, sig1, sig2, 100)

fig, ax = plt.subplots(layout='constrained')
ax.plot(t[:len(sigtr1)], sigtr1, t[:len(sigtr2)], sigtr2)
ax.legend(['Signal 1', 'Signal 2'])
ax.set_title('Shift: {:1.3f} s'.format(tau))
ax.set_xlabel('Time [s]')

(Source code, png, hires.png, pdf)

_images/dtk-18.png