Received November 11, 2015, accepted November 30, 2015, date of publication January 5, 2016, date of current version March 2, 2016. Digital Object Identifier 10.1109/ACCESS.2016.2514978

Magnetic Field Control for Haptic Display: System Design and Simulation QI ZHANG, HAIWEI DONG, (Member, IEEE), AND ABDULMOTALEB EL SADDIK, (Fellow, IEEE) Multimedia Computing Research Laboratory (MCRLab), School of Electrical Engineering and Computer Science, University of Ottawa, Ottawa, ON K1N 6N5, Canada

Corresponding author: H. Dong ([email protected])

ABSTRACT We present magnetic rendering, a new haptic display method applying an electromagnet array to produce magnetic field in mid-air where the force field can be felt as magnetic repulsive force exerted on the hand through the attached magnet disks. The magnetic field is generated by a specifically designed electromagnet array driven by direct current. By attaching small magnet disks on the hand, the tactile sensation can be perceived by the user. This method can provide a strong tactile force and avoid cumbersome attachments with wires, and thus, it is suitable for a colocated visual and haptic display. In this paper, we introduce the detailed design of the electromagnet array for haptic rendering purposes, which is modeled and tested using finite-element method simulations. We characterize the model mathematically, and apply recursive least squares adaptive control algorithm for controlling the magnetic field. We evaluate the performance of our simulated model in terms of force strength, operation distance, and force stiffness. We then implement and test the control algorithm, which results in fast and accurate convergence. We conclude with simulations on a 15-by-15 model to generate a haptic human face, which results in a smooth force field and accurate force exertion on the control points. INDEX TERMS Mid-air haptic display, 3D volumetric shape, magnetic force, recursive least squares. I. INTRODUCTION

In realistic interactions with virtual or remote objects, haptics has garnered enormous interest as the core of realizing manual interaction with environment [1]. Thrilling possibilities are opened up by introducing haptics to various disciplines involving virtual reality [2], medical training [3] and many other applications, where haptic rendering [4] plays a key role in generating force feedback. In the context of human interacting with virtual objects, rendering is labeled as the process of imposing certain stimulation on the user to convey the properties of the virtual object [5]. The properties may contain detailed perceptive information ranging from shape to texture, elasticity and so on. To explore the virtual environment, one of the most common approaches is to apply a proxy robotic device which has been adopted in many haptic devices such as the Geomagic Touch X (formerly Sensable Phantom Desktop) by 3D Systems [6]. Geomagic Touch X enables users to interact with virtual objects via a robotic pen with an articulated arm attached. The force feedback is exerted when the proxy pen contacts a virtual object in form of resistance for the users to perceive. A similar approach can be seen in Lorentz force magnetic levitation devices [7] where the articulated arm is replaced by magnetic force. This method provides VOLUME 4, 2016

point-wise force feedback, which is not an intuitive way for human perception using hands. Besides, due to the separation of the visual image and haptic interface, this method is not suitable for co-located 3D visual and haptic display. A solution to this problem is mid-air tactile sensation generation, which is addressed by wearable haptic devices [8]–[12], air jet [13], air vortex [14], [15], ultrasound [16], [17], and magnetic force [18], [19]. The former two approaches are known for long operation distance but with very weak force strength and can only generate discrete focal force points. Magnetic repelling force in the other hand, is much stronger within short distance and can form smooth continuous force field in space. Hence we introduce the idea of Magnetic Rendering, which refers to a new haptic display method applying an electromagnet array to produce magnetic field in the air where the force field can be felt as magnetic repulsive force exerted on the hand through the magnet disks attached. In this paper, we present a magnetic rendering system to render volumetric shapes [20], [21] by applying magnetic repelling force to produce tactile sensation in space. By attaching small magnet disks on the hands, users can feel the haptic shapes just by moving their hands through the air. The proposed system allows users to perceive continuous

2169-3536 2016 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

299

Q. Zhang et al.: Magnetic Field Control for Haptic Display

FIGURE 1. Magnetic Rendering concept.

force field as shown in Fig. 1. The magnetic field is generated by carefully designed electromagnet array whose advantage is proved by Finite Element Method (FEM) simulations and is controlled using recursive least squares (RLS) algorithm. The main contributions of our work are: •

•

•

We propose a Magnetic Rendering system for users to feel volumetric shapes with little attachment on the hand by using magnetic repelling force as tactile feedback approach. We design a new electromagnet array to produce high magnetic field strength outside one side of the array, which is evaluated by FEM simulations. For a single electromagnet model, the design of a magnetic concentrator lifts its magnetic field up and the mechanics behind which is derived mathematically. The introduction of a magnetic isolator shields the magnetic field between each two elements and avoids severe interference. The final configuration of an array is obtained by optimized design in terms of magnetic field interference and smoothness of the force field. We prove the linearity characteristics of the system and introduce recursive least squares (RLS) algorithm for controlling the electromagnet array to generate the desired magnetic field for users to perceive. Simulation shows that the algorithm converges accurately and fast.

II. BACKGROUND AND RELATED WORK A. CONTACTLESS TACTILE SENSATION GENERATION

A variety of methods have been studied to create contactless tactile sensation, however there is still no device that can render a continuous volume shape with strong force feedback for users to perceive. One of the most popular methods is wearable devices, which aims to attach force vector stimulator directly on the hands and impose feedback such as vibration [8], [9] and cutaneous forces [10], [11]. These methods all suffer from the fact that the attachments can be cumbersome and they can not produce actual pressure on skin thus are not suitable for haptic volumetric shape rendering purposes. Air jet is chosen intuitively in [13] to produce contactless haptic feedback with low accuracy. Unlike the standard jets of air, air vortices which are rings of air are also used in devices such as AirWave [15] and AIREAL [14]. However both these methods have low fidelity due to the non-focusing stimulating area. Ultrasound radiation force is another candidate applied to generate discrete focal force 300

points in 3D space [16], [17]. The drawback of this method is that the ultrasound acoustic force is very weak and only limited focal points can be generated. Magnetic force is widely used in magnetic levitation for that it is sufficient to counteract gravity or other force imposed, and in recent years it is adopted by a number of devices [18], [19] as a mid-air force feedback actuation method. Although the high strength of magnetic force is a desired property in haptic rendering applications, it still has drawbacks such as the force between magnets decreases rapidly with the distance and it is hard to control due to its nature of non-linearity. We apply magnetic force as the mid-air tactile feedback in our design for that it is strong and contactless in nature. The problem of short operation distance is addressed by a new design of electromagnet array. which proves to have near-linear characteristics that facilitate the design of control method. B. MAGNETIC FORCE FOR CONTACTLESS HAPTIC INTERACTION

Magnetic force between two magnets obeys the complex electromagnetism rules. When one of them is an electromagnet, the exerted force can be controlled by the amount of electric current running in the winding which is much stronger compared to the pressure of air vortex or ultrasound acoustic force. Another merit of magnetic force is that it can form continuous force volume comparing with other methods, which can only produce discrete force points. The limitation of using magnetic force for haptic rendering is that it can not render sharp shapes which origins from the fact that magnetic field is oval-shaped in nature. For electromagnets used for haptic force actuation purpose, ferrite cores are desired for that they have high magnetic permeability which can help generate strong magnetic field, whereas result in non-linearity of the produced force in respect to the stimulating current. Moreover, the produced magnetic force is highly nonlinear with the distance between the magnets and in fact, it drops with cube of the distance. When it comes to an array of electromagnets, the coils can not be controlled individually due to the superposition of magnetic field. All these properties make precise control a challenging task. FingerFlux [18] is a haptic device generating near-surface force feedback for tabletops using attractive and repulsive magnetic force and vibration to guide the users when approaching the screen. They also apply an array of electromagnet as force feedback actuator and magnets attached on the hand as force receiver. However since the device is not specifically designed for haptic volume shape display, the force produced is not strong enough for our design, and no numerical analysis of the force generated by the array is done, which leaves no chance for us to derive a control method. Co-located 3D graphic and haptic display [19] is another device using force between electromagnet arrays and magnets for haptic rendering. For the force receiver, a pen with VOLUME 4, 2016

Q. Zhang et al.: Magnetic Field Control for Haptic Display

magnets attached is designed for a user to hold so that the force feedback can be felt. The limitation of this method is that the force feedback is single-point, that is, users can only feel the force on one point at a certain time. For haptic rendering purposes, the force feedback not only needs to be strong enough for user to perceive from a distance, but it also needs to cover several parts of the user’s hand simultaneously for the user to feel multiple force points at the same time. The electromagnet array model in our design is optimized in terms of layout and space between every two elements. Besides that, physical characteristics of the force generated are determined in order to derive a control method.

electromagnet array whose size can be scaled as needed and each one of the electromagnets is driven by direct current (DC) output from power source. For a given set of 3D volumetric data, preprocessing is required to fit the data with the size and the force range of the electromagnet array, after which the data is fed into the controller and transformed into suitable current values based on the control algorithm. The final configuration of the current values is then imposed on the array to generate the magnetic field and when placing hands with magnets attached above the array, the pressure sensation of the virtual object can be perceived in the midair.

C. ELECTROMAGNET ARRAY DESIGN AND CONTROL

Electromagnet arrays are mostly used for the purpose of magnetic levitation as planer motors [22], [23], however, some ideas involving boosting the magnetic field strength can be tapped out to work for haptic rendering. To augment the magnetic field above the electromagnet, one can increase the number of turns of the coil or raise the stimulating current value, which all lead to the increase of spread magnetic field by boosting the total magnetic field. The current running in the wires has strict restraint due to material and manufacturing limitations such as the current capacity of the wire and the heating effect of the coil. Moreover, small size electromagnets are desired for that the force field generated can be felt with detail, which renders large number of wire turns impossible. Therefore the main consideration falls into concentrating the magnetic field to only one pole of the electromagnet, which can also increase the efficiency of energy usage. In fact, the model is ideally expected to be unipolar. However, since unipolar magnet does not exist in reality, one can only augment the magnetic field on one side by arranging permanent magnets in special patterns. One good example is Halbach array [24] which aligns permanent magnets in a special way where the magnetic field is reinforced on one side while cancelled on the other side to near zero. The Halbach effect is discovered in 1-dimensional arrays and is extended to 2-dimensional arrays [25] where nearly one-side magnetic field is obtained by the special layout arrangements. However, on the magnetic field strengthened side, the direction of the field is interchanging which is undesirable for haptic rendering where only one pole is needed. In our design of electromagnets suitable for haptic rendering, we develop a new model with a magnetic field concentrator inspired by Halbach effect which proves to have very good performance as discussed in Section III. III. DESIGN OF THE MAGNETIC RENDERING SYSTEM

The proposed magnetic rendering system is a haptic interface applying repulsive magnetic force as haptic feedback to render desired shapes haptically in mid-air for users to perceive. The haptic shape is rendered by a specially designed VOLUME 4, 2016

FIGURE 2. Magnetic Rendering System.

The architecture of the proposed magnetic rendering system is shown in Fig. 2 in which the four main components are: • FEM simulation, which is used for designing the electromagnet array model, and is responsible for analysing various electromagnetic properties of the model. • 3D surface data, which can be obtained from 3D scanner or other volumetric data generation tools. • Haptic rendering, which produces magnetic force against the user magnets under the control of a designed algorithm based on the properties of the model and the simulation results. Having the system structure determined, our main tasks falls into two parts: electromagnet array actuator design and control methods development. In this section, the design concerns and methods of the whole electromagnet array will be introduced. The control methods based on the properties of the actuator model are introduced in Section IV. For the following discussions, two sets of coordinates are used for the overall array and individual element respectively as shown in Fig. 3. The Cartesian coordinates (Ex , Ey, Ez) are used for describing the locations of the elements in the array, and the cylindrical coordinates (rE, θE, Ez) are used for magnetic field derivations for individual electromagnet. For haptic rendering devices, the repulsive magnetic force generated is the most important criterion when 301

Q. Zhang et al.: Magnetic Field Control for Haptic Display

TABLE 1. Summary of design parameters.

FIGURE 3. Global Cartesian coordinates and local cylindrical coordinates.

making decisions, besides which there is also stiffness which is the rigidity of the virtual object. In the simulation step, magnetic flux density B is solved by FEM simulation. B. CONFIGURATIONS OF THE INNER COIL

FIGURE 4. Structure of the designed model. (a) Core and coil parts. (b) Concentrator part. (c) Complete model.

A. MAGNETIC FIELD CONCENTRATED ELECTROMAGNET MODEL STRUCTURE AND DESIGN PARAMETERS

The components of the designed magnetic field concentrated electromagnet model is depicted in Fig. 4. We start with a basic electromagnet model shown in Fig. 4a which only consists of a metal core and a coil. To augment the power usage efficiency, we introduce a magnetic field concentrator which concentrates the magnetic field to the upper part of the model as shown in Fig. 4b. When constructing an array, the model is expected to be isolated so that the magnetic field generated by neighbouring elements would not cause serious interference. To realize this, an isolator is designed to shield the magnetic field and the complete structure of a single designed electromagnet model is depicted in Fig. 4c. From innermost to outermost the components of the electromagnet model are soft iron core, multi-turn coil which are for magnetic field generation, permanent magnet concentrator for magnetic field concentration and mu-metal shell for magnetic field isolation. A summary of the designed parameters is listed in Table 1 The detailed design is discussed in the following order: in Section III-B, we describe the configuration of the core and the coil of the electromagnet model. In Section III-C, we introduce the design and mechanics of the magnetic field concentrator. In Section III-D, we introduce the design of the isolator which shields the magnetic field from spreading widely. Finally, the layout of the array is discussed in Section III-E. 302

The electromagnet model is based on a solenoid structure which contains a coil of insulated wire tightly wound on a cylinder metal core. When stimulated by DC, magnetic field is exerted by the inner coil due to Ampere’s law and is magnified by the core in the center. The parameters of the inner coil are chosen based on industrial standards and manufacturing limits. The cylinder core of the coil has a height of 2 cm and a diameter of 1.4 cm and is designed to be made of soft iron for that it can withstand strong magnetic field while not reaching the saturation point which makes it a preferred material as electromagnetic core. However, soft iron is a non-linear magnetic material, whose permeability is not a constant but depends on the magnetic field strength. The dimensions of the coil is set based on a coil-winding manufacturer, where the outer diameter is set as 2.2 cm. The windings has a total cross-sectional area of 0.8 cm2 . To generate a powerful magnetic field, the number of turns is set to be 1000, thus the cross-sectional area for the wire is set to be 0.08mm2 which is very close to gauge 28 of the standard wire sizes according to American Society for Testing and Materials (ASTM) Standard B258-02 [26] whose current carrying capacity is about 1 A. Therefore the current range of the electromagnet is from −1 A to 1 A. C. MAGNETIC FIELD CONCENTRATION

Inspired by Halbach array [24] and the magnetic arrays with increased magnetic flux [25] introduced in Section II-C where the magnetic field is increased on one pole by special arrangement, in our design, the coil is encased in a cylinder ring permanent magnet acting as a magnetic field concentrator to enhance the magnetic flux on the upward side. The directions of magnetization and corresponding current flow imposed are depicted in Fig. 5. When imposing counter-clockwise current into the coil in the vertical view, the magnetic field lines generated are pointing upward according to the right-hand rule, thus the coil can be viewed as a magnet with north pole on the top. VOLUME 4, 2016

Q. Zhang et al.: Magnetic Field Control for Haptic Display

bottom JSM = (−Mc rE) × (−Ez) = −Mc θE(A/m) E E × (−r) E = 0(A/m) J inner = (−Mc r) SM

outer E E × rE = 0(A/m) JSM = (−Mc r)

(7) (8) (9)

top

FIGURE 5. Current and magnetization directions settings. (a) Magnetization direction of the concentrator and the magnet disk. (b) Mechanism of the magnetic field concentrator.

The concentrator, which is a permanent magnet cylinder ring depicted in Fig. 5a is radially magnetized with north pole pointing to the axis. The mechanism of the concentrator is that it contributes the current density in the form of magnetization currents that flow both within the volume and on the surface of the concentrator. The detailed statement and proof is shown as follows. Theorem 1: Given a cylinder ring uniformly radially magnetized with the north pole pointing to its central axis, the volume magnetization can be replaced by surface current on the top and bottom when calculating the magnetic field. Proof: The equivalent volume and surface currents of a volume magnetization under the cylindrical coordinates is shown in Fig. 3, which are calculated by JVM = ∇ × M (rE, θE, Ez) E θE, Ez) × nE JSM = M (r,

(1) (2)

where JVM is the volume magnetization current and JSM is the surface magnetization current [27]. M (rE, θE, Ez) is the magnetization vector field expressed in cylindrical coordinate system, and nE denotes the unit normal vector to the surface of the E θE, Ez) is the curl of the vector magnetized object. ∇ × M (r, field, which can be written explicitly in cylindrical coordinate system as ∂Mθ 1 ∂Mz − )rE ∇ × M (rE, θE, Ez) , ( r ∂θ ∂z ∂Mr ∂Mz 1 ∂(rMθ ) ∂Mr +( − )θE + ( − )Ez ∂z ∂r r ∂r ∂θ (3) By assuming the magnetization is uniform and radial, the magnetization vector field can be described as E θE, Ez) = −Mc rE M (r,

(4)

where Mc is a positive constant representing the magnitude of the magnetization. By substituting Eq. 4 and Eq. 3 into Eq. 1 and Eq. 2, we can obtain the volume magnetization current density as 2 E E = 0(A/m JVM = ∇ × (−Mc r) )

(5)

and the surface magnetization current density for the top, bottom, inner side and outer side surfaces as top

E JSM = (−Mc rE) × Ez = Mc θ(A/m) VOLUME 4, 2016

(6)

bottom , J inner and J outer are the magnetization where JSM , JSM SM SM surface current on the top, bottom, inner and outer side of the concentrator respectively. Remark 1: From Eq. 6 we can see that when looking down from above, the magnetic concentrator contributes to the top surface current with a counter-clockwise magnetization current and to the bottom surface current with a clockwise magnetization current according to Eq. 7. When visualizing the magnetization current in Fig. 5b, we can see that the directions of the magnetization current and the stimulating current in the coil are the same on the top surface of the concentrator while opposite on the bottom surface. Thus, the magnetic field is augmented on the top and reduced on the bottom. In reality, the magnetization value of a normal bar magnet can be 105 A/m while for iron it can be as large as 106 A/m. To ensure the electromagnet model and the magnet disks having magnetic field strength of approximately the same orders of magnitude, the magnetization of the concentrator is set to be 5 × 105 A/m and that of the disk magnets is 106 A/m.

D. MAGNETIC SHIELDING

For magnetic field isolation purpose, alloys with high permeabilities (µr  1) are suitable for shielding electromagnets [28]. The shield acts as a better conductor of magnetic flux thus provides a low reluctance path than air or other materials of the electromagnet. The magnetic field lines, which appear at the radius of the electromagnet, are then stay within the shell. Mu-metal is one of the most efficient materials for magnetic field shielding. It is a nickel-iron alloy composed of approximate 80% nickel and other materials including iron, copper and so on depending on the mix. The relative permeability of mu-metal can reach as high as 80, 000 − 100, 000 compared with several thousand for steel and 1 for air. The thickness of mu-metal shell is set as 2 mm which is one of the mu-metal sheet product specifications. For the shape of the mu-metal shell, we initially plan to make it a half-open cylindrical chamber closed only at one end while leaving the desired magnetic pole exposed. The reason is first to shield the radial magnetic field and second, to shield one pole of the model thereby make it a nearly unipolar electromagnet. However, although the magnetic field is shielded effectively by the mu-metal shield with single bottom, most of the electromagnetic energy is wasted in the bottom of the shell and only a minor part of the magnetic field lines generated are ejected from the top. Therefore, the mu-metal shell is designed to be an open-ended cylinder wrapping the concentrated model without a bottom cap. The magnetic flux density norm distribution and the values measured from 1 cm above 303

Q. Zhang et al.: Magnetic Field Control for Haptic Display

centre generated by the surrounding elements in both layouts as follows. 4 8 X 1 X square 1 square square + √ Bsurrounding = 3 Bi Bi 3 d ( 2d) i=5 i=1 square Bsurrounding

6 6 X hexagon = 3 Bj d

(10)

j=1

square

FIGURE 6. Comparison of magnetic field in two shell shapes. (a) Magneitc field of the model using mu-metal shell with bottom. (b) Magneitc field of the model using mu-metal shell without bottom.

where Bsurrounding is the magnetic field above the centre element generated by the surrounding elements in the hexagon square layout and Bsurrounding is that in the hexagon layout. By assuming the currents running in all the elements are of the same value, the magnetic field generated will be the same, square

the top of the electromagnet along the axis are shown in Fig. 6. Instead of concentrating at the bottom corners as shown in Fig. 6a, the magnetic field in the modified design concentrates at the upper part of the core as shown in Fig. 6b, which results in a stronger magnetic field above the electromagnet. Here the measurement is taken at 1 cm above along the central axes of both models and shows an increase of 37.56% in magnetic flux density norm.

hexagon

= Bs , (i = 1, 2...8, = Bj we denote it as Bs , and Bi j = 1, 2...6). Then from the Eq. 10 we can see the calculated interference from the surrounding elements is about 5.14 B d3 s in the square arrangement, which is less than d63 Bs of the hexagon layout, therefore square layout is chosen in our design.

E. LAYOUT OF THE ELECTROMAGNET ARRAY

To form an array using the electromagnet designed as a haptic display, the arrangement of the electromagnets is preferred to be compact while the details of the force field are not blurred by the magnetic field superposition. FIGURE 8. Experiments to determine the spacing between every two elements. (a) Layout of the two elements. (b) Force along x axis for two elements at 2cm above.

FIGURE 7. Two candidate layouts of electromagnet array. (a) Square layout. (b) Hexagon layout.

Two candidate arrangements of array have been considered, as shown in Fig. 7. We use the magnetic field interference as a criterion. When considering the influence of surrounding elements to the one in the center, we can show that the influence is smaller in the square layout. By assuming the distance between the nearest two elements in both layouts is d as shown in Fig. 7, we denote the magnetic field above the central element of the square layout by each of the 8 surroundsquare ing elements is Bi (i = 1, 2, 3...8) and that of the hexagon hexagon layout is Bj ( j = 1, 2, 3...6). Since the magnetic field obeys the superposition principle and decays with the cube of the distance, then we have the magnetic field above the 304

Having the layout of electromagnet array determined, the distance between each element is another parameter that needs concern as it is a trade-off factor between force field smoothness and controllability. As mentioned before, the superposition effect abates the precision of control while helps to form a smooth force field for haptic perception. Therefore the force field produced by two adjacent electromagnets with different spacings are simulated with the layout shown in Fig. 8a and the forces along the horizontal axis are measured at 2 cm above the top of the elements. The current values of the two elements are both set to be 1 A and Fig. 8b shows the differences of magnetic vertical force along X axis among three spacing strategies which are 0 cm, 0.5 cm and 1 cm spacing in between. The simulation result clearly shows the magnetic force superposition effect attenuates with the spacing at both vertical distances. All the three arrangements preserve the desired peak force along the axises of the two element and the force drops the least in between the two elements when there is no spacing. To summarize, the arrangement of the electromagnet array is decided to be of square layout and with no space between every two elements. VOLUME 4, 2016

Q. Zhang et al.: Magnetic Field Control for Haptic Display

IV. CONTROL OF THE MAGNETIC FIELD

with a linear transformation

To develop valid control methods for the newly designed electromagnet array model, we have to characterize the model regarding the current-to-force transformation pattern. In Section IV-A we analyse the model mathematically to find its linear property and verify it by FEM simulation. Recursive least squares (RLS) control method is introduced in Section IV-B to estimate the transformation matrix of the multiple-input-multiple-out linear system adaptively. A. LINEARITY OF THE SYSTEM CHARACTERISTICS

Having the M by K electromagnet array designed, control methods need to be developed to render certain 3D surface haptically. For the 3D surface data, no matter if it is static or dynamic, the generated magnetic field is either static or a slowly changing field which means in a short time slot, the magnetic field can be seen as static. The proposed magnetic rendering control processes for the N = M × K elements are shown in Fig. 9.

Fi,z =

The expected forces Fs are obtained by sampling and rescaling the 3D surface data S using certain function p(S) and the electromagnet array is seen as a current-to-force transformation process denoted by f (I) which depends on the electromagnet array model and the vertical distance between the magnet disk and the array. The control algorithm then needs to produce the current values of all the elements based on the expected forces and f (I). Note that the force is a vector, and here we use the vertical component of the magnetic force at the distance of 1 cm as the control points when developing the control algorithm. We denote the transformation from current values to force values as the following equation (11)

where N is the number of elements in the array, Fi,z (i = 1, 2, 3 . . . N ) denotes the force exerted on the magnet disk above the ith element and fi (i = 1, 2, 3 . . . N ) are N unknown functions depending on the system which need to be determined in order to model the system. The current-force transformation of the array proves to be a multi-input-multi-output linear system. The detailed proof is shown as follows. Theorem 2: Given the designed electromagnet array model described in Section. III and the stimulating current range [−1 A, 1 A], the vertical force on the disk above the ith element caused by all the electromagnets can be modelled VOLUME 4, 2016

N X

βij

(12)

j=1

where αij and βij are coefficients depending on the locations of source electromagnet and the magnet disk. Proof: The vertical magnetic force on a magnet disk is calculated as a volume integral [29] shown in Eq. 13. Z ∂Bz EzdV (13) Fi,z = Mdisk · ∂z  where Mdisk is the magnetization magnitude of the magnet disk which is along the z direction, Bz is the z component of the magnetic flux density and  is its volume. By assuming the N elements are stimulated by DC Ij , (j = 1, 2, 3 . . . N ), we define the magnetic flux density generated in the location of the ith magnet disk (i = 1, 2, 3 . . . N ) individually along the direction of z as i = 1, 2, 3 . . . N , j = 1, 2, 3 . . . N

(14)

where gi,j is a function depending on the designed electromagnet model. By substituting Eq. 14 into Eq. 13 we can obtain the magnetic force caused by the jth electromagnet on the ith magnet disk Z Z ∂gi (Ij ) − ∂Bi,j,z − → → z dV = Mdisk · z dV Fi,j,z = Mdisk · ∂z ∂z   (15)

FIGURE 9. Magnetic rendering control processes.

i = 1, 2, 3 . . . N

αij Ij +

j=1

Bi,j,z , gi,j (Ij ),

Fi,z , fi (I1 , I2 , I3 . . . IN ),

N X

When considering all the elements in the array, since the magnetic field superposition principle states that the total magnetic field generated by multiple separated sources is equal to the sum of the fields generated by each individual source as in Eq. 14, the total magnetic field in the location of magnet disk i along the direction of z generated can then be represented as Bi,z =

N X j=1

Bi,j,z =

N X

gi (Ij )

(16)

j=1

By substituting in Eq. 16 into Eq. 13 we can obtain the expression of the total magnetic force caused by all the electromagnet on the ith disk.  N Z X ∂gi (Ij ) − → Fi,z = Mdisk · z dV (17) ∂z  j=1

From the above analysis we can see that to obtain the total force caused by multiple current sources, the function in Eq. 15 needs to be determined. For this purpose, simulations are performed for one electromagnet stimulated by current from −20 A to 20 A and the corresponding forces are measured at 1 cm and 2 cm vertical distances. The force-current curves are shown in Fig. 10a. From the curves we can see that the force is nonlinear with the current in a large range due to the nonlinearity of the electromagnet core, however, for the current range [−4 A, 4 A] 305

Q. Zhang et al.: Magnetic Field Control for Haptic Display

where Mtrans is a N -by-N transformation matrix and Fz , I and b are N -by-1 column vectors. To validate this model, simulations with two electromagnets are performed and the repulsive magnetic force along the x-axis of both elements at the distance of 2 cm above are calculated using Eq. 13 which are depicted in Fig. 11.

FIGURE 10. Near linear interval of current-force curve. (a) Force-current curve for one electromagnet with current from −20A to 20A. (b) Forcecurrent curve fitting at 1cm above with current from −4A to 4A.

the curve is near linear. Note that the forces produced when there is no current are not zero due to the permanent magnet concentrator. To verify this assumption, polynomial curve fittings with a degree of one are performed using the method of least squares. The fitting results are shown together with the residuals in Fig. 10b. For the precision of the fittings, statistics of the fittings including the 95% confidence intervals, mean square errors are calculated and listed in Table 2. TABLE 2. Force-current curve fitting results.

FIGURE 11. Force-current curve of two elements at 2cm above. (a) F1 with I1 = 1A and I2 = ±1A. (b) F2 with I1 = 1A and I2 = ±1A.

The vertical forces at 2 cm above the electromagnets are simulated with I1 1 A and I2 ±1 A. The calculated forces are compared to the simulation of only one element. From Fig. 11 we can see each of the forces depend on both current values, and the resulting forces are a linear combination of both currents, which is consistent with our derivation. B. RECURSIVE LEAST SQUARES (RLS) ADAPTIVE CONTROL

From the fitting results we can see that the linear curve fitting error is negligible, thus the magnetic force is almost proportional to the current value within the range [−4 A, 4 A], and as is mentioned before, the current capacity of the wire is [−1 A, 1 A], which falls in the nearlinear range. Thus the relationship between the vertical force on the ith disk above the element caused by the jth electromagnet can be modeled with a linear function Fi,j,z = αij Ij + βij ,

i, j = 1, 2, 3 . . . N , Ij ∈ [−1 A, 1 A] (18)

where αij and βij are coefficients depending on the locations of source electromagnet and the magnet disk. Finally, by substituting Eq. 18 into Eq. 17 the total force caused by all the electromagnet on the ith magnet disk can be modeled as Fi,z =

N X j=1

(αij Ij + βij ) =

N X

αij Ij +

j=1

N X

βij

(19)

j=1

To estimate the parameters of the unknown multi-inputmulti-output linear system, adaptive filtering is adopted as it adjust its parameters according to an adaptive algorithm to recursively finds the filter coefficients. Recursive Least Squares (RLS) method [30] which aims to minimize a weighted least squares cost function is applied for its fast convergence. 1) PROBLEM FORMULATION

The detrended input currents and the output vertical forces are denoted as u and d which are N × 1 vectors. The system coefficients need to be estimated are the parameters of the ˆ trans which is of size N × N. The transformation matrix M RLS problem formulation in our work is • Given: observations u(n) and d(n), n = 0, 1 . . . r, where r is the total simulation steps ˆ trans that minimizes • Determine: M n X E(n) = λn−i ke(i)k2 (21) i=0

where 0 < λ ≤ 1 is called forgetting factor which gives older errors exponentially less weight and ˆ ˆ trans (n)u(n) is the error e(n) = d(n) − d(n) = d(n) − M vector at simulation step n. 2) RLS SOLUTION DERIVATION

Remark 2: The linear transformation in Eq. 12 can be further simplified using matrix expressions as Fz = Mtrans I + b 306

(20)

ˆ trans yields the solution Minimization of E with respect to M according to Wiener-Hopf equations [31] ˆ trans (n) = D(n) R(n)M (22) VOLUME 4, 2016

Q. Zhang et al.: Magnetic Field Control for Haptic Display

where R(n) is the deterministic correlation matrix for input: R(n) ,

r X

λn−i u(i)uT (i)

(23)

i=n

and D(n) is the cross-correlation vector between the input and desired output: D(n) ,

n X

λn−i u(i)dT (i)

(24)

i=0

Note that the update rules of R(n) and D(n) can be derived from Eq. 23 and Eq.24 as R(n) = λR(n − 1) + u(n)uT (n)

(25)

D(n) = λD(n − 1) + u(n)d (n)

(26)

T

Thus the solution

3) SUMMARY OF RLS ALGORITHM

ˆ trans (n − 1) = R−1 (n − 1)D(n − 1) M

(27)

can be updated as ˆ trans (n) = R−1 (n)D(n) M

(28)

ˆ trans (n), we check the old To find the update rule for M ˆ coefficients Mtrans (n − 1): ˆ trans (n−1) = D(n) − u(n)(d(n) − M ˆ T (n−1)u(n)) R(n)M trans (29) Thus we have ˆ trans (n) = M ˆ trans (n − 1) + R−1 (n)u(n)ξ (n) M

(30)

ˆ T (n − 1)u(n) is called the a priori where ξ (n) , d(n) − M trans ˆ trans (n − 1) is error which is the error we have at time n if M used. Define the gain vector as K(n) , P(n)u(n)

(31)

where P(n) , R−1 (n). Then Eq. 30 becomes ˆ trans (n) = M ˆ trans (n − 1) + K(n)ξ (n) M

(32)

By substituting Eq. 25 into Eq. 31 we can get K(n) =

λ−1 P(n − 1)u(n) 1 + λ−1 uT (n)P(n − 1)u(n)

(33)

To avoid matrix inversion, the Matrix Inversion Lemma [32] is applied: P(n) = λ−1 P(n − 1) − λ−1 K(n)uT (n)P(n − 1)

(34)

For initialization, normally P(0) is set to be a diagonal matrix with the non-zero entries being δ  1. VOLUME 4, 2016

Algorithm 1 Recursive Least Squares (RLS) Algorithm Input: u(n), d(n) ˆ trans Output: M 1: Initialization: δ ← 100, λ ← 0.99999 ˆ trans (0) ← 0 2: P(0) ← δI, M 3: for n ← 1 to maximum iteration do −1 4: K(n) ← 1+λλ−1 uP(n−1)u(n) T (n)P(n−1)u(n) ˆ T (n − 1)u(n) 5: ξ (n) ← d(n) − M trans ˆ trans (n) ← M ˆ trans (n − 1) + K(n)ξ (n) 6: M 7: P(n) ← λ−1 P(n − 1) − λ−1 K(n)uT (n)P(n − 1) ˆ T (n)u(n) 8: e(n) ← d(n) − M trans 9: if e(n) < threshold then 10: break 11: end if 12: end for

The RLS adaptive control algorithm is summarized in Algorithm 1. The algorithm is known for its fast convergence while its computational complexity is relatively high. To accelerate the convergence of the algorithm, a good first guess is important. Here for initialization we use the fitting result for one element in Fig. 10b where the interference is not considered in, which means the first guess for the transformation matrix is set to be an diagonal matrix whose entries along the diagonal are all 1.065. V. SIMULATION AND RESULTS A. SIMULATION SETUP

The Finite Element Method simulations are performed with COMSOL Multiphysics 4.4 which is a powerful FEM simulation software for modeling and solving engineering problems with coupled physics. In our simulations, two coupled physics are used which are Magnetic Field module for solving the magnetic field potential and Coefficient Form PDE module for obtaining the partial derivatives of the magnetic flux density used in force calculation. In addition, MATLAB scripting environment is integrated with Multiphysics via the LiveLink, which provides additional flexibility in generating, running and processing the model in MATLAB while applying its functions. Therefore the whole modeling and simulating process can be done in MATLAB from scratch. Our models are built directly in 3D and is under magnetostatic analysis framework where the current is considered as stable. The 3D model is meshed with tetrahedra elements. To solve the spatial magnetic field, all the field equations are set. B. PERFORMANCE OF THE DESIGNED ELECTROMAGNET MODEL 1) CONCENTRATION OF MAGNETIC FIELD

The magnetic field distributions of the electromagnet models depicted in Fig. 4 when stimulated by 1 A DC are shown in Fig. 12. 307

Q. Zhang et al.: Magnetic Field Control for Haptic Display

and is depicted in Fig. 14. The degree of stiffness is based on the study of human hand stiffness perception [33] where objects with stiffness higher than 650 N mm−1 are perceived as very firm. When having stiffness between 650 N mm−1 and 350 N mm−1 it feels like touching a piece of wood lying on foam which is considered firm. Stiffness lower than 350 N mm−1 is perceived as a soft object. FIGURE 12. Magnetic field distribution of the three models. (a) Magnetic field corresponds to Fig.4(a). (b) Magnetic field corresponds to Fig.4(b). (c) Magnetic field corresponds to Fig.4(c).

The basic electromagnet model produces uniformly distributed magnetic field as in Fig. 12a. The concentrating effect of the concentrator and the isolating effect of the isolator can be seen in Fig. 12b and Fig. 12c respectively. 2) STRENGTH OF FORCE

The vertical force strength exerted is compared among the basic model, the concentrated model and the complete model shown in Fig. 4. All of the three models have the same settings for the coil, i.e., the dimensions, the turns of wire and the current value. Simulations are done which scan the force generated along the x-axis at 1 cm and 2 cm above the array model. The current value of the three models is set as 1 A.

FIGURE 14. Axial stiffness of the designed model.

From Fig. 14 we can see that within the height of 0.9 cm, the magnetic force generated can be felt as very firm, and between 0.9 cm and 1.4 cm above the model a firm object sensation is produced, beyond which the force is considered soft. 4) OPERATION RANGE

FIGURE 13. Force produced by three models along x axis. (a) 1cm above. (b) 2cm above.

As can be seen in Fig. 13, the magnetic force forms a bell shaped field above the element. The force generated is augmented significantly by the introduction of the concentrator and the isolator in the complete model. The shielding effect can be observed in Fig. 13a where the bell generated by the designed model is thinner compared to the other two which indicates the magnetic field is bounded better. However the shielding effect becomes weaker with the vertical distance as can be seen in Fig. 13b where the bell shapes become wider compared to that in Fig. 13a which means the magnetic force field becomes smoother. The range of vertical force generated by a designed model is between [1N, 3.2N] at a vertical distance of 1 cm and [0.25N, 0.7N] at 2 cm above. 3) STIFFNESS OF FORCE

The stiffness of the axial magnetic force Kz can then be calculated using ∂ Kz = − Fz (35) ∂z 308

To determine the range of vertical distance the force can be perceived, data regarding human hand thresholds for static force are required. According to researched literature, for fingertips, the minimum force that can be perceived is 0.8mN, and for the palm, the threshold is 1.5mN [34]. Thus the vertical force generated by the designed model with a 1 A current at different vertical distances from 0 cm to 5.2 cm with a 0.1 cm interval are simulated and depicted in Fig. 15 with the results in the range of [4 cm, 5.2 cm] zoomed in. Magnetic force decays with the cube of the distance, which can be seen in Fig. 15. Based on the human hand threshold for static force, we find the farthest vertical distance where the force can be felt by both palm and fingertip is about 5.1 cm. C. ELECTROMAGNET ARRAY CONTROL EXPERIMENTS ON A 15-BY-15 MODEL 1) TRANSFORMATION MATRIX ESTIMATION BY RLS ALGORITHM

An experiment on a 15-by-15 electromagnet array model aiming to produce the force field of a human face is done. To first estimate the transformation matrix of this array, RLS algorithm summarized in Algorithm 1 is implemented and 450 iterations are done with the 15-by-15 model. In the adaptive control process, for each iteration we feed a set VOLUME 4, 2016

Q. Zhang et al.: Magnetic Field Control for Haptic Display

FIGURE 15. Force produced by designed model 1A current at different vertical distances above.

of detrended current value u(n) into the model and extract the force values d(n) produced. The algorithm updates the estimated transformation matrix to reduce the error between ˆ ˆ trans (n)u(n) and the force the force calculated by d(n) = M generated d(n) using FEM simulation. The first 226 current sets are showing in Table 3. TABLE 3. First 226 current sets used in RLS for 15-by-15 model.

FIGURE 16. Convergence analysis for a 15-by-15 model. (a) Convergence of current values for a 15-by-15 model. (b) Convergence of force values for a 15-by-15 model.

2) HUMAN FACE RENDERING

Note that we need to detrend the current and force record matrices by removing the sample mean values before performing the RLS algorithm. Since the current range is set to be [−1 A, 1 A], the mean value vector for each current is then a zero vector. The mean values for each force measurement are estimated by averaging the force measurement Fi,z (i = 1, 2 . . . 225) over the cases where Ii = ±1 A. The updated transformation matrix results in the convergence of current and force values are shown in Fig. 16. In Fig. 16a and Fig. 16b, the horizontal axis is the iteration number and the vertical axis indicates the error between the desired values and the estimated values of the present iteration. We can see from the convergence curves that the errors at the beginning are small owning to the good first guess, and after 225 iterations they converge and stay stable afterwards. VOLUME 4, 2016

Now having the estimated transformation matrix, we test it with human face volumetric data which is collected from 3D scanner. The original data is stored in the form of point cloud data and is saved as a PLY file, then some preprocessing techniques are applied on the data to clean it up and down sample the large scale data into a STL file which can be converted in to a text file for MATLAB to read in. The processed face surface data is shown in Fig. 17a. Since the original data set is not gridded, quantization is done first to convert the (X,Y) coordinates into a 15-by-15 gridded mesh, which results in multiple repeated data with the same (X,Y) coordinate, removal of repeated data with the same (X,Y) coordinate are done by selecting the data with the lowest Z component. Finally, the Z components of the pre-processed face surface data are scaled to the range of [0.8, 2.6] to fit the produced force range and all the spots in the grid without data are filled with 0.8 which is the minimum force value. Thus the desired force matrix is obtained and visualized in Fig. 17b. The corresponding estimated current configuration matrix is then calculated and is depicted using colour map in Fig. 17c. The mean square error (MSE) of the force values at the control points is 0.00001 which is negligible. The final force field generated at different distance above the array are shown in Fig. 17d, Fig. 17e and Fig. 17f from which we can see that the force field becomes weaker but smoother with the distance increasing. The force produced from 1 cm above ranges from 0.8N to 2.6N while from 2 cm above it is from 0.11N to 0.49N and from 3 cm above it is from 0.02N to 0.18N. The force field generated from 1 cm is full of bumps and holes 309

Q. Zhang et al.: Magnetic Field Control for Haptic Display

FIGURE 17. Human face rendering test by 15-by-15 electromagnet array. (a) Original scattered face data. (b) Desired force matrix. (c) Current configuration matrix. (d) Force field from 1cm above. (e) Force field from 2cm above. (f) Force field from 3cm above.

while the fluctuation is alleviated in the force field from 2 cm above, and from 3 cm above the force field becomes smooth. VI. CONCLUSION

In this paper we propose a magnetic rendering system for rendering volumetric shapes in mid-air for users to perceive with attached magnets on the hands. In order to generate strong enough magnetic field for haptic rendering, we design a new electromagnet model with a magnetic field concentrator and an isolator. The model demonstrates desirable performance in terms of the magnitude of the vertical force produced, the magnetic field isolation property, and stiffness with high fidelity via FEM simulations. The electromagnet array is configured and its linearity is characterized, which facilitates control in adaptive approach. FEM simulation experiments are performed on 15-by-15 model. A set of pre-processed 3D scanned face data is used to test the force field generation performance using RLS algorithm, which converges to results with insignificant error on the control points, and the force field produced is smooth. REFERENCES [1] A. El Saddik, M. Orozco, M. Eid, and J. Cha, Haptics Technologies: Bringing Touch to Multimedia. New York, NY, USA: Springer, 2011. [2] R. J. Adams and B. Hannaford, ‘‘Stable haptic interaction with virtual environments,’’ IEEE Trans. Robot. Autom., vol. 15, no. 3, pp. 465–474, Jun. 1999. 310

[3] T. R. Coles, D. Meglan, and N. John, ‘‘The role of haptics in medical training simulators: A survey of the state of the art,’’ IEEE Trans. Haptics, vol. 4, no. 1, pp. 51–66, Jan./Feb. 2011. [4] A. El Saddik, ‘‘The potential of haptics technologies,’’ IEEE Instrum. Meas. Mag., vol. 10, no. 1, pp. 10–17, Feb. 2007. [5] H. Dong, Y. Gao, H. Al Osman, and A. El Saddik, ‘‘Development of a Web-based haptic authoring tool for multimedia applications,’’ in Proc. IEEE Int. Symp. Multimedia, 2015. [6] Geomagic Touch X Haptic Device. [Online]. Available: http:// www.geomagic.com/en/products/phantom-desktop/overview, accessed Jan. 9, 2016. [7] P. Berkelman, ‘‘A novel coil configuration to extend the motion range of Lorentz force magnetic levitation devices for haptic interaction,’’ in Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst., Oct./Nov. 2007, pp. 2107–2112. [8] R. Traylor and H. Z. Tan, ‘‘Development of a wearable haptic display for situation awareness in altered-gravity environment: Some initial findings,’’ in Proc. 10th IEEE Symp. Haptic Int. Virtual Environ. Teleoper. Syst., Mar. 2002, pp. 159–164. [9] H. Kim, C. Seo, J. Lee, J. Ryu, S. Yu, and S. Lee, ‘‘Vibrotactile display for driving safety information,’’ in Proc. IEEE Intell. Transp. Syst. Conf., Sep. 2006, pp. 573–577. [10] K. Minamizawa, D. Prattichizzo, and S. Tachi, ‘‘Simplified design of haptic display by extending one-point kinesthetic feedback to multipoint tactile feedback,’’ in Proc. IEEE Haptics Symp., Mar. 2010, pp. 257–260. [11] D. Prattichizzo, F. Chinello, C. Pacchierotti, and M. Malvezzi, ‘‘Towards wearability in fingertip haptics: A 3-DoF wearable device for cutaneous force feedback,’’ IEEE Trans. Haptics, vol. 6, no. 4, pp. 506–516, Oct./Dec. 2013. [12] F. Arafsha, K. M. Alam, and A. El Saddik, ‘‘Design and development of a user centric affective haptic jacket,’’ Multimedia Tools Appl., vol. 74, no. 9, pp. 3035–3052, 2013. [13] Y. Suzuki and M. Kobayashi, ‘‘Air jet driven force feedback in virtual reality,’’ IEEE Trans. Comput. Graph. Appl., vol. 25, no. 1, pp. 44–47, Jan./Feb. 2005. VOLUME 4, 2016

Q. Zhang et al.: Magnetic Field Control for Haptic Display

[14] R. Sodhi, I. Poupyrev, M. Glisson, and A. Israr, ‘‘Aireal: Interactive tactile experiences in free air,’’ ACM Trans. Graph., vol. 32, no. 4, 2013, Art. ID 134. [15] S. Gupta, D. Morris, S. N. Patel, and D. Tan, ‘‘AirWave: Non-contact haptic feedback using air vortex rings,’’ in Proc. ACM Int. Joint Conf. Pervasive Ubiquitous Comput., 2013, pp. 419–428. [16] T. Hoshi, M. Takahashi, T. Iwamoto, and H. Shinoda, ‘‘Noncontact tactile display based on radiation pressure of airborne ultrasound,’’ IEEE Trans. Haptics, vol. 3, no. 3, pp. 155–165, Jul./Sep. 2010. [17] B. Long, S. A. Seah, T. Carter, and S. Subramanian, ‘‘Rendering volumetric haptic shapes in mid-air using ultrasound,’’ ACM Trans. Graph., vol. 33, no. 6, 2014, Art. ID 181. [18] M. Weiss, C. Wacharamanotham, S. Voelker, and J. Borchers, ‘‘FingerFlux: Near-surface haptic feedback on tabletops,’’ in Proc. 24th ACM Annu. Symp. User Interface Softw. Technol., 2011, pp. 615–620. [19] P. Berkelman, M. Miyasaka, and J. Anderson, ‘‘Co-located 3D graphic and haptic display using electromagnetic levitation,’’ in Proc. IEEE Haptics Symp., Mar. 2012, pp. 77–81. [20] L. Zhang, H. Dong, and A. El Saddik, ‘‘From 3D sensing to printing: A survey,’’ ACM Trans. Multimedia Comput., Commun. Appl., vol. 12, no. 2, 2015, Art. ID 27. [21] L. Yang, H. Dong, A. Alelaiwi, and A. El Saddik, ‘‘See in 3D: State of the art of 3D display technologies,’’ Multimedia Tools Appl., Oct. 2015. [22] J. W. Jansen, C. M. M. van Lierop, E. A. Lomonova, and A. J. A. Vandenput, ‘‘Modeling of magnetically levitated planar actuators with moving magnets,’’ IEEE Trans. Magn., vol. 43, no. 1, pp. 15–25, Jan. 2007. [23] J. C. Compter, ‘‘A planar motor with electro-dynamic propulsion and levitation under 6-DoF control,’’ in Proc. 4th Int. Symp. Linear Drives Ind. Appl., 2003, pp. 1–4. [24] J. C. Mallinson, ‘‘One-sided fluxes—A magnetic curiosity?’’ IEEE Trans. Magn., vol. 9, no. 4, pp. 678–682, Dec. 1973. [25] P. Sankar, ‘‘Magnetic arrays with increased magnetic flux,’’ U.S. Patent 8 514 047, Aug. 20, 2013. [26] Standard Specification for Standard Nominal Diameters and CrossSectional Areas of AWG Sizes of Solid Round Wires Used as Electrical Conductors, ASTM Standard B258-14, 2008. [27] S. V. Vonsovskiy, Magnetizm (Magnetism). Moscow, Russia: Nauka, 1971. [28] S. Celozzi, G. Lovat, and R. Araneo, Electromagnetic Shielding. New York, NY, USA: Wiley, 2008. [29] S. Sanz, L. Garcia-Tabares, I. Moya, D. Obradors, and F. Toral, ‘‘Evaluation of magnetic forces in permanent magnets,’’ IEEE Trans. Appl. Supercond., vol. 20, no. 3, pp. 846–850, Jun. 2010. [30] S. Haykin, Adaptive Filter Theory. New York, NY, USA: Pearson Education, 2007. [31] A. Leonard and T. W. Mullikin, ‘‘Integral equations with difference kernels on finite intervals,’’ Trans. Amer. Math. Soc., vol. 116, pp. 465–473, Apr. 1965. [32] A. S. Householder, The Theory of Matrices in Numerical Analysis. New York, NY, USA: Blaisdell, 1964. [33] K. J. Kuchenbecker, J. Fiene, and G. Niemeyer, ‘‘Improving contact realism through event-based haptic feedback,’’ IEEE Trans. Vis. Comput. Graphics, vol. 12, no. 2, pp. 219–230, Mar. 2006. [34] G. C. Burdea, Force and Touch Feedback for Virtual Reality. New York, NY, USA: Wiley, 1996.

VOLUME 4, 2016

QI ZHANG received the B.Eng. degree in communications engineering and the B.A. degree in advertising from Xiamen University, China, in 2012, and the M.Sc. degree in electrical and computer engineering from the University of Ottawa, in 2015. She is currently with Cisco Systems, Inc. She also has a fondness for a finite element method modeling and adaptive signal processing. Her academic interests include haptic interface design, multimedia computing, and control.

HAIWEI DONG (M’12) received the M.Eng. degree in control theory and control engineering from Shanghai Jiao Tong University, China, in 2008, and the D.Eng. degree in computer science and systems engineering from Kobe University, Japan, in 2010. He is currently with the University of Ottawa. Before that, he was appointed as a Post-Doctoral Fellow with New York University Abu Dhabi, a Research Associate with the University of Toronto, a Research Fellow (PD) of the Japan Society for the Promotion of Science, a Science Technology Researcher with Kobe University, and a Science Promotion Researcher with the Kobe Biotechnology Research and Human Resource Development Center. His research interests include robotics, haptics, control, and multimedia. He is a member of ACM.

ABDULMOTALEB EL SADDIK (M’01–SM’04– F’09) is the University Research Chair and Distinguished Professor with the School of Electrical Engineering and Computer Science, University of Ottawa. He is an Internationally-Recognized Scholar who has made strong contributions to the knowledge and understanding of multimedia computing, communications, and applications. He has authored or co-authored four books and over 450 publications. He chaired more than 40 conferences and workshops, and has received research grants and contracts totaling more than U.S. $18 million. He has supervised over 100 researchers. He received several international awards among others are the IEEE Canada Computer Medal. He is an ACM Distinguished Scientist and a fellow of the Engineering Institute of Canada and the Canadian Academy of Engineers.

311