WO2013009264A1

WO2013009264A1 - A method of estimating time difference of arrival

Info

Publication number: WO2013009264A1
Application number: PCT/SG2012/000250
Authority: WO
Inventors: Wai Hoong Andy KHONG; Vaninirappuputhenpurayil Gopalan REJU; Amir Bin SULAIMAN; Kattukandy Rajan ARUN; Woon Seng Gan
Original assignee: Nanyang Technological University
Priority date: 2011-07-14
Filing date: 2012-07-16
Publication date: 2013-01-17

Abstract

A method of determining the origin of non-standardised surface acoustic waves in a substantially homogenous and/or quasi-heterogeneous surface comprising passively detecting the acoustic waves in the surface, Short Time Fourier transforming the acoustic waves into time frequency domain data, and determining the time difference of arrival (TDOA) of the acoustic waves based on the data.

Description

A Method of Estimating Time Difference of Arrival FIELD OF THE INVENTION

The present invention relates to a method of estimating time difference of arrival.

BACKGROUND

Touch screens are a popular human-computer interface (HCI) because they are easy to use. Various technologies have been proposed, but each has disadvantages for certain applications.

Capacitive or resistive surfaces [2] are common in mobile phones, but may become uneconomical as the size of the surface increases. They may also not be suitable for high impacts e.g. foot-steps or bullets will permanently damage the surface.

Optical-based or camera-based approaches [1] are commonly used for larger surfaces, but require extensive calibration, and may be sensitive to different lighting conditions and/or occlusions. Microphones are also used to perform source localization through air, but this may be inaccurate in most outdoor environments due to high ambient noise.

An improved technique for localization in solid surfaces is to use to measure the waves travelling through the solid surface directly. However this is difficult due to velocity dispersion [3][4] and frequency variation [5]. An approach to overcome this is using a location template matching (LTM) [6] [7]. LTM requires extensive training taps at predetermined locations. As a result LTM is vulnerable to "intra-location" dissimilarity (e.g. tapping inconsistencies due to a different style of tapping, or due to a deviation from the test conditions used during calibration). Additionally, the number of training tap locations may increase with surface area.

A more promising localization technique is estimation of the time difference of arrival (TDOA). However with existing TDOA approaches, the propagation speed is assumed before-hand. In solid materials, the propagation speeds of different frequency components varies, and any assumption of the propagation speed made prior may no longer hold true. TDOA estimation by cross-correlation [8] may also be inaccurate because of velocity dispersion. TDOA estimation using time-frequency approaches e.g. wavelet transforms [9] [10] typically requires a source signal of known characteristics and may not work for non-standardised source signals (i.e. where the tap on the surface is variable). Also, when performing source localization, a set of non-linear equations may have to be solved, and significant effort may be required for calibration. Another approach is using Surface acoustic wave (SAW) based transducers, with active emitters. However, liquid or solid contaminants moving on the surface may cause false touches or result in non-touch sensitive areas; dirt and water on the touch surface may reduce the sensitivity of the surface. Also this approach may have higher power dissipation. SUMMARY OF THE INVENTION

In general terms the present invention proposes two main embodiments for estimating the location of an impact made on a solid surface. The TDOA of the surface acoustic waves or vibrations in the surface originating from an unspecified source location may be estimated without extensive training, extensive calibration and/or high processing cost. The received signal may be analysed in the time-frequency (TiF) domain after which a difference-averaging (DA) process may be performed. The TiF-DA method may be suitable for estimating the TDOA when the rate of change of amplitude of the received signal at the onset is high, such as for a signal generated using a metal stylus. For the case when such amplitude change is not significant, such as a finger tap, the Hermitian angle (HA) is used in the TiF-HA method for estimating the onset of the signal. TiF-HA may also be usable where the signal onset is high. Other variations and/or combinations of these main embodiments are also contemplated.

Unlike optical and camera-based approaches, touch surfaces based on vibration sensing may be used under different lighting conditions or when occlusion occurs. Other existing technologies such as resistive and capacitive touch surfaces, which observe the changes in electrical current and capacitance respectively, require specially developed and fixed touch-sensitive hardware which limits the size of the interface after fabrication. These technologies are also more expensive, especially when larger surfaces are required. In some applications the impact due to "touch" can be very hard, for example in foot-step detection. Embodiments may be used for applications where the size of the surface needs to be varied or where the source of impact damages the surface but not the surface-mounted sensors such as during firearm training. Embodiments may be scalable in size during on-site applications which may be attractive for large-scale commercial or large-scale personal applications. Active emitters are not required, and embodiments may be portable, scalable and lower cost. Sensors used in the embodiments may be clipped on to the surface with sufficient bonding between the surface and the sensor. There may be no need for the sensors to be located at specific positions as long as the sensor positions are known. Also, the use of sensors may not restrict the size of the surface as long as the sensors are able to detect vibration. Embodiments may use a substantially homogenous surface and/or a quasi-heterogeneous surface. Embodiments may also be robust towards variations in the signal shape, amplitude and/or noise of a tap or impact, and may not need the selection of a threshold value for determining time difference of arrival (TDOA).

In a first specific expression of the invention there is provided a method of estimating TDOA of a stylus impact made on a solid surface as claimed in claim 1. In a second specific expression of the invention there is provided a method of estimating TDOA of a finger impact made on a solid surface as claimed in claim 2.

In a third specific expression of the invention there is provided a method of determining the origin of non-standardised surface acoustic waves in a substantially homogenous and/or quasi-heterogeneous surface as claimed in claims 3 to 19.

In a fourth specific expression of the invention there is provided a system as claimed in claim 20.

In a fifth specific expression of the invention there is provided a touch detection apparatus as claimed in claim 21 . BRIEF DESCRIPTION OF THE DRAWINGS

One or more example embodiments of the invention will now be described, with reference to the following figures, in which:

Figure 1 is a graph of signals received by a sensors placed at x, y positions (0.1, 0.1) m and (1.1, 0.1) m due to an impact at (0.6, 0.7) m on a glass surface of dimension 1 .2 m x 1 .0 m x 5 mm. (a) and (b) using stylus (c) and (d) using finger;

Figure 2 is a flow chart of the TiF-DA method;

Figure 3 is a flow chart of the TiF-HA method;

Figure 4 is a graph of absolute value of the DFT coefficients of the sensor output, j , (n) shown in Fig. 1 (a), in the time-frequency plane, \ x_l (k, m) \ . Only the first

1 024 samples of x_x (n) are used;

Figure 5 is a graph of absolute value of the DFT coefficients of the sensor output, χ_λ { ) shown in Fig. 1 (a), in the time-frequency plane represented in logarithmic scale, log₁₀ \ x_l (k, m) \ . Only the first 1 024 samples of x, (n) are used; Figure 6 is a plot of (a) χ.(42) , (b) d,.(42) , (c) ¾(42) and (d) r. . ;

Figure 7 is a graph of Hermitian angles e_H (k, m) calculated for two different signals (a) for signal generated due to stylus impact, shown in Fig. 1.i (a) and (b) for signal generated due to finger tap, shown in Fig. 1 .ii (a).

Figure 8 is a graph of the first stage of TiF-HA (a) sensor output, ,.(«) (b)

Histogram of <r.(m) (c) standard deviation of the Hermitian angles, cr.(m) (d)

Cumulative histogram (e) detected transition points and (f) histogram of the cumulative histogram;

Figure 9 is a graph of the second stage of TiF-HA (a) sensor output, x, (n) (b) Histogram of a, (m ) (c) standard deviation of the Hermitian angles, σ_: (ηι) (d)

Cumulative histogram (e) detected transition points and (f) histogram of the cumulative histogram; and

Figure 10 is a graph of the signal picked up by a sensor due to a metal tap and the standard deviation of the Hermitian angles, (a) Sensor output and (b) standard deviation of the Hermitian angles.

DETAILED DESCRIPTION

A set of low-cost surface-mounted sensors is used in example embodiments for localizing an impact made on a solid surface. In practical situations the impact can be generated via various circumstances, commonly using a finger or a stylus. However the signals generated by a finger may be different from that of a stylus. Consider an illustrative example where two Murata PKS1 -4A10 shock sensors are placed at the corners of a glass surface. Figures 1 (a)-(b) and 1 (c)-(d) show the signals received by the sensors for an impact made by a metal stylus and a finger tap, respectively. We note that the change in amplitude at the onset of the wave arrival is less pronounced for signals generated by the finger compared to that of the stylus. In addition, unlike that of signals generated by the stylus, no single dominant peak exists for signals generated by the finger tap.

.

In view of the above, according to a first embodiment a TiF-DA method is proposed for the metal stylus, which estimates the onset of wave arrival by extracting the time instant corresponding to the change in power level of the signal in the time-frequency domain via a difference-averaging process. According to a second embodiment a TiF-HA method is proposed for a signal which increases at a modest rate during the onset, such as that generated by a finger. The TiF-HA method can also be applied for an impulsive source. The time difference of arrival (TDOA) between the signals picked up by the sensors is estimated via time-frequency analysis. Due to velocity dispersion, only a few frequency bins in the time-frequency domain are used for the analysis of the signals which leads to the utilization of sliding DFT for efficient computation. This enables a surface such as glass walls or table tops to convert to a low cost touch surface for the next generation human-computer interface applications. The validity in impact localization is tested on a prototype using the signals recorded on plate surfaces of two different materials namely aluminium and glass.

The TiF-DA method 200 is shown in Figure 2. The signal from each sensor is transformed into the time-frequency domain by a sliding window Discrete Fourier Transform (DFT) 201. Selected frequency bins are differentiated with respect to time 202. The differentiated value for the selected bins are difference-averaged to obtain a corresponding average value of each time sample 203. The TDOA is then estimated by performing a correlation between the difference-averaged values of two or more sensors 204. The source location is then estimated using the TDOA 205.

The TiF-HA method 300 is shown in Figure 3. For each sensor, the DC offset is removed 301 and the signal is converted into the time-frequency domain using a sliding window DFT 302. A reference noise matrix is constructed to approximate the noise component 303. The Hermitian angle is found using the sensor data and the reference matrix 304; the angle is found between a reference vector and another vector computed using elements of the reference matrix and the sensor data at the same time frequency point. The standard deviation of the Hermitian angle across the frequency bins is then histogramed 305. The histogram is then K-means clustered to obtain two cluster centroids 306. Optionally, each centroid is smoothened using a moving average filter 307. The steps are iterated 308 until there is a single transition. The TDOA is then estimated for the Hermitian Angle transition between two or more sensors 309. Then the TDOA is used to estimate the source location. The TiF-DA method

The proposed TiF-DA method 200 is suitable for the case where the rate of change in amplitude at the onset of each sensor signal is high such as when an impact is excited by a metal stylus. Let sin) be the source signal generated on a plate surface by an impact where j ,.(n) and are the received signals from two sensors placed at different locations on the surface. Defining /¾.(«) as the impulse response between the point of impact and the i th sensor, the corresponding sensor output can be expressed as in Equation 1 :

x_j (n) = h_i (n) * s(n), (1 ) where * denotes linear convolution. Depending on the distance between the impact and sensor locations, the length of the "leading zeros" corresponding to the delay of each impulse response will be different. As an illustrative example, Fig. 1 show the received signals where an impact is made on a piece of glass plate of dimension 1.2m xl .0mx5 mm at x, y position (0.6, 0.7) m while the sensors are located at (0.1, 0.1) m and (1.1, 0.1) m. In this work we employ vibration sensors that are manufactured by Murata (Model PKS1 -4A10). We note that estimation of the TDOA via onset sample index of the time-domain waveform poses a challenge since, as can be seen from Fig. 1 , the amplitude at the onset is not significantly high. In addition, due to dispersion, the waveforms are vastly different and as such, the delay between the two signals cannot be accurately estimated using the cross-correlation method such as that described in [8]. In such a scenario, since the wave propagates with a constant speed for each frequency, the TDOA can be estimated by cross-correlating the magnitude of the discrete Fourier transform (DFT) coefficients between two received signals at the same frequency [14]. Therefore by enhancing the DFT coefficients at the time frame corresponding to the onset due to the direct path in the impulse response compared to that with other time frames, the influence of multipath can be reduced and hence the TDOA estimation accuracy can be improved. We first evaluate the short-time Fourier transform (STFT) of the received signal (n) by considering a frame of length N in Equation 2:

x,( ) = [x_;(m),...,x,.^ + N-l)]^r (2) constitutes a data frame where m = Q,\,...,M is the frame index and ^T denotes transposition operator. Here, we note that the formulation of Equation (2) implies that the STFT is performed using a hop size of one. Defining as the KxK discrete Fourier matrix, K≥N , the STFT of x,.(m) is then given by the Kxl vector in Equations 3 and 4:

x,(m) = F_Kx,(m) (3) = [x_i(0,m),---,x_i(k,m),--,x_i(K-l,m)]^T , (4) where k = 0,l,---,K-l is the frequency bin index. Signals in the frequency domain are distinguished from its time domain counterpart using an underline. ^•Defining 0_lx(ir__W) as the Ix(K-N) null vector, the variable x/(m) = [x,^r(m) 0_lx(K__N)f is a sequence obtained from x_;.(m) after padding it with K-N zeros. Assuming that s(n) is an impulse and that a tap is made closer to Sensor i compared to Sensor j , we can therefore express the STFT of *,.(«) , for the kt frequency bin, as an Mxl vector in Equation 5:

, (5)

where /i,(^,m) is the fcth frequency bin of the impulse response STFT from source to the i th sensor while L_i is the length of the "leading zeros" corresponding to the propagation delay from the source to this ί th sensor at this frequency. We therefore note that the subvector z,(£) corresponds to the STFT coefficients of the "leading zeros" while subvector h_',-(fc) corresponds to the STFT of the non-zero coefficients of the impulse response. In a similar manner, the STFT of the received signal χ.(η) can now be expressed in Equation 6:

= [h_j (k,0), ...,h_j (k,L_i + T_:j - 1) , h . (k, L_f + T_ij),..., h_j (k,M- l)]^T , (6) where τ,._;. is defined as the TDOA between the received signals. The aim therefore is to estimate r_;/ using STFTs χ_.(£) and %_.(k).

To illustrate the above, Fig.4 shows the plot of I χ (fe) I for the data shown in Fig. 1 (a) according to Equation 7:

Ιχ.(*)Ι =[\h_i(k,0)\,...^h_i(k,L_i-l)\^h_ί(k,L_i)\,...^h_i(k,M-l)\]^T , k = 0,...,K-l.

(7) For clarity of presentation, only the first 1024 samples of x_{(n) are used for the computation of Ιχ (&)l. To extract the onset of the signal, we next proceed by taking the logarithm of Equation (7) for each frequency bin k -Q,...,K-\ in Equation 8: log I χ . (k) I = [log \h_i(k,0)\,..., log I A, (k, L, - 1) I , log I Λ,. (k, L, ) I, ... , log \h_i(k,M - 1) if (8) which is then plotted in Fig.5. As can be seen, due to the non-linear mapping of the logarithmic function, the disparity between the "leading zeros" and that of the non-zero samples of the impulse response is enhanced significantly. While this disparity provides an indication of the onset between the STFTs, the difference between the elemental positions of such onsets for two received signals can therefore be used to estimate τ,, .

To estimate _i consider an illustrative example where we investigate the STFT of signals received from two sensor index 1 and 5, each at frequency-bin index £ = 42, i.e., χ,(42) and χ₅(42) plotted as shown in Fig. 6 (a). Note that the frequency-bin index k may be selected to be in the range of resonance frequency of the transducer used e.g. using the Murata PKS1-4A10 sensor, the frequencies in the range 3kHz to 7kHz may be reliable. In this illustrative example, we have used = 1024 and N = 256. To enhance clarity of presentation, we have added a value of five to every element in ¾(42). The inter-element difference within the vector log I χ (Λ:) I can be obtained in Equation 9: d,. (A:) = log I χ. (A:) 1 -log I χ. (A:) I

= [d_i(k,0),---,d,(k,M-l)]^T ,

(9) where d,(£, -l) = 0 and C_M - . Similarly, for the jth channel, in

^•Ίχ(Μ-Ι)

Equation 10:

d,(*) = CVloglx.(*)l-loglx.(*)l

(10) We note that when computing the STFT, the data frame length must be selected such that N < L_min , where L_min is the length of z,.(fc) corresponding to the signal having the least number of "leading zeros". This constraint is imposed to ensure that the disparity between noise-only segments and the tap is included during the analysis process. In addition, due to the need for high frequency resolution, K is often large ( = 1024 ). In order to reduce the computational load, we have used N < K and each frame of x,(m) is padded with K -N zeros to make its length equal to K . For a sufficiently large value of N we can assume that the magnitude of the DFT coefficients between adjacent frames for a particular frequency bin to be almost equal, i.e. , log I h (k, m) \~\ h_j (k, m + l) \ , except at frame index m = L + τ^ -Ι . This effect is illustrated in Fig. 5 and 6 (a), using N = 256 , where there are two distinct energy levels across time for each of the frequency. Since the energies of is often less than that of h_:'(k) , the onset can clearly be observed. We therefore note that since z_;.(fc) < h__j(k) in Equation 1 1 and 12:

(1 1 )

This implies that, in theory, there exists a peak at rf,.(£,L, -l) and ^.(fc,L. -1) .

Therefore, cross-correlating the sequences d,.(fc) and d_;.(fc) will give a peak at ηι = -τ_ν corresponding to the TDOA between i th and jth sensor outputs.

In practical scenarios, however, spurious peaks may exist in d_;.(£) due to the small differences between \h_j(k,m)\ and I h_j(k,m + 1)1 when m≠L_j+T_ij-l. This effect can be observed from Fig. 6 (b), where we added a value of one to elements in d,(42) for clarity. Unlike the significant difference which can be found between h_j(k,m)\ and \h_j(k,m + \)\ for m = L+r_i]-l , these spurious peaks are not consistent over different frequency bins. Hence to improve the robustness of the method, magnitude of elements in vector d_;(fc) are summed over a few adjacent frequency bins giving an xl vector in Equation 13.

¾(£)= J d_j(k + b)\

= [_¾.(0),...,a.( -l)f (13) where, Δ£ is the number of adjacent frequency bins and \d_j(k+b)\=[\d_j(k+b,0)\,...,\d_j(k + b,M -l)\f ■ The vectors a,(42) and a₅(42) are as shown in Fig.6 (c). For this illustrative example, we have used Δ£ = 8. Therefore, the cross correlation sequence between ¾(fc) and a_;.(fc) is given in

Equation 14 and 15: r, . =[r_lt] (-M +!),...,/:, (0), ... , r._j (M - l)f (14)

>1.ί(^τ) + ρ), τ = -Μ +l,...,0,...,M -I- (15)

The estimated TDOA between channels i and j is then given in Equation 16: arg max r (16)

Fig.6 (d) illustrates elements of the cross correlation sequences r₁₅ while r_u is plotted as a reference. The elemental differences between the peaks of these sequences provide t .

The TiF-HA method

It is expected that the signals generated by a finger lacks sharp signal transients. The method 300 for detecting the signal onsets generated by a finger begins by removing the DC offset of x,(n) by subtracting E[v,.(n)] from x_;(n), where v,(n) of length L_v<L_min is a segment of the data taken from the "leading zeros" portion of the signal x_:(n). In order to detect the onset of *,_·(«), the method requires a reference noise, in the TiF domain, whose power must be proportional to the noise power of x_:.(n). In order. to reduce the computational load and maximize the standard deviation of the Hermitian angles during the "leading zeros" period, we construct a KxM reference noise matrix of complex numbers in Equation

where a is the proportionality constant which controls the reference noise level, £_start and £_end are the starting and ending indices of the frequency bins and w₍.(fc,m) is the (k,m)th element of the reference noise matrix. Note that &_start and k_end may be selected to be in the range of resonance frequency of the transducer used. The elements of the reference noise matrix have a uniform distribution in the interval -0.5 to 0.5.

We note that the sensor output is noisy in Equation 18: x_i(n) = k(n)*s(n) + v_i(n) (18) where v,.(«) is the noise picked up by the /th sensor. Converting Equation (18) in to the TiF domain, the signal and the reference noise at the fcth frequency bin can be constructed in matrix form given in Equation 19:

/¾(/,m)

i{k,m)

h_i(k,m) + _j(k,m) + W:(k,m) (19)

0

Assuming the source s(n) is an impulse, at all the points in the TiF plane s_(k) will be a constant, whereas h.(k,m) , v (&,m) and vv.(£,m) are time and frequency dependent. During the "leading zeros" period of x,.(«), the first term at the right hand side of Equation (19) will be zero. The Hermitian angle, 9_n{k,m) , between two complex vectors y and u_i =[x₁(k,m),w_i(k,m)]^T at any time frequency point (k,m) is given in Equation 20:

where

are constants, the Hermitian angle between these vectors and another reference vector , will remain the same even if they are multiplied by any complex scalars.

In Equation (19), these complex scalars are v.(£,m) and w,.(fc,m) . However, since v^{k,m) and w,(£,m) are non-zero at all the points in the TiF domain and the amplitude of vv,.(^,m) is almost in the same range of ^m), the Hermitian angle between the resultant vector

[ (£,m), w(k,m)f =[1, Of (k,m)+[0,l]^Tw_i(k,m) and y will vary within the maximum range of Hermitian angle given by 0 to π/2 rad. Note that we have chosen uniform distribution for the elements vv,.(fc, ) in the reference noise matrix to get maximum variation in Hermitian angle during the "leading zeros" period of *,.(«). During the period where h_i(k,m)≠0 , i.e., from the onset point onwards, all the three terms in Equation (19) will be non-zero. Hence, the Hermitian angle between the reference vector y and [ ,.(fc,m), ¼-,.(£, m)f will be different at different points in the TiF plane. However, the variation in angle will be smaller due to the lesser influence of second and third term on the first term.

This fact is illustrated in Fig.7 where the reference vector y = [1 + V^T, Of Now it is obvious that the standard deviations of the Hermitian angles given by Equation 21 :

calculated across the frequency bins during the "leading zeros" period will be much higher than those during the remaining period of x, (n) . Hence, the detected transition point in o, (m) calculated over the complete period of the signal corresponds to the onset point of the signal. For example, Fig. 8 (a) and (c) shows x_; (n) and σ_;(/η) respectively, calculated as explained above. An offset of 128 samples observed between the transition point in o,. (m) and the onset point in ,. (n) is due to the STFT analysis and it will be same for all the signals whose STFT analysis window length is same and hence it can be ignored. The main advantage of using σ,.(τη) instead of directly using *,.(«) for the detection of the onset point is that unlike in x. [n) , whose amplitude can vary from -∞ to +∞ , the variation of 0_H (k, m) and hence that of o,. (m) will be only between 0 and π / 2 rad. This makes the selection of the threshold very easy and robust.

For the detection of the transition point of c, (m) , we use the histogram method.

The problem of detection of the transition point is now translated to the problem of finding a threshold value as shown in Fig. 8 (c) and 9 (c) which will separate the samples a^m) into two groups, one corresponds to the points in the "leading zeros" region of the signal and the other corresponds to the remaining points of the signal. In this method, histogram method is used to estimate the threshold to group the samples of σ,.(τη) into two, where the histogram of o,.(m) is calculated first, using B bins. The larger values of cr.(m) in the "leading zeros" period compared to rest of the period will result in a valley point in the histogram as shown in Fig. 8 (b) and 9 (b). This corresponds to the threshold value that is to be applied on to σ,θ) . For the estimation of the valley point and hence the threshold value, we proceed with the following steps. First, the samples in the valleys at both ends of the histogram are removed. This is done by clustering <7, (m) into two clusters using k-means clustering algorithm and discarding the samples in the histogram which are before and after the bins whose centres are nearest to the cluster centroids obtained from the k-means algorithm. We next obtain the cumulative histogram using the remaining samples as shown in Fig. 8 (d) such that corresponding to the desired valley point in Fig. 8 (b) the slope of the curve in Fig. 8 (d) will be zero. The x -coordinate of this saddle point in Fig. 8 (d) corresponds to the threshold value. For the detection of this saddle point, again we compute the histogram of the cumulative histogram for number of bins equal to the last sample value in the cumulative histogram. Then the bin centre of the bin which contains maximum number of samples represents the y -coordinate of the saddle point in Fig. 8 (d) and the corresponding x -coordinate represents the threshold value. For example, for the case shown in Fig. 8, the threshold value is 7.33. In the second histogram, Fig. 8 (f), if there are more than one bin with equal number of maximum samples, the histogram calculation will be repeated using number of bins equals half of the number of bins used in the previous step. This process will be repeated until there is only a single maxima. The value cr. (m) corresponding to the peak position of the cumulative frequency will then be used as the threshold. Now apply the estimated threshold on and find out the transition point as shown in Fig. 8 (e). The detected transition points are the points where the plot of crosses the threshold level. However, in some cases, as shown in Fig. 8 (e), the threshold may give more than one transition points due to the spurious spikes present in σ,. (τη) . In such cases, the spikes from o,. (m) are eliminated in such a way that for all the points where the transitions are detected, the sample of a,. (m) at those points are replaced by the mean of its two adjacent samples, provided both of them are higher or lower than the sample at the point where the transition is detected. This is to keep the samples at the actual transition point of o,. (m) unaltered as the adjacent samples at this point will be at the two sided of the threshold level. However, if the false transition detection is due to two or more adjacent spikes in o, (m) , the above method will not remove all the spikes. To overcome this problem, in addition to the above method, o,. (m) is smoothened by

3-point moving average filter. After the spurious spike elimination stage, the threshold detection process is repeated. This process is repeated until a single transition point is detected or a certain number of cycles is reached. If the method fails to find out the threshold within a certain number of cycles, the corresponding sensor output will be discarded. Now the detected transition point will be considered as the onset of the sensor output and only those sensor outputs whose onset points were estimated successfully will be used for the TDOA estimation. In this work we have used a maximum of 10 cycles. Since we have a clear estimation of the onset points of the sensor outputs, the TDOA, in terms of number of samples, between two sensor outputs will be simply the difference between indices of the onset points. Unlike TiF-DA, here the cross correlation calculation is not required. It may be noted that when the SNR of the sensor output is high, such as impacts produced by a stylus, the chance for the spurious spikes in σ_;(0 will be very low, as shown in Fig. 10 (b).

Having known the TDOA of the source signal at different sensors with respect to the reference sensor, an algorithm may be used to determine the source localization. A first example algorithm uses pre-estimated velocity of the signal and another example algorithm is a closed form solution where both the source location and velocity of the signal are estimated simultaneously. In the first method, for pre-estimation of the signal velocity, source signals are generated by tapping at pre-defined locations on the plate surface and the signals are recorded using the sensors located at known positions. Let /,. be the distance between the source and the * sensor, c{f) is the velocity of the signal at frequency / and T_y is the estimated TDOA between the z^'th and the j^th sensors. The velocity of the signal at frequency / can then be estimated by minimizing the sum of the squared errors between TDOA calculated using the estimated velocity and r,₇ , in Equation 22:

by using any one of the optimization algorithms. The Levenberg-Marquardt algorithm [15] may be used. After estimation of c(f ) , the unknown source position, (x_s , y_s ) , can be estimated by minimizing the error function in Equation 22:

(x_s, jJ = argmin 2_J

c(f) (22) where U, , >', ) is the /^th sensor co-ordinate and c(f) is the estimated velocity.

Let us call this method for source localization as "calibration method" as it requires pre-calibration. In this work we use calibration method for all the experiments.

The ( -,., >',.) co-ordinates are converted into pixel co-ordinates for output to a computer. It is envisaged that the output may also contain the amplitude, type and duration of the touch. Additionally, it is envisaged that the output may be of a standard e.g. a USB output suitable for a USB mouse.

The sampling frequencies of the sensor(s) affect the resolution, as do the number of sensors, their positions and location of the touch. A higher sampling frequency may yield better accuracy. For example, if the signal velocity is 1500 m/s and sampling frequency 96000 Hz, a single sample delay may correspond to a distance of 1 .5 cm. By using optimization for source localization, a single sample error may result less than 1 .5 cm of error.

The potential applications include health care, home automation, contact-less switches, design of next generation musical instruments, foot step detection in stage show, point of impact localization in firearm training etc.

The surface for the methods 200 and 300 may be a substantially homogenous surface or may be a locally heterogeneous (i.e. "quasi-heterogeneous") surface. In a substantially homogenous surface, the velocity of a signal is substantially the same in all the directions due to uniform density of the material. In a locally heterogeneous surface, the material(s) of the surface is heterogeneous for an area, but the non-uniformity of the material is substantially consistent over the whole surface; the velocity of a signal in heterogeneous portions of the surface is different for different signal directions as the density of the medium is different at different points. An example of a locally heterogeneous material that is substantially globally homogenous is wood. Whilst exemplary embodiments of the invention have been described in detail, many variations are possible within the scope of the invention as claimed as will be clear to a skilled reader. It is envisaged that other transforms (e.g. wavelet transform) may be used in Steps 201 and/or 302 instead of sliding window DFT. Additionally, it is envisaged that any vibration sensor and/or piezo element may be used for the sensor(s), e.g. the piezo element ABT-441-RC from Multicomp

TM

References

[1 ] S. J. V. Nichols, "New interfaces at the touch of a fingertip," Computer, vol. 40, pp. 12 - 15, Aug. 2007. [2] P. I. S. Lei and A. K. Y. Wong, "The multiple-touch user interface revolution," IT Professional, vol. 1 1 , pp. 42 - 49, Jan. 2009.

[3] W. Sachse and Y. Pao, "On the determination of phase and group velocities of dispersive waves in solids," Journal of Applied Physics, vol. 49, pp. 4320 - 4327, Aug. 1978.

[4] A. Sulaiman, K. Poletkin, and A. W. H. Khong, "Source localization in the presence of dispersion for next generation touch interface," in Proc. Int'l Conf. on Cyberworlds, pp. 82 - 86, 2010.

[5] K. Poletkin, X. Yap, and A. W. H. Khong, "A touch interface exploiting the use of vibration theories and infinite impulse response filter modeling based localization algorithm," in Proc. IEEE Int'l Conf. on Multimedia and Expo, pp. 286 - 291 , 2010.

[6] R. K. Ing, N. Quieffin, S. Catheline, and M. Fink, "In solid localization of finger impacts using acoustic time-reversal process," Applied Physics Letters, vol. 87, pp. 1 - 3, Nov. 2005.

[7] D. T. Pham, Z. Wang, Z. Ji, M. Yang, M. Al-Kutubi, and S. Catheline, "Acoustic pattern registration for a new type of human-computer interface," in Proc. Virtual Conf., IPROMs, 2005.

[8] C. H. Knapp and G. C. Carter, "The generalized correlation method for estimation of time delay," IEEE Trans. Acoustics, Speech and Signal Processing, vol. 24, pp. 320 - 327, Aug. 1976.

[9] H. Jeong and Y. Jang, "Wavelet analysis of plate wave propagation in composite laminates," Composite Structures , vol. 49, pp. 443 - 450, Aug. 2000.

[10] H. Jeong and Y. Jang, "Fracture source location in thin plates using the wavelet transform of dispersive waves," IEEE Trans. Ultrasonics, Ferroelectrics and Frequency Control, vol. 47, pp. 612 - 619, Aug. 2002.

[1 1] Y. T. Chan, T. K. C. Lo, and H. C. So, "Passive range-difference estimation in a dispersive medium," IEEE Trans. Signal processing , vol. 58, pp. 1433 - 1439, Mar. 2010.

[12] J. Zheng, K. W. K. Lui, and H. C. So, "Accurate three-step algorithm for joint source position and propagation speed estimation," Signal Processing, vol. 87, pp. 3096 - 3100, Dec. 2007.

[13] C. W. Reed, R. Hudson, and K. Yao, "Direct joint source localization and propagation speed estimation," in Proc. IEEE Int'l Conf. on Acoustic Speech and Signal Processing, pp. 1 169 - 1 172, 1999.

[14] S. M. Ziola and M. R. Gorman, "Source location in thin plates using cross-correlation," J. Acoust. Soc. Am., vol. 90, pp. 2551 - 2556, Nov. 1991 .

[15] K. Levenberg, "A method for the solution of certain non-linear problems in least squares," The Quarterly of Applied Mathematics, vol. 2, pp. 164 - 168, 1944.

[16] M. Akay, S. Member, I. Marsic, A. Medl, and G. Bu, "A system for medical consultation and education using multimodal human/machine communication," IEEE Trans. Information Technology in Biomedicine, vol. 2, pp. 282 - 291 , Dec. 1998.

[17] A. Crevoisier and P. Polotti, "Tangible acoustic interfaces and their applications for the design of hew musical instruments," in Proc. of the Int'l Conf. on New Interfaces for Musical Expression, pp. 97 - 00, 2005.

Claims

1 . A method of estimating time difference of arrival (TDOA) of a stylus impact made on a solid surface comprising:

analysing a received signal in the time-frequency domain,

performing a difference- averaging process, and

estimating the TDOA when the rate of change of the difference- averaged signal at the onset is high.

2. A method of estimating time difference of arrival (TDOA) of a finger impact made on a solid surface comprising:

analysing a received signal, in the time-frequency domain,

performing a Hermitian Angle process, and

estimating the TDOA using a transition of the standard deviation of the Hermitian Angles.

3. A method of determining the origin of non-standardised surface acoustic waves in a substantially homogenous and/or quasi-heterogeneous surface comprising:

passively detecting the acoustic waves in the surface;

Short Time Fourier transforming the acoustic waves into time frequency domain data;

determining the time difference of arrival (TDOA) of the acoustic waves based on the data.

4. The method in claim 3 wherein the TDOA is used to estimate the origin location.

5. The method in claim 3 or 4 wherein the detecting comprises providing signals representative of the surface acoustic waves at a plurality of locations adjacent to a periphery of the surface.

6. .The method in any one of claims 3 to 5 wherein the Short Time Fourier transforming comprises a sliding Discrete Fourier Transform (DFT).

7. The method in any one of claims 3 to 6 wherein selected frequency bins from the time frequency domain data are used to determine the TDOA.

8. The method in claim 7 wherein the frequency bins are selected based on the resonance frequency of a transducer for passively detecting the acoustic waves.

9. The method in claim 8 further comprising differentiating the selected frequency bins with respect to time.

10. The method in claim 9 further comprising difference-averaging the differentiated values for the selected bins to obtain a corresponding average value of each time sample.

1 1. The method in claim 10 further comprising correlating the difference-averaged values of two or the plurality of locations to estimate the TDOA.

12. The method in any one of claims 3 to 1 1 wherein the unspecified surface acoustic waves correspond to a sharp impact.

13. The method in any one of claims 3 to 8 further comprising removing the DC offset.

14. The method in claim 13 further comprising determining a Hermitian angle using the sensor signal and a reference noise matrix.

15. The method in claim 14 further comprising histograming the standard deviation of the Hermitian angle.

16. The method in claim 15 further comprising K-means clustered the histogram to obtain cluster centroids.

17. The method in claim 16 further comprising smoothening or filtering each centroid.

18. The method in claim 17 further comprising iterating clustering and smoothening until there is a single transition in the histogram.

19. The method in any one of claims 13 to 18 wherein the unspecified surface acoustic waves correspond to a soft impact.

20. A system configured to store and execute the method of any preceding claim.

21. A touch detection apparatus comprising:

a plurality of passive transducers configured to detect signals representing surface acoustic waves in a substantially homogenous and/or quasi-heterogeneous surface; and

a processor configured to determine the origin of the acoustic waves on the surface based the time difference of arrival of the acoustic waves in the time frequency domain between signals for each sensor.