Introduction

Fig 1: Ultrasonic NDT system.

Modern ultrasonic and radar measurement systems are widely used in the field of non-destructive testing for a long time. The limitation of a currently available ultrasonic instruments hardly lies on the property of hardware but it may lie on the lack of sufficient signal processing techniques [1]. At present, the ultrasonic A-scan type instruments are most commonly used. It is believed that the received A-scan signal may carry a lot of information on material properties and defects, but information appears in various guises of noise, which is to be deciphered completely. A lot of research has been done on ultrasonic signal processing and still now is going on in search for more reliable and versatile signal processing techniques [2-6]. Generally, the flaw signals measured in ultrasonic NDT include the effects of the measurement system and are corrupted by different kind of noise. The highly complex interaction between the defect geometry and the back-scattered ultrasonic wave inside the test piece may not be assumed as a linear process. So, the signal processing techniques which require apriory knowledge of noise statistics, are subject to fail in many situations. Therefore the approach of signal processing should be involving the noisy signal itself in constructing the signal processing method.

Signal processing in ultrasonic NDT systems

Let us analyze a real time system in which a transmitter and a receiver are located at predetermined points. Several M-sequences (usually the Barker code), modulated by ultrasonic wave are used as the transmitted signal and receiver receives signals reflected from the target. Such a system is depicted in Fig. 1.

This system measures the thickness of moving object A. The reference signal x_ref consisting of a certain coded sequence is emitted by the ultrasonic transmitter at the time moment t₀. The propagating signal partially reflects from the front side of object A:

(1)

The other part of the emitted reference signal x_ref reflects from the rear side of object A:

(2)

where: k₁ and k₂ are the coefficients depending on a distance to the object, environment and object properties, Dt₁ and Dt₂ are the delay times directly proportional to the distance d and the thickness of object A:

(3)

(4)

Finally in the receiver we get the signal y:

(5)

The task of signal processing is to determine the time instances Dt₁ and Dt₂. Then values of control signals are calculated and transmitted to actuators. The signal processing time is restricted by properties of a technological process and the velocity of the object A. For determination of the time instances Dt₁ and Dt₂ usually it is used principle of obtaining the impulse response by a correlation process. Let’s consider x_ref to be the transmitted sequence, y(t) to be the received sequence and h(t) to be the impulse response of the composite system, which includes the test piece, the transceiver system and their associated electronics. The scattered ultrasound that is picked up by the receiver and the additive system noise n(t) constitute the received signal y(t). The signal of the output of the correlation filter can be represented by

(6)

If peaks corresponding to reflections from the targets were clearly identified in the cross correlation function (CCF), it would be easy to determine the time instances Dt₁ and Dt₂. In practice, however, it quite difficult to identify them because of suspicious peaks in the CCF due to a noise from the surrounding medium and it is essential to cancel out effects of noise. In order to reduce the effects of a noise during transmission and reception some measures have to be taken.

Band modification by moving average

Data processing is performed on both of the sampled received signals and the original M-sequence and the sampling frequency is such that there are j samples per unit pulse of the M-sequence. So, the expected peak on the CCF is supposed to consist of 2j samples. In order to minimize spurious peaks, widths of which are less than 2j samples, the sampled data are smoothed by performing the moving average of the data sequence:

(7)

(8)

This is equivalent to application of the Hanning spectral window. This method considerably minimizes the noise peaks of comparatively small width while keeping the expected peaks intact.

Averaging

The system emits the reference signal x_ref periodically

(9)

where N is the period of the reference signal. If the position of the object A during l periods changes a little, it is possible to average input signal l times:

(10)

The noise level is reduced times. This is effective, but time wasting method and is not used in the case of signal processing time restrictions in real time systems. This method does not allow eliminating peaks caused by surrounding environment.

Noise cancellation by subtraction

Some coherent peaks, additive to the expected peaks, appear on the CCF, which are confusing in regard to the clear distinction of a target. This is due to the surrounding structure or due to the effect of limitations of the measuring system. These clutter peaks appear irrespective of presence of any target. To perform the subtraction, first of all the data are collected from the test object without presence of any the target. Another data are taken with the presence of target. Coherent peaks are cancelled by taking the difference between the CCF of second data and that of the first data. This helps distinguishing the peaks corresponding to the reflections from the newly developed targets by removing the coherent noise of the system.

Inverse filtering

Passing the signal through an inverse filter can significantly reduce a random and clutter type noise. A major part of the long CCF is to be assumed as a noise except the portion corresponding to the direct signal and the signal reflected from the target. The inverse filtering operation [7-9] of a signal is described as:

(11)

where is the output of the filter and P is prediction or the model order. The inverse filter is designed calculating the coefficients {a_k} based on noisy data. Coefficients {a_k} are obtained by solving the equation [10]:

(12)

where R[i] is the autocorrelation function defined by:

(13)

This autocorrelation function is constructed with N samples of data from a suitable portion of the received signal, which presumably contains no expected peak. Such a filter is depicted in Fig. 2

Fig 2: Inverse filter.

This filter attempts to remove the predictable part of the signal and produce an output , which is completely unpredictable to the filter.

Wavelet transform based noise reduction

During the last time the wavelets have become a popular de-noising (or noise reduction) tool [11]. Donoho and Johnston [12] showed that this method has statistical optimality properties. Many algorithms define a criterion to divide wavelet transform coefficients into two groups. The first group contains the coefficients dominated by a noise, while other coefficients are rather clean. These algorithms eliminate all wavelet coefficients below a certain threshold, because these coefficients are dominated by a noise.

Let’s consider the following model of the received discrete noisy signal

(14)

or in a vector notation:

(15)

To reconstruct the original data, a wavelet representation is used. We use simple non-redundant orthogonal, discrete wavelet transforms. An orthogonal matrix W can be used to represent this operation. We consider the following transform:

(16)

These transforms localize the most important spatial and frequencies characteristics of a regular signal in a limited number of wavelet coefficients. On the other hand, it is easy to prove that an orthogonal transform of a stationary, white noise results in a stationary white noise. This means that the expected noise energy is the same in all coefficients. If this energy is not to large, the noise has a relatively small influence on the important large regular signal coefficients. These observations suggest replacing the small coefficients by zero, because they are dominated by noise and carry only a small amount of information.

The thresholding operations can be represented as

(17)

where

(18)

There are known two threshold methods – hard threshold and soft threshold (or shrinking function) [13-15].

In the case of the hard threshold the entries of the matrix D_d are

(19)

In the case of soft threshold the entries of the matrix D_d are

(20)

These threshold functions are shown in Fig. 3. A wavelet coefficient w between -d and d is set to zero, while others have the same value in the case of the hard threshold, or are shrunk in an absolute value in the case of the soft threshold.

Fig 3: Hard thresholding (a) and soft thresholding (b) functions.

A natural question arising from this procedure is how to chose the threshold. If y_d is the result of applying threshold procedure to the wavelet coefficients of signal y, and e_d = y_d-f is the noise of this result, then an often used criterion to measure the quality of this result is its signal to noise ratio (SNR(d)):

(21)

An optimal choice of d should maximize SNR(d). This is equivalent to minimizing the mean squared error R:

(22)

Because the wavelet transform is orthogonal, we can also compute R from the wavelet coefficients as:

(23)

w_d = We_d is the noise after operation in the wavelet domain.

However, because f is unknown, the function R_d is not computable and hence it cannot be used to find optimal d. For automatic spline smoothing it was suggested to adapt general cross validation (GCV). Our aim is to minimize error function based on an unknown exact signal. We therefore try to find a good compromise between a goodness of fit and smoothness. We assume that the original signal is regular to some extend, which means that the value f_i can be approximated by an linear GCV combination of its neighbors. So, by considering y_di a combination of f_i, not depending on f_i itself, we can eliminate the noise in this particular component. Since we replace it by a weighted average of its neighbors, the noise in these components is smoothed, and so we end up with a relatively clean, noise-independent value. Applied to the wavelet procedure this GCV should be a function of a threshold value using only known data and having approximately the same minimum as the residual function R(d).

We have a definition of general cross validation:

(24)

where .

Note that if iąj, then d_ij=0. For i=j we have

(25)

Thus, if Tr(D^') is the trace of D^',

(26)

The results of applying the threshold procedure on the reflected signal are depicted in Fig. 4. In this case only a fragment of the Barker code is used for formatting M-sequence.

a

b

c

Fig 4: The reflected signal: a - without applying de-noising procedure, b - after applying hard thresold procedure, c - after applying soft threshold procedure.

Wavelet transform based signal processing method for ultrasonic NDT system

In order to reduce computations, the reference signal x_ref and the reflected signal y_d (signal y after applying thresholding procedure) are transformed:

(27)

After applying this transform we get digital signals and with logical values "0" and "1". This transform is possible because the most important information is the pulse widths of M-sequence. Results of applying procedure (27) are shown in Fig. 5

Fig 5: The transformed reflected signal: a - the reflected signal without noise, b - the reflected noisy signal after applying hard threshold procedure, c - the reflected noisy signal after applying sof threshold procedure.

The soft threshold procedure allows achieving a better visual quality, than the hard threshold procedure as shown in Fig.4. But, when the noise level is high, applying the soft thresholding procedure more distorts pulse widths of the transformed signal as shown in Fig.5.

The results of these transforms depend on the threshold d₁ value, which optimal value varies in accordance to the noise level e. When the threshold is low, an additional pulses emerge in the transformed signal , as shown in Fig. 5b. When the threshold value is high, the pulse widths of the transformed signal are shrinking and pulses may be distorted as shown in Fig. 5c. The optimal threshold d₁ value was defined by computing maximum value of the correlation function:

(28)

at different threshold d₁ values.

The best results were achieved when d₁=1.4d, where d is the threshold value defined by computing GCV function. Note that GCV function computation in this case does not require any floating-point operation and may be computed by a hardware. These computations may be simplified in the case when a fixed M-sequence is used.

Such a procedure may be used for recovering a distance to the object from noisy data. It possesses three steps:

1. apply the interval adapted pyramidal algorithm of Cohen, Daubechies, Jawerth and Vial [11] to the measured data, obtaining empirical wavelet coefficients w_i
2. apply the soft threshold nonlinearity coordinatwise to the empirical wavelet coefficients with the specially chosen threshold d;
3. invert the pyramid filtering recovering ;
4. apply the transform (xx) to the reference signal x_ref and de-noised data y_d;
5. detect the argument of correlation function maximum value.

For a fast wavelet transform we need 2N2F flops, where F is the number of filter coefficients. For F=4 , we have 16N flops. To reconstruct the signal after operation with the optimal threshold d we need 16N flops.

Computation of GCV(d) can be performed completely in the wavelet domain. Because GCV(d) is an approximation itself it is not useful to compute its minimum very precisely. Moreover, in most cases this is not necessary to the curve of R(d) in the neighborhood of its minimum. A relative accuracy of 10^-3 is enough. Using a classic minimization procedure (such as Fibonacci) this requires approximately 15 function evaluations. The denominator N-Tr(D^')counts the number of coefficients that are set to zero. This does not require any floating-point operation. Computation of the nominator can be done with 2N floating point operations. So 15 function evaluations lead to some 30N floating-point operations.

Computation of the signal can be done with N floating point operations.

Computation of the correlation function does not require any floating-point operation.

So execution of the suggested signal processing algorithm leads to 63N operations. Execution of a classical signal processing algorithm leads to (L+2P)N operations, where N is the number of samples, L is the length of M-sequence, P is the model order of the inverse filter. The suggested algorithm requires less floating-point operations, when
(L+2P)>63.

Conclusions

Generally the flaw signals measured in ultrasonic NDT systems are spoiled by different kind of a noise. Therefore, the approach of signal processing should be involving the noisy signal itself in constructing the signal processing method. The noise in such systems is cancelled by band modification using moving average, signal averaging, inverse filtering and noise cancellation by subtraction. These methods are time consuming and due to signal processing time restrictions not always may be used in real time systems. During the last time the wavelets have become a popular de-noising (or noise reduction) tool and this method has statistical optimality properties. New data processing method based on the wavelet transform for real time systems is suggested. It is shown that the hard threshold algorithm is preferred to the soft threshold in such systems. Execution of this method leads to less amount of floating point operations than classical signal processing methods.

References

1. Sinclair A. An analysis of ultrasonic frequency response for flaw detection: a technique review. Materials evaluation. 1989. Vol. 43. P. 870-883.
2. Fomitchev M. I. et. al. Ultrasonic pulse shaping with optimal lag filters. Int. J. Imaging syst. technol. (USA). 1999. Vol. 10. P. 397-403.
3. Sallard J. et. al. Use of a priori information for the deconvolution of ultrasonic signals. Review of progress in quantitative nondestructive evaluation. 1998. P. 735-742.
4. Crawford D. C. et. al. Compensation for the signal processing characteristics of ultrasound B-mode scanners in adaptive speckle reduction. Ultrasound Med. Biol. (UK). 1993. Vol. 19. P. 469-485.
5. Andrade Lima E. et. al. Image processing techniques to remove depth bias effects in magnetic source images of deep cracks. Review of progress in quantitative nondestructive evaluation. Plenum press. New York. 1997. Vol.1. P. 797-803.
6. Jianzhong C. et. al. Noise analysis of digital ultrasonic system and elimination of pulse noise. Int. J. Pressure vessels piping. 1998. Vol.75. P. 887-890.
7. Bengt M. et. al. Weighted least squares pulse shaping filters with application to ultrasonic signals. IEEE Trans. UFFC. 1989. Vol. 36. P. 109-113.
8. Venkantraman S. et. al. Combining pulse compression and adaptive drive signal design to inverse filter the transducer system response and improve resolution in medical ultrasound. Med. Biomed. eng & comp. 1996. P/318-320.
9. Izguierdo M. A. G. et. al. Multipattern adaptive inverse filter for real time deconvolution of ultrasonic signals in scattering media. Sensors and actuators. 1999. Vol. 76. P. 26-31.
10. Lim J. S. et al. Advanced topics in signal processing. Prentice Hall. Tokyo. 1988. P. 1-55.
11. Daubechies.Ten Lectures on Wavelets. CBMS-NSF Regional Conf. Series in Appl. Math., Vol 61. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1992.
12. Donoho D. L. and I. M. Johnstone. Adapting to unknown smoothness via wavelet shrinkage. Journal of American Statistical Association, 1995, 1200-1224.
13. Mallat S. A theory for multiresolution signal decomposition: The wavelet representation.IEEE Transactions on Pattern Analysis and Machine Intelligence, 11, 1998, 674-693.
14. Donoho D. L.and I. M. Johnstone.Minimax estimation via Wavelet Shrinkage. The Annals of Statistics, 26, 1998, 879-921.
15. Neumann M. H. and Von Sachs R.Wavelet thresholding: beyond the Gaussian I.I.D. situation. In A. Antoniadis and G. Openheim Wavelets and Statistics, New York: Springer-Verlag 1995, 302-329.

E. Kazanavicius, A. Mikuckas, I. Mikuckiene, V. Kazanavicius

Iškraipytu signalu apdorojimas realaus laiko DSP sistemose
Reziume
Ultragarsinese sistemose daug demesio skiriama triukšmo sukeltu efektu pašalinimui. Daugelis triukšmo sumažinimo metodu nera pakankamai efektyvus, todel dažnai vienoje sistemoje naudojami keli metodai iš karto. Del to išauga skaiciavimu apimtis, o kartu ir signalu apdorojimo trukme. Šiuo metu placiai pradeti taikyti vilneliu transformacija pagristi triukšmu sumažinimo metodai, kurie yra statistiškai optimalus. Pasiulytas vilneliu transformacija pagristas signalo apdorojimo ultragarsinese nedestruktyvaus testavimo sistemose metodas. Parodyta, kad šiam signalo apdorojimo metodui igyvendinti reikia mažiau skaiciavimu ir pagereja labai svarbus realaus laiko sistemu parametras – signalo apdorojimo trukme

Pateikta spaudai 2003 02 27

© NDT.net - info@ndt.net

|Top|