A linearization-based method of simultaneous position and velocity measurement using ultrasonic waves

Sensors and Actuators A 233 (2015) 480–499

Contents lists available at ScienceDirect

Sensors and Actuators A: Physical journal homepage: www.elsevier.com/locate/sna

A linearization-based method of simultaneous position and velocity measurement using ultrasonic waves Natee Thong-un a,∗ , Shinnosuke Hirata b , Yuichiro Orino c , Minoru Kuribayashi Kurosawa a a

Department of Information Processing, Tokyo Institute of Technology, 4259 Nagatsuta, Midori, Yokohama 226-8502, Japan Department of Mechanical and Control Engineering, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro, Tokyo 152-8552, Japan c Collaborative Research Center, The University of Shiga Prefecture, 2500 Hassaka-cho, Hikone-shi, Shinga 522-8533, Japan b

a r t i c l e

i n f o

Article history: Received 20 January 2015 Received in revised form 15 July 2015 Accepted 26 July 2015 Available online 12 August 2015 Keywords: Echolocation Linear-period-modulated signal One-bit signal processing Three-dimensional position and velocity measurements

a b s t r a c t This paper proposes a novel system for three dimensional – position and velocity measurements with echolocation. It is improved from the former systems, based on low computational one-bit signal processing, to have wider detectability from one and two dimensions to three dimensions with the rearrangement of the ultrasonic receiver positions. The proposed method can measure the object location from nonlinear equations without using iterative methods. This method can concurrently measure position and velocity in three dimensional spaces with one-time measurement. The sound radiation and reflection studies using echolocation are also included in this paper. The additive noise: White Gaussian Noise, Colored Gaussian Noise, General Gaussian Noise, Laplacian Noise, and Gaussian Mixture Noise, is considered for the performance study. The system consists of one sound transmitter and three receivers. These devices are simple sound packages to support low-cost applications. In addition, to satisfy the requirement for a low-computation cost, an echo is converted into a one-bit signal by a three-channel delta–sigma-modulation board. Then, FPGA is used to determine recursive cross correlation based on one-bit computation. The object location is then defined with x–y–z coordinates. The object’s velocity is computed using linear-period-modulated signals. The validity and repeatability of experimental results are evaluated using statistics. © 2015 Elsevier B.V. All rights reserved.

1. Introduction An echolocation concept is a simple technique of object navigations by determining time-of-flight (TOF), which is calculated using an echo from a target. This is a remarkable ability of bats and dolphins to identify objects, that are in front of them, by transmitting a sound wave and receiving an echoed reflection by sensory receptors. Bats and dolphins can sense the distance from their own position to obstacles. Accordingly, echolocation is widely used in many ranging measurements applications, such as radar, sonar, and other acoustic systems [1–3]. Acoustic systems have the benefits of being low in cost and, small in size, and using relatively simple hardware compared to other systems; they also tend to be the most up-to-date and reliable. Robotics studies using acoustic systems for localization problems have been proposed [4,5]. The most

∗ Corresponding author. E-mail addresses: [email protected] (N. Thong-un), [email protected] (S. Hirata), orino.y@office.usp.ac.jp (Y. Orino), [email protected] (M.K. Kurosawa). http://dx.doi.org/10.1016/j.sna.2015.07.029 0924-4247/© 2015 Elsevier B.V. All rights reserved.

critical problem for robot navigation is to determine the location of obstacles. The determination of an obstacle’s position must be as reliable as possible to satisfy the self-localization procedure. Currently, advanced devices based on laser ranging and a vision finder must scan in a single plain, and are relatively expensive; moreover, certain targets cannot be located accurately [6]. Three dimensional range finders are also available, but their measurement procedures are time consuming [7,8]. Thus, it is not appropriate to apply laser range finders for real-time robot navigation because the laser range finders have very narrow beam when relatively compared with ultrasound when we use a laser pointer in 3D space, it needs a long-time scanning of the position identification to cover 3D space. It is too slow if used for a moving object. When considering vision systems, a primary weak point is computationally expensive methods and high cost [9]. Three-dimensional-positioning methods based on echolocation have already been proposed [10,11]; however, these research papers have simply studied the vision of bats and require significant computational time. In the case of concurrently position and velocity measurement, the ambiguity function method, which is well known in radar technology [1], is used to

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

detect the stationary and moving targets in various target detection scenarios. A delay and a Doppler shift of a chirp signal are computed by means of the cross ambiguity function. A peak can be located on both the delay axis and the Doppler shift axis. However, this method requires Fourier transform of lag products and consists of huge iterations of multiplications and accumulations when relatively compared with the general cross correlation [12]. In the field of acoustic, it can play a vital role in the positioning techniques in both two and three dimensions including sound-source localization, self-localization system, sensor networks, and robot navigation, for examples. Initially, 2-D simple measurement systems composing of two microphones mounted at the left and right ears for video conference and robotics have been realized [13,14]. 3-D sound source localization systems using four microphones have been developed [15,16]. Then, a 3-D localization system for simultaneous position and velocity measurements of a moving sound source has been described using four narrow-band ultrasonic receivers [17]. In the case of sensor networks, a system enabling both 2-D position and temperature measurements by means of speed difference between light and supersonic has been proposed [18] 2-D and 3-D communication architectures for under water acoustic sensor networks have been presented to support the robustness of the sensor network to node failures [19]. A grid-based method for the location of multiple sound sources in an acoustic sensor network containing a microphone array has also been proposed [20]. In addition, an acoustic ranging system for underwater localization of robotic fish involving simple hardware: a pair of monotone buzzer and microphone has been accomplished [21]. Several research groups as above papers are examples of noteworthy studies in acoustic positioning system. However, it seems that we lack a system-based echolocation for concurrent position and velocity measurements in three dimensional spaces. Therefore, it is of interest to explore a method to use echolocation for robot navigation with both simple computation and low-cost devices. One-bit stream signal processing is well known as a digital decoder of Super Audio CD (SACD) with Direct Stream Digital (DSD) Technology. This processing allows SACD to achieve high audio quality that is better than any other digital or analog technology [22] because quantization noise of Analog to Digital conversion (ADC) can be relatively reduced inside a signal band by delta–sigma modulation. When considering an oversampling technique, cross correlation with accuracy and resolution based on one-bit signal processing has been proposed for ultrasonic distance measurement [23]. This method can compute TOF with low computational costs using one-bit signal processing-based a linear-frequency-modulated (LFM) signal. However, the important problem of LFM signals cannot not exactly correlate with an echo of a moving object due to the Doppler effect. To solve this issue, a linear-period-modulated (LPM) signal, a period of which is linearly swept with time, has been proposed for ultrasonic systems identifying moving objects [24,25]. Although these methods can provide cross-correlation functionality to calculate the TOF from the reflected echo, they require significant computation time due to envelope-signal computations. Therefore, a low cost and highresolution method of ultrasonic distance measurement using pulse compression with two cycles of LPM signals and Doppler-shift compensation has been presented [26]. In the case of an ultrasonic two-dimensional system relied on one-bit signal processing, a solution for position and velocity measurements using two ultrasonic receivers has also been proposed [27]. However, this system uses expensive acoustical receivers from the B&K Company, which are often used in calibration applications; this is inconvenient for the development of larger systems. Instead, to support low-cost applications, a system using a silicon MEMS microphone as an acoustical receiver, which is much cheaper than the acoustical microphone in the previous version,

481

Fig. 1. Sound source used in the proposed system: (a) loudspeaker model PT-R4; and (b) radiation of a rectangular surface source.

is used in the proposed three-dimensional ultrasonic positioning system [28]. These systems can estimate an object’s position as a distance from the sound source, the azimuth angle, and the elevation angle; however, they cannot measure the object’s velocity. Therefore, this paper has been developed from the earlier works for a new three dimensional ultrasonic position and velocity measurement by reducing the number of the acoustical receivers and uses a general home stereo as a sound source to support low-cost application. As claimed, recursive cross correlation based on onebit signal processing and Doppler-shift compensation is applied to achieve low computational times. The object position is calculated as x–y–z coordinates. Then, a delta–sigma modulation board and FPGA are adapted to support one-bit signal processing. The validity of the proposed system is confirmed by the repeatability of the experimental results. 2. Sound source and acoustic receiver 2.1. Sound source In the proposed system, a home stereo is used as a sound source that radiates sound directly into the air. The loudspeaker used is a direct-radiator Pioneer PT-R4. This speaker is often used in small public address system. The principle advantage of a direct-radiator speaker is a satisfactory response across a comparatively wide frequency range. However, the primary disadvantages of this type are low efficiency, narrow directivity pattern and frequently irregular response at high frequencies. To fully describe the loudspeaker, we need to know the intensity levels it produces at all frequencies of interest. The loudspeaker used has a rectangular plane diaphragm [29], which is shown in Fig. 1(a). The sound-pressure direction of the rectangular plane surface at a point P, which is placed far from the diaphragm at a distance R in Fig. 1(b), can be expressed by the directivity index Dϕ,g [30,31] as: D, (dB) = 10log10 |pT |

(1)

where pT is the total sound pressure. Thus, on the Z–Y plane; DZ−Y (dB) = 10log10 |

iω0 uT 4absinc(ka cos )| 2R

(2)

and on the X–Y plane DX−Y (dB) = 10log10 |

iω0 uT 4absinc(kb cos )|. 2R

(3)

Fig. 2(a) and (b) shows the directivity intensity of the loudspeaker used for a rectangular plane with a = 25 mm. and b = 3 mm. In this state, uT is the total amplitude of small elements on the diaphragm area 4ab, the frequency  is 50 kHz, the density of the air 0 is 1.20 kg/m3 , the amplitude of vibration on the plane is 2 m, and the wave number k is /v, where the speed of sound propaga-

482

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

Fig. 2. Directivity beam of the loudspeaker at 50 kHz: (a) sound intensity on the X–Y plane; (b) sound intensity on the Z–Y plane (c) actual sound intensity on the X–Y plane (d) actual sound intensity on Z–Y plane.

tion v = 345 m/s. In Fig. 2(a) and (b), this analysis can describe the radiation of the loudspeaker and if it is sufficient for detecting an object. The loudspeaker at the normal position can radiate a sound beam to the object on the Z–Y plane, at approximately ±20◦ and on the X–Y plane, at ±45◦ about the +Y axis. Because pulse compression is used to sweep from 50 to 20 kHz, the area of vertical detection expands; this system can be designed to scan an object’s position by moving the loudspeaker position up and down along to the Z-axis. To realize the actual sound radiation, sound amplitude on ground at many points in front of a direct-radiator is evaluated experimentally. Fig. 2(c) and (d) shows the directivity intensity of the loudspeaker used in the controlled temperature room. 2.2. Acoustical receiver Acoustical receivers are electrostatic transducers for transforming acoustic energy to electric energy. In general, a standard acoustic receiver for performing accuracy measurements is expensive [24,26,27]. This is an important limitation of producing measurements that require using more sensors. An acoustical receiver must be cared for to avoid damage, and thus, it is not convenient for use outdoors or in robust applications. Therefore, an acoustical receiver, that has a low cost, a small size, high reliability and robustness is required. A Knowles Acoustic model SPM0204UD5 silicon microphone that is primary used in ultrasonic sensors is a potential candidate. This model has the ability to sense sound pressure or particle velocity in all directions (i.e., omnidirectional). Frequencies can be 10 kHz up to 100 kHz. Thus, the sensor is selected and embedded into a signal processing board with a low pass frequency circuit with a 60 kHz frequency cutoff and a preamplifier of 20 dB. Sensitivity has a minimum level at – 47 dB when the humidity does not exceed 70 R.H. Due to the pulse compression, an LPM signal, is applied to perform position measurements; frequency- distortion measured by the acoustic sensors must be examined with a frequency-distortion check to maintain system

reliability. In Fig. 3, a frequency-distortion check of the LPM signal at frequencies of 50 kHz to 20 kHz is shown. The results show that the transducer can measure the original signal, as referenced in [25]. 3. Proposed three-dimensional position and velocity measurements 3.1. Cross correlation using one-bit signal processing and Doppler—shift compensation to determine time-of-flight The recursive cross-correlation used in this paper is shown in Fig. 4. The transmitted signal that excites to the loudspeaker is a pair of two LPM signals. An echo from the moving object is incident on the receivers and converted into a one-bit stream by three delta–sigma modulators. Conversely, a reference signal utilizes a digital comparator as a signal converter. Then, both one-bit signals are correlated to compute TOF using the recursive cross-correlation function. This algorithm is described in [23] and is critical to reducing computational costs. When no Doppler effect is present, the TOF of a received LPM signal is generally estimated from the maximum peak time in the cross-correlation function. However, the cross-correlation function is shifted due to the Doppler effect in the case of a moving object. Accordingly, TOF is estimated by performing a Doppler-shift compensation between the maximum and minimum peaks of the cross-correlation function as shown in [26]. 3.2. Three-dimensional position measurement In Fig. 5, the designed system can detect an object with ranging measurements from the −X to +X axis and from the −Z to +Z axis. The object appears in front of the loudspeaker in the +Y axis. We assume that the position of the three microphones can be placed anywhere on the X–Z plane. Their position are (x1 , 0, z1 ), (x2 , 0, z2 ),

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

 A=

⎡

−x1 −x2 −x3

−z1 −z2 −z3

TOF21 v2 − x12 − z12

483



vTOF1 vTOF2 , b = ⎤ vTOF3

  x z d

, andc =

⎣ TOF22 v2 − x2 − z2 ⎦where b is an unknown vector. The matrix 2 2 TOF23 v2 − x32 − z32 A consists of any locations of the microphone position and the TOFs. In practice, the matrix A might be a deficient matrix from these values if its determinant is zero. Thus, this system should be designed to avoid the singular-matrix problem with det. ( A) = 0. A singularity of the matrix A will not be produced as long as the three microphones do not lie on a line. In addition, Eq. (5) is concerned with an attainable right-hand side c but also with the solutions b that attain them. With c = 0, the solution b = 0 is always valid; however, in this case, there may be infinitely many solutions; this is the null space of matrix A. However, because |v2 TOF2i | is always much more than | − xi2 − zi2 |, the null space of matrix A does not exist in the case c = 0. Next, assuming that the system 2 Ab = c is transformed to be 2 Tb = e, where T is an upper triangular matrix, then:

⎡ ⎢

−x1

T =⎣ 0 0

vTOF1

a1 z1 − z2

v (TOF2 − a1 TOF1 )

0

v ((a1 a3 − a2 ) TOF1 − a3 TOF2 + TOF3 )

and

⎡

⎥ ⎦,

⎤

c(1, 1)

e=⎣

⎤

−z1

⎦

c(2, 1) − a1 c(1, 1) (a1 a3 − a2 ) c (1, 1) − a3 c (2, 1) + c (3, 1)

where c(1,1), c(2,1), and c(3,1) are the first, second, and third rows of the vector c, respectively, and the coefficient numbers are a1 = x2 x1 , a2

=

x3 x1 , anda3

x

=

−z3 + x3 z1 1

x

−z2 + x2 z1

. The null space of matrix A is the

1

same as the null space of T [32]. Thus, the unknown parameters d, z, and x can be solved, respectively, as: d=

0.5e(3, 1)

v ((a1 a3 − a2 ) TOF1 − a3 TOF2 + TOF3 )

z= Fig. 3. LPM pulse compression signal from 50 kHz to 20 kHz measured by a Knowles Acoustic model SPM0204UD5 silicon microphone in (a) frequency and (b) time domains.

and (x3 , 0, z3 ), and the object has an unknown position as a small point on (x, y, z). d is a distance between from the original point to the object. A simple equation of the microphone position, the object position, and the TOF is:

 d+

(x − xi )2 + y2 + (z − zi )2 = vTOFi ,

(4)

where TOFi is the time-of-flight at microphone positions i, where i = 1,2, and 3, respectively. Thus, there are three nonlinear equations: we can reform (4) to be a linear form as: 2·A·b=c and:

(5)

x=

,

0.5e(2, 1) − vd(TOF2 − a1 TOF1 ) , and a1 z1 − z2

vdTOF1 − 0.5e(1, 1) − z1 z x1

(6)

where e(1,1), e(2,1), and e(3,1) are the first, second, and third rows of the vector e, respectively. Considering in the first pivot of the upper triangular matrix T . x1 of the first position is designed to not

be equal to zero. Then, y can be calculated as y = d2 − x 2 − z 2 . In the next step, to ensure matrix A is not singular and to measure a direct sound from the speaker, which has a strong radiation in a horizontal line from the speaker in Fig. 2(b), the first microphone position is set to be (x1 = 10 cm, y1 = 0 cm, z1 = −10 cm), the second microphone position is set to be (x2 = 0 cm, y2 = 0 cm, z2 = 10 cm), and the third microphone position is set to be (x3 = −10 cm, y3 = 0 cm, z3 = −10 cm), as shown as Fig. 6. det. A from the given microphone positions can be computed as 3.45TOF1 + 6.9TOF2 + 3.45TOF3 . det. A is a positive number in every case of the object position.

3.3. Three-dimensional velocity vector measurement For velocity measurements, moving-object velocity is estimated using the relative velocity measurement (vd ), the instantaneous

484

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

Fig. 4. Recursive cross-correlation of one-bit signal processing.

position (x, y, z) of the moving object, and the microphone positions (xi , yi , zi ). The relative velocity can be defined by the interval between the first peak and the second peak in the modulated crosscorrelation function with the Doppler-shift length (ld ), and those in the cross-correlation function of the reference LPM signal (l0 ), which are shown as [26]:

vd =

l0 − ld v. l0 + ld

(7)

We can computed ld and l0 from two peaks of the cross correlation function by directly using a pair of LPM signal, shown in Fig. 4. Based on the fundamentals of three-dimensional Doppler velocimetry [33], vector projection can be thought of as the dot product between the two vectors shown to Fig. 5. Hence, the results of vector projection between the moving-object direction (i.e., blue line) and the direction (i.e., red line) between the object and microphones are exactly the relative Doppler velocity estimation at each microphone. Assuming that an unknown vector of the moving object based on x, y, z co-ordinates is u = [ux , uy , uz ]T , and the relative velocity measurements at microphones 1, 2, and 3 are recorded as vd = [vd1 , vd2 , vd3 ]T , then, the projection is expressed as: vd = [

u T · p1 u T · p2 u T · p3 , , ] p1  p2  p3 

T

(8)

where pi at i = 1,2, and 3 is a vector [xi – x, −y, zi – z]T , that has a direction from the object position to the microphone position, and therefore, the unknown vector u can be expressed as: u = H −1 · W · g, where:

H=



x1 − x x2 − x x3 − x

−y −y −y

z1 − z z2 − z z3 − z



 ,W =

(9) 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3



,

andg = [vd1 · p1 , vd2 · p2 , vd3 · p3 ]T . The weighted averaging matrix W is inserted to reduce the variance of the relative velocity measurements between all microphones. 4. Evaluation of proposed system by computer simulation The proposed three dimensional positioning method obtained from one-bit signal processing was evaluated by computer simulations using MATLAB. The frequency of a pair of LPM signals varied from 50 to 20 kHz. The LPM signal range was 3.274 ms. An object

Fig. 5. Three-dimensional position and velocity measurements of a moving object.

position is assumed as a small point and no effect of the object’s surface attenuation. The propagation velocity of an ultrasonic wave in air was approximately 345.1 m/s at 22.4 ◦ C. The attenuation of ultrasonic wave in air was 2.11 dB/m [34]. The sampling frequencies of the one-bit delta–sigma modulated signal were 12.5 MHz. The cross-correlation function of the one-bit received signal and the one-bit reference signal was performed from the recursive cross correlation and the smoothing operation accomplished by the 141taptriangular weighted moving average filter. 4.1. Comparison with the Cramer-Rao lower bound In general, any unknown parameter estimators should be unbiased on the average, and the variance is smaller than that of the other estimators. The Cramer-Rao Lower Bound (CRLB) is the lowest possible variance that an unbiased estimator can obtain. First, we have to consider a model of the proposed method. A object in a three-dimension scenario used to determine the position is p = [x, y, z]T , and the microphone position si = [xi , 0, yi ]T i = 1, 2, and 3 is assumed known, as shown Fig. 5. As proposed in Eq. (4), we have three unknown parameter b = [x, z, d]T to be estimated and can consider a linear model of Eq. (4) as ri = |p-si | + d =



(p-si )T (p-si ) + d,

(10)

where ri is the distance between the speaker, object, and the microphone i. We will assume that the TOF measurements can be expressed by the additive noise model as r = ri + ri

(11)

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

485

Fig. 6. Distance errors from various noisy signals: (a) white Gaussian noise, (b) colored Gaussian noise, (c) general Gaussian noise,(d) Laplacian noise, (e) Gaussian mixture noise.

where r = [r1 , r2 , r3 ]T , and ri is noise from TOFi . It is zero mean  2 

1 0 0 and has covariance matrix Q = 0 22 0 . In the presence of 0 0 32 TOF noise, noisy values result in the equation error vector from Eq. (4)

 = c − Ab.

(12)

The least square (LS) solution of b that minimize  is b = (AT A)−1 AT c

(13)

The Probability Density Function (PDF) is viewed as a function of the unknown parameters as pi (ri ; b)



N−1

=

n=0

1 exp √ 2 2

−

 1 2 2

ri [n] −



(p − si )T (p − si ) − d

(14)

J = −E

∂ ln p(r; b) ∂b

T 

∂ ln p(r; b) ∂b



CRLB (b) = J

−1

=

∂r(b) ∂b

T Q

−1∂r(b)

−1

−1

(16)

∂b

The terms along the diagonal of inverse Fisher information matrix result in the CRLB bounds on the variances: var(ˆx) ≥

 

2 9 + 6A1 − 4A22 det(J)N

(17a)

var(ˆz ) ≥

 

2 9 + 6B1 − 4B22 det(J)N

(17b)

ˆ ≥ var(d)

 

2 4A1 B1 − 6(B1 + A1 + 2C1 ) − C12 det(J)N

(17c)

where N is all trials and

2 

Next, in the presence of noise, the variances cannot be less than those expressed by the CRLB. CRLB is equal to the inverse of the Fisher matrix defined as [35]



After taking natural logarithm and then differentiating this function, the CRLB for the underlying problem can be reduced to [36]

 (15)

A1 =



A2 =



B1 =



B2 =



C1 =



(z − z1 )2 (p-s1 )T (p-s1 ) (z − z1 ) (p-s1 )T (p-s1 ) (x − x1 )2

(p-s1 )T (p-s1 ) (x − x1 )

(p-s1 )T (p-s1 ) (z − z1 ) (x − x1 ) (p-s1 )T (p-s1 )

+



+

(z − z2 )2

+





+



+



+



+



+



+



+



(p-s2 )T (p-s2 ) (z − z2 ) (p-s2 )T (p-s2 ) (x − x2 )2

(p-s2 )T (p-s2 ) (x − x2 )

(p-s2 )T (p-s2 ) (z − z2 ) (x − x2 ) (p-s2 )T (p-s2 )

(z − z3 )2 (p-s3 )T (p-s3 ) (z − z3 ) (p-s3 )T (p-s3 ) (x − x3 )2

(p-s3 )T (p-s3 ) (x − x3 )

(p-s3 )T (p-s3 ) (z − z3 ) (x − x3 ) (p-s3 )T (p-s3 )

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

65.22 0.31 0.12 0.24 0.18

−41.51 65.14 −17.06 0.17

−41.43 0.11 −17.18 0.16 65.11 0.13 −41.47 0.13 −17.13 0.18 65.15 0.16 −41.44 0.08 −17.16 0.11 64.02 0.15 −41.28

0.07 0.05 0.13 −41.41 63.75 −16.88 −41.44 63.78 −16.84 Fourth position (x = −40 cm, y = 60 cm, z = −20 cm) −41.57 0.47 −40.94 1.00 x y 64.04 0.38 64.42 0.89 −16.99 0.35 −17.95 0.81 z

0.23 0.16 0.16

0.77 1.09 0.34 -49.65 75.48 −2.49 0.23 0.16 0.21 49.81 74.76 −2.76 0.10 0.06 0.14 0.23 −49.82 0.17 74.78 0.73 −2.73 -49.88 74.80 −2.99 0.23 0.15 0.11 −49.83 73.73 -2.72 0.12 0.08 0.14 −49.95 73.27 -2.71 −49.89 73.32 −2.68 Third position (x = −50 cm, y = 75 cm, z = 0 cm) −50.47 1.15 −49.96 0.48 x 74.85 1.13 73.42 0.38 y 1.23 1.51 −2.92 0.75 z

0.44 0.30 0.27

0.31 0.12 0.24 24.49 72.16 −25.66 1.09 0.46 0.18 25.94 73.30 −25.98 0.07 0.06 0.08 0.08 25.09 0.06 73.68 0.12 −25.73 0.08 25.10 0.03 73.60 0.08 −25.78 24.93 72.21 −26.82 0.85 0.50 0.13 25.12 72.41 −25.67

0.4 m/s AVG. STD. 0.3 m/s AVG. STD. 0.2 m/s AVG. STD. 0.1 m/s AVG. STD. 0 m/s AVG. STD. −0.1 m/s AVG. STD. −0.2 m/s AVG. STD. −0.3 m/s AVG. STD. −0.4 m/s AVG. STD.

Table 2 Averages (avg.) and standardeviation (std.) of determined positions (cm – unit).

In this study, we address the position error due to different noises when the measured LPM signal contains additive white Gaussian noise (WGN), colored Gaussian noise (CGN), general Gaussian noise (GGN), Laplacian noise (LN), and Gaussian mixture noise (GMN). To proceed further with the measurement issues, we need to define the necessary background of each noise. The WGN has the following important points: each sample has the identity distributed PDF and the PDF is Gaussian, each noise is independent of other samples and it cannot be predicted based upon observing the others, Power Spectral Density (PSD) is flat. We can typically use WGN in communication systems to model the ambient noise for virtually any signal processing [26]. The CGN can often be seen in clutter in radar [37] and reverberator in sonar [38]. PSD of colored Gaussian noise is decidedly non-flat in the band of interest. The useful model for a non-flat PSD is the autoregressive (AR) model, and the first order AR for the CGN model is performed in this study. In the previous two noises, Gaussian models for stationary noise data allow us to define a PDS but non-stationary noise data does not. The GGN, which is a type of non-stationary Gaussian random processes, is used for an additive noise. It can occur in practice due to abrupt changes. Then, two non-Gaussian PDFs, the LN and the GMN generating large noise spikes of noise samples, can be assumed to model such physical noise processes as man-made noise [39], acoustic transients [40], and geomagnetic noise [41]. Noise processes as previously demonstrated can be added in the LPM signal to account for noise effects in measurement with a SNR of 0 dB. We assumed that the object’s distance was varied to increase from 65 to 90 cm. Absolute errors due to noisy signals of the WGN, CGN, GGN, LN, and GMN are illustrated in Fig. 6. For the comparison, it is observed that the maximum absolute errors attain 0.7 mm at the Laplacian noisy signal, as illustrated in Fig. 6(d). This noisy signal yields the most variance when relatively compared with others. The absolute errors of other noisy signals oscillate in a range approximately from 0.675 to 0.695 mm. In Fig. 7, the Cumulative density function (CDF) of position errors at a reference position (x = −80 cm, y = 65 cm, z = −45 cm) is evaluated with different noisy signals at SNRs of 10 dB, 0 dB, and -10 dB. The simulation results are

0.35 0.17 0.26

4.2. Position error due to different kinds of noise

0.20 0.35 0.44

18.42 96.60 4.67

0.14 0.23 0.18

18.49 96.60 4.61

0.10 0.06 0.16

To perform estimation, we used a Monte-Carlo simulation to achieve the means and the variances against the CRLB. A position at x = 50 cm, y = 50 cm, z = 50 cm and d = 86.6025 cm was simulated. A zero mean WGN process with the variance 2 is included to the LPM signal. The means and variance are tabulated in Table 1 with SNR of 10, 0 and −10 dB. The Cramer-Rao bounds, computed using (17a)–(17c) are obtained under the variance for comparison. The result can be observed that the means closely achieve the assumed position, and the variance can reach the CRLB for SNR as low as 0 dB, as shown in Table 1.

25.08 73.49 −24.64

86.5692 0.02557 0.55801 × 10−3

Second position (x = 25 cm, y = 70 cm, z = −25 cm 24.93 0.75 25.36 0.78 x 68.80 0.40 73.31 0.68 y −28.15 0.97 −25.62 0.54 z

49.4945 0.16008 0.07425

0.17 0.39 0.17

51.8095 0.05256 0.02919

18.512 97.138 4.70

86.6711 0.21002 × 10−3 0.17619 × 10−3

0.13 0.07 0.24

50.0445 2.05674 × 10−3 1.50044 × 10−3

18.52 96.89 4.74

50.0392 0.898697 × 10−3 0.76897 × 10−3

18.32 95.96 4.65

86.6712 0.13606 × 10−3 0.14998 × 10−3

0.14 0.05 0.13

50.0456 0.92723 ×10−3 1.20002 × 10−3

18.35 95.59 4.78

50.0311 0.73231 × 10−3 0.72180 × 10−3

0.78 0.21 0.17

Distance (cm) 86.6025 cm

18.47 95.64 4.79

z-position (cm) 50 cm

First position (x = 19 cm, y = 92 cm, z = 7 cm) 18.02 1.43 18.48 0.16 x 96.39 0.51 96.06 0.56 y 5.15 0.24 5.12 0.28 z

10-dB SNR Mean Variance CRLB 0-dB SNR Mean Variance CRLB −10-dB SNR Mean Variance CRLB

x-position (cm) 50 cm

0.21

Table 1 Comparison of the crlb and the variances for snrs 100 trials – monte carlo simulation.

−16.99

486

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

487

Fig. 7. CDFs of x–y–z position errors from various noisy LPM signals at a coordinate in 3D x = −80 cm, y = 65, and z = −45 cm

observed that, at SNRs of 10 dB and 0 dB, the results can accuracy at the minimum errors of x–y–z parameters (x = 0.1 mm, y = 0.8 mm, and z = 0.1 mm) and at the maximum errors of x–y–z parameters (x = 0.5 mm, y = 1.4 mm, and z = 0.5 mm). Next, the proposed method was tested in the worst environment at SNR of −10 dB. The results in Fig. 7(c), (f), and (i) shows that the proposed method cannot obtain accuracy from almost noisy signals because more than 100 mm of errors are observed. However, the LPM signal, which is embedded in General Gaussian Noise (GGN), can provide robustness of SNR as low as −10 dB at the maximum error about 5 mm. This observation can imply that the cross correlation, which was used for TOF computation, is strong to GGN as low as SNR of −10 dB, and

is highly sensitive to WGN, CGN, LN, and GMN in the same condition. In the large picture, we can conclude that accuracy of the proposed method can be guaranteed from simulations of different noisy signals at the maximum error 1.5 mm of SNR as low as 0 dB. 4.3. Accuracy comparison with an iterative method The proposed method is realized using geometrical constraints of the nonlinear equations to the linear equations, and sometimes the nonlinear equations can be solved using iterative methods. To compare the performance of the proposed method, we have surveyed the different iterative methods for the nonlinear equations.

488

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

Fig. 8. Comparison of accuracy between the Newton Raphson method and the proposed method when positions swept from x = −80 to −79 cm, y = 64 to 65 cm, and z = −46 to −45 cm.

The simplex search method [42] (a first order method) utilizes a function evaluation by means of iterating the unknown parameters. This method has a serious problem that is significant iteration requirement. The Lavenberg–Marquardt method [43] (a second order) employs gradient as well as the Hessian matrix approximation for the next routine. Although this method get faster than the simplex method but it is weak to local minima and requires a good initial guess to achieve the optimal solution. A well-known iterative method based on only the first two terms of Taylor-series expansion is the Newton–Raphson method for unknown-parameters estimation. The Newton–Raphson method has been developed for fast computation to evaluate high accuracy model-based estimation [44]. To assess the accuracy of the proposed method, we utilize Newton–Raphson method to observe comparison from a reference position with x = −80 to −79 cm, y = 64–65 cm, and z = −46 to

−45 cm at SNR of 0 dB. In position comparison, as shown in Fig. 8, it is clear that the Newton–Raphson method for the position measurement can acquire very high accuracy. However, the proposed method provides less accuracy than the Newton–Raphson but it can get relatively close to the reference point. Accuracy of x-position is lower than 1 mm, as shown in Fig. 8(a). Results of y-position can be obtained not over than 2 mm of accuracy, as shown in Fig. 8(b). In z-position, it seems that the proposed method is possible to attain the Newton–Raphson method with high accuracy. Although the Newton–Raphson method can achieve the optimum but it takes significant iterations and a good initial guess does not guarantee the optimal solution [45]. For real-time navigation, very high accuracy is not a big problem, but computational-time issue is more necessary. In this paper, we selected to use a method without taking iterations to reduce a computational time.

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

489

Fig. 9. A plane wave traveling across a planar surface (a) reflection of a plane pressure wave incident on a fluid–solid interface. (b) reflection of different surfaces.

4.4. Reflection at a fluid–solid interface for echolocation In general, echolocation for robot’s navigation is an issue involving with the study of sound propagation in medium and the reflection of target. The sound waves generated by the vibrating diaphragm of a loudspeaker can move through air, water and solids as longitudinal waves and also as transverse waves in solids. Amplitude of sound propagation is dependent on attenuation as a function of distance through the compressible medium, and attenuation in air is typically constant at 22.1 dB/m [34]. Reflection is a change in direction of a wavefront at an interface between two different media to return from an origin that sound is generated. Consider a sound pressure wave in a fluid (air) incident on a plane fluid–solid boundary (air to object), as shown in Fig. 9(a). To satisfy the interface boundary situations, this incident P-wave (compression wave) generates both transmitted P-wave and S-wave (shear wave), and this condition is called mode conversion. Next, we assume that the plane of incidence is considered as the x–y plane. Interface conditions of incident, reflected, and transmitted waves can then reduce to a set of linear equations for the unknown amplitudes, that is [46]

⎡

cos p1

cos p2 cp2

− sin s2 cs2

2 cp2

2 cs2

⎤

⎡ ⎤ ⎢ 1 cp1 ω2 ⎥ Pr ⎢ ⎥ ⎢ ⎥ ⎢ 1 −2 ω2 cos 2 s2 2 ω2 sin 2 s2 ⎥ ⎢ ⎥ ⎣ At ⎦ ⎢ ⎥ sin 2 p2 ⎣ ⎦ Bt sin 2 s2 0

⎡ P cos p1 i ⎢ 1 cp1 ω2 =⎢ ⎣ −Pi

⎤

⎥ ⎥, ⎦

(18)

0 Fig. 10. Comparison of reflected wavefront (a) rectangular shape (b) circular shape (c) triangular shape.

where  is a frequency, and Pi , Pr , Ai , and Bi are amplitude of incident, reflected, transmitted P-, and transmitted

490

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

Fig. 11. Accuracy of object position due to microphone positions at SNR = 0 dB (a) direction of varying microphone positions (b) x-coordinate error (c) y-coordinate error (d) z-coordinate error.

S-waves. We can define reflection and transmission amplitude ratios:



Pr = Pi

At = Pi

Bt = Pi

2 cp2 cos p1 1 cp1 2 cp2 cos p1 1 cp1

 1 cp1

ω2

2 s2 +

c 2 sin 2 s2 sin 2 p2 s2

c2

p2

 cos2 2 s2 +

c 2 sin 2 s2 sin 2 p2 s2



− cos p2



c2

(19a) + cos p2

p2

2c

2 cp2 cos p1 1 cp1

 1 cp1 cp2 ω2

cos2

 p2

cos p1 cos 2 s2

cos2 2 s2 +

c 2 sin 2 s2 sin 2 p2 s2



c2 p2

2 cos p1 sin 2 p2 −2cs2 2 cp2 cos p1 1 cp1



cos2 2 s2 +

c 2 sin 2 s2 sin 2 p2 s2

c2

p2



(19b)

+ cos p2



 + cos p2

(19c)

Now, in our study (air to stainless steel interface), we consider the behavior of reflected and transmitted waves from the definition of Snell’s law. There are two critical angles present in the solid medium because speed wave in stainless steel (5790 m/s) is much more than that in air (345.1 m/s). The transmitted P- and S-waves become an inhomogeneous wave. These waves do not affect to the reflected echo sensed by microphones. Consider in Eq. (19a), this is directly proportional to amplitude of the reflected wave. It is observed that the reflection and transmission amplitude ratio is approximately equal to 1 when the incident angle ( p1 ) is 90◦ , as shown in a rectangular surface of Fig. 9(b). In earlier works involving with one-bit signal processing [26,27], a mirror-surface object was assumed as a target and it resulted in the high SNR of reflection. In this work, we have studied and developed a system to satisfy other surfaces over a mirror-surface object. Using Eq. (19a)–(19c), ABAQUS, finite element analysis, was used to study the intensity of

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

491

Fig. 12. PDOP values in the X–Z plane (a) at y = 50 cm and (b) at y = 1 m from the original point.

reflections from different surfaces when air and stainless steel densities of interface boundary were assumed as 7890 and 0.0012 g/m3 , respectively. The simulation results shown in Fig. 10 summarize the intensity of reflections generated by rectangular, circular, and triangular objects. The sound source was located as a small color point opposite to the object position. Wavefront is observed that the reflection of the rectangular surface, as illustrated in Fig. 10(a), is the strongest intensity when relatively compared with other shapes. The reflection of the triangular surface in Fig. 10(c) has dispersion moving toward other directions compared with before incidence. It is difficult for navigation-based echolocation because the received echo at microphones is very low SNR. For the circular shape in Fig. 10(b), although some wavefront of the echo can

be reflected to the same direction that a sound source generates but it is lower intensity than the rectangular shape. Therefore, the reflection of the spherical object as a target of the system-based echolocation is possible and fascinating to be used for the study. 4.5. Effect of microphone position As the proposed system, the object position expressed in Eq. (5) is solved from TOFs and microphone positions on X–Z plane. In this subsection, accuracy of the proposed method is studied when microphone positions are varied directly whereas the object position is constant . One of the microphone arrangements is shown in Fig. 11(a). The first microphone position was assumed as (x1 , 0,

492

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

Fig. 13. Experimental setup of three-dimensional ultrasonic position and velocity.

z1 ) and had a direction of moving on Z axis. The second and third microphones were at (x2 , 0, z2 ) and (x3 , 0, z3 ), respectively and their directions were changed on diagonal lines between X and Z axes. We assumed that microphone coordinates were varied with

the same magnitude in each line from 0 to 50 cm, x1 = x2 = −x3 and z1 = −z2 = −z3 . Simulation results were observed that x–y–z coordinates of the object had more variance when microphone positions were moving from 0 cm to approximately 10 cm as shown

Fig. 14. Cumulative density function of position error at x = 19 cm, y = 92 cm, and z = 7 cm when varying the velocity: (a) x-position; (b) y-position; and (c) z-position.

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

493

Fig. 15. Cumulative density function of position error at x = 25 cm, y = 70 cm, and z = −25 cm when varying the velocity: (a) x-position; (b) y-position; and (c) z-position.

in Fig. 11(b)–(d). After 10 cm, errors of x–y–z coordinates became quite constant at 0.05 cm, 0.2 cm, and 0.04 cm, respectively. For the reason why error of y-coordinate larger than that of x and z coordinates, since the proposed method estimates the object position with the parameters d, z, and x, earlier than the parameter y as defined in Eq. (6), therefore, error of the parameter y is accumulation of errors from the parameters d, z, x. We have used this observation for the microphone-position design at 0, 0, 10 cm, 10, 0, −10 cm, and −10, 0, −10cm for the experimental evaluation, as shown in Fig. 13.

4.6. Study of position dilution of precision (PDOP) PDOP is a task of selecting the appropriate evaluation units, which results in a better geometric configuration and more accurate estimation. PDOP involve with the uncertainty in the object position (x, y, z) and the uncertainty in distance (d). This factor can provide a ratio about errors in the distance directly effecting to errors in the object position. The small the PDOP is, the less sensitive to errors the proposed method becomes [47]. The definition for PDOP can be expressed in Equation (20), where x , y , and z

are the standard deviations for the object position and ␴d is the standard deviation of the distance.



PDOP =

x2 + y2 + z2

d

(20)

Simulation results explain the PDOP in the X–Z plane in front of the speaker position 0,0,0 cm, the microphone positions 0, 0, 10 cm, 10, 0, −10 cm, and −10, 0, −10 cm for an object position at y = 50 cm in Fig. 12(a) and at y = 100 cm in Fig. 12(b). PDOP is computed and expressed on the X-Z plane. From these results, when the object positions, at the same y-coordinate, were moved far away from a central line (x = 0 cm, z = 0 cm), errors of the proposed method became more sensitive. Moreover, if the object position was increased in the direction of +Y axis, PDOP was larger in proportion to the distance. It means that the more distance the object belongs to, the more sensitive to errors the proposed method obtains. 5. Evaluation of fluctuation in air for the proposed system 5.1. Experimental setup The experimental setup for the three-dimensional positioning system is proposed in Fig. 11. The frequency of the transmitted LPM

494

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

Fig. 16. Cumulative density function of position error at x = −50 cm, y = 75 cm, and z = 0 cm when varying the velocity: (a) x-position; (b) y-position; and (c) z-position.

signal was swept from 50 kHz to 20 kHz. The length of the transmitted LPM signal was 3.274 ms, and the length of two cycles of the LPM signal was 6.548 ms. The driving voltage of the function generator was 4 Vp-p. The LPM signal was distributed from the function generator and enlarged 10 times with an amplifier. The loudspeaker transmitted the LPM signal, and then echoes were sensed by the silicon MEMS microphone based on a development of hand held telecommunication instruments using Knowles Acoustic model SPM0204UD5 [28]. A stainless ball with a diameter 10 cm was used as a test object. According to the calculations, the first microphone was located at x1 10 cm, y1 = 0 cm, z1 = -10 cm, the second microphone was located at x2 =0 cm, y2 = 0 cm, z2 = 10 cm, the third microphone was located at x3 = −10 cm, y3 = 0 cm, z3 = −10 cm. The propagation velocity of an ultrasonic wave in the air was approximately 345 m/s at a temperature of 20◦ –25 ◦ C and a humidity of 40–50 R.H. The signals measured by the acoustical receivers were converted into one-bit signal by 7th order delta–sigma modulators AD7720. The sampling frequency of the delta–sigma modulator was 12.5 MHz. The cross correlation and the smoothing operation performed by 141 taps of weighted moving average filter were designed and embedded on an FPGA board model cyclone V 5CGXFC5C6F27C7. The logic utilization for the cross correlation using one-bit signal, programmed into the FPGA board were 2602 logic elements; the number of total registers was 5777; and the total block of memory used included 175,948 bits. The moving

object was driven using SIGMA KOKI SGMA46-300 motorized stage, which can drive a moving object in only the +Y and –Y direction with maximum and minimum velocities of ±0.4 m/s; velocities could be adjusted in ±0.1 m/s steps. A reference point of the object position is set up on the object surface and synchronized with the acoustic transmitter by using an optical sensor shown in Fig. 11. When the moving object passed through the reference point, instantaneously the speaker was trigged to generate the LPM signals. 5.2. Experimental results and discussion The three-dimensional position determined by the proposed method was evaluated by the experimental results. In the experiments, the velocity of the motorized stage was adjusted from -0.4 m/s to +0.4 m/s. The velocity of the motorized stage was determined by 50 experiments per position of the moving object. The cumulative density functions (CDFs) of the positions are illustrated in Figs. 14–17. Averages and standard deviations of the determined positions are also shown in Table 2. The examples of the determined positions in the proposed system are assumed to be located at four positions. To ensure that the real moving object positions were covered in every quadrant of detection, the position was assumed to be on the +X, −X, +Z, −Z, and +Y axes. The first position was 19, 94, 7 cm, the second was 25, 70, −25 cm, the third was −50, 75, 0 cm, and the fourth was −40, 60, −20 cm. The velocities determined by the

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

495

Fig. 17. Cumulative density function of position error at x = −40 cm, y = 60 cm, and z = −20 cm when varying the velocity: (a) x-position; (b) y-position; and (c) z-position.

proposed method of Doppler velocity estimation and Doppler-shift compensation are illustrated in Fig. 20 and Table 3. The velocity and position of the moving object can be concurrently determined by the proposed ultrasonic system based on one-bit signal processing. The experimental results at a reference point x = 19 cm in Fig. 14(a) shows that the x-position error measurement at −0.4–0.4 m/s is approximately a ranging performance from 0.5 to 1.5 cm and there are the 10% x-position error at −0.4 m/s provides higher than 2 cm. At a reference point y = 92 cm in Fig. 14(b), the errors vary from 3 to 5 cm but the 40% y-position error is more than 5 cm at the velocity 4 m/s. In Fig 14(c), the CDF of the z-position errors has a narrow bandwidth interval 1.5–2.5 cm. In the second position 25, 70, −25 cm of Fig. 15, the x-position errors are less than 0.25 cm at the velocity range −0.2–0.2 m/s while the rest of the x-position errors varies from 0 to 2 cm, the y-position errors are approximately in a range measurement 0– 4.5 cm, the z-position errors swing from 0 to 2 cm but at the error range (0–5 cm) of the velocity −0.4 m/s the CDF is wider distribution than other velocities. The errors of the third position −50, 75, 0 cm have deviation from x = 0–2 cm, y = 0–2.5 cm, and z = 0–3 cm as shown in Fig. 16. In the fourth position (−40, 60, −20), the x–y–z position errors are determined as the interval 0–2 cm, 0–2.5 cm, and 0–4 cm, respectively, as shown in Fig. 17. In the experiment, deviations between the real position and assumed reference position were found; this was likely caused by

the spherical object used in this paper, which had a diameter of 12 cm, while the proposed three-dimensions-positioning measurement assumed that the object was a very small point. Therefore, the sound that was incident on the object propagated from many points, which resulted in many different TOFs at the same acoustical receivers as shown in Fig. 18. This affected the efficiency from the echoed-sound propagation of the object, producing an uncertainty due to variance of the position and velocity measurements used in the proposed system. To evaluate the experimental results, precision from 50 times of observation is mentioned further . The reference object position at 19, 92, 7 cm in Fig. 14 was set up at a location near the central line of +Y axis. In this location, it is noticed that the sound pressure was moving toward the object and reflected to the receivers at an intense and constant sound-beam level because precision of the point (x–y–z coordinates) of contact on the object was higher than other positions far from the central line of +Y axis. The precision obtained from experiment is approximately x = 18.4 ± 0.3 cm, y = 96.3 ± 0.4 cm, and z = 4.9 ± 0.28 cm. Next, the second position was shifted to the negative direction of Z-axis at 25, 70, −25 cm. The experimental results in Fig. 15 show that when the object was at a location much lower than the central line of +Y axis, the precision had many points of contact on the object and the dispersion of standard deviation at the observed position. This was a position comparatively far away from the main sound beam expressed in Fig. 2(b). Therefore, the

496

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

precision of contact at the same location was less likely. In the third reference position at −50, 75, 0 cm, the object position was assumed in the direction of –X axis on the right hand side of the loudspeaker. Although the object was assumed on X–Y plane having the strong intensity of sound beam Fig 2(a), it was placed far away from the central line of +Y axis. The results in Fig. 16 from the third position were similar to the second position because the precision had several points of contact on the object and the dispersion of standard deviations. Finally, with the object positioned at a reference −40, 60, −20 cm, this position was on a quadrant of –X and –Z axes. The results in Fig. 17 still had many points of contact on the object when relatively compared with the first reference position at the approximate central point of the main sound beam. The experimental results can be concluded that the object at a location near the main sound beam had higher precision than that far away from the central line of the main sound beam due to more oblique contact with the main beam. The sound pressure moving toward the object was not constant. This agrees with our previous paper [48], which studied the position estimation of the non-moving rigid object in three dimensions. The precision is also consistent with the PDOP expressed in subsection 4.6 because the PDOP provides the greater value at a position far away from the central line of +Y axis. If carefully scrutinizing in the repeatability of x, y, and z coordinates, we observed that z-coordinate has the highest precision when relatively compared with others. This can be explained from Eq. (6). x parameter was computed using the values of d and z, and then y parameter was calculated from those of d, z, and x. It means that the uncertainty of x and y was accumulated from the earlier parameters (d, z). In addition, the experimental results at each velocity indicate that at high velocity (−0.3, 0.3, −0.4, and 0.4 m/s) the object position had higher variance than at low velocity (0, −0.1, 0.1, −0.2, and 0.2 m/s). The motorized stage for driving the object was controlled by an automatic system to evaluate repeatability. Then, a rod was connected between the object and the motorized stage. When the moving object was accelerated to a high velocity in each turn, the object did not completely stop because the rod swung with higher vibration. The object assumed as a rigid body thus was not at a static condition during the next measurements. We can explain this condition by experimentally measuring vibration velocity versus various velocities. This effect was evaluated by a Polytech OFV-302 vibrometer during vibration

Fig. 18. Possibility of reflection area on the spherical object.

testing. Relatability to high velocity measurements of a moving part in the Surface Acoustic Wave (SAW) motor [49] was found along with vibration velocity of the titanium-based hydrothermal lead zirconate titanate thick film transducer. The vibrometer was placed 50 cm away from a measurement point. At that point, the object was accelerated to −0.4, −0.3, −0.2, −0.1, 0, 0.1, 0.2, 0.3 and 0.4 m/s, and the rod, that was connected to the moving object and the motorized stage, had a length 40 cm. Experiment results regarding the vibration effect of the moving object is shown in Fig. 19; higher velocities tend to show higher vibration velocities. Regarding object position at 0, 0.1, 0.2, −0.1, and −0.2 m/s, they are shown to have satisfactory repeatability. A moving object had a lower vibration than the objects moving at 0.3, 0.4, −0.3, and −0.4 m/s. However, vibration velocity did not proportionally affect to the contact point on the object, for instance, at the reference position (−40, 60, −20) cm the position of −0.3 m/s provides more dispersion than that of 0.3, 0.4, and −0.4 m/s while at the reference point (−50, 75, 0) cm the position of −0.4 m/s is the worst case of precision. In an overall picture, the higher vibration velocity on surface due to the moving object directly affects to the repeatability. We will study the relationship between vibration velocity and reflection in the future work. In Fig. 20, the determined velocity measurement in the ±Y direction was presented alone. The component uy was analyzed because it is in the same direction as the motorized stage movement. The component vectors of ux and uz were estimated to be nearly zero.

Fig. 19. Vibration velocity of the moving object at given velocities.

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

497

Fig. 20. Cumulative density function of velocity error at various positions: (a) first position; (b) second position; (c) third position; and (d) fourth position.

Table 3 Averages (avg.) and standardeviation (std.) of determined velocities (m/s – unit).

0.4 m/s 0.3 m/s 0.2 m/s 0.1 m/s 0 m/s −0.1 m/s −0.2 m/s −0.3 m/s −0.4 m/s

First position (19,92,7) AVG. STD.

Second position (25,70,−25) AVG. STD.

Third position (−50,75,0) AVG. STD.

Fourth position (−40,60,−20) AVG. STD.

0.4017 0.0253 0.3075 0.0156 0.1994 0.0133 0.0925 0.0173 2.222 × 10−4 0.0103 −0.1008 0.0146 −0.2034 0.0170 −0.3047 0.0204 −0.4022 0.0157

0.4324 0.0101 0.3239 0.0125 0.2132 0.0087 0.1100 0.0076 −5.637 × 10−3 0.011 −0.1079 0.0116 −0.2132 0.0103 −0.3339 0.0184 −0.4481 0.0132

0.4290 0.0247 0.3311 0.0087 0.2201 0.0084 0.1027 0.0074 −1.65 x −3 0.0044 −0.1144 0.0102 −0.2215 0.0121 −0.3384 0.0148 −0.4254 0.0082

0.4212 0.0065 0.2932 0.0096 0.2118 0.0076 0.0864 0.0095 −1.40 × 10−3 0.0052 −0.0921 0.0091 −0.2080 0.0107 −0.3341 0.0101 −0.4332 0.0112

At the first positon in Fig. 20(a), the velocity errors have more satisfaction than other positions and 10% of the velocity errors over 0.05 m/s. The velocity-measurement range, from −0.2 to 0.2 m/s, of the position (x = 25, y = 70, z = −25) cm provides the velocity errors smaller than 0.04 m/s but for other velocities 10 % of their results have error at nearly 0.07 m/s, as shown in Fig. 20(b). In Fig. 20(c), the velocity errors are less than 0.03 m/s at the velocity measurement from −0.1 to 0.1 m/s, and the velocity errors of −0.4, −0.3, 0.3, and 0.4 m/s have the interval from 0 to 0.08 m/s. In Fig. 17 (d), the most

velocity error varies in the range of 0 to 0.045 m/s, except −0.4, −0.3, and 0.4 m/s the results are inside the error range of 0.02–0.065 m/s. Measurement was shown to have the highest error and dispersion at velocities of 0.3, 0.4, −0.3, and −0.4 m/s, which was likely because three-dimensional velocity vector measurement relied on the position computation and thus directly varied according to the accuracy of the determined position. The performance of the proposed system form the observations has the precision for the maximum ranging: 49.85 ± 0.23 cm

498

(x = 50 cm), (z = 25 cm).

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499

95.96 ± 0.17 cm

(y = 92 cm),

and

26.82 ± 0.12 cm

6. Conclusion An over-sampling signal processing method using ultrasonic position and velocity measurements for a three-dimensional system using a pair of LPM signals has been proposed. The proposed system by redesigning the receiver arrangement uses cross-correlation one-bit signal processing technology, low-computation-cost Doppler-shift compensation, and low-cost devices. This is the development from the previous space measurement of one-bit signal processing to simultaneously measure position and velocity in three-dimensional space. Through MonteCarlo computer simulation of 100 trials, the variances of estimated parameters can achieve the analytical CRLB for SNR as low as 0 dB. Accuracy of the proposed method was tested by the five noisy signals. The results can provide accuracy less than 0.1 mm at a SNR of 0 dB. Since the position-measurement model is a form of nonlinear systems, it is then improved into a linear model. Accuracy of the proposed method was compared with an iterative method, which is well known for dealing with the nonlinear problems. The simulation results show that the performance of the proposed method can attain accuracy of the Newton–Raphson method. In addition, the simulation results show the wavefront reflection from rectangular, circular, and triangular surfaces. The ultrasound receiver position also affects to accuracy of the object position. PDOP is noticed that the object position far away from the central line of +Y axis provide the high value. A delta–sigma modulation board was used to convert analog signals into one-bit digital signals. Then, cross-correlation and Doppler-shift were computed with a FPGA. The object used to test the system was a spherical stainless object that was driven by a motorized stage during velocity-measurement evaluation. In this paper, three-dimensional—positioning measurement was determined by three acoustical receivers and one sound source. We proposed a method to design the optimal receivers-positions. The velocity of the moving object was calculated using the concept of vector projection. In the experiments, four-example positions in each area were evaluated with velocities between −0.4 m/s and +0.4 m/s by in steps of 0.1 m/s. The repeatability of the proposed system was determined via 50 measurements at each position and velocity. From the experimental results, it is observed that accuracy of the position measurement compared with a reference point can be achieved between 0.2 and 5 cm, and accuracy of the velocity measurement is between 0.02 and 0.07 m/s. Moreover, we notice that the moving-object’s position showed the most variance at high velocities due to the vibration velocity on the object surface. For future work, multiple moving objects in the three dimensional area and the vibration velocity of the moving object will be studied. References [1] B.R. Mahafza, Radar Systems Analysis and Design using MATLAB, CRC Press, Florida, 2000, pp. 5–8. [2] D.C. Finfer, P.R. White, G.H. Chua, T.G. Leighton, Review of the occurrence of multiple pulse echolocation clicks in recordings from small odontocetes, IET Radar Sonar Nav. 6 (6) (2012) 545–555. [3] B.G. Ferguson, B.G. Quinn, Application of the short-time fourier transform and the Wigner–Ville distribution to the acoustic localization of aircraft, J. Acoust. Soc. Am. 96 (August (2)) (1994) 821–827. [4] J.M. Martin, A.R. Jimenez, F. Seco, L. Calderon, J.L. Pons, R. Ceres, Estimating the 3D-position from time delay data of US-waves: experimental analysis and a new processing algorithm, Sens. Actuat. A-Phys. 101 (2002) 311–321. [5] J. Jiminez, M. Mazo, J. Urena, A. Hernandez, F. Alvarez, J. Garcia, E. Santiso, Using PCA in time-of-flight vectors for reflector recognition and 3-D localization, IEEE Trans. Robot. 21 (October (5)) (2005) 909–924.

[6] B. Kreczmer, (2010), Objects Localization and Differentiation Using Ultrasonic Sensor, Robot Localization and Map Building, Hanafiah Yussof (Ed), InTech. [online] http:// www. Intechopen.com/books/robot-localization -and -and -map -buiding/object -localization -and differentiation -using -ultrasonic -sensor. [7] D. Klimentjew, N. Hendrich, J. Zang, Multi sensor fusion of camera and 3D laser range finder for object recognition, In Proc. IEEE Int. Conf. MF-IIS (2010) 236–241. [8] S. Hussmann, T. Edeler, Pseudo-four-phase-shift algorithm for performance enhancement of 3D-TOF vision systems, IEEE Trans. Instrum. Meas. 59 (May 5) (2010) 1175–1181. [9] J. Smit, R., Kleihorst, A., Abbo, J. Meuleman, and G. van Willigenburg, Real time depth mapping performed on an autonomous stereo vision module pp. 306–310. [online] http://www.researchgate.net/publication/40122241 Realtime depth mapping performed on an autonomous stereo vision module [10] M. Yano, I. Matsuo, J. Tani, Echolocation of multiple targets in 3-D space from a single emission, J. Biol. Phys. 28 (2002) 509–525. [11] I. Matsuo, J. Tani, M. Yano, A model of echolocation of multiple target in 3D space from a single emission, J. Acoust. Soc. Am. 110 (July (1)) (2001) 607–624. [12] L. Christopher, Yatrakis Computing the Cross Ambiguity Function—A Review Master Thesis, Binghamton University, State University of New York, 2001. [13] L. Bastista, E. Schena, G. Schiavone, S.A. Sciuto, S. Silvestri, Calibration and uncertainly evaluation using monte carlo method of a simple 2D sound localization system, IEEE Sens. J. vol. 13 (September (9)) (2013) 3312–3318. [14] J.C. Murray, H. Erwin, S. Wermter, Robotic sound-source localization and tracking using interaural time difference and cross-correlation, In Proc. AI workshop NeuroBot. (2004). [15] A. Pourmohammad, S.M. Ahadi, Real time high accuracy 3-D PHAT-based sound source localization using a simple 4-microphone arrangement, IEEE Syst. J. 6 (September (3)) (2012) 455–468. [16] B. Kwon, G. Kim, Y. Park, Sound source localization methods with considering of microphone placement in robot platform, In Proc. IEEE RO-MAN 26–29 (August) (2007) 127–130. [17] S. Nakamura, T. Sato, M. Sugimoto, H. Hashizume, An accurate technique for simultaneous measurement of 3D position and velocity of a moving object using a single ultrasonic receiver unit, In Proc. IPIN (September) (2010) 15–17. [18] S. Ohyama, A.H. Alasiry, J. Takayama, A. Kobayashi, Determining 2D position of sensor network nodes for temperature distribution measurement, Sens. Actuators A-Phys. 135 (2007) 203–208. [19] D. Pompili, T. Melodia, I.F. Akyildiz, Three-dimensional and two-dimensional deployment analysis for underwater acoustic sensor networks, Ad Hoc Networks 7 (2009) 778–790. [20] A. Griffin, A. Alexandridis, D. Pavlidi, Y. Mastorakis, A. Mouchtaris, Localizing multiple audio sources in a wireless acoustic sensor network, Signal Process. 107 (2015) 54–67. [21] S. Shatara, X. Tan, An efficient, time-of-flight-based underwater acoustic ranging system for small robotic fish, IEEE J. Oceanic Eng. 35 (October (4)) (2010). [22] Super Audio CD. Production Using Direct Stream Digital Technology, Electronics and Sony Corporation copyright, 2000 and 2001. [online] http:// www.canadapromedia.com [23] S. Hirata, M.K. Kurosawa, and T. Katagiri Accuracy and resolution of ultrasonic distance measurement with high-time-resolution cross-correlation function obtained by one-bit signal processing, Acoust. Sci. Tech. 30 (6) (2009) 429–438. [24] R.A. Altes, Sonar velocity resolution with a linear-period-modulated pulse, J. Acoust. Soc. Am. 61 (S1) (1977) 1019–1030. [25] J.J. Kroszczynski, Pulse compression by means of linear period modulation, In Proc. IEEE 5 (1969) 1260–1266. [26] S. Hirata, M.K. Kurosawa, Ultrasonic distance and velocity measurement using a pair of LPM signals for cross-correlation methodimprovement of Doppler-shift compensation of Doppler velocity estimation, Ultrasonics 52 (February) (2012) 873–879. [27] S. Saito, M.K. Kurosawa, Y. Orino, S. Hirata, Airborne ultrasonic position and velocity measurement using two cycles of linear-period-modulated signal, in: In Proc. TAROS 2011, Sheffield, UK, August 31–September 2, 2011, pp. 46–53. [28] S. Thong-un, M.K. Kurosawa, Y. Orino, Three-dimensional-positioning measurements based on echolocation using linear-period-modulated ultrasonic signal, in: In Proc. IEEE IUS, Chicago, September 2014, 2014, pp. 2478–2481. [29] K.B. Ocheltree, L.A. Frizzell, Sound field calculation for rectangular sources, IEEE Trans. Ultras0n. Ferroelectr. Freq. Control 36 (March (2)) (1989) 242–248. [30] D. Royer, E. Dieulesant, Elastic Waves in Solids I Free and Guided Propagation, Springer-Verlage, Berlin Heidelberg, 2000, pp. 46–49. [31] L.L. Beranek, Acoustics Newyork, The Acoustical Society of America, 1996, pp. 109–115. [32] G. Strang, Linear Algebra and Its applications, Thomson, Belmont, 2006, pp. 71–73. [33] M.D. Fox, W.M. Gardiner, Three-dimensional Doppler velocimetry of flow jets, IEEE Trans. Biomed. Eng. 35 (October (10)) (1988) 834–840. [34] E.J. Evans, E.N. Bazley, The absorption of sound in air at audio frequencies, Acoustica 6 (1956) 238–245. [35] S.M. Kay, Fundamentals of Statistical Signal Processing, Prentice Hall, 1993.

N. Thong-un et al. / Sensors and Actuators A 233 (2015) 480–499 [36] K.C. Ho, W. Xu, An accurate algebraic solution for moving source location using TDOA and FDOA measurements, IEEE Trans. Signal Process. 52 (September (9)) (2004) 2453–2463. [37] F.E. Nathanson, Radar Design Principles, SciTech Pubs., NJ, 1999. [38] W.S. Burdic, Underwater Acoustic Systems Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1984. [39] D. Middleton, A.D. Spaulding, Element of Weak Signal Detection in Non Gaussian Environment, Advance in Statistical Signal processing, 2, JAI Press, Greenwich, CT, 1993, pp. 142–215. [40] R. Dwyer, A technique for improving detection and estimation of signals contaminated by under ice noise, J. Acoust. Soc. Am. 74 (July (1)) (1983) 124–130. [41] J. Schweiger, EValuation of geomagnetic activity in the MAD frequency band (0.04 to 0.6 Hz), in: Master Thesis, October, Naval Postgraduate School, 1982. [42] J.A. Nelder, R. Mead, A Simplex method for function minimization, Comput. J. 7 (1967) 308–313. [43] K.P. Chong, S.H. Zak, Introduction to optimization, John Wiley & Sons, Inc., New York, 2001. [44] R. Demirli, J. Saniie, Model-based estimation of ultrasonic echoes part I: analysis and algorithms, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 48 (May (3)) (2001). [45] W.H. Foy, Position-location solutions by Taylor-series estimation, IEEE Trans. Aerosp. Electron. Syst. AES-12 (1976) 187–193. [46] L.W. Schmerr Jr, Fundamentals of Ultrasonic Nondestructive Evaluation, Plenum Press, New York, 1998. [47] A.J. Martin, A.H. Alonso, D. Ruiz, I. Gude, C.D. Marziani, M.C. Perez, F.J. Alvarez, C. Gutierrez, J. Urena, EMFi-based ultrasonic sensory array for 3D localization of reflectors using positioning algorithms, IEEE Sens. J. 15 (May (5)) (2015) 2951–2962. [48] N. Thong-un, S. Hirata, M.K. Kurosawa, Improvement in airborne position measurements based on an ultrasonic linear-period-modulated wave by 1-bit signal processing, Jpn. J. Appl. Phys. 54 (7S1) (2015) 07HC06, 6 pages. [49] T. Shigematsu, M.K. Kurosawa, Friction drive of an SAW motor. Part I: measurements, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 55 (September (9)) (2008) 2005–2015.

Biographies Natee Thong-un was born in Samutprakarn, Thailand. He received the B.Eng. degree in Instrumentation system engineering from King Mongkut University of Technology North Bangkok, Thailand, in 2004. He received the M.Eng. degree in

499

Instrumentation Engineering from King Mongkut Institute of Technology Ladkrabang, Thailand, in 2006. He is currently working toward the Ph.D. degree in the Department of Information Processing, Interdisciplinary Graduate School of Science and Engineering, Tokyo Institute of Technology, Japan. He is a student member of the Institute of Electronics Information and Communication Engineers, the Acoustical Society of Japan, and IEEE. Shinnosuke Hirata received his Dr. Eng. degree in information processing from Tokyo Institute of Technology in 2009. Since 2009, he has been an assistant professor in the Department of Mechanical Engineering and Intelligent Systems, The University of Electro-Communications. Since 2012, he has been an assistant professor in the Department of Mechanical and Control Engineering, Tokyo Institute of Technology. His current research interests include ultrasound in medicine, measurement using airborne ultrasound, signal processing for ultrasound, and numerical analysis of ultrasonic propagation. He is a member of the Institute of Electronics Information and Communication Engineers, the Acoustical Society of Japan, IEEE, and the Japan Society of Ultrasonics in Medicine. Yuichiro Orino received the B.Eng. degree in physical electronics, and the M.Eng. degree in advanced applied electronics from Tokyo Institute of Technology, Tokyo, in 2002 and 2004, respectively. From 2004 to 2007 he was a Ph.D. student in the Interdisciplinary Graduate School of Engineering, Tokyo Institute of Technology, and from 2007 to 2012 he worked as a part-time researcher at Tokyo Institute of Technology. He is currently working as a Junior Researcher at the Research Center for Socioecological Systems, University of Shiga Prefecture. His research interests include digital signal processing and acoustics. Minoru Kuribayashi Kurosawa (formerly Kuribayashi) (M’95) was born in Nagano, Japan, on April 24, 1959. He received the B.Eng. Degree in electrical and electronic engineering, and the M.Eng. and D.Eng. degrees from Tokyo Institute of Technology, Tokyo, in 1982, 1984, and 1990, respectively. Beginning in 1984, he was a research associate at the Precision and Intelligence Laboratory, Tokyo Institute of Technology, Yokohama, Japan. From 1992, he was an associate professor at the Department of Precision Machinery Engineering, Graduate School of Engineering, The University of Tokyo, Tokyo, Japan. Since 1999, he has been an associate professor in the Department of Advanced Applied Electronics, then Department of Information Processing, Interdisciplinary Graduate School of Engineering, Tokyo Institute of Technology, Yokohama, Japan. His current research interests include ultrasonic motors, micro actuators, PZT film and its application to ultrasonic transducers, SAW actuators, and single bit digital signal processing and its application. Dr. Kurosawa is a member of the Institute of Electronics Information and Communication Engineers, the Acoustical Society of Japan, IEEE, the Institute of Electrical Engineers of Japan, and the Japan Society for Precision Engineering.