Ultrasonics 55 (2015) 15–25
Contents lists available at ScienceDirect
Ultrasonics journal homepage: www.elsevier.com/locate/ultras
Wave separation: Application for arrival time detection in ultrasonic signals Patrick Avanesians ⇑, Moe Momayez 1 Department of Mining and Geological Engineering, University of Arizona, United States
a r t i c l e
i n f o
Article history: Received 19 August 2013 Received in revised form 24 June 2014 Accepted 14 August 2014 Available online 26 August 2014 Keywords: Wave separation Nonlinear Decomposition Time-frequency analysis Arrival time
a b s t r a c t A method to detect and accurately measure the arrival time of wave packets in ultrasonic signals using a nonlinear decomposition technique is presented. We specifically address the problem of extracting events that are not well separated in the time, space and frequency domains. Analysis of complex ultrasonic signals, even those composed of poorly separated echoes, provided exceptional estimates of the desired time of arrival, from the media under investigation. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction The ultrasonic pitch-catch technique is one of the common tests for inspecting the quality of materials. Two piezoelectric transducers separated by a certain distance are placed on a free surface of the test object, where one transducer is used to transmit an ultrasonic pulse inside the material while the other sensor captures reflections due to differences in the acoustic impedance between two media. One goal of conducting ultrasonic tests is to measure the thickness of the material under investigation or to locate defects/objects inside the material. Typically, the equation d = (m t)/2, where m is the ultrasound speed and t is the arrival time is used to determine the distance from the transducer to the target. Thus, distance is related directly to arrival time if velocity is known. Fig. 1 shows two signals recorded on the same test object. Tracing the signal in time, where the wave typically travels through the material, the recorded amplitude is zero or close to zero. The direct wave (one that travels near the surface of the material when the transmitting and receiving transducers are place close to each other) is usually the first one to reach the receiver as shown by the change in the amplitude in the recorded signal (in the case of Fig. 1 the first arrival time is around 45 ls). However, in ultrasonic
⇑ Corresponding author. Tel.: +1 818 331 3161. E-mail addresses: [email protected] (P. Avanesians), moe.momayez @arizona.edu (M. Momayez). 1 Co-author. Tel.: +1 520 621 6580. http://dx.doi.org/10.1016/j.ultras.2014.08.019 0041-624X/Ó 2014 Elsevier B.V. All rights reserved.
tests, this arrival time is occasionally used to get an estimate of the wave velocity since the distance between the transducers is known. The typical ultrasonic signal is composed of multiple modes and reflections with different frequencies. In testing concrete, since the reflections from defects inside the slab are rarely well separated in time, it becomes highly problematic to identify the arrival time from the desired targets. The detection of the arrival times from multiple scatterers inside the medium has been the subject of intense research over many decades. Different methods have been proposed to decompose the recorded ultrasonic waves and measure the associated arrival times. The following is an inventory of the main techniques that have become popular among the nondestructive testing practitioners: 1. Application of the traditional quadrature-demodulation scheme to suitably extract the envelope of the main echo and locate its onset [1]. 2. Correlation between input and output for high signal-tonoise ratio waveforms and L2-norm with low noise level [2]. 3. Application of the square-root unscented Kalman filter (SRUKF) to identify the shape parameters of an ultrasonic echo envelope [3]. 4. Application of Synthetic Aperture Focusing Technique (SAFT) where a constructive interference procedure enhances the lateral spatial resolution and signal-to-noise ratio. The arrival time is estimated from enhanced images of the interior of the medium [4].
16
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
71
10
(a)
4500
6
4000
4
3500
Frequency [ Hz ]
Amplitude
8
2 0 -2
3000 2500 2000
-4
1500
-6
1000 500
-8 -10
Model Space
5000
0 0
50
100
150
200
250
300
350
0
0.1
0.2
0.3
0.6
0.7
0.8
0.9
1
0.5
(b)
6
0.45 0.4
4
0.35
Frequency [Hz]
2
Amplitude
0.5
Time-Frequency Distribution
8
0 -2 -4
0.3 0.25 0.2 0.15 0.1
-6
0.05
-8
0
-10
0.4
Time [ seconds ]
Time (us)
0
50
100
150
200
250
300
350
0.1
0.2
0.3
400
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Time [seconds]
Time (us) Fig. 2. Model space and the associated Wigner–Ville transform. Fig. 1. (a) Well separated ultrasonic signal; (b) mixed ultrasonic signal.
5. Separation of non-linear and non-stationary ultrasonic signals by means of the empirical mode decomposition method to obtain the intrinsic mode functions used for signal reconstruction. The arrival time is then obtained from the reconstructed signal [5]. 6. Combination of both the improved self-interference driving technique and the optional optimization signal processing algorithms for measuring the arrival time [6]. In some ultrasonic signals, picking the first reflection (arrival time) from the desired target is an easy task since the direct wave and the main reflection are well separated in time (with little or no interference from other sources of ultrasonic scattering), and therefore can be discriminated visually. Fig. 1(a) shows a signal in which the direct wave and the main reflection are relatively well separated and the arrival time was correctly determined by visual inspection. Fig. 1(b) on the other hand, depicts a commonly encountered situation in practice where multiple waves are superimposed, turning the detection of the main reflection into a challenging problem. In these cases, signal decomposition has shown promise and is used in this study to separate the signal into different modes in order to extract the desired reflected wave pack. More specifically, the synchrosqueezed transform (SST), a new signal decomposition method, is shown to provide useful results.
Conventional time–frequency representations, such as the short-time Fourier transform (STFT), the Wigner–Ville transform (WVT), the Wavelet Transform (WT), and the Empirical Mode Decomposition (EMD) method, perform poorly when individual components of the recorded signal are not sufficiently separated in the time–frequency domain. The synchro-squeezed transform builds on the philosophy of the EMD method using a more robust theoretical foundation and employs a different approach for the construction of a signal’s constituent components. In fact, synchro-squeezing can be overlaid with many commonly used time– frequency methods [7–10]. SST is a combination of the wavelet transform and reassignment technique [11]. 2. Theoretical background The simplest form of a time–frequency representation is a spectrogram created using a short-time Fourier transform (STFT). This transform employs a short-time window that limits the evaluation of the Fourier transform to some specified region along the time index. It is well-known that STFT spectrograms suffer from the time–frequency resolution trade-off, that is, a shorter window provides better resolution in the time domain, and the larger bandwidth of the window’s spectrum affords a poorer resolution in the frequency domain. The Wigner–Ville transform (WVT) takes advantage of the powerful concept of match filtering, effectively employing the time-reversed copy of the signal as the window to carry out a short-time window Fourier transform. A
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
17
Fig. 3. (a) Concrete slab with a 15° bottom surface; (b) rectangular grid on concrete slab; (c) ultrasonic transmitter and receiver; (d) ultrasonic transducers while conducting the experiment; (e) typical recorded signal.
number of properties make the WVT valuable for decomposing complex signals. For example, the analysis is performed using only the transform’s ability to be entirely localized in the case of frequency modulated signals. The trade-off in using WVT is the manifestation of prominent interference (oscillations for the most part) between any two interacting components of the signal [12]. The effect is illustrated in Fig. 2 below. Smoothing somewhat enhances the readability of the transformed signal through the application of the pseudo-Wigner–Ville distribution which is essentially a windowed version of the WVT. Here, the cross-terms are attenuated at the cost of spreading signal energy along the frequency axis, adversely affecting the time–frequency localization of the signal content. To optimize the tradeoff between attenuating oscillations and focusing the signal energy in the time–frequency domain, the coefficients (t, w) of the WV distribution are reassigned to new points ~ which are the center of gravity of the signal energy ð~t; wÞ
distribution around (t, w). The reassignment method was first introduced by Kodera et al. [13] to overcome the spectrogram’s drawbacks such as biased estimation of the signals instantaneous frequency, group delay, and the Gabor–Heisenberg inequality which controls the locality in either the time or frequency domain, but never both [14]. The concept of reassignment can also be applied to improve the spectral resolution obtained from wavelet transforms. In this context, the process has been dubbed synchro-squeezing, and was first introduced to analyze speech signals [15]. Subsequent studies [8– 11,16] showed that synchro-squeezing performance is superior and could be considered as a viable alternative to the EMD method. In this paper, we show that the synchro-squeezing method could also be used to decompose ultrasonic signals into their constituent wave packets, providing the means to detect and determine the arrival time of the desired echoed pulses from inside the medium under investigation. This is especially valuable since
18
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
Original Signal 10
26
Time - Frequency - Amplitude
200
162
(a)
8
(a)
8
180 6
6 160
Frequency (KHz)
4
Amplitude
2 0 -2 -4
Time index for maximum amplitude 121.6
-6
4 2
140
0
120
-2
100 -4
80
-6
60
-8
-8 -10 0
50
100
150
200
250
300
350
400
50
Time (us)
150
200
250
300
350
400
Time (us)
Truncated and Padded signal
Synchrosqueezed Wavelet Transform
10
(b)
8
(b)
180
4.5
160
6
4
140
4 2
Frequency (KHz)
Amplitude
100
0 -2 -4
3.5
120
3
100
2.5 2
80
1.5
-6
1
-8 -10 0
60 50
100
150
200
250
300
350
Central Frequency = 0.8
0.5
400 0
Time (us)
50
100
150
200
250
300
350
400
Time (us) Fig. 4. (a) Original signal with lower and upper time interval and time index for maximum amplitude; (b) truncated and padded signal.
the shape of the mother wavelet is tuned to closely match the ultrasonic pulse waveform used in the inspection process. The following is a brief description of the synchro-squeezing transform (SST). If s(t) is the desired signal, its continuous wavelet transform (CWT), Ws(a, b) is given by:
1 W s ða; bÞ ¼ pffiffiffi a
Z
1
1
sðtÞw
tb dt a
ð1Þ
where w⁄ is the complex conjugate of the mother wavelet (a continuous function in both the time and frequency domains), and b is the time shift applied to the mother wavelet which is scaled by a. We can think of the CWT as the cross-correlation of the signal s(t) with pffiffiffi the scaled and time shifted wavelet wððt bÞ=aÞ= a. This cross-correlation provides a measure of similarity between the signal and the scaled shifted wavelet. The wavelet coefficients Ws(a, b), however, frequently spread out along the scale dimension a, producing a blurred image in the time-scale domain. The next step in the process involves the transformation of the points in the time-scale domain into the time–frequency domain for all Ws(a, b) – 0. Different expressions may be used to handle the transformation. We adopted the following relationship [16]:
1 d xs ða; bÞ ¼ arg½W s ða; bÞ 2p db
ð2Þ
Fig. 5. (a) TFA distribution; (b) SST distribution at central frequency of 0.8.
where arg () denotes the unwrapped phase of the complex coefficients and the 1/2p factor is introduced to convert from circular to normal frequencies. Eq. (2) is the expression for instantaneous frequency (derivative of the unwrapped phase of the continuous wavelet transform). When the coefficients xs(a, b) are determined, we select the desired frequencies xk and form the frequency interval [wi–Dw/2, wi + Dw/2] where Dw = wi–wi1. The synchrosqueezing transform Ts(wi, b) is computed only at the centers wi of the frequency range as:
T s ðwi ; bÞ ¼
1 Cw a
X
3=2
W s ðaj ; bÞ aj
Daj
ð3Þ
j:jwðaj ;bÞwi j6Dw=2
where Daj are the distances between two adjacent scales aj, and R1 ^ eÞ de, with wð ^ eÞ being the Fourier transform of the mother C w B 12 0 wð e wavelet. Eq. (3) indicates that synchro-squeezing takes place along the frequency/scale axis, producing a more accurate and detailed time–frequency representation. The synchro-squeezed transform is invertible and the original time domain signal can be reconstructed (with no loss of information) using [11]:
! X X X sðtÞ ¼ Re T s ðwi ; bÞ ¼ T s ðwi ; bÞ cos arg T s ðwi ; bÞ i i i
ð4Þ
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
19
form of the wavelet because, compared to other mother wavelets, it represents most accurately the actual ultrasonic pulse that probed the interior of the concrete material in this study:
Input Signal 50
t2 2 1 wðtÞ ¼ Bf 0 ei2pf 0 t e2ð2pf 0 Þ e 2
0
ð5Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffi where Bf 0 is a normalization factor equal to 2=pf 0 and f0 is the socalled central frequency, which determines the relationship between the time and the scale (or frequency) resolutions of the wavelet transform.
Amplitude
-50
-100
3. Proposed technique -150
-200
0
50
100
150
200
250
300
350
400
Time (us)
Time-Frequency-Amplitude 300
44.4 23.9
Frequency (KHz)
250
3.52 -16.9
200
-37.3 150
-57.8 -78.2
100
-98.6 -119
50
-139 0
-160 0
50
100
150
200
250
300
350
400
Time (us) Fig. 6. Input signal and its TFA distribution.
Synchrosqueezed Wavelet Transform 180KHz
180
Frequency (KHz)
160
140KHz
140
120KHz
120
100
90KHz 50
100
150
Time (us) Fig. 7. Selection of frequency bands.
The choice of the mother wavelet is based on many factors, and should reflect the type of features present in the signal. A commonly used function is the Morlet wavelet, consisting of a plane wave modulated by a Gaussian function. The Morlet wavelet is non-orthogonal, complex and offers a good tradeoff between detecting oscillations and peaks or discontinuities making it especially suited for analyzing ultrasonic signals. We used the following
We describe a method to unscramble closely spaced events commonly encountered in ultrasonic testing of materials. This condition occurs when the direct wave and the reflection from a scatterer superimpose and create an interference in the recorded signal because the events are not well separated in time, space, and frequency domains, or in other words, the length of the ultrasonic pulse is larger than the time of flight between the echoes. We specifically tackle the case where the detection and estimation of the arrival time of the first reflection is needed. Our approach uses both the synchro-squeezed transform and the measurement of instantaneous frequencies to decompose a signal into its components. The arrival time is then determined for the desired reflection. First, the arrival time associated with the highest peak is determined from the time domain signal. Let us call this tmax. Since both the input and output signals have a Gaussian shaped envelope, we know that the actual arrival time should be before tmax. Next, a certain range is selected to truncate the signal and make the decomposition procedure go faster. The lower limit is chosen at some point before the onset of the direct wave arrival, while the upper limit is a point after tmax where the change in the shape of the envelope of the signal is discernable, or tmax + tmax/3 in extremely mixed waveforms. The signal is then truncated and padded with zeroes to the same length as the original signal. To decompose the signal and measure the arrival time, the following steps are performed. Calculate the instantaneous frequencies and plot the time–frequency–amplitude (TFA) distribution. Note that the instantaneous frequency is the derivative of the unwrapped instantaneous phase, which is obtained by performing a Hilbert transform on the truncated s(t) signal. The amplitude values in the TFA distribution are taken from the original time domain signal and color coded for better readability. Next a synchrosqueezed transform is carried out on the truncated signal using an appropriate central frequency. The choice of the central frequency is crucial and will affect the outcome of the analysis. The rule we propose for the selection of the central frequency is based on visual interpretation of both the SST and time–frequency– amplitude distributions. The instantaneous frequency plot is compared with the SST and the optimum central frequency value is found from the closest match between the two distributions. The best match between the TFA and SST distributions is obtained by changing the value for the central frequency in the SST and comparing the resultant maximum amplitudes (values in the third dimension) in the TFA and SST plots. The time and frequency values corresponding to the maximum amplitudes should be the same or very close for the optimal central frequency. In almost all cases, the central frequency is found to be between 0.5 and 1. We continue the process by defining the frequency bands needed for signal decomposition. Since the first reflection normally has the strongest echo, we look for the frequency range in both the TFA and SST plots where the highest amplitudes (absolute value) are located.
20
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
Decomposed Signal : Fmin=50000 Fmax=90000
Decomposed Signal : Fmin=90000 Fmax=120000
0.2
2
A
B
1.5
0.15
1
0.1
0.5
Amplitude
Amplitude
0.05 0 -0.05
0 -0.5 -1
-0.1
-1.5
-0.15 -0.2
-2
0
50
100
150
200
250
300
350
-2.5
400
0
50
100
150
Time (us)
200
250
300
350
400
Time (us)
Decomposed Signal : Fmin=120000 Fmax=140000
Decomposed Signal : Fmin=140000 Fmax=180000
8
3
C
6
D 2
4 1
Amplitude
Amplitude
2 0 -2
107
0
-1
-4 -2
-6 -8
0
50
100
150
200
250
300
350
400
-3
0
50
100
Time (us)
150
200
250
300
350
400
Time (us)
Fig. 8. Decomposed signal and associated waveforms from defined frequency bands.
Time-Frequency-Amplitude 9.59
300 7.67
Sudden Change in Frequency Value
250
5.75 3.83
Frequency (KHz)
4 200
1.91
150
1
-0.00497
3
2
-1.92
100
-3.84 -5.76
50 -7.68
0 40
60
80
100
possible that the main frequency band is selected incorrectly and the resultant waveform does not correspond to the actual reflection. To mitigate those special circumstances, the two methodologies described below can be used to validate the results.
120
140
160
-9.6
Method 1: In the time–frequency–amplitude plot, abrupt changes in frequency are indicative of different wave types detected by the sensor. The last peak before the cluster of high amplitude values corresponds to the time index of the main reflected wave. Method 2: In the time–frequency–amplitude plot, select the frequency range where the maximum amplitude values are located (the cluster of red and blue dots). Use this frequency range to reconstruct a time-domain waveform from the original time-domain signal. In other words, amplitudes outside of this range are removed from the original signal s(t).
Time (us) Fig. 9. Sudden changes in the time–frequency–amplitude distribution.
Once the frequency bands are selected, we extract the waveform from each band by performing an inverse synchro-squeezed transform of the SST coefficients. The reconstructed time-domain signal with the highest energy is the waveform associated with the main reflection from the discontinuity. In some cases, it is
4. Experimental setup The proposed decomposition method was tested by analyzing signals recorded on three concrete slabs, each with a surface area of 1 1 m2 and a sub-horizontal inclined bottom surface of 5°, 10°, and 15° respectively. Fig. 3(a) shows the slab with a 15° inclined bottom surface. The concrete used in this work had a max-
21
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
Reconstructed Signal : Fmin=90000 Fmax=180000
Time-Frequency-Amplitude 10
300
8
106.3
6 4
200
2
Amplitude
Frequency (KHz)
250
150
100
0 -2 -4
50 -6
107
-8
0 40
60
80
100
120
140
160 -10
Time (us)
0
20
40
60
Decomposed Signal : Fmin=120000 Fmax=140000 8
100
120
140
160
Fig. 12. Signal in time domain after removing amplitudes value outside certain frequency range.
C
6
80
Time (us)
Reconstructed Signal : Fmin=90000 Fmax=180000
4
10
Amplitude
2
8
0
6
-2
4
107
2
Amplitude
-4 -6 -8
0
50
100
150
200
250
300
350
0 -2 -4
400
Time (us)
-6
Fig. 10. Comparison of frequency changes with reconstructed waveform for validating Method 1.
107
-8 -10
0
20
40
60
80
100
120
140
160
Time (us) Decomposed Signal : Fmin=120000 Fmax=140000 8
Time-Frequency-Amplitude plot 250
8
6
6
4
180KHz
200
2 150 0
90KHz
100
-2 -4 -6
50
2
Amplitude
Frequency (KHz)
4
C
0 -2
107
-4 -6
-8
-8 20
40
60
80
100
120
140
160
Time (us) Fig. 11. Selection of the frequency band in the TFA distribution for validating Method 2.
0
50
100
150
200
250
300
350
400
Time (us) Fig. 13. Validation of the main reflection wave packet arrival time using Method 2 for signal S15.
22
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
15
Time-Frequency-Amplitude
122 Original Signal
26
200 180
10
160
5
5
Frequency (KHz)
Amplitude
10
0
-5
Time index for maximum amplitude 94.68
-10
-15
140 0
120
-5
100 80
-10
60 -20 0
50
100
150
200
250
300
350
400
0
50
100
Time (us)
150
200
250
300
350
400
-15
Time (us)
Tuncated and Padded Signal
Synchrosqueezed Wavelet Transform
15
200
7
180
10
160
6
140 5
Frequency (KHz)
Amplitude
5
0
-5
120 4
100 3
80 2
-10 1
60
Central Frequency = 0.6
-15 0
-20
0
50
100
150
200
250
300
350
50
100
250
300
350
0
Fig. 15. (a) TFA distribution; (b) SST distribution at central frequency of 0.6.
Fig. 14. (a) Original signal with lower and upper time interval; (b) truncated and padded signal.
Synchrosqueezed Wavelet Transform 160
150KHz 140
Frequency K(Hz)
imum grain size of 20 mm and compressive strength of 32 MPa. It was made by a local company with the same properties as the normal concrete used in constructions to simulate a ‘‘real world’’ concrete structure. The incline surface at the bottom of the slabs was built using a layer of wood with a thickness of 5 mm. The slabs simulate wide inclined cracks in large concrete structures such as dams where water can seep through and expand the cracks resulting in significant damage. Data were collected from rectangular grids mapped in the middle of slabs with 10 mm spacing between the points as shown in Fig. 3(b). The size of the grid was different for each slab and varied from 21 to 27 cm in both x and y directions, corresponding to a depth of 100–300 mm, as shown in Fig. 3(a). The tests were conducted in pitch–catch mode using two ultrasonic transducers, one as a transmitter and the other as a receiver. Fig. 3(c) shows ultrasonic transducers used b this experiment and Fig. 3(d) demonstrates ultrasonic transducers conducting experiment. The transducers had a nominal frequency of 140 kHz and an actual bandwidth of 400 kHz. A pair of flat acrylic plates were attached to the surface of the transducers as a protection measure against scratch and damage. The thickness of the plates are 5 cm, and together, they introduce an additional delay of 26 ls in the signal.
200
Time (us)
400
Time (us)
150
130KHz 120
100KHz
100
40
60
80
100
120
140
Time (us) Fig. 16. Selection of frequency bands.
A two-channel data acquisition and function generator was used to input a 140 kHz square pulse and to capture the reflected waveforms. The data were collected at a high sampling frequency of 25 MHz using 12-bit resolution. For each waveform, 10,000 samples with a peak-to-peak amplitude of 12 V were collected. The signal captured from the test medium was passed through a signal
23
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
Decomposed Signal : Fmin=10000 Fmax=130000
Time-Frequency-Amplitude
15
250
10
69 200
Frequency (KHz)
Amplitude
5
0
-5
68
150
100
50
-10
0 -15
20 0
50
100
150
200
40
60
80
100
120
140
Time (us)
250
Time (us)
Decomposed Signal : Fmin=10000 Fmax=130000 15
Decomposed Signal : Fmin=130000 Fmax=150000 2 10 1.5 1
5
Amplitude
Amplitude
0.5 0
0
-5
-0.5
68
-1
-10
-1.5 -15 -2
0 0
50
100
150
200
250
300
350
50
100
Time (us) Fig. 17. Decomposed signal and associated waveforms from defined frequency bands.
150
200
250
Time (us)
400
Fig. 18. Comparison of frequency change with reconstructed waveforms for validating Methods 1.
Time-Frequency-Amplitude 200
5. Data analysis For the sake of brevity, the analysis results for two signals are presented in this section. The first signal is from the 15-degree slab (see Fig. 3 above) at a known depth of 16.69 cm. This gives an actual arrival time of 81.41 ls, using the velocity of 4100 m/s. As explained in the experimental setup section, there is a delay of 26 ls due to the thickness of the guides, and therefore the total time we should get after signal decomposition is 107.41 ls (81.41 + 26). We name this signal S15.The second signal, S05, is taken on the 5-degree slab at a known depth of 8.5 cm. The arrival time for this signal is 67.47 ls (41.47 + 26).The analysis steps are explained in detail for the S15 signal only. For signal S05, we only present the results.
180
10
160
5
Frequency (KHz)
conditioner with the appropriate low-pass filter before feeding it to the data acquisition system. Fig. 3(f) shows a screenshot of the software used in this work displaying a typical waveform obtained from the test slab and recorded by the data acquisition system. Calibration tests on five cylindrical concrete samples provided an average velocity of 4100 m/s. This value is used to calculate the depth of the inclined bottom surface.
140
130KHz 0
120
100KHz
100 80
-5
-10
60
-15 20
40
60
80
100
120
140
Time (us) Fig. 19. Selection of the frequency band from TFA distribution for validating Method 2.
The first step is the selection of tmax which is 121.6 ls for S15. We take as the lower limit, the time index of 26 ls and as upper limit, the time index 162 ls. Fig. 4(a) shows the signal and the
24
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25 15
Reconstructed Signal : Fmin=10000 Fmax=130000
10
10
5
Amplitude
Amplitude
5
0
-5
0
-5
68 -10
-10
-15
-15 40
60
80
100
120
93.93
0
50
100
150
200
250
300
350
400
250
300
350
400
Time (us)
140
Time (us) 5
Decomposed Signal : Fmin=10000 Fmax=130000
4
15
3 2
10
Amplitude
1
Amplitude
5
0
0 -1 -2 -3
-5
102.8
-4
68
-5
-10
0
50
100
150
200
Time (us) -15
0
50
100
150
200
250
Time (us) Fig. 20. Validation of the main reflection wave packet arrival time using Methods 2 for signal S05.
selected time interval. The signal is then truncated and padded to zeroes as shown in Fig. 4(b). Next, the instantaneous frequencies, and the synchro-squeezed transform (at different central frequencies) are computed, and the time–frequency–amplitude (TFA) and the SST distributions are plotted. The optimal value of 0.8 for the central frequency is found by comparing the maximum values in the TFA and SST distributions, and the best match is selected. Fig. 5(a) shows the TFA distribution and Fig. 5(b) presents the SST at the optimum central frequency for signal S15. Let is take a closer look at the SST plot and select the frequency band that encompass the largest SST coefficients and call it the main reflection. Also, the selected frequency band should include the main input signal frequencies (80–140 kHz in our case). Fig. 6 shows the input signal and its associated time–frequency– amplitude plot. Using the SST distribution shown in Fig. 7, we identify the main frequency range from 120 kHz to 140 kHz, and we define the following frequency bands for signal decomposition: 50–90 kHz, 90–120 kHz, 120–140 kHz, and 140–180 kHz. Once the waveforms are reconstructed in the time-domain for each of the defined frequency band, the waveform with the largest energy (sum of the amplitudes squared) is used to determine the arrival time. Fig. 8
Fig. 21. Checking accuracy of arrival time by subtracting tmax of two consecutive signals.
presents the reconstructed waveforms for each of the defined different frequency bands. From the previous analysis, we conclude that the time corresponding to the wave reaching the sensor 107 ls with a frequency bandwidth from 120 kHz to 140 kHz. The last step involved in determining the time of arrival is to verify that the reconstructed waveform carries the correct information as discussed in the proposed method section. 5.1. Method 1 We consult the time–frequency–amplitude distribution, shown in Fig. 9, and locate the frequency peak (peak #4) before the cluster of high amplitude values in the signal. The time index corresponding to the peak #4 equals to 106.3 ls which is very close to the actual arrival time of 107.41 ls. Fig. 10 shows the comparison of the sudden change in the time–frequency–plot with the decomposed signal (Fig. 8(c)) in order to confirm that the correct reflection wave packet is selected. 5.2. Method 2 Let is take a look at the time–frequency–amplitude distribution and define the frequency range encompassing the highest amplitude values. In the case of signal S15 illustrated in Fig. 11, the lower
25
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
and upper limit of the frequency range are 90 kHz and 180 kHz. Then, amplitude values outside this range are removed from the signal s(t), and the desired arrival time is selected from the reconstructed waveform as shown in Fig. 12. A comparison of the reconstructed waveform in the band 120– 140 kHz (Fig. 4C) and the reconstructed waveform shown in Fig. 13, validates the results of the analysis and provides the correct time of arrival. The second signal from the 5-degree slab has an arrival time 41.47 ls corresponding to a depth of 8.5 cm. The actual total time is 67.47 ls (includes the delay from the acrylic guides). As in the previous example, we start with the time–frequency–amplitude distribution plot and compare this with the SST distribution in order to find the optimal central frequency. This step is followed by the selection of frequency bands, extraction of the SST coefficients and reconstruction of the waveform for each frequency band. The waveform with the highest energy content which also contains the main input frequency range is used to determine the time of arrival. Validation is performed by comparing the change in the TFA distribution with the reconstructed signal from the main frequency band (Method 1), and by removing the samples with amplitudes outside the defined frequency range from the original signal s(t) (Method 2). The reconstructed waveforms are plotted and compared to obtain the arrival time. Figs. 14–20 illustrate the preceding steps for S05.
yields
Dta1 a2 ¼ jta2 t a1 j ¼ j48:63 39:73j ! Dt a1 a2 ¼ 8:9 ls yields
Dt1max 2max ¼ jt 2max t 1max j ¼ j102:8 93:93j ! Dt 1max 2max ¼ 8:87 ls Dta1 a2 ffi Dt1max 2max 7. Conclusions The synchro-squeezed transform is an effective tool for decomposing many types of noisy, time-varying frequency, and nonlinear signals. In this paper, we presented a new technique for determining the arrival time of ultrasonic signals used in nondestructive testing of geomaterials. The method uses the synchro-squeezed transform to decompose the signal into constituent waveforms reconstructed from suitably defined frequency bands. We have shown that the time of arrival of the primary reflection (from a discontinuity inside the material) can be accurately measured from signals where the events are poorly separated in time and frequency. This technique could also be used to detect the arrival time of well separated time, frequency and space domain signals. It was successfully applied in the processing of lower frequency elastic waves from earthquake events and seismic exploration. References
6. Discussion Over 500 signals were analyzed in this study and arrival times were successfully determined. The results were verified with known depths of the bottom discontinuity in the concrete slabs. In the process, we observed that the choice of the central frequency affects the outcome of the analysis to some extent, especially with respect to the creation of the synchro-squeezed transform plot. As mentioned above, at present, the selection of the central frequency relies on a visual inspection and comparison of both the SST and TFA distributions. The best match between the instantaneous frequency and SST defines the optimal value for the central frequency. Care should also be exercised in the selection of the main frequency interval and the corresponding frequency bands used for the extraction of the SST coefficients and the subsequent reconstruction of waveforms in the time-domain. The suggested approach where the size of the frequency bands is taken to be equal to the width of the main frequency interval (obtained from the SST distribution) is a good initial guess and provided accurate estimates of the time of arrival for all analyzed signals. To obtain the arrival time for the first reflection, the original signal s(t), the time-instantaneous frequency–amplitude and the decomposed waveforms using the synchro-squeezed transform are analyzed simultaneously. The accuracy of the measured values can further be investigated by comparing two closely-spaced signals. For example, if time indexes corresponding to the maximum amplitude for the first and second time-domain signals are t 1max and t 2max , and the actual arrival times are ta1 and ta2 , then difference between t 1max and t2max must be equal to the difference between t a1 and t a2 . From Fig. 21 below, we read off the following values:
t a1 ¼ 39:73 ls; ta2 ¼ 48:63 ls; t 1max ¼ 93:93 ls; t2max ¼ 102:8 ls We then form the differences for the maximum amplitude time indexes and the measured arrival times:
[1] L. Angrisani, R.D.L. Moriello, Estimating ultrasonic time-of-flight through quadrature demodulation, IEEE Trans. Instrum. Meas. 55 (2006) 1. [2] M.A.H. Qayyum, High accuracy time of flight measurement using digital signal processing techniques for subsea applications, J. Sig. Inf. Process. 2 (2011) 330– 335. [3] Z.J. Yao, Q.H. Meng, M. Zeng, Improvement in the accuracy of estimating the time of-flight in an ultrasonic ranging system using multiple square-root unscented Kalman filters, Am. Ins. Phys. Rev. Sci. Instrum. (2010) 81. [4] Z. Hosseini, M. Momayez, F. Hassani, D. Lévesque, Detection of inclined cracks inside concrete structures by ultrasonic SAFT, AIP Conf. Proc. 975 (2008) 1298– 1304. [5] T.L. Chen, P. Que, Q. Zhang, Q.K. Liu, Ultrasonic signal identification by empirical mode decomposition and Hilbert transform, Am. Inst. Phys. Rev. Sci. Instrum. (2006) 76. [6] X.F. Wang, Z.A. Tang, A novel method for digital ultrasonic time-of-flight measurement, Am. Inst. Phys. Rev. Sci. Instrum. (2010) 81. [7] Daubechies I. Ten lectures on wavelets. CBMS-NSF Regional Conference Series in Applied Mathematics, 1992. [8] C. Li, M. Liang, A generalized synchro squeezing transform for enhancing signal time–frequency representation, Signal Process. 92 (2012) 2264–2274. [9] C. Li, M. Liang, Time-frequency analysis for gearbox fault diagnosis using a generalized synchrosqueezing transform, Mech. Syst. Sig. Process. 26 (2012) 205–217. [10] G. Thakur, H.T. Wu, Synchrosqueezing-based recovery of instantaneous frequency from nonuniform samples, SIAM J. Math. Anal. 43 (5) (2011) 2078–2095. [11] I. Daubechies, J. Lu, H. Wu, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl. Comput. Harmon. Anal. 30 (2) (2011) 243–261. [12] G. Matz, F. Hlawatsch, Wigner distributions (nearly) everywhere: timefrequency analysis of signals, systems, random processes, signal spaces, and frames, Sig. Process. 83 (2003) 1355–1378. [13] K. Kodera, R. Gendrin, C. Villedary, Analysis of time-varying signals with small BT values, IEEE Trans Acoust. Speech Sig. Process. 26 (1978) 64–76. [14] D. Gabor, Theory of communication, Inst. Electron. Eng. 93 (11) (1946) 429– 457. [15] Daubechies I, Maes S. A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models. Wavelets Med. Biol., CRC Press 1996, pp. 527–546. [16] D. Iatsenko, A. Stefanovska, P.V.E. McClintock, Nonlinear mode decomposition: a noise-robust, adaptive, decomposition method based on the synchrosqueezed wavelet transform, Appl. Comput. Harmon. Anal. (2012).
Contents lists available at ScienceDirect
Ultrasonics journal homepage: www.elsevier.com/locate/ultras
Wave separation: Application for arrival time detection in ultrasonic signals Patrick Avanesians ⇑, Moe Momayez 1 Department of Mining and Geological Engineering, University of Arizona, United States
a r t i c l e
i n f o
Article history: Received 19 August 2013 Received in revised form 24 June 2014 Accepted 14 August 2014 Available online 26 August 2014 Keywords: Wave separation Nonlinear Decomposition Time-frequency analysis Arrival time
a b s t r a c t A method to detect and accurately measure the arrival time of wave packets in ultrasonic signals using a nonlinear decomposition technique is presented. We specifically address the problem of extracting events that are not well separated in the time, space and frequency domains. Analysis of complex ultrasonic signals, even those composed of poorly separated echoes, provided exceptional estimates of the desired time of arrival, from the media under investigation. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction The ultrasonic pitch-catch technique is one of the common tests for inspecting the quality of materials. Two piezoelectric transducers separated by a certain distance are placed on a free surface of the test object, where one transducer is used to transmit an ultrasonic pulse inside the material while the other sensor captures reflections due to differences in the acoustic impedance between two media. One goal of conducting ultrasonic tests is to measure the thickness of the material under investigation or to locate defects/objects inside the material. Typically, the equation d = (m t)/2, where m is the ultrasound speed and t is the arrival time is used to determine the distance from the transducer to the target. Thus, distance is related directly to arrival time if velocity is known. Fig. 1 shows two signals recorded on the same test object. Tracing the signal in time, where the wave typically travels through the material, the recorded amplitude is zero or close to zero. The direct wave (one that travels near the surface of the material when the transmitting and receiving transducers are place close to each other) is usually the first one to reach the receiver as shown by the change in the amplitude in the recorded signal (in the case of Fig. 1 the first arrival time is around 45 ls). However, in ultrasonic
⇑ Corresponding author. Tel.: +1 818 331 3161. E-mail addresses: [email protected] (P. Avanesians), moe.momayez @arizona.edu (M. Momayez). 1 Co-author. Tel.: +1 520 621 6580. http://dx.doi.org/10.1016/j.ultras.2014.08.019 0041-624X/Ó 2014 Elsevier B.V. All rights reserved.
tests, this arrival time is occasionally used to get an estimate of the wave velocity since the distance between the transducers is known. The typical ultrasonic signal is composed of multiple modes and reflections with different frequencies. In testing concrete, since the reflections from defects inside the slab are rarely well separated in time, it becomes highly problematic to identify the arrival time from the desired targets. The detection of the arrival times from multiple scatterers inside the medium has been the subject of intense research over many decades. Different methods have been proposed to decompose the recorded ultrasonic waves and measure the associated arrival times. The following is an inventory of the main techniques that have become popular among the nondestructive testing practitioners: 1. Application of the traditional quadrature-demodulation scheme to suitably extract the envelope of the main echo and locate its onset [1]. 2. Correlation between input and output for high signal-tonoise ratio waveforms and L2-norm with low noise level [2]. 3. Application of the square-root unscented Kalman filter (SRUKF) to identify the shape parameters of an ultrasonic echo envelope [3]. 4. Application of Synthetic Aperture Focusing Technique (SAFT) where a constructive interference procedure enhances the lateral spatial resolution and signal-to-noise ratio. The arrival time is estimated from enhanced images of the interior of the medium [4].
16
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
71
10
(a)
4500
6
4000
4
3500
Frequency [ Hz ]
Amplitude
8
2 0 -2
3000 2500 2000
-4
1500
-6
1000 500
-8 -10
Model Space
5000
0 0
50
100
150
200
250
300
350
0
0.1
0.2
0.3
0.6
0.7
0.8
0.9
1
0.5
(b)
6
0.45 0.4
4
0.35
Frequency [Hz]
2
Amplitude
0.5
Time-Frequency Distribution
8
0 -2 -4
0.3 0.25 0.2 0.15 0.1
-6
0.05
-8
0
-10
0.4
Time [ seconds ]
Time (us)
0
50
100
150
200
250
300
350
0.1
0.2
0.3
400
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Time [seconds]
Time (us) Fig. 2. Model space and the associated Wigner–Ville transform. Fig. 1. (a) Well separated ultrasonic signal; (b) mixed ultrasonic signal.
5. Separation of non-linear and non-stationary ultrasonic signals by means of the empirical mode decomposition method to obtain the intrinsic mode functions used for signal reconstruction. The arrival time is then obtained from the reconstructed signal [5]. 6. Combination of both the improved self-interference driving technique and the optional optimization signal processing algorithms for measuring the arrival time [6]. In some ultrasonic signals, picking the first reflection (arrival time) from the desired target is an easy task since the direct wave and the main reflection are well separated in time (with little or no interference from other sources of ultrasonic scattering), and therefore can be discriminated visually. Fig. 1(a) shows a signal in which the direct wave and the main reflection are relatively well separated and the arrival time was correctly determined by visual inspection. Fig. 1(b) on the other hand, depicts a commonly encountered situation in practice where multiple waves are superimposed, turning the detection of the main reflection into a challenging problem. In these cases, signal decomposition has shown promise and is used in this study to separate the signal into different modes in order to extract the desired reflected wave pack. More specifically, the synchrosqueezed transform (SST), a new signal decomposition method, is shown to provide useful results.
Conventional time–frequency representations, such as the short-time Fourier transform (STFT), the Wigner–Ville transform (WVT), the Wavelet Transform (WT), and the Empirical Mode Decomposition (EMD) method, perform poorly when individual components of the recorded signal are not sufficiently separated in the time–frequency domain. The synchro-squeezed transform builds on the philosophy of the EMD method using a more robust theoretical foundation and employs a different approach for the construction of a signal’s constituent components. In fact, synchro-squeezing can be overlaid with many commonly used time– frequency methods [7–10]. SST is a combination of the wavelet transform and reassignment technique [11]. 2. Theoretical background The simplest form of a time–frequency representation is a spectrogram created using a short-time Fourier transform (STFT). This transform employs a short-time window that limits the evaluation of the Fourier transform to some specified region along the time index. It is well-known that STFT spectrograms suffer from the time–frequency resolution trade-off, that is, a shorter window provides better resolution in the time domain, and the larger bandwidth of the window’s spectrum affords a poorer resolution in the frequency domain. The Wigner–Ville transform (WVT) takes advantage of the powerful concept of match filtering, effectively employing the time-reversed copy of the signal as the window to carry out a short-time window Fourier transform. A
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
17
Fig. 3. (a) Concrete slab with a 15° bottom surface; (b) rectangular grid on concrete slab; (c) ultrasonic transmitter and receiver; (d) ultrasonic transducers while conducting the experiment; (e) typical recorded signal.
number of properties make the WVT valuable for decomposing complex signals. For example, the analysis is performed using only the transform’s ability to be entirely localized in the case of frequency modulated signals. The trade-off in using WVT is the manifestation of prominent interference (oscillations for the most part) between any two interacting components of the signal [12]. The effect is illustrated in Fig. 2 below. Smoothing somewhat enhances the readability of the transformed signal through the application of the pseudo-Wigner–Ville distribution which is essentially a windowed version of the WVT. Here, the cross-terms are attenuated at the cost of spreading signal energy along the frequency axis, adversely affecting the time–frequency localization of the signal content. To optimize the tradeoff between attenuating oscillations and focusing the signal energy in the time–frequency domain, the coefficients (t, w) of the WV distribution are reassigned to new points ~ which are the center of gravity of the signal energy ð~t; wÞ
distribution around (t, w). The reassignment method was first introduced by Kodera et al. [13] to overcome the spectrogram’s drawbacks such as biased estimation of the signals instantaneous frequency, group delay, and the Gabor–Heisenberg inequality which controls the locality in either the time or frequency domain, but never both [14]. The concept of reassignment can also be applied to improve the spectral resolution obtained from wavelet transforms. In this context, the process has been dubbed synchro-squeezing, and was first introduced to analyze speech signals [15]. Subsequent studies [8– 11,16] showed that synchro-squeezing performance is superior and could be considered as a viable alternative to the EMD method. In this paper, we show that the synchro-squeezing method could also be used to decompose ultrasonic signals into their constituent wave packets, providing the means to detect and determine the arrival time of the desired echoed pulses from inside the medium under investigation. This is especially valuable since
18
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
Original Signal 10
26
Time - Frequency - Amplitude
200
162
(a)
8
(a)
8
180 6
6 160
Frequency (KHz)
4
Amplitude
2 0 -2 -4
Time index for maximum amplitude 121.6
-6
4 2
140
0
120
-2
100 -4
80
-6
60
-8
-8 -10 0
50
100
150
200
250
300
350
400
50
Time (us)
150
200
250
300
350
400
Time (us)
Truncated and Padded signal
Synchrosqueezed Wavelet Transform
10
(b)
8
(b)
180
4.5
160
6
4
140
4 2
Frequency (KHz)
Amplitude
100
0 -2 -4
3.5
120
3
100
2.5 2
80
1.5
-6
1
-8 -10 0
60 50
100
150
200
250
300
350
Central Frequency = 0.8
0.5
400 0
Time (us)
50
100
150
200
250
300
350
400
Time (us) Fig. 4. (a) Original signal with lower and upper time interval and time index for maximum amplitude; (b) truncated and padded signal.
the shape of the mother wavelet is tuned to closely match the ultrasonic pulse waveform used in the inspection process. The following is a brief description of the synchro-squeezing transform (SST). If s(t) is the desired signal, its continuous wavelet transform (CWT), Ws(a, b) is given by:
1 W s ða; bÞ ¼ pffiffiffi a
Z
1
1
sðtÞw
tb dt a
ð1Þ
where w⁄ is the complex conjugate of the mother wavelet (a continuous function in both the time and frequency domains), and b is the time shift applied to the mother wavelet which is scaled by a. We can think of the CWT as the cross-correlation of the signal s(t) with pffiffiffi the scaled and time shifted wavelet wððt bÞ=aÞ= a. This cross-correlation provides a measure of similarity between the signal and the scaled shifted wavelet. The wavelet coefficients Ws(a, b), however, frequently spread out along the scale dimension a, producing a blurred image in the time-scale domain. The next step in the process involves the transformation of the points in the time-scale domain into the time–frequency domain for all Ws(a, b) – 0. Different expressions may be used to handle the transformation. We adopted the following relationship [16]:
1 d xs ða; bÞ ¼ arg½W s ða; bÞ 2p db
ð2Þ
Fig. 5. (a) TFA distribution; (b) SST distribution at central frequency of 0.8.
where arg () denotes the unwrapped phase of the complex coefficients and the 1/2p factor is introduced to convert from circular to normal frequencies. Eq. (2) is the expression for instantaneous frequency (derivative of the unwrapped phase of the continuous wavelet transform). When the coefficients xs(a, b) are determined, we select the desired frequencies xk and form the frequency interval [wi–Dw/2, wi + Dw/2] where Dw = wi–wi1. The synchrosqueezing transform Ts(wi, b) is computed only at the centers wi of the frequency range as:
T s ðwi ; bÞ ¼
1 Cw a
X
3=2
W s ðaj ; bÞ aj
Daj
ð3Þ
j:jwðaj ;bÞwi j6Dw=2
where Daj are the distances between two adjacent scales aj, and R1 ^ eÞ de, with wð ^ eÞ being the Fourier transform of the mother C w B 12 0 wð e wavelet. Eq. (3) indicates that synchro-squeezing takes place along the frequency/scale axis, producing a more accurate and detailed time–frequency representation. The synchro-squeezed transform is invertible and the original time domain signal can be reconstructed (with no loss of information) using [11]:
! X X X sðtÞ ¼ Re T s ðwi ; bÞ ¼ T s ðwi ; bÞ cos arg T s ðwi ; bÞ i i i
ð4Þ
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
19
form of the wavelet because, compared to other mother wavelets, it represents most accurately the actual ultrasonic pulse that probed the interior of the concrete material in this study:
Input Signal 50
t2 2 1 wðtÞ ¼ Bf 0 ei2pf 0 t e2ð2pf 0 Þ e 2
0
ð5Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffi where Bf 0 is a normalization factor equal to 2=pf 0 and f0 is the socalled central frequency, which determines the relationship between the time and the scale (or frequency) resolutions of the wavelet transform.
Amplitude
-50
-100
3. Proposed technique -150
-200
0
50
100
150
200
250
300
350
400
Time (us)
Time-Frequency-Amplitude 300
44.4 23.9
Frequency (KHz)
250
3.52 -16.9
200
-37.3 150
-57.8 -78.2
100
-98.6 -119
50
-139 0
-160 0
50
100
150
200
250
300
350
400
Time (us) Fig. 6. Input signal and its TFA distribution.
Synchrosqueezed Wavelet Transform 180KHz
180
Frequency (KHz)
160
140KHz
140
120KHz
120
100
90KHz 50
100
150
Time (us) Fig. 7. Selection of frequency bands.
The choice of the mother wavelet is based on many factors, and should reflect the type of features present in the signal. A commonly used function is the Morlet wavelet, consisting of a plane wave modulated by a Gaussian function. The Morlet wavelet is non-orthogonal, complex and offers a good tradeoff between detecting oscillations and peaks or discontinuities making it especially suited for analyzing ultrasonic signals. We used the following
We describe a method to unscramble closely spaced events commonly encountered in ultrasonic testing of materials. This condition occurs when the direct wave and the reflection from a scatterer superimpose and create an interference in the recorded signal because the events are not well separated in time, space, and frequency domains, or in other words, the length of the ultrasonic pulse is larger than the time of flight between the echoes. We specifically tackle the case where the detection and estimation of the arrival time of the first reflection is needed. Our approach uses both the synchro-squeezed transform and the measurement of instantaneous frequencies to decompose a signal into its components. The arrival time is then determined for the desired reflection. First, the arrival time associated with the highest peak is determined from the time domain signal. Let us call this tmax. Since both the input and output signals have a Gaussian shaped envelope, we know that the actual arrival time should be before tmax. Next, a certain range is selected to truncate the signal and make the decomposition procedure go faster. The lower limit is chosen at some point before the onset of the direct wave arrival, while the upper limit is a point after tmax where the change in the shape of the envelope of the signal is discernable, or tmax + tmax/3 in extremely mixed waveforms. The signal is then truncated and padded with zeroes to the same length as the original signal. To decompose the signal and measure the arrival time, the following steps are performed. Calculate the instantaneous frequencies and plot the time–frequency–amplitude (TFA) distribution. Note that the instantaneous frequency is the derivative of the unwrapped instantaneous phase, which is obtained by performing a Hilbert transform on the truncated s(t) signal. The amplitude values in the TFA distribution are taken from the original time domain signal and color coded for better readability. Next a synchrosqueezed transform is carried out on the truncated signal using an appropriate central frequency. The choice of the central frequency is crucial and will affect the outcome of the analysis. The rule we propose for the selection of the central frequency is based on visual interpretation of both the SST and time–frequency– amplitude distributions. The instantaneous frequency plot is compared with the SST and the optimum central frequency value is found from the closest match between the two distributions. The best match between the TFA and SST distributions is obtained by changing the value for the central frequency in the SST and comparing the resultant maximum amplitudes (values in the third dimension) in the TFA and SST plots. The time and frequency values corresponding to the maximum amplitudes should be the same or very close for the optimal central frequency. In almost all cases, the central frequency is found to be between 0.5 and 1. We continue the process by defining the frequency bands needed for signal decomposition. Since the first reflection normally has the strongest echo, we look for the frequency range in both the TFA and SST plots where the highest amplitudes (absolute value) are located.
20
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
Decomposed Signal : Fmin=50000 Fmax=90000
Decomposed Signal : Fmin=90000 Fmax=120000
0.2
2
A
B
1.5
0.15
1
0.1
0.5
Amplitude
Amplitude
0.05 0 -0.05
0 -0.5 -1
-0.1
-1.5
-0.15 -0.2
-2
0
50
100
150
200
250
300
350
-2.5
400
0
50
100
150
Time (us)
200
250
300
350
400
Time (us)
Decomposed Signal : Fmin=120000 Fmax=140000
Decomposed Signal : Fmin=140000 Fmax=180000
8
3
C
6
D 2
4 1
Amplitude
Amplitude
2 0 -2
107
0
-1
-4 -2
-6 -8
0
50
100
150
200
250
300
350
400
-3
0
50
100
Time (us)
150
200
250
300
350
400
Time (us)
Fig. 8. Decomposed signal and associated waveforms from defined frequency bands.
Time-Frequency-Amplitude 9.59
300 7.67
Sudden Change in Frequency Value
250
5.75 3.83
Frequency (KHz)
4 200
1.91
150
1
-0.00497
3
2
-1.92
100
-3.84 -5.76
50 -7.68
0 40
60
80
100
possible that the main frequency band is selected incorrectly and the resultant waveform does not correspond to the actual reflection. To mitigate those special circumstances, the two methodologies described below can be used to validate the results.
120
140
160
-9.6
Method 1: In the time–frequency–amplitude plot, abrupt changes in frequency are indicative of different wave types detected by the sensor. The last peak before the cluster of high amplitude values corresponds to the time index of the main reflected wave. Method 2: In the time–frequency–amplitude plot, select the frequency range where the maximum amplitude values are located (the cluster of red and blue dots). Use this frequency range to reconstruct a time-domain waveform from the original time-domain signal. In other words, amplitudes outside of this range are removed from the original signal s(t).
Time (us) Fig. 9. Sudden changes in the time–frequency–amplitude distribution.
Once the frequency bands are selected, we extract the waveform from each band by performing an inverse synchro-squeezed transform of the SST coefficients. The reconstructed time-domain signal with the highest energy is the waveform associated with the main reflection from the discontinuity. In some cases, it is
4. Experimental setup The proposed decomposition method was tested by analyzing signals recorded on three concrete slabs, each with a surface area of 1 1 m2 and a sub-horizontal inclined bottom surface of 5°, 10°, and 15° respectively. Fig. 3(a) shows the slab with a 15° inclined bottom surface. The concrete used in this work had a max-
21
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
Reconstructed Signal : Fmin=90000 Fmax=180000
Time-Frequency-Amplitude 10
300
8
106.3
6 4
200
2
Amplitude
Frequency (KHz)
250
150
100
0 -2 -4
50 -6
107
-8
0 40
60
80
100
120
140
160 -10
Time (us)
0
20
40
60
Decomposed Signal : Fmin=120000 Fmax=140000 8
100
120
140
160
Fig. 12. Signal in time domain after removing amplitudes value outside certain frequency range.
C
6
80
Time (us)
Reconstructed Signal : Fmin=90000 Fmax=180000
4
10
Amplitude
2
8
0
6
-2
4
107
2
Amplitude
-4 -6 -8
0
50
100
150
200
250
300
350
0 -2 -4
400
Time (us)
-6
Fig. 10. Comparison of frequency changes with reconstructed waveform for validating Method 1.
107
-8 -10
0
20
40
60
80
100
120
140
160
Time (us) Decomposed Signal : Fmin=120000 Fmax=140000 8
Time-Frequency-Amplitude plot 250
8
6
6
4
180KHz
200
2 150 0
90KHz
100
-2 -4 -6
50
2
Amplitude
Frequency (KHz)
4
C
0 -2
107
-4 -6
-8
-8 20
40
60
80
100
120
140
160
Time (us) Fig. 11. Selection of the frequency band in the TFA distribution for validating Method 2.
0
50
100
150
200
250
300
350
400
Time (us) Fig. 13. Validation of the main reflection wave packet arrival time using Method 2 for signal S15.
22
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
15
Time-Frequency-Amplitude
122 Original Signal
26
200 180
10
160
5
5
Frequency (KHz)
Amplitude
10
0
-5
Time index for maximum amplitude 94.68
-10
-15
140 0
120
-5
100 80
-10
60 -20 0
50
100
150
200
250
300
350
400
0
50
100
Time (us)
150
200
250
300
350
400
-15
Time (us)
Tuncated and Padded Signal
Synchrosqueezed Wavelet Transform
15
200
7
180
10
160
6
140 5
Frequency (KHz)
Amplitude
5
0
-5
120 4
100 3
80 2
-10 1
60
Central Frequency = 0.6
-15 0
-20
0
50
100
150
200
250
300
350
50
100
250
300
350
0
Fig. 15. (a) TFA distribution; (b) SST distribution at central frequency of 0.6.
Fig. 14. (a) Original signal with lower and upper time interval; (b) truncated and padded signal.
Synchrosqueezed Wavelet Transform 160
150KHz 140
Frequency K(Hz)
imum grain size of 20 mm and compressive strength of 32 MPa. It was made by a local company with the same properties as the normal concrete used in constructions to simulate a ‘‘real world’’ concrete structure. The incline surface at the bottom of the slabs was built using a layer of wood with a thickness of 5 mm. The slabs simulate wide inclined cracks in large concrete structures such as dams where water can seep through and expand the cracks resulting in significant damage. Data were collected from rectangular grids mapped in the middle of slabs with 10 mm spacing between the points as shown in Fig. 3(b). The size of the grid was different for each slab and varied from 21 to 27 cm in both x and y directions, corresponding to a depth of 100–300 mm, as shown in Fig. 3(a). The tests were conducted in pitch–catch mode using two ultrasonic transducers, one as a transmitter and the other as a receiver. Fig. 3(c) shows ultrasonic transducers used b this experiment and Fig. 3(d) demonstrates ultrasonic transducers conducting experiment. The transducers had a nominal frequency of 140 kHz and an actual bandwidth of 400 kHz. A pair of flat acrylic plates were attached to the surface of the transducers as a protection measure against scratch and damage. The thickness of the plates are 5 cm, and together, they introduce an additional delay of 26 ls in the signal.
200
Time (us)
400
Time (us)
150
130KHz 120
100KHz
100
40
60
80
100
120
140
Time (us) Fig. 16. Selection of frequency bands.
A two-channel data acquisition and function generator was used to input a 140 kHz square pulse and to capture the reflected waveforms. The data were collected at a high sampling frequency of 25 MHz using 12-bit resolution. For each waveform, 10,000 samples with a peak-to-peak amplitude of 12 V were collected. The signal captured from the test medium was passed through a signal
23
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
Decomposed Signal : Fmin=10000 Fmax=130000
Time-Frequency-Amplitude
15
250
10
69 200
Frequency (KHz)
Amplitude
5
0
-5
68
150
100
50
-10
0 -15
20 0
50
100
150
200
40
60
80
100
120
140
Time (us)
250
Time (us)
Decomposed Signal : Fmin=10000 Fmax=130000 15
Decomposed Signal : Fmin=130000 Fmax=150000 2 10 1.5 1
5
Amplitude
Amplitude
0.5 0
0
-5
-0.5
68
-1
-10
-1.5 -15 -2
0 0
50
100
150
200
250
300
350
50
100
Time (us) Fig. 17. Decomposed signal and associated waveforms from defined frequency bands.
150
200
250
Time (us)
400
Fig. 18. Comparison of frequency change with reconstructed waveforms for validating Methods 1.
Time-Frequency-Amplitude 200
5. Data analysis For the sake of brevity, the analysis results for two signals are presented in this section. The first signal is from the 15-degree slab (see Fig. 3 above) at a known depth of 16.69 cm. This gives an actual arrival time of 81.41 ls, using the velocity of 4100 m/s. As explained in the experimental setup section, there is a delay of 26 ls due to the thickness of the guides, and therefore the total time we should get after signal decomposition is 107.41 ls (81.41 + 26). We name this signal S15.The second signal, S05, is taken on the 5-degree slab at a known depth of 8.5 cm. The arrival time for this signal is 67.47 ls (41.47 + 26).The analysis steps are explained in detail for the S15 signal only. For signal S05, we only present the results.
180
10
160
5
Frequency (KHz)
conditioner with the appropriate low-pass filter before feeding it to the data acquisition system. Fig. 3(f) shows a screenshot of the software used in this work displaying a typical waveform obtained from the test slab and recorded by the data acquisition system. Calibration tests on five cylindrical concrete samples provided an average velocity of 4100 m/s. This value is used to calculate the depth of the inclined bottom surface.
140
130KHz 0
120
100KHz
100 80
-5
-10
60
-15 20
40
60
80
100
120
140
Time (us) Fig. 19. Selection of the frequency band from TFA distribution for validating Method 2.
The first step is the selection of tmax which is 121.6 ls for S15. We take as the lower limit, the time index of 26 ls and as upper limit, the time index 162 ls. Fig. 4(a) shows the signal and the
24
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25 15
Reconstructed Signal : Fmin=10000 Fmax=130000
10
10
5
Amplitude
Amplitude
5
0
-5
0
-5
68 -10
-10
-15
-15 40
60
80
100
120
93.93
0
50
100
150
200
250
300
350
400
250
300
350
400
Time (us)
140
Time (us) 5
Decomposed Signal : Fmin=10000 Fmax=130000
4
15
3 2
10
Amplitude
1
Amplitude
5
0
0 -1 -2 -3
-5
102.8
-4
68
-5
-10
0
50
100
150
200
Time (us) -15
0
50
100
150
200
250
Time (us) Fig. 20. Validation of the main reflection wave packet arrival time using Methods 2 for signal S05.
selected time interval. The signal is then truncated and padded to zeroes as shown in Fig. 4(b). Next, the instantaneous frequencies, and the synchro-squeezed transform (at different central frequencies) are computed, and the time–frequency–amplitude (TFA) and the SST distributions are plotted. The optimal value of 0.8 for the central frequency is found by comparing the maximum values in the TFA and SST distributions, and the best match is selected. Fig. 5(a) shows the TFA distribution and Fig. 5(b) presents the SST at the optimum central frequency for signal S15. Let is take a closer look at the SST plot and select the frequency band that encompass the largest SST coefficients and call it the main reflection. Also, the selected frequency band should include the main input signal frequencies (80–140 kHz in our case). Fig. 6 shows the input signal and its associated time–frequency– amplitude plot. Using the SST distribution shown in Fig. 7, we identify the main frequency range from 120 kHz to 140 kHz, and we define the following frequency bands for signal decomposition: 50–90 kHz, 90–120 kHz, 120–140 kHz, and 140–180 kHz. Once the waveforms are reconstructed in the time-domain for each of the defined frequency band, the waveform with the largest energy (sum of the amplitudes squared) is used to determine the arrival time. Fig. 8
Fig. 21. Checking accuracy of arrival time by subtracting tmax of two consecutive signals.
presents the reconstructed waveforms for each of the defined different frequency bands. From the previous analysis, we conclude that the time corresponding to the wave reaching the sensor 107 ls with a frequency bandwidth from 120 kHz to 140 kHz. The last step involved in determining the time of arrival is to verify that the reconstructed waveform carries the correct information as discussed in the proposed method section. 5.1. Method 1 We consult the time–frequency–amplitude distribution, shown in Fig. 9, and locate the frequency peak (peak #4) before the cluster of high amplitude values in the signal. The time index corresponding to the peak #4 equals to 106.3 ls which is very close to the actual arrival time of 107.41 ls. Fig. 10 shows the comparison of the sudden change in the time–frequency–plot with the decomposed signal (Fig. 8(c)) in order to confirm that the correct reflection wave packet is selected. 5.2. Method 2 Let is take a look at the time–frequency–amplitude distribution and define the frequency range encompassing the highest amplitude values. In the case of signal S15 illustrated in Fig. 11, the lower
25
P. Avanesians, M. Momayez / Ultrasonics 55 (2015) 15–25
and upper limit of the frequency range are 90 kHz and 180 kHz. Then, amplitude values outside this range are removed from the signal s(t), and the desired arrival time is selected from the reconstructed waveform as shown in Fig. 12. A comparison of the reconstructed waveform in the band 120– 140 kHz (Fig. 4C) and the reconstructed waveform shown in Fig. 13, validates the results of the analysis and provides the correct time of arrival. The second signal from the 5-degree slab has an arrival time 41.47 ls corresponding to a depth of 8.5 cm. The actual total time is 67.47 ls (includes the delay from the acrylic guides). As in the previous example, we start with the time–frequency–amplitude distribution plot and compare this with the SST distribution in order to find the optimal central frequency. This step is followed by the selection of frequency bands, extraction of the SST coefficients and reconstruction of the waveform for each frequency band. The waveform with the highest energy content which also contains the main input frequency range is used to determine the time of arrival. Validation is performed by comparing the change in the TFA distribution with the reconstructed signal from the main frequency band (Method 1), and by removing the samples with amplitudes outside the defined frequency range from the original signal s(t) (Method 2). The reconstructed waveforms are plotted and compared to obtain the arrival time. Figs. 14–20 illustrate the preceding steps for S05.
yields
Dta1 a2 ¼ jta2 t a1 j ¼ j48:63 39:73j ! Dt a1 a2 ¼ 8:9 ls yields
Dt1max 2max ¼ jt 2max t 1max j ¼ j102:8 93:93j ! Dt 1max 2max ¼ 8:87 ls Dta1 a2 ffi Dt1max 2max 7. Conclusions The synchro-squeezed transform is an effective tool for decomposing many types of noisy, time-varying frequency, and nonlinear signals. In this paper, we presented a new technique for determining the arrival time of ultrasonic signals used in nondestructive testing of geomaterials. The method uses the synchro-squeezed transform to decompose the signal into constituent waveforms reconstructed from suitably defined frequency bands. We have shown that the time of arrival of the primary reflection (from a discontinuity inside the material) can be accurately measured from signals where the events are poorly separated in time and frequency. This technique could also be used to detect the arrival time of well separated time, frequency and space domain signals. It was successfully applied in the processing of lower frequency elastic waves from earthquake events and seismic exploration. References
6. Discussion Over 500 signals were analyzed in this study and arrival times were successfully determined. The results were verified with known depths of the bottom discontinuity in the concrete slabs. In the process, we observed that the choice of the central frequency affects the outcome of the analysis to some extent, especially with respect to the creation of the synchro-squeezed transform plot. As mentioned above, at present, the selection of the central frequency relies on a visual inspection and comparison of both the SST and TFA distributions. The best match between the instantaneous frequency and SST defines the optimal value for the central frequency. Care should also be exercised in the selection of the main frequency interval and the corresponding frequency bands used for the extraction of the SST coefficients and the subsequent reconstruction of waveforms in the time-domain. The suggested approach where the size of the frequency bands is taken to be equal to the width of the main frequency interval (obtained from the SST distribution) is a good initial guess and provided accurate estimates of the time of arrival for all analyzed signals. To obtain the arrival time for the first reflection, the original signal s(t), the time-instantaneous frequency–amplitude and the decomposed waveforms using the synchro-squeezed transform are analyzed simultaneously. The accuracy of the measured values can further be investigated by comparing two closely-spaced signals. For example, if time indexes corresponding to the maximum amplitude for the first and second time-domain signals are t 1max and t 2max , and the actual arrival times are ta1 and ta2 , then difference between t 1max and t2max must be equal to the difference between t a1 and t a2 . From Fig. 21 below, we read off the following values:
t a1 ¼ 39:73 ls; ta2 ¼ 48:63 ls; t 1max ¼ 93:93 ls; t2max ¼ 102:8 ls We then form the differences for the maximum amplitude time indexes and the measured arrival times:
[1] L. Angrisani, R.D.L. Moriello, Estimating ultrasonic time-of-flight through quadrature demodulation, IEEE Trans. Instrum. Meas. 55 (2006) 1. [2] M.A.H. Qayyum, High accuracy time of flight measurement using digital signal processing techniques for subsea applications, J. Sig. Inf. Process. 2 (2011) 330– 335. [3] Z.J. Yao, Q.H. Meng, M. Zeng, Improvement in the accuracy of estimating the time of-flight in an ultrasonic ranging system using multiple square-root unscented Kalman filters, Am. Ins. Phys. Rev. Sci. Instrum. (2010) 81. [4] Z. Hosseini, M. Momayez, F. Hassani, D. Lévesque, Detection of inclined cracks inside concrete structures by ultrasonic SAFT, AIP Conf. Proc. 975 (2008) 1298– 1304. [5] T.L. Chen, P. Que, Q. Zhang, Q.K. Liu, Ultrasonic signal identification by empirical mode decomposition and Hilbert transform, Am. Inst. Phys. Rev. Sci. Instrum. (2006) 76. [6] X.F. Wang, Z.A. Tang, A novel method for digital ultrasonic time-of-flight measurement, Am. Inst. Phys. Rev. Sci. Instrum. (2010) 81. [7] Daubechies I. Ten lectures on wavelets. CBMS-NSF Regional Conference Series in Applied Mathematics, 1992. [8] C. Li, M. Liang, A generalized synchro squeezing transform for enhancing signal time–frequency representation, Signal Process. 92 (2012) 2264–2274. [9] C. Li, M. Liang, Time-frequency analysis for gearbox fault diagnosis using a generalized synchrosqueezing transform, Mech. Syst. Sig. Process. 26 (2012) 205–217. [10] G. Thakur, H.T. Wu, Synchrosqueezing-based recovery of instantaneous frequency from nonuniform samples, SIAM J. Math. Anal. 43 (5) (2011) 2078–2095. [11] I. Daubechies, J. Lu, H. Wu, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl. Comput. Harmon. Anal. 30 (2) (2011) 243–261. [12] G. Matz, F. Hlawatsch, Wigner distributions (nearly) everywhere: timefrequency analysis of signals, systems, random processes, signal spaces, and frames, Sig. Process. 83 (2003) 1355–1378. [13] K. Kodera, R. Gendrin, C. Villedary, Analysis of time-varying signals with small BT values, IEEE Trans Acoust. Speech Sig. Process. 26 (1978) 64–76. [14] D. Gabor, Theory of communication, Inst. Electron. Eng. 93 (11) (1946) 429– 457. [15] Daubechies I, Maes S. A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models. Wavelets Med. Biol., CRC Press 1996, pp. 527–546. [16] D. Iatsenko, A. Stefanovska, P.V.E. McClintock, Nonlinear mode decomposition: a noise-robust, adaptive, decomposition method based on the synchrosqueezed wavelet transform, Appl. Comput. Harmon. Anal. (2012).
