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
)
- 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()
- 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
andlinestyle
attributes on this object and passes the rest through toplot_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)
- 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
)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
)
- 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])
- 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
)
- 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
)
- 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
, and3
where1
is about the x axis of the first reference frame,2
is about the y axis of the this new frame and3
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
)
- 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
)
- 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
)
- 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
)
- 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 andm
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 tom
, 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
)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
)
- 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
)
- 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
)
- 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
)
- 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
)