DETAILED DESCRIPTION OF THE INVENTION

[0127] The present invention provides a method for quantifying muscle tone and an apparatus that carries out the method.

[0128] In general, the device in accordance with the present invention is an automated tone assessment device that non-invasively and properly quantifies tone. Quantification of tone allows a clinician to decide on the nature and degree of intervention to administer. A pre and post-evaluation can also determine the effectiveness of the intervention as well as track the progress of the patient. The system software can easily be modified to output the muscle stiffness of the patient to the clinician in real time. The implications of this are in the operating room where a neuro-surgeuon makes temporary lesions in the brain while monitoring the patient's muscle tone in real time. Once the surgeon finds the proper location in the brain that alleviates the degree of muscle tone with out too much loss in function, he/she can make the lesion permanent.

[0129] In an exemplary embodiment the device and method of the present invention are utilized to quantify the forearm muscle tone, in particular, the wrist in a spastic patient. However, it is to be understood that the present method and device are not limited to quantification of forearm muscle tone or to the spastic patient.

[0130] When used to quantify muscle tone in the spastic patient, upon determination of the muscle tone, the degree of spasticity that the patient has can be determined, thereby allowing a clinician to decide what form of intervention to take and to what degree. The muscle tone can further be monitored over time to allow the clinician to determine whether the intervention is providing benefits to the patient.

[0131] More particularly, the present invention utilizes equation 1 to determine the stiffness, viscosity and inertial parameters. 9$\begin{array}{cc}{\tau }_s\left(t\right)=K_H\theta \left(t\right)+B_H\stackrel{.}{\theta }\left(t\right)+J_T\ddot{\theta }\left(t\right)+{\tau }_{\mathrm{off}}& \text{[Eq. 1]}\end{array}$

[0132] While equation 1 has been used in the past, the methods and devices that were used to apply equation 1 utilized sinusoidal displacements of the ankle or finger joint at various frequencies. Using sinusoidal displacements, as set out above, leads to an indeterminate system, and it is impossible to isolate the stiffness values from the inertia values.

[0133] In accordance with the method of the present invention, a desired trajectory that is non-sinusoidal and not a ramp is first determined. The desired trajectory is then input into the device of the present invention, which is an electromechanical system that utilizes a feedback controller to track the desired displacement trajectory.

[0134] In accordance with the present method, a desired displacement of the wrist is first determined using the following set of conditions.

[0135] (1) The trajectory is ‘smooth’. In other words, when perturbing an appendage about its joint of rotation, there should be “no sudden impulsive movements that might synchronize a large number of sensory receptors in an artificial or unphysiological way.” (See P. M. H. Rack, H. F. Ross and T. I. H. Brown: Reflex Response during Sinusoidal Movements of Human Limbs. Prog. Clin. Neurophysiol., vol. 4, Ed. J. E. Desmedt, pp 216-228, 1978). As used herein, a “smooth discrete trajectory” is mathematically defined as follows:

[0136] θ_d(t_k) is a discrete signal consisting of set of N (1000) points that could be sampled from a continuous, differentiable function, θ_d(t) for all time t such that:

[0137] a. The sampling frequency (fhz_mc) of the discrete signal must be near to or greater than or equal to 100 Hz, particularly if the maximum frequency required to represent the discrete signal approaches or exceeds 10 Hz.

[0138] b. In the frequency domain, the maximum frequency needed to represent the desired trajectory must be less than the smaller of one-sixth frequency of the sampling rate of the motor controller, fhz_mc, or 15 Hz.

[0139] (2) At time, t=0, θ_d(0)=0, the motor does not jump rapidly to the initial point.

[0140] (3) The trajectory θ_d(t) is periodic (repeatable) and smooth for all time t>0. In other words θ_d(t_N+1) equals θ_d(t₁).

[0141] (4) The trajectory has a finite range of motion with a controllable, defined distribution.

[0142] (5) The matrix Ψ has linearly independent columns, which ensures that the matrix Ψ^TΨ is nonsingular.

[0143] Filtered random trajectories are then generated using, for example, MATLAB in the following way. There exist three discrete signals, r_a(t_k), r_b(t_k) and r_c(t_k) in the time domain which correspond to R_a(ω), R_b(ω), and R_c(ω) respectively in the frequency domain. The time between two consecutive points given to the controller is dt_mc seconds apart.

[0144] (1) A set of N(1000) uniformly distributed numbers is generated from −0.5 to 0.5. The frequency of the motor controller, fhz_mc, is 128 Hz, so consecutive points are dt_mc=0.0078 seconds apart in time. This discrete signal is called r_a(t_k).

[0145] (2) The discrete Fourier Transform of r_a(t_k) is taken, R_a(ω) to obtain the magnitude, |R_a(ω)| and phase, φ(R_a(ω)) of the signal in the frequency domain.

[0146] (3) Now, R_b(ω) is generated. For all ω, φ(R_b(ω))=φ(R_a(ω)). An ω_max is selected such that |R_b(0)|=0, for 0<ω≦ω_max, |R_b(ω)|=|R_a(ω)|, and for all ω>ω_max, |R_b(ω)|=0. In FIG. 2, ω_max=2π*4 Hz.

[0147] (4) The Inverse Fourier Transform of R_b(ω) is taken to obtain r_b(t). Note that the mean value of r_b(t_k)=0.

[0148] (5) There is no guarantee that r_b(0₁)=0, however, the trajectory does cross θ=0 at several instances in time. The index, min_i is found where r_b(tmin_i) has the lowest absolute value which is approximately 0. Note that r_b(t_k) is a continuous periodic function for all time with a period of T_θ=dt_mc*(N)=7.8125 s.

[0149] (6) r_c(t_k)=α(r_b(t_k+t_min—i)−r_b(t_min—i)). By shifting the signal r_b(t_k) back tmin_i seconds in time, r_c(0)=0. α is a scalar number which scales the range of motion of the trajectory. α was selected so that the root mean squared value of r_c(t_k) is 7000/{square root}{square root over (2)} counts. Thus, α=7000*RMS(r_b(t_k))/{square root}{square root over (2)}

[0150] (7) Finally, each data point in r_c(t_k) is rounded to the nearest integer which gives the desired trajectory, θ_d(t_k). This is done because the encoder (position sensor) and the controller deal with an integer number of digital counts.

[0151] The discrete set of finite data points defining the trajectory is then inputted into the controller, and a motor shaft tracks these points at a certain sampling rate. The desired trajectories generated in this fashion satisfy all conditions listed above. Faster desired trajectories can be further generated that cover the same distribution of angular positions from the same random number set. For example, as shown in FIG. 4, two random number sets were used with maximum frequency components of 4, 5, 6, 7, and 8 Hz to give 10 trials.

[0152] Generally, the trajectory is determined for testing of the right hand. However, if the left hand is to be tested, these trajectories for the right hand can simply be flipped or negated about the bias position, or relative origin, 0, to ensure symmetry across contralateral limbs.

[0153] The device in accordance with the present invention is an electromechanical system that is designed to drive the wrist with the filtered random trajectories described above. Because these trajectories have complex variable position and velocity profiles, an automated feedback controller, or robot, is used which constantly varies the torque, speed and direction of the actuator (motor). The desired trajectory is given as input, and the actuator moves back and fourth tracking the desired trajectory. Thus, design of the system in accordance with the present invention includes a mechanical design, the use of sensors and actuators, feedback control, data acquisition and programming.

[0154] One embodiment of the electromechanical system is shown in FIG. 5. As shown, the device includes a housing 1. The materials used in fabricating the housing are not particularly limited provided they supply the required structural support for holding a patient's arm on the top surface 2 of the housing 1 during testing. One preferred material is a metal such as aluminum. The housing 1 is shown as box-like in shape, however, the shape of the housing 1 is not limited and can come in a variety of shapes. Preferably, the top surface 2 of the housing 1, where the patient's arm rests during testing, is flat.

[0155] On the top surface 2 of the housing 1 is a hand securing mechanism 3 where the patent places his hand during the test. The hand securing mechanism 3 is designed such that the patient's hand can be placed sideways within the hand securing mechanism 3 with the palm facing the left (when the right hand is tested) and the thumb facing upwards. To accommodate various sized hands, the hand securing mechanism 3 is preferably adjustable to varying sizes.

[0156] In one preferred embodiment, as shown in FIG. 5, the hand securing mechanism 3 is formed of two vertical bars 6 that extend upwards from the top surface 2 of the housing 1 so that the patient can place his or her hand between the two bars as shown in FIG. 7. The two vertical bars 6 preferably have padding or cushions 7 on their outer surfaces to provide comfort, as shown in FIG. 11. The two vertical bars 6 may accommodate various sized hands by, for example, slidably attaching the two vertical bars 6 to the housing 1 so that they can be moved inwards and outwards and locked into a desired position. The padding or cushions 7 can also be removable and can come in a variety of thicknesses to provide smaller and larger openings between the two vertical bars 6, thereby accommodating various sized hands. In yet another embodiment, the padding or cushions 7 are compressible and expandable such that the padding or cushions 7 compress or expand as required to hold hands of various sizes in place.

[0157] A hand rest 4 and an arm rest 5 may further be mounted on top surface 2 of the housing 1 to assist in providing proper positioning of the arm and hand during testing. Preferably, the arm rest 5 includes an arm securing mechanism 10 such as, for example, one or more Velcro straps, as shown in FIG. 11, to secure the arm in place during testing.

[0158] In a preferred embodiment, a motor shaft 11 drives the hand about the wrist joint back and forth to follow the desired trajectories. In this embodiment, a direct drive system is preferably used wherein the hand rest 4 is directly connected to an extension of the motor shaft 11, called the output shaft 12. Preferably, there are no gears in between the motor shaft 11 and the output shaft 12 so that there is no jitteriness or backlash when the motor shaft 11 moves or changes direction. During use, the hand is positioned on the hand rest 4, which may be in the form of a small metal platform as shown in FIG. 6. The hand is secured as it rotates about the wrist joint.

[0159] The hand securing mechanism 3, e.g. two vertical bars 6, secure the area of the hand directly proximal to the first joint between the metacarpals and phalanges so that the knuckles lie past the hand securing mechanism 3, as shown in FIG. 11. The hand is secured so that the patient being tested does not grasp the hand securing mechanism 3. In fact, the patient should be completely relaxed so that all muscle contractions are completely involuntary due to the stretch reflex arc. Thus, the patient should not grasp the hand securing mechanism 3 during the test, as there are muscle groups that control flexion/extension of the wrist as well as the fingers.

[0160] In a preferred embodiment, the hand rest 4 lies on a sliding forked wrench 22 that is screwed into the output shaft 12. The hand rest 4 is then secured to a position on a sliding track, such that the axis of rotation of the wrist joint lies directly above the axis of rotation of the output shaft 12. During the perturbation of the hand, the remainder of the arm is secured to the arm rest 5 with the arm securing mechanism 10 as set out above. The arm rest 5 is preferably made out of low temperature thermoplastic that is covered with a thin layer of foam for comfort. In a preferred embodiment, the arm rest 5 is on an adjustable track, e.g. a friction screw track, that allows for minor adjustments in the positioning of the arm and hand if necessary.

[0161] The arm rest 5 and associated arm securing mechanism 10 hold the arm securely in place with quick setup time and is sized to accommodate various arm sizes. Removable cushions or pads of varying sizes can, for example, be placed in-between the arm rest 4 and the arm to accommodate various sized arms. In one embodiment, a compressible and expandable cushion (not shown) is mounted on the arm rest 5 such that the cushion fits about varying sized arms by compressing or expanding as necessary. Further, the arm rest 5 and arm securing mechanism 10 is preferably adjustable to varying sized arms. For example, in one embodiment, Velcro straps, which can be quickly and easily secured tighter or looser, are used as the arm rest arm securing mechanism 10. In another embodiment, one or more straps having multiple snaps along the length of the straps can be used for varying sized fits. Still further, one or more straps similar to a belt with multiple holes can be used to provide varying sized fits.

[0162] A torque transducer 20 (sensor) is placed between the motor shaft 11 and the output shaft 12. The torque transducer 20 may, for example, be placed between the motor shaft 11 and the output shaft via two machined couplings secured with screws. Preferably, a top coupling has a flat surface for a flat end screw to securely hold the output shaft 12 onto the coupling. During fabrication of the device, before the top surface 2 of the housing 1 is put in place, an incremental optical encoder 24 (position sensor) is placed around the output shaft 12 as shown in FIG. 5. After the top surface 2 is secured on top of the housing 1, the inner rotating ring 26 of the encoder 24 is secured to the output shaft 12 with, for example, a clamp adapter ring positioned between the optical encoder 24 and output shaft 12. In a preferred embodiment, the rotating ring 26 is secured to the output shaft 12 via a plurality of hex-screws with the clamp adapter ring in between the inner rotating ring 26 of the optical encoder 24 and the output shaft 12. A rotation of the output shaft 12 corresponds to the same amount of rotation of the code wheel disc of the encoder 24.

[0163] The incremental optical encoder 24 is a device that converts motion into a series of digital pulses that can be measured. A preferred incremental optical encoder 24 in accordance with the present invention is shown in FIG. 6 and consists of a disc with a series of equally spaced slits or transparent segments 32 on it. Separating these slits or segments 32 are opaque radial lines 38 of equal arc thickness. On either side of the disk are two LED emitters 34 paired with two phototransistors 36. The typical circuitry of an LED/phototransistor pair is shown in FIG. 7.

[0164] Typically the forward voltage drop across the LED, V_AK is rated at around 1V at 20 mA. Since V_S is 5V, the resistor, R₁ is needed to dissipate 4V. Thus by Ohm's law,

R₁=4V/20 mA=200Ω

[0165] When there is no opaque object between the emitter and detector, the light from the emitter supplies a base current to the transistor. By Kirchoff's Voltage Law,

V_C=V_S=5V

[0166] The saturated collector emitter voltage, V_CE is rated at 0.2V, thus

V_out=5V−0.2V=4.8V

[0167] If R₂=1 MΩ,

I_E=V_out/R₂=4.8V/1 MΩ=4.8 μA

[0168] By Kirchoff's Current Law, inside the transistor the emitter current is equal to the sum of the base and collector currents,

I_E=I_C+I_B

[0169] Also in the transistor,

I_C=h*I_B

[0170] where h is some positive number. In summary, when there is no opaque object between the LED emitter 34 and the phototransistor 36, there is some finite base current, IB from the LED emitter 34, and V_out has a reading of almost 5V (TTL high). When the opaque radial line blocks the light between the LED emitter 34 and the phototransistor 36, I_B≈0, and I_C≈0, thus V_out≈0V. As the disc 30 rotates between the transparent slits 32 and the radial lines 38, the output voltage shifts between TTL high and TTL low respectively, resulting in a square wave output. Typically, the number of square wave cycles (ncypr) generated by an LED/phototransistor pair 34, 36 is equal to the number of radial lines 38 or transparent slits 32 on the disc 30. However if more sophisticated circuitry and interpolation techniques is used, the number of square wave cycles may be greater than the number of radial lines 38 on the disc. The preferred incremental optical encoder used in the present device is the Gurley hollow shaft encoder, which has 11,250 radial lines on the disc and generates 90,000 square wave cycles per revolution.

[0171] The more radial lines 38 (square wave cycles) that there are on the disc 30, the higher the count resolution. The simplest way to obtain angular position from a pair of LED emitters 34 and phototransistors 36 is to count the number of cycles that go by, either by counting the rising or falling edges of one square wave output (CH A) as depicted in FIG. 8, (c) and (d). In this case, the angular resolution is equal to the number of cycles per revolution [ncypr] counts/rev or [360/ncypr] degrees. As shown in FIG. 8, (e) and (f), both the rising and falling edges of CHA are counted, the resolution can be doubled to [2*ncypr] counts/rev or [720/ncypr] degrees.

[0172] If only one channel of square wave output is used, it is impossible to tell the direction of rotation of the disc 30 as the counting pattern of the clockwise and counterclockwise direction are the same (see FIGS. 9(c) and 9(d)). To overcome this problem, an additional LED/phototransistor pair 34, 36 is physically placed around the disk 30 such that it's square wave output (CH B) is a quarter cycle out of phase with that of CH A. Thus, the square wave cycles from both channels are called quadrature signals. When the disk 30 rotates in one direction (i.e. clockwise), CH A leads CH B by a quarter of a cycle, but when the disk 30 rotates in the opposite direction, CH A lags CH B by a quarter of a cycle. Another advantage of using quadrature signals is that the count resolution of the optical encoder 24 can be increased to [4*ncypr] counts/rev by counting the rising and falling edges of both CH S A and CH B. This is done using a quadrature decoder chip (for example, the US-Digital PC-6-84-4). As shown in FIG. 8, (g) and 5(h), the count patterns in the clockwise and counterclockwise directions are different—they are reversed. The quadrature decoder converts the signals from CH A and CH B into a direction (up/down pulse) and count pulses, where a count pulse is generated every time the TTL state of either channel changes. Since the Gurley hollow shaft encoder, preferably used in the present device, generates 90,000 cycles per revolution in each channel, after quadrature edge detection, the encoder gives 360,000 counts per revolution, which is an amazing resolution of one thousandth (10⁻³) of a degree.

[0173] In the setup of the present invention, the outputs of the optical encoder 24 split and go to both the motion controller for feedback control and to a National Instruments data acquisition card for acquiring data. Theoretically, position information can be obtained by interrogating the motion controller via software, in which case, data acquisition of the position signal is not necessary. However, when the motion controller is interrogated at a high sampling rate, its operation and control of the motor shaft 11 is impeded. The controller software cannot both operate the feedback loop to follow a desired trajectory and output data to the user at 500 Hz.

[0174] The optical encoder 24 outputs three additional signals relative to ground, not just A and B. These outputs are the index signals, {overscore (A)}, and {overscore (B)}. Once activated, the optical encoder 24 keeps track of angular position, in counts, relative to the initial angular position where it started counting. The index signal provides an absolute reference so that the absolute position of the hand rest 4 is known in space. In addition to the many equally spaced slits on the disk 30, there is one additional slit that is at a different radial distance away from its center, where the index signal is at TTL 1 only in one absolute position using an additional LED/phototransistor pair 34, 36 as described above. The counter value determined by channels A and B is reset to zero where the TTL value of the index is 1 and, thus, the index serves as an origin in space. Once this position is found, the system knows where the absolute position of the hand rest 4 is in real space.

[0175] As used herein, the axis of rotation of the output shaft 12 is defined as pointing upward, so that a counter clockwise rotation is positive change in direction, and a clockwise rotation is negative change in direction. The motion controller however, defines a clockwise rotation as positive, and a counterclockwise rotation as negative.

[0176] Preferably, in the device of the present invention, the index signal is not used. Rather, in one preferred embodiment, the reference point used to determine absolute position is the left mechanical stop 26. A right mechanical stop 27 may further be located opposite the left mechanical stop 26. The motor shaft 11 moves counterclockwise at a slow speed until the forked wrench of the hand rest 4 hits the left mechanical stop 26 at which point the motor stops. In this configuration, the mid-axis of the forked wrench is 90 degrees (90940 counts) away from the midline, the counter value with the motion controller is set to −90940 which is −90.940 degrees.

[0177] Now the counter value reads, and will continue to read, the correct absolute position until the computer that powers the optical encoder 24 is turned off.

[0178] The output pulse width determines the minimum separation between consecutive A and B pulses that can be detected by the decoder chip. If the separation between consecutive A and B pulses is less than the output pulse width of the decoder, then two output pulses overlap and only one count is registered instead of two, resulting in an erroneous position reading. The width of the clock pulse is influenced by the value of the resistor between pins 1 and 3 of the LS7084 quadrature decoder chip embedded in the US-Digital PC-6-84-4 as shown in FIG. 127. In an exemplary embodiment, an added resister across pins 1 and 3 is used to decrease the decoder pulse width to less than 750 ns.

[0179] An actuator is a device that converts electromechanical force into motion. This electromechanical force can be explained by Lorentz's Law:

{right arrow over (F)}={right arrow over (I)}×{right arrow over (B)}

[0180] When a current in a conductor, {right arrow over (I)}, is exposed to a magnetic field, {right arrow over (B)}, a force {right arrow over (F)} is produced in the direction perpendicular to both {right arrow over (I)} and {right arrow over (B)} as described by the right hand rule.

[0181] Since the hand is perturbed about the wrist joint during testing, a motor 40 was used as the actuator for this task. Generally, the stationary outer housing of a motor, the stator, is magnetized either by permanent magnets or by currents flowing through coiled wires wrapped around iron cores. The rotor or rotating shaft is also magnetized either by permanent magnets or winded field coils. The rotor and its windings are collectively known as armature. If the rotor is magnetized via coils, then it has commutator segments where the current is supplied. Brushes provide stationary electrical contact between the power source and the moving commutator and rotor. FIG. 9 shows the theory of operation of a simple two-pole DC motor in view of the long axis of the rotor, wherein “DC” denotes Direct Current. The setup shown consists of a stator that is a permanent magnet and a rotor that is magnetized via coils. In the starting position, (i), the right brush touches commutator segment A and the left brush touches segment B creating a current that flows through the rotor winding. This current magnetizes the rotor into N and S poles as shown. Since like poles repel and opposite poles attract, the rotor rotates in the clockwise direction as shown in (ii) and (iii). In (iv), for a brief instant, the commutators lose contact with the brushes, no current flows through the armature and the rotor is demagnetized. However, the rotor continues to turn clockwise due to momentum. In (v), the commutator contacts switch brushes and the armature current now flows in the opposite direction, so the poles of the rotor are switched. Thus, the motor continues rotating as the armature current constantly switches direction, maintaining torque in the clockwise direction. Like most motors, the motor shown in FIG. 9 is a brush motor. Its stator is a permanent magnet and its rotor armature has a constantly switching current provided by an external current source through the brushes. The brushes of a DC motor have several limitations: brush life, brush residue, maximum speed, and electrical noise. The currents produced in the armature produce a large amount of heat in the rotor that is difficult to dissipate.

[0182] The brushless DC motor operates by the same principle as the brush motor except that the rotor consists of a permanent magnet and the stator is magnetized via rotating fields in the stator. It is like a brush motor turned inside out. Indicative of its name, the brushless motor eliminates the need for brushes and commutators. Brushless DC motors are potentially cleaner, faster, more efficient, less noisy, more reliable, produce no sparks and the rotor does not heat up. These motors are also known as permanent magnet motors. Motors that use feedback for a position or trajectory control applications are called servomotors.

[0183] Preferably, the motor according to the present invention is a brushless DC servomotor, although other types of motors are compatible for use in the present invention. A particularly preferred motor for use in the present invention is the B-202-B-3 1 brushless servomotor from Kollmorgen as its speed range, torque-speed curves and maximum stall torque are acceptable for use in the present invention. The speed torque curve shows the maximum torque the motor can produce at a certain speed using the amplifier of the present invention (see FIG. 10). For the present invention, this is the torque applied by the hand in response to the perturbations given, and, more importantly, the inertias of the hand and the hand rest. As shown in FIG. 10, the maximum possible operating speed of the motor is 3800 RPM. The motor can generate torques up to 3.5 Nm while maintaining a speed fairly close to the maximum value. If the motor is to generate torques above 3.5 N, the maximum speed it can maintain drops very sharply until the motor finally stops moving at the maximum torque of 4.7 Nm. This maximum torque value where the motor shaft is stationary is called the stall torque.

[0184] When the load torque on the motor is equal and opposite to the toque produced by the motor while the shaft is stationary, the motor is said to be stalling. This value is limited by the maximum amount of current that the amplifier can supply to the motor. What is considered a fast rotation in most industrial applications is at least an order of magnitude larger than what one would consider a rapid angular joint velocity. In industrial applications, motors may move on the order of 10³ RPM. In the present invention, the wrist joint is not expected to be moved to speeds above 25 rads/s or 238 RPM. As shown on the torque-speed curve, the motor can produce close to the maximal stall torque through the entire range of speeds that the present invention is operating at and, thus, the torque-speed relation is not of concern. Rather the concern is that the maximum toque of 4.8 Nm is sufficient for the purpose of the present invention. To determine this, the magnitudes of torques that the motor has to oppose must be analyzed. By Newton's second law, 10$\sum \tau ={\tau }_m+{\tau }_1+{\tau }_f=J_t\ddot{\theta }$

[0185] Rearranging to: 11${\tau }_m=J_t\ddot{\theta }-{\tau }_1-{\tau }_f$

[0186] wherein τ_m is the torque generated by the motor, τ_l is the load torque applied by the passive and active viscoelastic properties of the muscles and τ_f is the negligible friction torque of the mechanical apperatus. J_t is the total inertia of the mechanical apperatus (arm rest) J_m plus the inertia of the human hand, J_l. The signs of τ_l and τ_f are opposite that of τ_m since they oppose the motor torque and rotation direction, thus: 12${\tau }_m=J_t\ddot{\theta }-{\tau }_1-{\tau }_f$

[0187] The majority of the torque that the motor is required to produce is to overcome the total inertia of the system, J_t. The inertia of the typical hand was crudely estimated to be around 0.0015 kg m². The inertia of the armrest is around 0.005 kg m². J_t=J_m+J_l≈0.0065 kg m². Since it is not expected that the wrist joint will be perturbed at accelerations that exceed ±500 rads/s², the maximum torque needed to overcome inertia is 13$J_t{\ddot{\theta }}_{\max }=\left(0.0065{\mathrm{kgm}}^2\right)\left(500\mathrm{rads}\text{/}s^2\right)=3.25\mathrm{Nm}$

[0188] With the load torque applied in response to the perturbations, the motor selected is sufficient for the present invention because the maximum possible torque that the motor can produce is not expected to be exceeded. Further, based on these specifications, other suitable motors may be readily determined by one of skill in the art. The torque capibility of the motor is expressed in terms of the peak stall torque and the continuous stall torque. The peak stall torque is the absolute maximum torque that the motor can produce. This peak torque can only be maintained for a few seconds before the torque settles down to to the continuous stall torque.

[0189] The back emf, V_emf is an induced voltage by the rotating armature that opposes the applied voltage that is determined by: 14$V_{\mathrm{emf}}=K_e\stackrel{.}{\theta }\mathrm{or}V_{\mathrm{emf}}\left(s\right)=K_e\left(s\right)\stackrel{.}{\Theta }$

[0190] in the Laplace domain. [eq. 8]

[0191] K_e is the emf constant of the motor and θ is the angular velocity of the rotor. If the power supplied to the motor is a voltage source and the voltage supplied at a particular time is V_m, then: 15$\begin{array}{cc}{V}_m=L{\stackrel{.}{I}}_m+{\mathrm{RI}}_m+V_{\mathrm{emf}}=L{\stackrel{.}{I}}_m+{\mathrm{RI}}_m+K_e\stackrel{.}{\theta }\mathrm{or}V_m\left(s\right)=\left(\mathrm{Ls}+R\right)I_m\left(s\right)+V_{\mathrm{emf}}\left(s\right)\mathrm{in}\mathrm{the}\mathrm{Laplace}\mathrm{domain}& \left[\mathrm{eq}.9\right]\end{array}$

[0192] I_m is the resulting armature current that can be derived in terms of V_m, L and R. L is the inductance, and R is the motor resistance. I is the rate of change of the armature current. The torque produced by the motor is proportional to the current supplied to it:

τ_m=K_tI_m

[0193] or T_m(s)=K_tI_m(s) in the Laplace domain [eq. 10]

[0194] where K_t is the torque constant of the motor.

[0195] Equation 10 can be rewritten as: 16${\tau }_m+{\tau }_1=B_m\stackrel{.}{\theta }+J_t\ddot{\theta }$ or T_m(s)+T_l(s)=(B_ms+J_ts²)Θ(s) in the Laplace domain. [eq.11]

[0196] The amplifier is the variable power supply to the motor. It amplifies the command voltage from the controller into either a voltage or a current to the motor 40. The input to the amplifier is the command voltage given by the Galil motion controller. There are three modes of operation of an amplifier: voltage source mode, current source mode, and velocity loop mode, all of which are described below:

[0197] Voltage Source Mode:

[0198] In this mode, the output to the motor is a voltage, V_m, that is some unitless gain, K_V, times the command voltage, V_cmd.

V_m=K_VV_cmd [eq. 12]

[0199] The resulting armature current can be derived from equation 9. The result is quite cumbersome since equation 9 is a first order differential equation in terms of I_m and also depends on the time varying velocity and the back emf constant, as well as the impedance and resistance of the motor. The resulting motor toque, τ_m is simply proportional to the current I_m by the torque constant. It is not easy to derive the torque output of the motor given the amplifier output of V_m because of the complexity of equation 9. The transfer function that is the ratio of motor position to voltage command V_cmd is: 17$\begin{array}{cc}\frac{\Theta \left(s\right)}{V_{\mathrm{cmd}}\left(s\right)}=\frac{K_vK_t}{J_t{\mathrm{Ls}}^3+\left(B_mL_m+J_tR_m\right)s^2+\left(B_mR_m+K_tK_e\right)}& \left[\mathrm{eq}.13\right]\end{array}$

[0200] Current Source Mode:

[0201] The relationship between the motor voltage, V_m and motor torque, τ_m is a complicated differential equation with respect to many variables. In current source mode, advantage can be taken of the simple linear relationship between the motor current and torque output as described by the torque constant, K_t as shown in equation 10. Here the amplifier takes the command voltage, V_cmd from the Galil controller and amplifies it to the current I_m the gain K_a:

I_m=K_aV_cmd [eq. 14]

[0202] The units of K_a are [Amps*RMS/Volt]. Here, the transfer function between the motor position and voltage command is: 18$\begin{array}{cc}\frac{\Theta \left(s\right)}{V_{\mathrm{cmd}}\left(s\right)}=\frac{K_aK_t}{J_ts^2+B_ms}& \left[\mathrm{eq}.15\right]\end{array}$

[0203] Velocity Loop Mode:

[0204] The velocity loop mode is the same as the current mode with an added feedback component. A tachometer that senses velocity converts the velocity into a voltage by the gain, Ku, and sends it back to the amplifier. The units of K_u are [V*s/count]. In this mode the transfer function between the motor position and command voltage is: 19$\begin{array}{cc}\frac{\Theta \left(s\right)}{V_{\mathrm{cmd}}\left(s\right)}=\frac{K_aK_t}{J_ts^2+B_ms+K_uK_aK_t}& \left[\mathrm{eq}.16\right]\end{array}$

[0205] A preferred amplifier for use in the present invention is the Kollmorgen BDS4-1033-0101-202B2, which matches the desired motor characteristics. In an illustrative embodiment, the amplifier is set up in the current mode. This mode makes the analysis of the control system fairly simple. The preferred power supply used for the amplifier is the PSR4/5-112-0102, which can supply a peak current to the above-described motor of about 6 Amps RMS, that lasts for around 2 seconds, and a maximum continuous current of about 3 Amps RMS, as shown by the dashed line in FIG. 12.

[0206] After setting up the system, prior to operating the motor 40, the gain of the amplifier and other adjustments are set via tuning pots located on the chassis of the amplifier. The command scale pot is used to adjust the velocity loop gain, K_u. However, since the velocity loop is disabled, this is not necessary. Pin 18 on Connector C1 of the amplifier is used to monitor the current output relative to common (pin 17). The output is a voltage proportional to the current output where 1V 0.75 Amps RMS. This feature is used to monitor the current supplied by the amplifier. The balance pot is used to make sure that when 0V is commanded from the Galil controller, the current output of the amplifier is 0 Amps RMS and the motor produces no torque. Equation 14 above is really of the form:

I_m=K_aV_cmd+α [eq. 14]

[0207] where α is the offset current that is adjustable with the balance pot. With the motor enabled, the input to the amplifier, V_cmd is set to 0V and the balance pot is rotated such that I_m read 0 Amps to ensure no offset current. The current limit pot is turned all the way clockwise to maximize the maximum amount of current that the amplifier can supply (up to ±6Amps RMS). The most important adjustment is the stability or gain adjustment. This determines the magnitude of K_a. This adjustable amplifier gain is set experimentally so that the maximum possible value of V_cmd (10 V) would correspond to the maximum current limit (approximately 6A) of the amplifier to take advantage of the full range of currents. Thus, the stability pot is adjusted so that value of K_a is set to 0.6 Amps RMS/V.

[0208] An experiment was conducted where the voltage output from the Galil controller was set from 0V to 10V at 0.5V increments. To measure the torque output of the motor, τ_m, the motor shaft was kept stationary using the mechanical stops. In this state, 20$\stackrel{.}{\theta }=0\mathrm{and}\ddot{\theta }=0$

[0209] and

τ_m=−τ_m from equation 11.

[0210] This torque was measured with the torque transducer. The motor torque and amplifier current was measured at each command voltage to produce the plots shown in FIG. 12.

[0211] For each command voltage, the correspontind amplifier currents and motor torques were plotted against each other and K_t was determined by a linear fit through the origin as shown in FIG. 13.

[0212] The maximum torque the motor can produce is limited by the maximum current supplied by the amplifier as defined by equation 10. Thus the peak transient torque the motor can produce is 4.86 Nm, which lasts for about 2 seconds before it drops to about 3 Nm.

[0213] The controller is the brain of the feedback system, which consists of hardware and software. It controls the command voltage, V_cmd that is sent to the amplifier that controls the output torque of the motor shaft since τ_m=K_tK_aV_cmd and K_a and K_t are known constants. The controller has two inputs at any given discrete time instant, t_i: the desired position of the motor shaft, θ_d(t_k) and the actual position of the motor shaft, θ(t_i). The encoder 24 gives the actual position of the motor shaft 11 with a resolution accuracy of 0.001 degrees, as there are 360,000 counts in one revolution of the motor shaft 11. The controller deals with positions in [counts] since that is the digital output of the encoder/decoder. The command voltage it sends to the amplifier is some function of the position error, θ_e where:

θ_e(t_i)=θ_d(t_i)−θ(t_i)

[0214] The desired trajectory, θ(t_k) is given as a set of digital points dt_mc seconds apart. This is the period of one cycle of the feed back loop. A new command voltage is sent to the amplifier every dt_mc seconds. Thus the frequency of the desired trajectory fed into the controller in Hz is 21${\mathrm{fhz}}_{mc}=\frac{1}{{\mathrm{dt}}_{mc}}$

[0215] The points in the desired trajectory are either generated offline or online. If they are generated offline, the exact points of the trajectory have to predetermined and fed into the controller as an array. This is what is done, for example, when the desired trajectories are the filtered random trajectories that are generated with MATLAB. Since the desired trajectory is periodic, only one period of the trajectory needs to be inputted to the controller. If the points in the desired trajectory are dependent on some output state such as torque or continuously change with time, then they have to be generated online within the loop. When the desired trajectory is a sinusoid, each point of the desired trajectory is generated online based on the time elapsed from the start and the user defined amplitude and frequency. If the points are to be generated within the loop, then everything from sensing the necessary output states to computing the next desired position or torque has to take place in less that dt_mc seconds. This is done so the sampling frequency of the controller is maintained at 128 Hz, otherwise the frequency of the controller needs to be decreased to increase the amount of time per loop. At substantially low sampling frequencies, the trajectory will become “choppy” and will not be smooth as defined earlier.

[0216] Preferred components of the controller selected for the present invention can be obtained from Galil Motion Control. Several high-level motion controllers can be considered. Preferably, the controllers can be interfaced and, thus, controlled by LabVIEW software, which is the preferred code used in the present device. Preferably, the device is entirely automated and, thus, intervention of the user is minimal as he/she only gives certain desired inputs such is file identifier and the desired trajectories for the device to follow. The components of the Galil controller include the following:

[0217] the DMC-1410 single-axis motion controller for the ISA bus.

[0218] the ICM-1460 interconnect module.

[0219] the WSDK-32 Servo Design Kit for Windows 95, 98

[0220] the Setup-32 software for Windows 95, 98

[0221] a 37 pin connector cable from the DMC-1410 to the ICM-1460

[0222] The controller uses a high level code, which is attached hereto as Appendix A. Of course, controllers having similar components and that function as required by the present invention could also be used.

[0223] FIG. 15 shows one preferred embodiment of a block diagram of components of the PD closed loop control system used to follow the desired trajectories.

[0224] As shown in FIG. 15, a PD Compensator is utilized. A preferred PD Compensator for use with the present invention has the following specifications:

[0225] Manufacturer: Galil Motion Control

[0226] Input: position error, θ_e [counts]

[0227] Output: digital number, N_d from −32767 to 32767

[0228] Description: motion control software downloaded to PLC takes desired trajectory input and feedback from encoder to control motor output.

[0229] Transfer Function: 22$\frac{N_d}{{\theta }_e}=\left(K_P+K_DS\right)$

[0230] Constants: K_P=0.75 (chosen) K_D=0.1563 s⁻¹ (chosen)

[0231] Further shown in FIG. 15 is a Digital to Analog converter (D/A). A preferred D/A for use with the present invention has the following specifications:

[0232] Manufacturer: Galil Motion Control

[0233] Input: digital number, N_d from −32767 to 32767.

[0234] Output: command voltage, V_cmd from −10 V to 10 V.

[0235] Description: converts the digital number from −32767 to 32767 to an analog voltage output from −10 V to 10 V.

[0236] Transfer Function: 23$\frac{V_{\mathrm{cmd}}}{{\theta }_e}=G_{\mathrm{dac}}$

[0237] Constants: G_dac=20/65535 (fixed, set by PWM & hardware)

[0238] Further shown in FIG. 15 is an amplifier. A preferred amplifier for use with the present invention has the following specifications:

[0239] Manufacturer: Kollmorgen Industrial Drives

[0240] Input: command voltage, V_cmd from −10 V to 10 V.

[0241] Output: motor armature current, i_m, from −6.16 Amps RMS to 6.16 Amps RMS.

[0242] Description: produces a current proportional to the input voltage. The adjustable amplifier gain was set experimentally so that the maximum possible value of V_cmd would correspond to the maximum current limit of the amplifier to take advantage of the full range of currents.

[0243] Transfer Function: 24$\frac{I_m}{V_{\mathrm{cmd}}}=K_a$

[0244] Constants: K_a=0.603 Amps RMS/V (experimentally set)

[0245] Further shown in FIG. 15 is a servo motor. A preferred servo motor for use with the present invention has the following specifications:

[0246] Manufacturer: Kollmorgen Industrial Drives

[0247] Input: motor armature current, i_m, from −6.16 Amps RMS to 6.16 Amps RMS.

[0248] Output: torque produced by motor from −5.19 Nm to 5.19 Nm.

[0249] Description: produces a torque proportional to the input current.

[0250] Transfer Function: 25$\frac{T_m}{I_m}=K_t$

[0251] Constants: K_t=0.843 Nm/Amps RMS (fixed, experimentally determined)

[0252] Further shown in FIG. 15 is a plant. A preferred plant for use with the present invention has the following specifications:

[0253] Input: motor torque, T_m plus external load torque, T_l applied by human [Nm].

[0254] Output: actual motor shaft position, 0 [counts].

[0255] Description: the resulting motor shaft position is determined by the input torques, combined viscosities and inertia of the motor shaft, all mechanical loads attached to the motor shaft including the arm rest B_m and J_m. The position is also affected by the inertia of the hand, J_l if attached.

[0256] Transfer Function: 26$\frac{\theta }{\left(T_l+T_m\right)}=\frac{1}{B_ms+\left(J_l+J_m\right)s^2}·\frac{4N_{\mathrm{cypr}}}{2\pi }$

[0257] Note: all units in the control system are expressed in SI units except for the angular positions which are not in radians but in counts. (4N_cypr counts=2π radians)

[0258] Constants: B_m=0.003 Nm s/rad experimentally determined,

[0259] J_m≈0.0048 kg m². The exact value slightly varies depending on the position of the armrest on the sliding bar.

[0260] J_l≈0.0023 kg m². The exact value varies depending on the mass distribution, and geometry of the hand of the subject being tested.

[0261] N_cypr=90,000. This is the number of square wave cycles per revolution of the incremental encoder. With quadrature edge detection, there are 4N_cypr counts per revolution.

[0262] The entire closed loop PD system preferably has the following specifications:

[0263] Input: desired motor shaft positions (trajectory), θ_d [counts].

[0264] Output: actual motor shaft position, θ [counts].

[0265] Description: the system implemented is a Linear Time Invariant second order model of the system.

[0266] Transfer Function: 27$\begin{array}{c}\frac{\theta }{{\theta }_d}=\frac{\left(K_DG_{\mathrm{dac}}K_aK_t\right)s+K_PG_{\mathrm{dac}}K_aK_t}{\frac{\left(J_l+J_m\right)\pi }{2N_{\mathrm{cypr}}}s^2+\left(K_DG_{\mathrm{dac}}K_aK_t+\frac{B_m\pi }{2N_{\mathrm{cypr}}}\right)s+K_PG_{\mathrm{dac}}K_aK_t}\\ =\frac{24.256s+465.705}{0.1239s^2+24.308s+465.705}\end{array}$

[0267] It is important that the system characteristic equation or denominator of R(s), gives poles at: 28$\frac{-\left(K_DG_{\mathrm{dac}}K_aK_t+\frac{B_m\pi }{2N_{\mathrm{cypr}}}\right)\pm \sqrt{{\left(K_DG_{\mathrm{dac}}K_aK_t+\frac{B_m\pi }{2N_{\mathrm{cypr}}}\right)}^2-\frac{2\pi \left(J_l+J_m\right)\left(K_PG_{\mathrm{dac}}K_aK_t\right)}{N_{\mathrm{cypr}}}}}{\frac{\left(J_l+J_m\right)\pi }{N_{\mathrm{cypr}}}}$

[0268] −87.56 and −10.73 after substitution and simplification. Since the roots have negative real parts, it can be said that the linear time-invariant system is bounded-input, bounded-output and asymptotically stable. The fact that roots are unequal and have no imaginary parts says that the system is overdamped (ζ>1). In fact, for the inertia constants used above, ζ=1.6. For smaller inertias, ζ increases slightly and for larger hand inertias, ζ decreases slightly.

[0269] Using typical values for constants J_m, and J_l shown above and substituting jω for s, the analytical magnitude and phase of the frequency response is shown in FIG. 16. The analytical step response is also shown in FIG. 16.

[0270] The tuning parameters, K_p and K_d, are chosen so that the system will have a fast response time without making them too large. The torque saturates when small position and velocity error gives rise to a large digital output N_d and a command voltage, V_cmd close to +/−10V. Preferably, from 0 to 8 Hz, the magnitude of frequency remains close to 1. The actual trajectory of the shaft was checked for all types of desired trajectories both analytically using the MATLAB function mypid.m (see Appendix A) as well as experimentally by measuring the actual output position of the shaft.

[0271] A preferred test setup and procedure in accordance with the present invention is as follows. The patient being tested with the device of the present invention is preferably seated in front of the device. Two surface electrodes are placed on the patient's arm to monitor activity of the flexors and extensors. The ground electrode is placed on an inactive or neutral site. After the absolute position of the motor shaft is calibrated, the seat is adjusted such that the patent's left or right arm can be placed in the arm rest 5. At this point, the patient's elbow is preferably at about 145 to 160 degrees. The patient's hand is placed on the hand rest 4 in the neutral supination/pronation position such that his or her knuckles are slightly passed between the padded vertical bars 6. The vertical bars 7 are moved closer together and tightened such that the hand is comfortably secured in place. The hand rest 4 is then slid forward or backward so that the wrist joint is lined up with the motor shaft 11. Now the subject is simply told to relax and not to resist the robot's movement as the robot moves to the initial bias position and follows the desired trajectories as described above. The robot automatically cycles through the desired trajectories (trials 1-10) while collecting torque, position and EMG data. If an additional test is to be performed in a different bias position, the patient's arm is preferably removed from the apparatus for 2 or 3 minutes before the procedure is repeated at the new initial bias position.

[0272] The bias position is the absolute initial angular position of the wrist joint that the motor shaft 11 moves to before each perturbation trajectory or trial is carried out. A bias position of 0 degrees corresponds to the neutral flexion/extension position. It is also the median angular position of the trajectory. The user selects the bias or initial position, and the motor shaft 11 automatically moves to that position before the desired trajectory is carried out. In testing control subjects, selecting a 0-degree bias position is not a problem as control subjects have a wide range of motion. Almost all patients with spasticity and others with neurological problems, however, have a limited range of motion. Their resting position is in a hyper-flexed position and it is painful or impossible for them to extend their hand beyond the 0-degree absolute position due to the contracted state of the flexor muscles. Due to this constraint, subjects in the patient population were tested in the 30 degree flexed position. Control subjects were tested at both 0 and 30 degree flexed bias position. Since the left or right arm may be tested, the four possible combinations of bias positions are shown in FIG. 17.

[0273] For more detailed instructions about the procedure, please refer to the user's manual attached hereto as Appendix B.

[0274] Next, the raw data is taken and transformed into the torque offset, elastic stiffness, viscosity and inertia parameters described above and represented by the variables τ_off, K_H, C_H, and I_T respectively. τ_off is the offset torque, K_H is the combined angular elastic stiffness of all the muscles acting on the wrist joint, C_H is the combined angular viscosity of all the muscles acting on the wrist joint, and I_T is the moment of inertia of the moving hand as well as the moving mechanical parts of the apparatus (hand rest).

[0275] Data obtained from subject JF (a 28 year old, male control subject) is identified as experiment DMJF. In this experiment, the 10 trial random trajectories were run on the right hand at a bias of 0 degrees. The entire data analysis procedure will be described from obtaining the raw data until the computation of impendence parameters in equation 1.

[0276] As mentioned above, for each perturbation trajectory given, five measurements are collected and stored using LabVIEW and a National instrument PCI-E series data acquisition board at 500 Hz. A reading from each sensor is synchronously collected every 2 ms. These raw data points are stored as an n×5 array in ASCII format. Data is collected for 20 seconds (n=10,000) per trial. The first data column is just a timestamp of when data from each sensor is collected (i.e. 0 ms, 2 ms, 4 ms . . . ). The second data column is a set of integer counts from the encoder 24, which make up the angular position signal of the motor shaft 11 and the wrist joint. The third data column is a list of torque readings from the torque transducer in Nm. The fourth and fifth data columns are flexor and extensor EMG readings respectively, in mV. Since each experiment consists of ten filtered random trajectories, there are ten arrays of raw data stored as DMJFI1t1, DMJFI1t2, . . . , DMJFI1t10.

[0277] The raw signals of each of the 10 trials in experiment DMJF are shown in FIGS. 18-27.

[0278] In one preferred embodiment, the device utilizes a filter. One preferred filter is a Savitsky-Golay filter (Savitsky & Golay, Press et. al: Numerical Recipes). A filter can provide a number of advantages. The discretization of position gives a finite resolution of position based on the resolution of the encoder. In the present invention, this resolution is 0.001 degrees. Filtering enables displacement positions to fall between these increments. Since the resolution is very high, this is particularly important when using low-resolution encoders, with a low ncypr. Further, filters can assist in the timing of data acquisition. Still further, filters are useful in obtaining smooth first and higher order derivatives of a discrete signal. A crude method to obtain the derivative (slope or tangent) at a point in a discrete signal is to take the difference between neighboring points and then divide it by the sampling period. This is based on the assumption that the two neighboring points are sufficiently close enough such that the points in between lie on a straight-line between the two. Sharp discontinuities in the derivative signal are observed using this technique. As in most smooth signals, the present displacement signal between any two points is actually a curve. Thus, the filter utilizes a better technique to obtain derivatives by fitting a polynomial through a window of neighboring points. The Savitsky-Golay filter, for example, is a smoothing filter that is useful when determining derivatives such as velocity and acceleration offline when the raw displacement signal has already been collected and stored. This smoothing technique fits a polynomial to a moving window of data using least squares fitting. Derivatives, as well as the function value, can be determined at the center point of each window. Since the coefficients of the fitted polynomial are linear in the data values, the filter coefficients can be pre-computed. The filter coefficients are then convoluted with the raw digital signal to obtain the actual smoothed version of the raw waveform as well as well as its derivatives at each center point.

[0279] The Savitsky-Golay algorithm operates as follows with the following variables:

[0280] d_raw: the column vector of actual displacement data points

[0281] N: total number of points in the raw digital displacement signal (rads)

[0282] ord: the order of the polynomial.

[0283] wins: the number of equally spaced points to be fit by the polynomial, the size of the window. This number is always odd.

[0284] y_i: a sub vector of d_raw of length wins with i being the index of the center point.

[0285] hw=(wins−1)/2: half width of the window.

[0286] The matrix X of dimension wins by (ord+1) is taken as 29$X_{\mathrm{winsXord}+1}=\begin{array}{ccccc}-{\mathrm{hw}}^0& -{\mathrm{hw}}^1& \cdots & -{\mathrm{hw}}^{\mathrm{ord}-1}& -{\mathrm{hw}}^{\mathrm{ord}}\\ -{\left(\mathrm{hw}-1\right)}^0& -{\left(\mathrm{hw}-1\right)}^1& \cdots & -{\left(\mathrm{hw}-1\right)}^{\mathrm{ord}-1}& -{\left(\mathrm{hw}-1\right)}^{\mathrm{ord}}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& 0& 0\\ {\left(\mathrm{hw}-1\right)}^0& {\left(\mathrm{hw}-1\right)}^1& \cdots & {\left(\mathrm{hw}-1\right)}^{\mathrm{ord}-1}& {\left(\mathrm{hw}-1\right)}^{\mathrm{ord}}\\ -{\mathrm{hw}}^0& -{\mathrm{hw}}^1& \cdots & -{\mathrm{hw}}^{\mathrm{ord}-1}& -{\mathrm{hw}}^{\mathrm{ord}}\end{array}$

[0287] The sum of the squared error, S, is minimized between the elements of y_i and y_i, an estimate of y_i:

y_i=X*F*y_i

[0288] where F is a (ord+1) by wins matrix. The sum of the squared errors can be written as:

S=(y_i−y_i)^T*(y_i−y_i)=(y_i−X*F*y_i)^T(y_i−X*F*y_i)

[0289] A matrix F is evaluated such that the least possible sum-squared error, S, is obtained. To do this, the derivative of S is set with respect to F to 0, dS/dF=0. This yields the pseudo-inverse solution,

F=(X^T*X)⁻¹*X^T

[0290] The columns of F are then reversed and the k^th derivative is then given by k! multiplied with the (k+1)^th row of F divided by (dt_daq)^k convolved with the raw displacement data. Before applying the Savitsky-Golay filter, the encoder count signal is converted into radians. Since the resolution of the encoder gives us 1000 counts for a degree of angle, the raw encoder count signal is multiplied by 30$\frac{\pi }{180000}$

[0291] to give the raw displacement data in radians in the column vector, rawrads. In equation form, the i^th element of the k^th derivative using the Savitsky-Golay filter is given as: 31$\frac{{}^k\left(rawrads\right)}{t^k}\left(i\right)=\left(\frac{1}{\mathrm{t\_}aq^k}\right)*k!*\sum _{j=-hw}^{hw}\left(rawrads_i*F_{k+1,i+1-j}\right)$

[0292] A third order polynomial (ords=3) is used with a window size of 21 points. Thus, F is a 4 by 21 matrix. From the equation above, the i^th element of the filtered displacement signal is given as: 32$disp\left(i\right)=\sum _{j=-10}^{10}\left(rawrads_i*F_{1,i+1-j}\right)$

[0293] the velocity signal is given as: 33$\mathrm{vel}\left(i\right)=\frac{\left(rawrads\right)}{t}\left(i\right)=\left(\frac{1}{0.002}\right)*\sum _{j=-10}^{10}\left(rawrads_i*F_{2,i+1-j}\right)$

[0294] and the acceleration signal is given as: 34$\mathrm{acc}\left(i\right)=\frac{{}^2\left(rawrads\right)}{t^2}\left(i\right)=\left(\frac{1}{4*{10}^{-6}}\right)*2*\sum _{j=-10}^{10}\left(rawrads_i*F_{3,i+1-j}\right).$

[0295] Since the units of the elements in the raw rads vector are in rads or radians, elements in the displacement, velocity, and acceleration array are in rads, rads/s, and rads/s² respectively. The beginning and the end of each vector are truncated by a number of points equal to the half width of the window. The frequency response of the Savitsky-Golay filter using a third order polynomial with a window size of 21 points (0.04 s) in the frequency domain can be seen using the freqz MATLAB function which shows the Z-transform digital filter frequency response as follows:

[0296] >>sav=sav_golay(3,21);

[0297] >>[h,f]=freqz(sav(1,:),1);

[0298] >>plot(f*500/(2*pi),abs(h),‘k’);

[0299] This produces the frequency response plot shown in FIG. 28. There is no phase lag since the window is symmetric with respect to the point at which the polynomial is evaluated (center point).

[0300] Unlike digital signals, such as that from the encoder 24, analog signals from torque transducers have high frequency noise components that do not give a true representation of the actual torque exerted. Low pass digital filters are often used offline to attenuate the high frequencies in the torque signal. A zero lag 8^th order digital Butterworth filter can be used to greatly attenuate these high frequencies.

[0301] In MATLAB, [b,a]=butter(n, Wn) designs an order n lowpass digital Butterworth filter with cutoff frequency Wn. It returns the filter coefficients of length n+1 in row vectors b and a, with coefficients in descending powers of z: 35$H\left(z\right)=\frac{B\left(z\right)}{A\left(z\right)}=\frac{b\left(1\right)+b\left(2\right)z^{-1}+\dots +b\left(n+1\right)z^{-n}}{1+a\left(2\right)z^{-1}+\dots +a\left(n+1\right)z^{-n}}$

[0302] The cutoff frequency is that frequency where the magnitude response of the filter is {square root}{square root over ({fraction (1/2)})}. For example, for MATLAB's butter function, the cutoff frequency, Wn must be a number between 0 and 1, where 1 corresponds to half the sampling frequency (the Nyquist frequency). A preferred embodiment of the present invention uses a Wn of 0.075, which corresponds to a cutoff frequency of 18.75 Hz (see getstatesgbw.m in Appendix A). Normally, Butterworth filters have a phase response where there is a lag between the filtered and the raw signal as a function of frequency of the raw signal. In other words the filtered torque signal is nonlinearly shifted in time with respect to the original signal as a function of frequency. Thus, the filtered torque signal also lags behind the encoder signal. These phase lags are particularly significant in high order filters. This lag is problematic because it is important for the torque waveform to be synchronized with the position velocity and acceleration waveforms in order to evaluate the parameters in the viscoelastic model described herein. However, the filtfilt function in MATLAB is prefereably used. y=filtfilt(b,a,x) performs zero-phase digital filtering by processing the input data in both the forward and reverse directions.

[0303] After filtering in the forward direction, the filtered sequence is reversed and it runs it back through the filter. The resulting sequence has precisely zero-phase distortion and double the filter order (A. V. Oppenheim, and R. W. Schafer. Discrete-Time Signal Processing. Englewood Cliffs, N.J.: Prentice Hall, 1989. Pgs. 311-312).

[0304] Stiffness, viscosity, and inertial parameters are obtained from position and torque signals obtained from the optical encoder and torque transducer. Velocity and acceleration signals are also required, as shown in equation 1. These signals are derived inherently from the position signal via the Savitsky-Golay filter described above. Filtering of the flexor and extensor EMG signals is not crucial since these signals are not necessary for post processing. However, monitoring live EMG information can be used to monitor possible voluntary contractions. Ideally the muscle should be fully relaxed before perturbing the wrist, since spasticity is being measured, which is an involuntary reflex. Thus, the EMG information can be used live to ensure proper integrity of the data collected. Since flexor and extensor EMG activity is monitored, the EMG data can also be collected and stored. Post processing on EMG signals can be useful to answer more basic questions, such as the investigating reflex delay times from the onset of single or repetitive muscle perturbations to muscle contraction. Thus, the total time it takes for the entire reflex loop can be measured. Researchers have shown that two separate response times can be measured, one due to spinal pathways to the spinal cord and back, and one due to supra-spinal pathways which involve upper motor neurons from the brain (See J. C. Houk: Participation of Reflex Mechanisms and Reaction-Time Processes in the Compensatory Adjustments to Mechanical Disturbances. Cerebral Motor Control in Man: Long Loop mechanisms. Prog. clin. Neurophysiol., vol 4, Ed. J. E. Desmedt, pp 193-215, Krager Basel, 1978; W. G. Tafton, P. Bawa, I. C. Bruce, R. G. Lee: Long Loop Reflexes in Monkeys: An Interpretative Base for Human Reflexes. Cerebral Motor Control in Man: Long Loop mechanisms. Prog. clin. Neurophysiol., vol 4, Ed. J. E. Desmedt, pp 229-245, Krager Basel, 1978). Naturally, the supra-spinal response occurs after the spinal response.

[0305] EMG signals are extremely noisy, particularly if surface electrodes are used. EMG signals measured with surface electrodes are influenced by hair, oil, lotions and dead skin. The actual action potentials sent to the muscles are internally amplified by 1000, typically resulting in a signal from 0 to ±5 V. Preferably, internal to the Delsys 2-channel EMG unit is a 20-450 Hz band pass filter. Once the data is collected through the data acquisition board and stored, further filtering is then performed offline in MATLAB as follows.

[0306] First a moving root-mean-square (RMS) window is applied to the raw output of the EMG signal where the size of the window is 15 points or 28 milliseconds. Thus, the k^th sample point is of the RMS EMG signal is equal to 36$EMG_{RMS}\left(k\right)=\sqrt{\frac{1}{15}\sum _{i=k-7}^{k+7}{\left(EMG_{RAW}\left(k\right)\right)}^2}$

[0307] After the RMS EMG signal is obtained, it is passed through a low pass butterworth filter as applied to the torque signal. The result gives the final filtered EMG signal. The cutoff frequency used in this butterworth filter is 37.5 Hz.

[0308] The displacement, velocity and acceleration signals using the Savitsky-Golay filter on the encoder signal for each of the 10 trials are shown in FIGS. 30-29. The filtered torque and EMG signals are also shown. Filtered data obtained from two cycles or periods of each random displacement signal are used for further processing.

[0309] Now all the proper state measurements 37$\left({\tau }_s\left(t_i\right),\theta \left(t_i\right),\stackrel{.}{\theta }\left(t_i\right),\mathrm{and}\ddot{\theta }\left(t_i\right)\right)$

[0310] are known after the signal processing described above, and the data can now be fit to the model of the present invention: 38${\tau }_s\left(t_i\right)={\tau }_{\mathrm{off}}+K_H*\theta \left(t_i\right)+B_H*\stackrel{.}{\theta }\left(t_i\right)+J_T*\ddot{\theta }\left(t_i\right)$

[0311] where τ_S is the filtered torque signal from the transducer. K_H and B_H are the angular stiffness and viscosity of the combined flexor and extensor muscle groups that act on the wrist joint. J_T is the combined inertia of oscillating appendage, J_H, and the rotating components of the apparatus, J_A. J_T=J_H+J_A. τ_off, the offset torque is equal to −K_H*θ_RP as described above. θ is the angular displacement of the system. 39$\stackrel{.}{\theta }\mathrm{and}\ddot{\theta }$

[0312] are velocity and acceleration, respectively which are the first and second derivatives of displacement. The data was fit using a linear least squares fit pseudo-inverse (See A. Bjorck, Numerical Methods for Least Squares Problems, SIAM, Philadelphia, 1996; H. Anton, C. Rorres, Elementary Linear Algebra, Applications version, sixth edition, Wiley & Sons, New York, 1991), as follows. For a set of discrete samples, the following matrices are constructed: 40$Ts=|\begin{array}{c}\tau s\left(t_1\right)\\ \tau s\left(t_2\right)\\ \tau s\left(t_3\right)\\ ⋮\\ \tau s\left(t_n\right)\end{array}{|}_{n\times 1}\mathrm{and}\Psi =|\begin{array}{cccc}1& \theta \left(t_1\right)& \stackrel{.}{\theta }\left(t_1\right)& \ddot{\theta }\left(t_1\right)\\ 1& \theta \left(t_2\right)& \stackrel{.}{\theta }\left(t_2\right)& \ddot{\theta }\left(t_2\right)\\ 1& \theta \left(t_3\right)& \stackrel{.}{\theta }\left(t_3\right)& \ddot{\theta }\left(t_3\right)\\ ⋮& ⋮& ⋮& ⋮\\ 1& \theta \left(t_n\right)& \stackrel{.}{\theta }\left(t_n\right)& \ddot{\theta }\left(t_n\right)\end{array}{|}_{n\times 4}$

[0313] It is possible to obtain the set of three parameters, 41$\tilde{N}={\begin{array}{c}\stackrel{{\widetilde{^}}_{\mathrm{off}}}{o}\\ {\tilde{K}}_H\\ {\tilde{B}}_H\\ {\tilde{J}}_T\end{array}}_{4\times 1}$

[0314] that best fits the present model in equation 1: 42${\tau }_S\left(t_i\right)={\tau }_{\mathrm{off}}+K_H*\theta \left(t_i\right)+B_H*\stackrel{.}{\theta }\left(t_i\right)+J_T*\ddot{\theta }\left(t_i\right)$

[0315] using the pseudo-inverse which minimizes the sum squared error, a positive scalar value. In matrix equation form the sum-squared error can be written as:

SSE=(Ô_S−Ø*Ñ)^T*(Ô_S—Ø*Ñ)

[0316] The matrix P that minimizes SSE is evaluated as:

Ñ=pinv(Ø)*Ô_S=(Ø^T*Ø)⁻¹*Ø^T*Ô_S [eq. 15]

[0317] First the pseudo-inverse is performed to obtain one set of parameters utilizing all data over 2 whole periods (15.625 s or n=7815) that construct one large Ô_S (7815×1) vector of torque measurements and one large Ø (7815×4) state matrix. Let us denote the four evaluated parameters as {tilde over (τ)}_off, {tilde over (K)}_H, {tilde over (B)}_H, and {tilde over (J)}_T. These evaluated parameters for the experiment DMJF whose data is shown below. 2

[0318] To determine whether the model used in the present invention to obtain overall parameters for each trial is a reasonable one, the residual torque, τ_R is plotted. The residual torque is the amount of torque that is unaccounted for by the model and for a particular time sample. In equation form, it is defined as: 43${\tau }_R\left(t_i\right)={\tau }_s\left(t_i\right)-\left[{\tilde{\tau }}_{\mathrm{off}}+{\tilde{K}}_H*\theta \left(t_i\right)+{\tilde{B}}_H*\stackrel{.}{\theta }\left(t_i\right)+{\tilde{J}}_T*\ddot{\theta }\left(t_i\right)\right]$

[0319] From the residual torque plots in FIG. 40, it is shown that the amount of measured torque, õ_S(t_i) that is not accounted by the model (the residual, τ_R(t)) is only about 4% on average. This demonstrates that the present model described by equation 1 fits to the data well, even when estimating only one set of parameters, {tilde over (τ)}_off, {tilde over (K)}_H, {tilde over (B)}_H and {tilde over (J)}_T over the entire data set in each trial. If, however the residuals are large, this would indicate that the equation of the present model is inadequate, or that a single set of time invariant, four parameters applied to the state data 44$\left(\theta \left(t\right),\stackrel{.}{\theta }\left(t\right),\mathrm{and}\ddot{\theta }\left(t\right)\right)$

[0320] over the whole trial is insufficient to describe the torque output. In other words, if the residual is large using only one set of time invariant parameters, the model may be good, but the impedance parameters of torque offset, stiffness and viscosity may be changing in time, within a given trial. To estimate how the impedance parameters change within a given trial, rather than have a single T_S, P_S and Ψ matrix over the whole trial as shown above, multiple T_S, P_S and Ψ matrices can be used by applying a smaller, one second moving window (501 samples) over the data as follows. First it is assumed that the total inertia does not change within each trial. This is a reasonable assumption since the distribution of mass of the moving parts about the center of the motor shaft does not change and the subject's hand is secured in place with the vertical bars 6. For each trial, {tilde over (J)}_T is used, which is evaluated above as the total inertia of the system. It is possible for the other three parameters to change within a trial. For a particular point in time, these now time varying parameters are denoted with subscripts as τ_off(t_i), K_H(t_i), B_H(t_i). Using these time varying parameters and subtracting the torque due to inertia results in a slightly modified equation of the model: 45${\tau }_S\left(t_i\right)-J_T*\ddot{\theta }\left(t_i\right)={\tau }_{\mathrm{off}}\left(t_i\right)+K_H*\theta \left(t_i\right)+B_H*\stackrel{.}{\theta }\left(t_i\right)$

[0321] The torque data set with out the contribution of inertia is now fit as follows: The data is fit using a linear least squares fit pseudo-inverse similar to above, but this is performed multiple times, once for each window as follows: 46$T_{\mathrm{si}}-{\tilde{J}}_T*{\ddot{\Theta }}_i={\begin{array}{c}{\tau }_s\left(t_{i-250}\right)-{\tilde{J}}_T*\ddot{\theta }\left(t_{i-250}\right)\\ ⋮\\ {\tau }_s\left(t_i\right)-{\tilde{J}}_T*\ddot{\theta }\left(t_i\right)\\ ⋮\\ {\tau }_s\left(t_{i+250}\right)-{\tilde{J}}_T*\ddot{\theta }\left(t_{i+250}\right)\end{array}}_{501\times 1}\mathrm{and}$${\Psi }_i={\begin{array}{ccc}1& \theta \left(t_{i-250}\right)& \stackrel{.}{\theta }\left(t_{i-250}\right)\\ ⋮& ⋮& ⋮\\ 1& \theta \left(t_i\right)& \stackrel{.}{\theta }\left(t_i\right)\\ ⋮& ⋮& ⋮\\ 1& \theta \left(t_{i+250}\right)& \stackrel{.}{\theta }\left(t_{i+250}\right)\end{array}}_{501\times 3}$

[0322] where i goes form 251 to 8064 that cover two periods of the random trajectory or 15.63 seconds). For each i, the matrices 47$\left[T_{\mathrm{si}}-{\tilde{J}}_T*{\ddot{\Theta }}_i\right]$

[0323] and Ø_i are obtained and P_i is computed for each i as follows: 48${\tilde{N}}_i={\begin{array}{c}{\hat{\tau }}_{\mathrm{OFF}}\left(t_i\right)\\ {\hat{K}}_H\left(t_i\right)\\ {\hat{B}}_H\left(t_i\right)\end{array}}_{3\times 1}$

[0324] where P_i is calculated similar to P above as follows:

Ñ_i=pinv(Ø_i)*Ô_Si=(Ø_i^T*Ø_i)⁻¹*Ø_i^T*Ô_Si

[0325] Since I goes from 251 to 8064, there are 7813 {circumflex over (τ)}_off, {circumflex over (K)}_H, and {circumflex over (B)}_H values per trial. The mean and variance or standard deviation of these values for each trial is determined. If the standard deviation of these values: std{{circumflex over (τ)}_off(t_i)}, std{{circumflex over (K)}_H(t_i)}, and std{{circumflex over (B)}_H(t_i)} are small, this indicates that the stiffness and viscosity of the subject's flexors and extensors remain fairy constant within each trial and the torque residuals of the time invariant model would be small as shown in FIG. 40. Note also that the mean of the time-variant impedance parameters values (denoted with the {circumflex over (x)}) should be very close to the evaluated time-invariant parameters (denoted with the {tilde over (x)}), that is:

mean{{circumflex over (τ)}_off(t_i)}≈{tilde over (τ)}_off;

mean{{circumflex over (K)}_H(t_i)}≈{tilde over (K)}_H; and

mean{{circumflex over (B)}_H(t_i)}≈{tilde over (B)}_H;

[0326] This is shown in FIG. 41.

[0327] The first row of bar plots in FIG. 41 show the results for the torque offset, stifffiess, and viscosity parameters obtained for each of the ten trials in experiment DMJF. The asterisk ‘*’ plots the values of {tilde over (τ)}_off, {tilde over (K)}_H, and {tilde over (B)}_H where one set of best-fit parameters over the whole trial as described above. The location of the tip of each bar graph represents mean{{circumflex over (τ)}_off(t_i)}, mean{{circumflex over (K)}_H(t_i)}, and mean{{circumflex over (B)}_H(t_i)} where a set of {circumflex over (τ)}_off, {circumflex over (K)}_H, and {circumflex over (B)}_H is obtained for each moving window of 501 data samples as described above. The height of the error bar indicates the standard deviation of these values, std{{circumflex over (τ)}_off(t_i)} std{{circumflex over (K)}_H(t_i)}, and std{{circumflex over (B)}_H(t_i)}.

[0328] The eleventh bar shows a summary of all ten trials. The asterisk on the eleventh bar shows the average of the ten {tilde over (τ)}_off, {tilde over (K)}_H, or {tilde over (B)}_H values, or in equation form: 49$\frac{1}{10}\sum _{k=1}^{10}\left({\tilde{\tau }}_{\mathrm{off}}\mathrm{for}\mathrm{trial}k\right),\frac{1}{10}\sum _{k=1}^{10}\left({\tilde{K}}_H\mathrm{for}\mathrm{trial}k\right),\mathrm{and}$$\frac{1}{10}\sum _{k=1}^{10}\left({\tilde{B}}_H\mathrm{for}\mathrm{trial}k\right)$

[0329] The location of the tip of the eleventh bar is the average of the height of the ten bar graphs, or in equation form: 50$\frac{1}{10}\sum _{k=1}^{10}\left(\mathrm{mean}\left\{{\hat{\tau }}_{\mathrm{off}}\right\}\mathrm{for}\mathrm{trial}k\right),\frac{1}{10}\sum _{k=1}^{10}\left(\mathrm{mean}\left\{{\tilde{K}}_H\right\}\mathrm{for}\mathrm{trial}k\right),\mathrm{and}\frac{1}{10}\sum _{k=1}^{10}\left(\mathrm{mean}\left\{{\hat{B}}_H\right\}\mathrm{for}\mathrm{trial}k\right)$

[0330] As mentioned earlier, there are 7813 {circumflex over (τ)}_off, {circumflex over (K)}_H, and {circumflex over (B)}_H values per trial. If all ten trials are grouped together, then there are 78130 {circumflex over (τ)}_off, {circumflex over (K)}_H, and {circumflex over (B)}_H values. The height of the error bar of the eleventh bar plot is the standard deviation of all 78130 values.

[0331] Experiment DMJF was performed on the right arm of a control adult subject tested at a bias of 0 degrees. Each control subject was also tested at a bias position of degrees flexion. The bottom row of graphs in FIG. 41 shows the results of parameters for two sets of experiments, DMJF, and DNJF, as labeled on the left of the graphs. The ‘M’ in the parentheses, indicates that the subject JF is a male. The bar graphs for experiment DMJF are plotted in white (bias of 0 degrees) and bar graphs for experiment DNJF are plotted in grey (bias of 30 degrees flexion) as indicated by the legend.

[0332] The results of subject JF show that the stiffness values when tested at a bias of degrees are markedly higher than when tested at a bias of 0 degrees in all ten trials. The average {tilde over (K)}_H across all ten trials are 0.876 N*m/rad and 1.30 N*m/rad at a bias of 0 and 30 degrees respectively. The asterisks on the eleventh bars show these values. Thus, the cumulative stiffness of all the muscles acting on the wrist joint is 48% higher at a bias of 30 degrees flexion.

[0333] Similarly, the viscosity values are markedly higher when tested at a bias of 30 degrees flexion than at a bias of 0 degrees for all ten trials. The average {tilde over (C)}_H across all ten trials are 0.0460 N*m*s/rad and 0.0714 N*m*s/rad at a bias of 0 and 30 degrees respectively and the cumulative stiffness of all the muscles acting on the wrist joint is 55% higher at a bias of 30 degrees flexion.

[0334] The resulting torque offsets, {tilde over (τ)}_off tell us something very interesting. The first noticeable thing about these values is not only that they are different in magnitude but also opposite in sign. τ_off is an important variable that tells us the angular position the hand prefers to be relative to the selected bias position or origin, θ₀. If the bias position was chosen to be the resting position of the wrist joint, then in the equation below θ_RP=0,

τ_off=−K_H*θ_RP

[0335] and ideally, τ_off would be 0 Nm. For a given angular stiffness, the further the bias position is selected from the preferred resting position of the wrist joint, the greater the torque offset, τ_off will be. The preferred resting position of the wrist joint is determined by the relative passive stiffness of flexor muscle group to the extensor muscle group. It should also be noted that τ_OFF is also dependent on the stiffness K_H. Since the active stiffness of muscles are influenced when they are perturbed, τ_OFF is influenced by these wrist perturbations or oscillations. If, however, the bias position is very close to the resting position of the wristjoint, τ_off will be very close to 0 Nm regardless of the magnitude of K_H since θ_RP≈0. A positive torque offset indicates that this contribution of torque is applied in the clockwise direction by the subject's arm and a negative torque offset indicates that this contribution of toque is applied in the counter clockwise direction. This information is important when interpreting the sign of {tilde over (τ)}off. The subject JF's right arm was tested in experiments DMJF and DNJF. As shown in the plots in FIG. 41, when tested at a bias of 0 degrees, the torque offset values are very slightly negative. The mean overall torque offset, {tilde over (τ)}_off for all ten trials shown by the asterisk on the eleventh bar is equal to −0.0538 Nm at a bias of 0 degrees, and the contribution of torque by {tilde over (τ)}_off is in the counter clockwise direction. At a bias of 30 degrees flexion, this value is equal to 0.110 Nm and the contribution of torque by {tilde over (τ)}_off in this case is in the clockwise direction. This means that there is less resistance in the wrist flexors at a bias of 0 than of 30 degrees flexion and the resting position of the hand of the control subject is close to the neutral flexion-extension position. The magnitude of these numbers is small as compared to the torque offset values in the spastic population as will be shown. It is important to note that if contra-lateral left arm of the subject was tested instead, then the sign of τ_off is expected to switch by virtue of symmetry. Thus, in order to compare the left hand to the right hand, if the left hand of a subject is tested, then the sign of τ_off is switched. Hence, the sign of τ_off no longer represents clockwise or counterclockwise torque. Rather, a positive torque offset indicates that the contribution of τ_off is in the direction towards the extension direction of the hand being tested relative to the selected bias position. Conversely, a negative torque offset indicates that the contribution of τ_off is in the direction towards the flexion direction of the hand being tested relative to the selected bias position.

[0336] Now, the results of each individual control adult will be shown. Either the left or right hand of 31 control adults was tested, at a bias position of both 0 and 30 degrees flexion. 10 of the subjects were male and 21 were female. The average age of control subjects tested was 31.29, the youngest was 18 years of age and the eldest was 60 years of age. The plots of the results for each individual subject are shown in FIGS. 42-52. All graphs that are plots of results of the same variable are plotted on the same scale so that visual comparisons can be made.

[0337] As shown in FIGS. 42-52, when looking at the average of {tilde over (τ)}_off, {tilde over (K)}_H and {tilde over (B)}_H over all ten trials for each experiment (denoted by the asterisk over the eleventh bar labeled trial 11), the following conclusions can be made. At a bias position of 0, the magnitude of {tilde over (τ)}off across all ten trials, for all 31 subjects is negligible, some being positive in value, some negative. This indicates that the resting position of a control subject's hand is very close to the neutral flexion/extension where muscles acting on the wrist exhibit the least amount of tension. At a bias of 30 degrees the torque-offset values are significantly higher and positive in sign. This indicates that the hand exhibits more elastic tension at this bias position. The positive sign indicates that the hand prefers to be at the neutral 0 degree bias position. All 31 subjects had a higher stiffness at a 30 degree flexion bias than at a 0 degree bias. 30 subjects had a higher mean viscosity value at a 30 degree flexion bias than at a 0 degree bias. Only 1 subject (MT) had a lower mean viscosity value at a 30 degree flexion bias than at a 0 degree bias. The average and standard deviation (denoted by the height of the error bars) of {tilde over (τ)}_off, {tilde over (K)}_H and {tilde over (B)}_H values across all control adults for each trail are shown in FIG. 53. The eleventh bar shows the overall average and standard deviation across all control adults, across all ten trials.

[0338] Null-hypothesis significance tests can be used to quantitatively test how significant these different observations are between testing at the bias of 30 degrees flexion versus the bias of 0 degrees. For each experiment, the mean {tilde over (τ)}_off, {tilde over (K)}_H, and {tilde over (B)}_H are taken over the ten trials and these parameters are separated into two groups according to bias position. Thus, there are 31 mean{{tilde over (τ)}_off}, mean{{tilde over (K)}_H}, and mean{{tilde over (B)}_H} parameters for each of the two bias positions. The null hypotheses test indicates that the mean torque offsets, mean stiffnesses and mean viscosities are not significantly higher at a bias of 30 degrees than they are at a bias of 0 in control subjects. From the null hypothesis test, the p-value is obtained, wherein the p-value is the probability of observing the given result by chance given that the null hypothesis is true. Small values of p cast doubt on the validity of the null hypothesis. The probability of the null hypothesis stating that the average torque offset values are higher at 0 than at 30 degrees is true is p=9.2*10⁻⁴ The probability of the null hypothesis stating that the average stiffness values are higher at 0 than at 30 degrees is true is p=2*10⁻⁵ The probability of the null hypothesis stating that the average viscosity values are higher at 0 than at 30 degrees is true is p=4.2*10⁻⁵. Since a resulting p-value of 5*10⁻³ is considered to be very significant, there is no question that the increase in all three parameters {tilde over (τ)}_off, {tilde over (K)}_H, and {tilde over (B)}_H is extremely significant when testing at a bias position of 30 degrees flexion rather than 0 in the control adult subject population.

[0339] Table 3 below shows the average and standard deviation of {tilde over (τ)}_off, {tilde over (K)}_H and {tilde over (B)}_H values across all ten trials across all control adult subjects which are indicated by the height of the 11^th bars shown in FIG. 53. Table 4 shows the range of average {tilde over (τ)}_off, {tilde over (K)}_H and {tilde over (B)}_H values across all ten trials in the control adult group. 3

[0340] 4

[0341] An additional observation from FIG. 53 is that when testing at a bias of 0 degrees, the elastic stiffness increases from trial 1 to trial 10. As stated previously, trials 1 and 2 have a maximum frequency component of 4 Hz, trials 3 and 4 have a maximum frequency component of 5 Hz, trials 5 and 6 have a maximum frequency component of 6 Hz, and trials 7 and 8 have a maximum frequency component of 8 Hz. One hypothesis that can explain this observation is that because velocities of the perturbations increase with trial number, the primary afferents of the muscle spindle sense these greater rates of muscle stretch in turn activating the stretch reflex causing muscles to contract and become stiffer. Thus, the hypothesis says that the combined active stiffness of the muscles acting on the wrist may be increasing with trial number. An observed increase in EMG activity at the greater trial numbers supports this hypothesis. This trend, however, is not shown at the bias of 30. This is likely because extensor muscles are already active in all 10 trials as they are stretched to a length further from the neutral position.

[0342] Many subjects with spasticity cannot be tested at a neutral bias position of 0 degrees as they have difficulty extending their wrist joint to that position. The resting configuration of the hand of someone with spasticity is in the hyper-flexed position due to hyperactivity of the flexor muscle group. Four adult patients with hemiplegic spasticity were tested, AM, JS, SH, and BY, 3 males and 1 female. A hemiplegic patient is one where one side of his or her body is affected by spasticity whereas the contra-lateral side is practically unaffected or less affected. Because the contra-lateral arm is unaffected or less affected, it can serve as a control to the spastic arm when testing a patient with hemiplegic spasticity. JS is a patient with Parkinson's disease, SH suffers from traumatic brain injury due to a motorcycle accident, and BY is a patient who severely lost function of his left hand after a stroke. AM's affected arm was tested at a bias of 30 degrees flexion. For subjects JS, SH, and BY the affected arms and unaffected contra-lateral arms were tested, both at a bias of 30 degrees. The results of these experiments are shown in FIGS. 54 and 55. Unlike the figures showing results of the control subjects, all plots are not on the same scale due to the large range in values. The white bar graphs show results for the unaffected arm and the grey bar graphs show results for the contra-lateral affected arm. Both arms were tested at a bias of 30 degrees flexion.

[0343] Looking at the average of {tilde over (τ)}_off, {tilde over (K)}_H and {tilde over (B)}_H over all ten trials for each experiment denoted by the asterisk over the eleventh bar labeled trial I 1, the following conclusions can be made. In the experiment with the unaffected hand at the bias of 30 degrees flexion, the average {tilde over (τ)}_off across all ten trials are comparable to what is seen in control subjects when tested at the same bias. As in the control subjects, the torque offset values are positive indicating that the unaffected hand prefers to be at the neutral 0 degree bias position as in the control subjects. However, when testing the affected, spastic side at the 30-degree bias, the {tilde over (τ)}_off values are negative in sign and greater in magnitude. This indicates that the preferred resting position of the subject's spastic wrist joint is at a flexed position even greater than 30 degrees, as the hand is pushing towards an even further flexed position. This phenomenon is what is called hyper-flexion. This observation of increase in magnitude and change in sign in the average {tilde over (τ)}_off across the ten trials when testing the affected hand versus the unaffected hand is significant (p=0.0177). All 3 subjects tested on both arms have a higher average stiffness (mean of {tilde over (K)}_H across the ten trials) in their spastic arm as compared to their contra-lateral arm (p=0.0743). When comparing the average stiffness between contralateral arms, in subjects JS, and SH, the difference is large. Also, in general, the average viscosity, {tilde over (B)}_H, across the ten trials is greater in the affected hand when compared to the unaffected hand, but the difference seen in viscosity values is not as great as the difference seen in the stiffness values (p=0.199). Here, the difference between the affected spastic arm is compared with the unaffected contra lateral arm within the spastic group.

[0344] The average and standard deviation (denoted by the height of the error bars) of {tilde over (τ)}_off, {tilde over (K)}_H and {tilde over (B)}_H values across all control adults for each trail are shown in FIG. 56. The eleventh bar shows the overall average and standard deviation across all control adults, across all ten trials.

[0345] Table 5 below shows the average and standard deviation of {tilde over (τ)}_off, {tilde over (K)}_H and {tilde over (B)}_H values across all ten FIG. 56. Table 6 below shows the range of average {tilde over (τ)}_off, {tilde over (K)}_H and {tilde over (B)}_H values across all ten trials in the group of adults with hemiplegic spasticity. 5

[0346] 6

[0347] Tables 5 and 6 show that at the bias position of 30 degrees flexion, the average parameter values of the unaffected or less affected arms of the hemiplegic population are still somewhat higher than when compared to the control population. This demonstrates that the tone of the less affected side of the patient group in the is still greater on average than what would be seen in the control group, though this difference is of course greater when comparing the severely affected side of the hemiplegic population with the control population. This is a fairer comparison since subjects in the control group have no known neurological disorders. When comparing the hemiplegic subjects to the controls, the observation of the negation of sign and increase in magnitude of the average torque offset values, {tilde over (τ)}_off across the ten trials shown by the 11^th bars are extremely significant (p=3.0*10⁻¹⁰). The increase in average stiffness, {tilde over (K)}_H and viscosity, {tilde over (B)}_H values across the ten trials when compared across populations are also extremely significant (p=1.5*10⁻⁸ and p=6.5*10⁻⁵ respectively).

[0348] On average, the stiffness and viscosity values across all trials of the control group were 1.305 N*m/rad and 0.0512 N*m*s/rad respectively. The corresponding stiffness and viscosity values of the severely affected side of the hemiplegic population were 3.858 N*m/rad and 0.0992 N*m*s/rad respectively. These elevated values when compared to those of the control group are extremely significant (p=1.5*10⁻⁸ for stiffness and p=6.5*10⁻⁵ for viscosity).

[0349] While the device and method have been described in detail for use in quantifying muscle tone, particularly in the wrist, in a spastic patient, the invention is not so limited. Rather, the device and method are useful on both the upper and lower extremities including, for example, the ankle. Thus, for example, the device could be modified to suit the leg and ankle accordingly. Further, the method and device can be used to measure muscle tone in the non-spastic patient. For example, the method and device can be used to measure rigidity in a patient. Further, the method and device could be used in analyzing a patient's muscle tone in line with, for example, physical therapy that the patient is undergoing to regain control of the muscles in the legs after a spinal accident.

[0350] Although a preferred embodiment of the invention has been described using specific terms, such description is for illustrative purposes only, and it is to be understood that changes and variations may be made without departing from the spirit or scope of the following claims.