Analysis of split-spectrum algorithms in an automatic detection framework

Signal Processing 92 (2012) 2293–2307

Contents lists available at SciVerse ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Analysis of split-spectrum algorithms in an automatic detection framework A. Rodrı´guez a,n, A. Salazar b,1, L. Vergara b,1 a b

´ndez de Elche, 03203, Spain Departamento de Ingenierı´a de Comunicaciones, Universidad Miguel Herna Instituto de Telecomunicaciones y Aplicaciones Multimedia, Universidad Polite´cnica de Valencia, 46022 Valencia, Spain

a r t i c l e i n f o

abstract

Article history: Received 7 November 2011 Received in revised form 23 January 2012 Accepted 6 March 2012 Available online 15 March 2012

In this paper we study the problem of automatic detection of ultrasonic echo pulses in a grain noise background considering split-spectrum (SS) algorithms as sub-optimum solutions. First, SS algorithms are reformulated following an algebraic approach which is more appropriate from the perspective of automatic detection. Then, recombination methods will be modified according to the previous reformulation. We will consider some of the popular methods based in the phase observation (Polarity Thresholding and Scaled Polarity Thresholding) and in the order statistics (Minimization, Normalized Minimization and Frequency Multiplication). Different experiments with simulated and real data will support our theoretical analysis, and will show the advantages of the Frequency Multiplication method. Derivation of the formulas of the probability of detection and the probability of false alarm in every detector are included in the paper. & 2012 Elsevier B.V. All rights reserved.

Keywords: Non-destructive testing Automatic detection Split-spectrum

1. Introduction Split-spectrum (SS) algorithms are recognized as suboptimal but practical alternatives to the problem of detecting the presence of an ultrasonic pulse in a background grain noise [1–7]. They are based on the general assumption that the superposition of the many echoes coming from the scattering centers of the material under analysis is frequency sensitive, i.e., the grain noise response is expected to be very different at different frequency channels. On the contrary, the echo from a possible defect is expected to be similar in different bands. Hence, most of the approaches rely on implementing a filter bank followed by a (usually non-linear) combination of the different filter outputs [8–10]. The signal-to-noise ratio gain (SNRG), which is the quotient between the signal-to-noise ratio (SNR) at the output and the SNR at the input of the SS processor, is usually considered as the performance measure [11,12].

In spite of the many works existing on SS algorithms [8,13,14], there is still a lack of analysis which can help in establishing theoretical comparisons among the expected performances of the proposed SS variations, as well as giving criteria for the selection of the appropriate parameters in a specific implementation. The main purpose of this work is to contribute to a better understanding of the expected comparative performance of five SS algorithms, selected as appropriately representative of all the existing ones. Moreover the terms of the comparison will afford insights into the best choice of the involved parameters. To facilitate this task, the work has been developed inside the detection theory framework, which basically implies:

 Filter bank implementation has been changed by 

n

Corresponding author. Tel.: þ34 96 665 8825; fax: þ 34 96 665 8903. E-mail address: [email protected] (A. Rodrı´guez). 1 Tel.: þ34 96 3877308; fax: þ 3496 3877919.

0165-1684/$ - see front matter & 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2012.03.005



discrete Fourier transform (DFT) of the moving sample segment under analysis. The non-linear combination of the filter (now DFT) outputs is formulated in terms of computing an appropriate statistics and comparing it with a threshold. Performance comparison is made by the Receiver Operating Characteristics (ROC). This implies that the

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

2294



most important task to face is the computation of the probability of false alarm (PFA) and the probability of detection (PD) for every analyzed SS algorithm. The duration of the observation vector is assumed to be approximately equal to the ultrasonic pulse duration. That is done to improve the SNR, to avoid the problems of SS algorithms when dealing with multiple echoes and to improve the echo localization.

In Section 2 we will introduce a new algebraic formulation from the perspective of detection for the SS algorithm and for the recombination methods. In Section 3, the statistics (PD and PFA) for all recombination methods will be computed, which will be analyzed in Section 4. Then, theoretical results will be compared with those obtained from experiments with simulated and real data in Sections 5 and 6, respectively. Section 7 summarized the conclusions of this work. The derivation of the PFA and the PD for the different recombination methods is included in Appendices A–E.

2.2. Detection framework We will start our approach to the SS algorithms from the detection perspective replacing the usual filter bank by the DFT of the segment under analysis. This algebraic implementation is equivalent to the SS algorithm filtering process if the impulse response of filter k matches the DFT vector tuned to the normalized frequency k/N [20], with N the size of the DFT. This is a simple form of implementing the frequency channels with the advantage of the matrix formulation, more appropriate for a detection approach. Let us call r ¼[r(ns),y, r(ns þN  1)]T the vector of samples corresponding to the segment where the possible presence of an echo is to be decided. First of all we compute the DFT of r defining ei as the vector of the DFT at frequency i as: i 1 h ð4Þ ei ¼ pffiffiffiffi 1ejð2p=NÞi ejð2p=NÞi2    ejð2p=NÞiðN1Þ N For 0 rirN 1, so we can write the DFT of r as:

2. SS algorithm from the perspective of detection

v ¼ DFTðrÞ-vi ¼ rH  ei

2.1. Detection problem statement

where vi is the sample of the DFT at frequency i. Now we select an appropriate band to work with, thus obtaining a new vector z defined as:

The detection problem is set when we wish to detect the presence of a possible ultrasonic echo pulse p(n) in an interval of the sampled ultrasonic signal r(n). Thus, we have two possible hypothesis: H1 (echo from a defect) or H0 (grain noise), defined as: H0 : rðnÞ ¼ w0 ðnÞ H1 : rðnÞ ¼ pðnÞ þw1 ðnÞ

n ¼ ns ,:::,ns þ N1

ð1Þ

where ns and ns þ N 1 are respectively the starting and the final sample number in the analyzed interval (i.e., N is the segment length), and wi (n) are the grain noise samples under hypothesis i. We will assume that there are no multiple defects inside the interval, what we achieve selecting N approximately equal to the transmitted pulse length. This has the additional benefit of improving the input SNR, as using an observation interval much longer than the pulse duration implies that more amount of noise is present at the detector input [15]. Moreover using shorter observation segments improves the location of the echo and so the defect. This constraint on the selection of N has some significant implications in the practical application of SS algorithms as we will explain in Section 3.1. Determining the presence of pulse p(n) implies some processing f{  } on the segment: zðns Þ ¼ f frðnÞg

ð2Þ

and then a comparison with a detection threshold t so that:   ( if zðns Þ 4 t decides H1   ð3Þ if zðns Þ o t decides H0 Unfortunately, if we do not know the statistical model under H0 and/or H1, it is impossible to establish an optimal detector that solves the problem, but we can use the SS algorithm as a sub-optimal solution since it does not require any specific model.

z ¼ ½z0 ,. . .,zL1 T ¼ ½viL ,. . .,viU T

ð5Þ

0 r iL r ir iU rN1

ð6Þ

where iL and iU are respectively the index of the lower and upper bands of the equivalent SS filter bank and L¼iU  iL is the number of analyzed bands. The different SS algorithms make different processing on z to compute a number where possible presence of an echo from a defect should be enhanced with respect to the presence of grain noise only. As indicated, in the detection framework this can be formulated by considering two hypotheses: H1 (presence of defect echo) and H0 (only grain noise present), deciding between both hypotheses by comparing that number with an appropriate threshold. It is recognized that SS algorithms do not lead to optimal detectors (see for example [1,2]), but they can deal with very general signal and noise models. This is because the only implicit assumption in all the SS variations is that presence of defect echo bias all the different channel outputs in a rather uniform manner, meanwhile grain noise (due to the mentioned tuning frequency sensitivity) affects the outputs in a more random manner. Let us formulate the five selected SS algorithms inside the detection framework. 2.3. Minimization (MIN) In this case the minimum [14,6,16] of the magnitude square of the elements of z is to be compared with a threshold t so that: ( 2 if min9zi 9 4 t decides H1 ð7Þ otherwise decides H0 This aims to exploit the fact that the presence of a defect echo must contribute significantly to all the elements of vector z so that an increasing of the minimum

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

magnitude is expected in comparison with the case of only grain noise presence. 2.4. Normalized Minimization (NMIN) This method is similar to MIN, but now the statistics are normalized [17] by the maximum of the square magnitudes of the elements of z: 8 2 < if minjzi j 2 4t decides H1 maxjzi j ð8Þ : otherwise decides H0 This normalization tries to exploit the fact that it is the variability among the different channels, rather than the magnitude level, what should actually indicate the presence or not of the defect. 2.5. Polarity Thresholding (PT) In this case [16,17] the sign of the real part of the elements of z is evaluated so that: 8 > < if Re½zi  4 0 8i decides H1 if Re½zi  o 0 8i decides H1 ð9Þ > : otherwise decides H 0

PT assumes that the presence of the echo from a defect bias the output in all channels in a uniform manner (always positive or always negative). 2.6. Scaled Polarity Thresholding (SPT) This SS algorithm is a combination of MIN and PT [18]. The minimum of the square magnitudes is weighted by the difference between the number of outputs having positive (L þ ) real part and the number of outputs having negative real part (L  ): (  2 if 9L þ L 9 minzi  4 t decides H1 ð10Þ otherwise decides H0 2.7. Frequency Multiplication (FM) In this variation of SS algorithms [6,16,19] the square magnitudes of the outputs are multiplied and compared with a threshold. This variation of SS algorithms is not very usual, but we consider interesting to include it as it is substantially different from the others. 8 LY 1   > > < if zi 2 4t decides H1 ð11Þ i ¼ 0 > > : otherwise decides H

2295

consider that the real and imaginary parts of the elements of vector z are zero-mean unit-variance independent and identically distributed (i  i  d  ) Gaussian random variables. Identical, independent and variance normalized distributions is a reasonable assumption if we assume that a prewhitening transformation [20] or a (more usual in SS context) spectrum equalization is implemented. Gaussianity is also a reasonable assumption even when grain noise could be non-Gaussian as vector z is the result of a linear transformation on vector r (see [2]) and the central limit theorem may be claimed. Under H1 we will rely on the main, though rather loose, assumptions implicit in all SS algorithms: Presence of the echo from a defect uniformly bias all the channel outputs, hence we will simply consider that under H1 a constant real value s is to be added to every element of vector z, since these DFT samples correspond to the central part of the spectrum of the transmitted pulse, which will be approximately flat, spectral regularization is not needed. Moreover an interesting aspect to be included in the model is that, under H1, the grain noise level should reduce, as the presence of a defect will decrease the number of scatters in the vicinity of it. In consequence the assumed model is: H 0 : z ¼ wg

wg : Nð0,2IÞ

H1 : z ¼ awg þ s1

ð12Þ

where wg is the grain noise component of the DFT of r, I is the identity matrix, 0 is a vector of 0’s, 1 is a vector of 1’s, 0 r a r1 is a constant modeling the possible grain noise reduction in presence of a defect, due to the small number of scatterers expected in the vicinity of it because of the volume occupied by the own defect, and s is a constant that models the presence of a defect and is proportional to its amplitude, which will be referred as defect amplitude within the text. 3.2. Computing the statistics The derivation corresponding to the PFA and the PD for the different methods is included in Appendices A–E for MIN, NMIN, PT, SPT and FM methods, respectively. The final relevant formulas are showed all together in Table 1. These equations will be used to compare the different methods in the next sections. In Table 1, I0(  ) is the modified Bessel function of zero order, Q is the Marqum Q-function, erfc is the complementary error function, Ei is the exponential integral function and d, p and q are variables defined in the corresponding Appendices.

0

4. Theoretical analysis 3. Computing PFA and PD 4.1. Previous considerations 3.1. Assumed model To compute the PFA and the PD we have to assume some statistical model of vector z under both H0 and H1. To make the analysis viable we have to assume simple but reasonably realistic models. Hence, under H0 we will

The Receiver Operating Curve (ROC) will be used as the figure of merit to make the comparison between the detectors. The different ROC will be calculated using the analytical expressions of the PFA and PD presented in the previous section.

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

2296

Table 1 Statistics of recombination methods. PFA MIN

NMIN

" LU

L1 X

ð1Þ

k



L1



1 1  t þ 1 þ k 2þ k

k

PT

2U "

L 1 X L

2 FM

Q L1

e

k¼0

SPT

PD

Lt 2

d odd

"

a

  L   L 1 s 1 s erfc  pffiffiffi þ 1 erfc  pffiffiffi 2 2 a 2 a 2 #L

!#

L þ

rffiffiffiffiffiffiffiffiffi! s2 y ½Q 3 Q 2@y 2

 Z 1 L1  L1 LX 2 2 Q k2 eððs =a Þ þ yÞ=2 I0 ð1Þk 2k¼0 k 0

 L 1 2

!

L ðLdÞ=2

 #L

ðL þ dÞ=2

"

Uet=2d

L X d odd

! pðLdÞ=2 qLððLdÞ=2Þ þ

L

!

ðLþ dÞ=2

#L

pðL þ dÞ=2 qLððL þ dÞ=2Þ Q 4

1    C Bln t L ln s2 Ei s2 C B ð2a2 ÞL 1 2a2 2a2 C sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi erfcB   C B 2 2 A @ 9zi 9 Lvar ln 2a2 =H1 0

  1 lnt0:116L erfc pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2U1:484L

  1 2 erfc  pffiffiffi q ¼ ð1pÞ 2 a 2  s pffiffiffiffiffi  s pffiffiffi , y , ty Q4 ¼ Q Q2 ¼ Q Q3 ¼ Q

  d ¼ L þ L   Q1 ¼ Q

L ðLdÞ=2

s

,

pffiffi t

a a

p¼

In all cases the procedure will be the same, we will establish a vector of values for the PFA between 0 and 1 and then calculate the threshold for each of these values. Once obtained the vector of thresholds, we can calculate the associated PD vector to finally represent the ROC. The only exception to this procedure will be the PT method, because (see Table 1) it is not involved in any threshold calculation, as both the PFA and the PD only depend on the number of bands used for the analysis. The above procedure will be repeated for different number of bands L and for different values of the defect amplitude s, restricted to L¼{2, 5, 10, 15, 20} and s¼ {1, 1.5, 2, 2.5, 3}, respectively, as these values are representative enough of all cases and prevent graphs to be unclear. Note that parameter a is set to 1 in all cases because its effect can be included in s.

a

a

s

a

,

pffiffiffiffiffiffiffi! t=d

a

 The FM method is very sensitive to both the SNR and the number of bands.

 For very low SNR, PT and SPT methods are better than 

FM, but when the SNR increases slightly the FM method is clearly superior to the rest. NMIN and MIN methods are not very sensitive to the number of bands and slightly sensitive to the SNR. In any case, these methods provide the worst results.

In summary, results evidence the superiority of the FM method in terms of detection capability and efficiency, as it achieves the best PFA/PD ratios with the lowest number of bands, as can be seen in Figs. 4 and 5 where we compare the different methods. 5. Analysis of simulation 5.1. Simulation models

4.2. Comparing the detectors Figs. 1 and 2 show ROC’s obtained for MIN, NMIN, SPT and FM as a function of the number of bands and the defect amplitude, respectively. Fig. 3 shows ROC’s for PT method, where as expected we get a series of discrete values for the PFA and PD. After the analysis of these figures and assuming the direct relationship between the amplitude of the defect and the input SNR, we can say that:

 PT and SPT methods are very sensitive to SNR and moderately sensitive to the number of bands, with SPT slightly better for both parameters.

In order to generalize the study, different stochastic models of grain noise will be used to generate experimental data. We have broadly classified the materials as low or high dispersive. Therefore, we are going to consider two propagation models of the ultrasonic pulse. The first one will be called low dispersive model (LDM) and the second one high dispersive model (HDM). In LDM it is assumed that the ultrasonic pulse remains essentially the same (except for a possible constant attenuation) when it propagates inside the material, thus a linear time invariant system model is possible, where the reflectivity of the material is to be convolved with the sent pulse. From a stochastic perspective, LDM implies a stationary model. The basic equation to simulate LDM

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

1

1

L=20

0.9

0.9

0.8

0.8

0.7

0.7

L=2

L=20

0.6

PD

PD

0.6 0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

2297

L=2

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

PFA 1

1

0.9

0.9

0.8

0.8 0.7

L=2

PD

PD

1

0.6

0.8

1

L=2

0.6

0.6 0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0.8

L=20

L=20

0.7

0.6

PFA

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

PFA

PFA

Fig. 1. ROC’s as a function of the number of bands L for defect amplitude s ¼2 (a) MIN (b) NMIN (c) SPT (d) FM.

is given by [9,21,22]: ( ) Kd Ks X X rðtÞ ¼ xðtÞn rsk Udðttsk Þ þ rdk Udðttdk Þ k¼1

ð13Þ

k¼1

with Ks and Kd the number of scatterers and defects, respectively, rsk and rdk the scatterers and defects reflection coefficients, respectively, tsk ¼2zsk/cl and tdk ¼2zdk/cl the delays corresponding to the scatterers and defects respectively (zsk and zdk are the scatterers and defects location and cl the propagation velocity of the ultrasonic wavefront in the material. Finally, x(t) is the transducer impulse response. In flaw detection approach, the reflection coefficients of the scatterers are modeled as a white random process, assuming constant the gap between them [23]. In HDM it is assumed that the absorption, scattering and reflection effects of the material microstructure are frequency dependent, hence the pulse suffer significant deformations as it propagates inside the material. Thus a

non-linear variant system model is required and nonstationarity appears. The basic equation to simulate HDM is given, in the frequency domain, by [1,9,24]: RðoÞ ¼ XðoÞU ( Ks X o2

b

k¼1

zsk

4

UeaR 2zsk o Ue2jozsk =cl þ

Kd X

) 4

rdk UeaR 2zdk o Ue2jozdk =cl

k¼1

ð14Þ 2 ¼QV/cl

is a characteristic constant of the material where b which depends on the characteristic geometry of the reflectors Q, the analyzed volume V and the propagation velocity of the ultrasonic wavefront in the material [25]. On the other hand aR is the fraction of the dispersion coefficient in the Rayleigh zone which can be also considered characteristic of the material, and rdk are the reflection coefficients of the defects. Finally, X(o) is the Fourier transform of the transducer impulse response in the frequency domain.

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

2298

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

PD

PD

s=3 1

0.5 0.4

0.3

0.2

0.2

0.1

0.1

0

0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

0.6

0.8

1

PFA 1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

PD

s=1

0.5

0.4

0.3

0.3

0.2

0.2

0.1

0.1 0.2

0.4

0.6

0.8

PFA

1

s=1

0.5

0.4

0

0

s=3

1

0

s=1

PFA

s=3

PD

0.5 0.4

s=1

0.3

s=3

0

0

0.2

0.4

PFA

Fig. 2. ROC’s as a function of the amplitude s for L ¼5 bands (a) MIN (b) NMIN (c) SPT (d) FM.

Regarding the transducer, models will be based both on deterministic Gaussian envelope or decreasing exponential signals, depending on the transducer to be considered. In the first case, the response of the transducer can be modeled by using a band-pass signal with Gaussian envelope [26], and in the second one by using a growing potential term combined with a decreasing exponential term, modulated to the desired central frequency [27]. In both cases, the desired waveforms can be achieved by modifying the appropriate parameters. Beamforming techniques [26] will be used to have better control on the transducer frequency response. Additive white Gaussian noise will be considered to simulate incoherent noise present in real data. To validate the results obtained for LDM, we simulate the behavior of a duraluminum specimen analyzed with a

2 MHz transducer, which will also be studied in the real data experiment section. At these frequencies, the aluminum alloy used behaves as a low dispersive material. The specimen has a hole in one of its faces acting as a defect. To simulate this environment, the transducer model used was the Gaussian envelope model centered at 2 MHz and with 50% of bandwidth. Scans were generated with a single delta at defect location with variable reflection coefficient, all immersed in a dispersive environment modeled by the stationary model. Fig. 6(c) shows an example of A-Scan obtained with this model which can be compared with a real A-Scan obtained from the duraluminum specimen shown in Fig. 6(a). For the HDM, we simulate the behavior of a block of cement with a defect in the middle of the material analyzed with a 5 MHz transducer. At these frequencies,

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

should be noted that the sampling frequency does not influence the results if N (defect size in number of samples) is chosen properly and calculations are consistent in every stage of the process for each case. Note that Gaussian noise hypothesis is only necessary for the Real and Imaginary samples of the DFT, as commented in Section 3.1. More evidences of the validity of these models can be found in [21,23,24,28].

the selected cement specimens (cement 32.5 30% humidity) are very dispersive. To simulate this environment, the transducer model used was a growing potential term combined with a decreasing exponential centered at 5 MHz and with 50% of bandwidth. Scans were generated with a single defect with variable reflection coefficient, all immersed in a dispersive environment modeled by the non-stationary model with attenuation parameter aR ¼4.86  10  30 and constant b ¼5  10  12 and with 5000 scatterers. Fig. 6(d) shows an example of A-Scan obtained with this model which can be compared with a real A-Scan obtained from the cement specimen shown in Fig. 6(b). The difference in the sampling frequency in the two experiments is not based on any specific reason other than the availability of the records in our database, but may serve to illustrate the fact that the validity of the results is independent of the oversampling factor. It

s=3

L=5 L=4 L=3

5.2. Previous considerations Before starting the experiments, it would be necessary to make certain comments on the limitations that the new algebraic implementation of the SS algorithm has as a method of processing. In Section 3.1, we saw that the set of spectral samples zi on which we apply the detectors must be independent, which will be true only when the size of the DFT equals N. Besides, as N is defined by the duration of the pulse, the number of independent spectral samples that fall within the bandwidth of the pulse cannot be chosen arbitrarily. For this reason, we will limit the experimental study to the effect of the defect amplitude s on the detection, fixing the number of bands at a value high enough to take advantage of the SNR gain but not so high as to compromise the independence of the spectral samples. Thus, we have restricted the number of bands from 2 to 10, obtaining in all cases very satisfactory results We only show results for 5 bands as this value is representative enough of all cases and prevent graphs to be unclear. For the PT method, as it only depicts discrete values, we can include more values without affect its interpretation.

L=2L=1

1 0.8

0.6

s=1

0.4

2299

0.2

5.3. Comparing the detectors 0

0

0.2

0.4

0.6

0.8

As previously discussed, two signal models will be used – stationary and non-stationary – to generate the records in which the amplitude of the defect will be modified with the same values as in the theoretical analysis. Then, PFA and PD will be calculated modifying

1

PFA Fig. 3. ROC’s of method PT as a function of the number of bands and defect amplitude.

1

0.8

0.8

0.6

0.6

PD

PD

1

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

PFA

0.6 PFA

FM

SPT

MIN

NMIN

Fig. 4. ROC’s of FM, SPT, MIN and NMIN (a) s ¼1 L¼ 5 (b) s ¼ 2 L¼ 5.

0.8

1

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

2300

L=5 L=4 L=3

L=2

L=5 L=4 L=3 1

0.8

0.8

0.6

0.6

PD

PD

1

L=2

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

PFA

0.6

0.8

1

PFA FM

PT

SPT

MIN

NORM

Fig. 5. ROC’s of FM, SPT, MIN, NMIN and PT as a function of L with the PFA calculated for the PT method. (a) s ¼1 (b) s¼ 2.

Amplitude

Amplitude

flaw echo

2

4

6

8

0

10

80

time (s)

160

240

time (s)

Amplitude

Amplitude

flaw echo

2

4

6

8

10

time (s)

0

80

160

240

time (s)

Fig. 6. ROC’s Example of real and simulated A-Scans. (a) Real duraluminum A-Scan (b) Real cement A-Scan (c) Simulated A-Scan with the LDM (d) Simulated A-Scan with the HDM.

the detection threshold within the limits calculated theoretically and counting the number of false alarms and detections in each case. This process will be repeated in series according to Monte Carlo method with 1000 iterations per series. Fig. 7 shows the obtained ROC for the different methods (except for the PT method that should be studied separately) using 5 bands and for amplitudes s¼1, 2, 3. As can be seen, in all cases results show a large coincidence between theoretical and experimental values, regardless of the model of noise used. Regarding the PT method, in Fig. 8 we can see the ROC obtained for 2 and 5 bands and amplitudes s ¼1 and 3 for the theoretical analysis, the stationary model and the

non-stationary model. We only show these results to avoid confusion in the graph due to a high density of points, but the trends are the same for intermediate values of the variables. Also in this case simulation results match expectations, especially regarding the values of the PD, being the PFA slightly better than expected, and repeating the theoretical pattern according to which for a given number of bands, the PFA is the same regardless the defect amplitude. In any case, it is clear the consistency of the statistical models we have developed for the study of the SS algorithm from the detection perspective, showing again the advantages of using the FM method compared to conventional options used to date.

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

1 0.8

s=2

0.6

s=1

0.4 0.2 0

0

0.2

0.4

0.6 s=1

0.4

0.6

0.8

0

1

0

PFA

s=3

s=1

PD

0.4

0.4

0.6

0.8

1

0.6

0.8

1

PFA

s=2

0.8

0.6

0.2 s=3

1 s=2

0.8

PD

s=2

0.2

1

0.6

s=1

0.4 0.2

0.2 0

s=3

0.8

PD

PD

1

s=3

2301

0

0.2

0.4

0.6

0.8

0

1

PFA

0

0.2

0.4

PFA Stationary

Non Stationary

Fig. 7. Comparison of ROC. Theoretical, Stationary model and Non-Stationary model for L¼ 5 bands. (a) MIN (b) NMIN (c) SPT (d) FM.

L=5

L=2

s=3

s=3

s=1

s=1

Fig. 8. PT method. Comparison of ROC.

6. A real data experiment To finish the tests, we will apply the different detectors to a series of real records obtained in the laboratory from two different sorts of materials. To analyze the records, windows of size equal to the transmitted pulse will be chosen, which will depend on the transducer and the sampling frequency used in each case. Then, we will move

the analysis window throughout the records applying the detector to each segment. In both experiments, detectors has been designed with 5 bands in the recombination stage and thresholds calculated for a PFA ¼0.01, except for the PT method, whose PFA is set to 0.0625 if we set the number of bands to 5. In the first experiment we analyze a specimen of duraluminum with a hole in one of its faces using a 2 MHz transducer (GE-B2S 2 MHz General Electric). Scans were acquired from the upper face along its longitudinal axis at 40 equally spaced scanning points, thus obtaining 40 scans. Data were collected with Dasel Ultrascope Usb hardware and software acquisition system at a sampling frequency of 250 MHz. Fig. 9(a) shows an image of the specimen, Fig. 9(b) an example of A-Scan obtained at the defect location and Fig. 9(c) the B-Scan obtained with all the acquired scans. In Fig. 9(b) and (c), we can see multiple back-surface echoes. Defect location is marked with a dotted circle. Dotted square in Fig. 9(c) indicates the data analyzed with the detectors. Fig. 10 shows the results obtained after processing the selected area. Black areas denote positive, but not necessary true, detection. Note that the first and the last two detection segments correspond to echoes coming from the back surface of the material. In all cases, true defect location is marked with a black circle (Segment 7 for A-Scans 13, 14 and 15). Table 2 shows PD and PFA obtained for all methods counting the true and false detections and dividing by the number of possible events in each case.

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

2302

5

Defect

10

70 mm

10 mm 60 mm

A Scan

Amplitude

70 mm

15 20

Defect

25 30

30 mm

35

210 mm

40 0

4000

8000

1000

2000

3000

Samples

4000

5000

6000

7000

8000

Samples

Fig. 9. Analysis of the duraluminum specimen. (a) Analyzed specimen (b) A-Scan (c) B-Scan.

10

A Scan

A Scan

10 20 30

20 30

40

40 500

1000

1500

2000

2500

3000

1

2

3

4

Samples

7

8

9

6

7

8

9

6

7

8

9

10

A Scan

A Scan

6

Segments

10 20 30

20 30

40

40 1

2

3

4

5

6

7

8

9

1

2

3

4

5

Segments

Segments 10

10

A Scan

A Scan

5

20 30

20 30

40

40 1

2

3

4

5

6

7

8

9

1

2

3

4

Segments

5

Segments

Fig. 10. Analysis of the duraluminum specimen. Outputs after detection for L¼ 5 bands. Black areas denote positive detection. In all cases, defect location is marked with a black circle. (a) Analyzed B-Scan (b) FM (c) MIN (d) NORM (e) SPT (f) PT.

For the second experiment we chose a block of cement with a linear defect in the middle of the material. The transducer used to collect scans was a 5 MHz transducer with exponential envelope (KBA-5 MHz General Electric). At these frequencies, the selected cement specimens (cement 32.5 30% humidity) are very dispersive. 10 Scans were acquired from the upper face along its longitudinal

Table 2 PD–PFA for the duraluminum specimen.

PD PFA

FM

MIN

NMIN

SPT

PT

0.9877 0.005

0.9754 0.06

0.8711 0.115

0.9693 0.055

0.8466 0.7

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

Scan

FM

MIN

2303

NORM

1

1

1

0

0

1

1

0

0

1

1

0

0

PT

SPT

1

1

0

0

0

1

1

1

0

0

0

1

1

1

0

0

0

0

0

0 1

1

1

1

1

0

0

0

0

0

0 1

1

1

1

1

0

0

0

0

0

0 1

1

1

1

1

0

0

0

0

0

1

1

1

1

1

0

0

0

0

0

1

1

1

1

1

0

0

0

0

0

1

1

1

1

1

0

0

0

0

0

1

1

1

1

1

0

0

0

0

0 1000

2000

0

0

0 0

5

10

15

5

10

15

0

0 5

10

15

5

10

15

5

10

15

Fig. 11. Analysis of the cement specimen. From up to dawn results for each analyzed scan. (a) A-Scan’s (b) FM (c) MIN (d) NMIN (e) SPT (f) PT. Table 3 PD–PFA for the cement specimen.

PD PFA

FM

MIN

NMIN

PT

SPT

0.9 0.007

0.8 0.021

0.5 0.078

0.6 0.085

0.8 0.014

axis at equally spaced scanning points, using the same Dasel Ultrascope Usb hardware and software acquisition system at a sampling frequency of 80 MHz. Fig. 11 shows results obtained after the analysis of the 10 scans from the cement specimen. In all cases, dotted box shows the location of the defect. Table 3 shows PD and PFA obtained for all methods counting the true and false detections and dividing by the number of possible events in each case. Again, results are consistent with those obtained from the theoretical and simulation analysis, showing the FM method the best performance in all cases.

ROC’s have been analyzed in Section 4 as a function of those variables. In Sections 5 and 6, we have tested the detectors using simulated and real scans, confirming the expectations generated from the results of the theoretical analysis. DFT-based formulation is constrained by the limited size of the observation vector dimension which is chosen approximately equal to the pulse duration. However, we may conclude that in spite of that limitation, the detector based on the FM method provides the best performance, which makes it the ideal method to be used as automatic detector of defects in high dispersive materials. Following this detector would come, in this order, those based on SPT, MIN and PT, which, despite being able to provide fairly good results, do not reach the levels of reliability achieved by the FM method. Finally, it should be mentioned that the automatic detection framework considered in this paper may has also practical interest to implement automatic feature extraction in pattern recognition algorithms [32,33].

7. Conclusions and further extensions Acknowledgment In this paper we have tackled the problem of defect detection using the SS algorithm. We have started in Section 2 with an algebraic reformulation of it based on the DFT, which has led us to reformulate also the different recombination methods according to this new approach. In Section 3 we have calculated analytically the statistics of the detectors, i.e., the probability of false alarm and probability of detection, summarized in the ROC curves that allow us to assess the performance of the different recombination methods in terms of the number of bands and the input SNR or the defect amplitude. Theoretical

This work has been supported by the Ministerio de Fomento (Spain), the Ministerio de Ciencia e Innovacio´n (Spain) and FEDER funds under the projects T39/2006, TEC2008-06728. Also partially funded by Generalitat Valenciana (Spain) under program PROMETEO 2010/040. Appendix A. PFA and PD of MIN 2

The probability that min9zi 9 is above the threshold t is the same that the probability that 9zi92 is above the

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

2304

threshold for all possible values of i. On the other hand, under H0, vector z is simply the grain noise component w, hence 9zi92 will be a chi-squared random variable with 2 degrees of freedom [29]. Hence h i Z 1 2 1 ðx=2Þ prob 9zi 9 ¼ x 4t=H0 ¼ @x ¼ eðt=2Þ ðA1Þ 2e t

and considering the L channels the PFA of the method can be written as 2

PFAmin ¼ prob½9zi 9 4t=H0 ; 2

i ¼ 0,. . .,L1 ¼ eðLt=2Þ

ðA2Þ

2

Under H1, 9zi9 /a will have a non-central chi-squared distribution with 2 degrees of freedom with non-centrality parameter equal to 9si92/a2 [4], hence we can write: " # 2 h i 9zi 9 t 2 ¼ x 4 =H prob 9zi 9 4t=H1 ¼ prob 1 2 2

a

Z ¼ Z ¼

a

1

1 ðððs2 =a2 Þ þ xÞ=2Þ e I0 t=a2 2 1

2

neðððs pffi t =a

=a2 Þ þ n2 Þ=2Þ

I0

rffiffiffiffiffiffiffiffiffi! s2 x @x 2

a



ðA3Þ

a a

f X ðx=H0 Þ ¼ 12eðx=2Þ UuðxÞ

where I0(  ) is the modified Bessel function of order zero and Q(a,b) is the Marcum Q-function defined as: Z 1 2 2 Q ða,bÞ ¼ vUeðða þ v Þ=2Þ I0 ðaUvÞdv ðA4Þ b

which can be found tabulated elsewhere. Thus considering the L channels the PD of the method will be  pffiffi h i t s 2 PDMIN ¼ prob 9zi 9 4 t=H1 ; i ¼ 0,. . .,L1 ¼ Q L ,

a a

ðA5Þ

ðB5Þ

and F X ðx=H0 Þ ¼ ð1eðx=2Þ ÞUuðxÞ

ðB6Þ

where u(x) is the unit-step function. Substituting in (24) and using the fact that  L1  X L1 ð1eðy=2Þ Þ ¼ ðB7Þ ðeðy=2Þ Þk k k¼0 we can write the probability density function of v, under H0 in the form f V ðv=H0 Þ ¼ " #  Z 1 L1 L1 L X U ð1Þk yUeðy=2Þðv þ 1Þ Ueðy=2Þk @y UuðnÞ 4 k¼0 k 0 "

a

rffiffiffiffiffiffi ! s2 Un @x 2

pffiffi t s , ¼Q

Under H0, x is a chi-squared random variable with 2 degrees of freedom, therefore

¼ LU

L1 X

ð1Þk

k¼0



L1 k

#  ðv þ 1 þkÞ2 UuðnÞ

Finally, integrating (note that v o1) and considering independence among the different channels, we may find the expressions of the PFA corresponding to NMIN. " # Z 2 1 9zi 9 prob ¼ v 4 t=H0 ¼ f V ðv=H0 Þdv 2 t max9zi 9  Z 1 L1 X L1 ¼ LU ð1Þk ðz þ 1þ kÞ2 @z k t k¼0    L1 X L1 1 1  ðB9Þ ð1Þk ¼ LU t þ 1þ k 2 þk k k¼0 And for the L channels " #L   L1 X L1 1 1  PFANMIN ¼ LU ð1Þk t þ 1 þk 2 þ k k k¼0

Appendix B. PFA and PD of NMIN

ðB10Þ

Let us call v¼x/y with x ¼9zi92 and y¼max9zi92 for i¼ 0yL 1. We can express the probability density function of v in the form: Z 1 f V ðnÞ ¼ f V ðn=yÞf Y ðyÞ@y ðB1Þ 1

But as for a given value of y, v is simply a scaled version of x f V ðn=yÞ ¼ 9y9f X ðynÞ Therefore we can write Z 1 9y9f X ðynÞf Y ðyÞ@y f V ðnÞ ¼

2

2

Now again, under H1, 9zi9 /a will be a non-central chisquared random variable with 2 degrees of freedom and non-centrality parameter equal to s2/a2. Let us make a normalization in the definitions of x ¼9zi92/a2 and y¼max9zi92/a2, which does not affect the quotient n ¼x/ y. Now we have that rffiffiffiffiffiffiffiffiffi! 1 s2 2 2 x UuðxÞ f X ðx=H1 Þ ¼ eðððs =a Þ þ xÞ=2Þ UI0 ðB11Þ 2 a2

ðB2Þ and hence

ðB3Þ

1

where we have assumed that x and y are independent, which can be considered a reasonable assumption for large values of L. Now we may consider the well-known expression [30] for the distributions of the maximum of a set of L i  i  d. random variables, having probability density function fX(x) and cumulative distribution function FX(x), namely f Y ðyÞ ¼ LUF L1 X ðyÞUf X ðyÞ

ðB8Þ

ðB4Þ

"Z

x

1 ðððs2 =a2 Þ þ wÞ=2Þ e I0 F X ðx=H1 Þ ¼ 0 2 h  s pffiffiffi i , x UuðxÞ UuðxÞ ¼ 1Q

a

rffiffiffiffiffiffiffiffiffiffi! # s2 w @w 2

a

ðB12Þ

so that f Y ðy=H1 Þ ¼ LUF L1 X ðy=H 1 Þf X ðy=H 1 Þ rffiffiffiffiffiffiffiffiffi! h  s pffiffiffi iL1 L s2 2 2 ¼ , y eðððs =a Þ þ yÞ=2Þ I0 y uðyÞ 1Q a 2 a2 ðB13Þ

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

and for the L channels, the PFA of the method will be

 PFAPT ¼ prob Re½zi  4 0=H0 ; i ¼ 0,. . .,L1

 L ðC2Þ þ prob Re½zi  o 0=H0 ; i ¼ 0,. . .,L1 ¼ 2U 12

and finally Z

f V ðv=H1 Þ ¼

1

0

1 2 2 y eðððs =a Þ þ yvÞ=2Þ I0 2

 s pffiffiffi iL1 Lh 2 2 1Q , y eðððs =a Þ þ yÞ=2Þ I0 2 a

rffiffiffiffiffiffiffiffiffiffiffiffi! s2 yv 2

a

rffiffiffiffiffiffiffiffiffi! s2 y @y 2

ðB14Þ

a

To achieve the PD of the method we will have to integrate "  2 # Z 1 zi  f V ðv=H1 Þ@v PDNMIN ¼ prob  2 ¼ v 4 t=H1 ¼   t max zi ! rffiffiffiffiffiffiffiffiffiffiffiffi Z 1Z 1  s pffiffiffi iL1 1 s2 Lh 2 2 1Q ¼ y eðððs =a Þ þ yvÞ=2Þ I0 yv , y 2 2 a a2 t 0 2

2

eðððs =a L 4 "Z 

Z

1

¼

Þ þ yÞ=2Þ

a

h  s pffiffiffi iL1 2 2 y 1Q , y eðððs =a Þ þ yÞ=2Þ I0

a

0 1

rffiffiffiffiffiffiffiffiffi! s2 y @y@v 2

I0

2 2 eðððs =a Þ þ yUvÞ=2Þ I0

t

L 4 "Z 

Z

h

1

y

¼

2 2 eðððs =a Þ þ yvÞ=2Þ I0

t

Z 

1

2

2

eðððs =a

1

Þ þ yvÞ=2Þ

I0

a

a

a

0 1

rffiffiffiffiffiffiffiffiffi! s2 y 2

rffiffiffiffiffiffiffiffiffiffiffiffiffiffi! # s2 yUv @v @y 2

 s pffiffiffi iL1 2 2 1Q , y eðððs =a Þ þ yÞ=2Þ I0

rffiffiffiffiffiffiffiffiffi! s2 y 2

a

rffiffiffiffiffiffiffiffiffiffiffiffi! s2 yv @v 2

a

rffiffiffiffiffiffiffiffiffiffiffiffi! # s2 yv @v @y 2

a

rffiffiffiffiffiffiffiffiffi! L s2 y 4 0 a a2 "Z rffiffiffiffiffiffiffiffiffiffi! 1 2 s2 2 2 w @w  pffiffiffiffi weðððs =a Þ þ wÞ=2Þ I0 a2 ty y rffiffiffiffiffiffiffiffiffiffi! # Z 1 2 s2 2 2 w @w @y  pffi weðððs =a Þ þ wÞ=2Þ I0 a2 t y rffiffiffiffiffiffiffiffiffi! Z L 1 s2 2 2 ½1Q 1 L1 eðððs =a Þ þ yÞ=2Þ I0 y ½Q 2 Q 1 @y ¼ 2 0 a2 Z

¼

1

h  s pffiffiffi iL1 2 2 y 1Q , y eðððs =a Þ þ yÞ=2Þ I0

ðB15Þ pffiffiffiffiffi pffiffiffi with Q 1 ¼ Q ðs=a, yÞ and Q 2 ¼ Q ðs=a, tyÞ, so the PD of the method, considering the L channels will be  L1  L1 LX PDNMIN ¼ ð1Þk 2k¼0 k rffiffiffiffiffiffiffiffiffi! Z 1 s2 2 2 ðB16Þ Q k1 eðððs =a Þ þ yÞ=2Þ I0 y ½Q 2 Q 1 @y 2 0

2305

a

that must be solved by integral calculus for each value of the threshold t. Appendix C. PFA and PD of PT This is a much simpler case. Under H0, Re[zi] is a zeromean unit variance Gaussian random variable, hence we can write



 ðC1Þ prob Re½zi  4 0=H0 ¼ prob Re½zi  o 0=H0 ¼ 12

On the other hand, under H1, Re[zi] is a Gaussian random variable having mean s and variance a2, hence for the positive values of Re[zi] we can write Z 1

 2 1 2 eðððxsÞ Þ=2a Þ @x prob Re½zi  ¼ x 4 0=H1 ¼ pffiffiffiffiffiffi a 2p 0   1 s ðC3Þ ¼ erfc  pffiffiffi 2 a 2 And for its negative values

 

 1 s prob Re½zi  ¼ x o 0=H1 ¼ 1 erfc  pffiffiffi 2 a 2

ðC4Þ

Thus, the PD of the method for the L channels will be

 PDPT ¼ prob Re½zi  40=H1 ; i ¼ 0,. . .,L1

 þ prob Re½zi  o 0=H1 ; i ¼ 0,. . .,L1   L   L 1 s 1 s erfc  pffiffiffi ¼ þ 1 erfc  pffiffiffi ðC5Þ 2 2 a 2 a 2 where erfc(x) is the complementary error function defined as: Z 1 2 2 et @t ðC6Þ erfcðxÞ ¼ pffiffiffiffi

p

x

Note that either the PD or the PFA of the method are not dependent on any threshold t as it happens in the rest of the methods. This implies some lack of flexibility for fitting those probabilities, which can be only modified by changing the number of channels L. Appendix D. PFA and PD of SPT 2

In SPT the probability that 9L þ -L  9  min9zi 9 is above the threshold t is equivalent to the probability that 9L þ 2 L  9  9zi 9 is above the threshold t for all possible values of i. Let us define d ¼9L þ -L  9. We consider at this point that L is an odd number, although a similar derivation could be done for an even value. Random variable d is a discrete random variable having values in the set of odd numbers d ¼{1,3,y,L}. Firstly, we have to calculate the probability of every possible value of d. A value d is obtained when L  ¼L þ þ d or L þ ¼L  þd. On the other hand we have that L  ¼L-L þ , so that substituting in the two foregoing possibilities we deduce that a value d is obtained when L þ ¼(L-d)/2 or L þ ¼(Lþd)/2. But L þ is a binomial discrete random variable with a different parameter p under every hypothesis, that is, under H0

 p ¼ prob Re½zi  4 0=H0 ¼ 12 ðD1Þ and under H1

 

 1 s p ¼ prob Re½zi  4 0=H1 ¼ erfc  pffiffiffi 2 a 2

ðD2Þ

Therefore we can express the probability of every possible value d in the form ! ! L L LððLdÞ=2Þ ðLdÞ=2 pðdÞ ¼ Up ð1pÞ þ ðLdÞ=2 ðLþ dÞ=2

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

2306

UpðL þ dÞ=2 ð1pÞLððL þ dÞ=2Þ

ðD3Þ

for d¼1,y,L. Now let us call n ¼9zi92 and x¼d  n. We may write f X ðxÞ ¼

L X

f X ðx=dÞUpðdÞ

ðD4Þ

d odd

But f X ðx=dÞ ¼

1 x f d V d

ðD5Þ

because for every value d, x is a scaled version of v. The probability density function of v is known under both hypotheses (central and non-central chi-square, respectively). Then, under H0 (notice that p ¼0.5 in (42)), we may write for the L channels (integration and sum are interchanged) " ! L h i L 1 X 1 2 U prob dU9zi 9 ¼ x 4 t=H0 ¼ L ðLdÞ=2 2 d odd d !# Z 1 L 1 ðx=2dÞ U þ e dx ðL þ dÞ=2 2 t " ! !# L L L 1 X ðD6Þ ¼ L þ Ueðt=2dÞ ðLdÞ=2 ðL þ dÞ=2 2 d odd

So considering all channels, the PFA of the method will

having mean L  E{ln9zi92} and variance L  var{ln9zi92}. Notice that, considering that 9zi92 is a central chi-squared random variable, we may write h i Z 1 1 2 lnxU Ueðx=2Þ @x E ln9zi 9 ¼ lnx=H0 ¼ 2 0 Z 1 ln2vUev @v ¼ Z 1 Z0 1 ln2Uev @v þ lnvUev @v ¼ ln2g ¼ m ðE2Þ ¼ 0

2

2

¼ ln x=H0 g þ Efln9zi 9 ¼ lnx=H0 g

In a similar manner, under H1 we may write h i 2 prob dU9zi 9 ¼ x 4 t=H1 ! L X L ¼ UpðLdÞ=2 Uð1pÞLððLdÞ=2Þ ðLdÞ=2 d odd ! pffiffiffiffiffiffiffi! L s t=d þ UpðL þ dÞ=2 Uð1pÞLððL þ dÞ=2Þ Q , ðLdÞ=2 a a

0

ðD9Þ

Noting that

i¼0

2

9zi 9 ¼

L1 X

ln9zi 9

2

6

¼Z

ðE4Þ

ðE5Þ

Hence, under H0, the logarithm of the statistic corresponding to FM method is a Gaussian random variable having mean equal to m  L¼ 0.116  L and variance s2  LE1.484  L. Therefore, the PFA corresponding to FM will be given by " # LY 1 2 PFAFM ¼ prob 9zi 9 4 t=H0 ; i ¼ 0,. . .,L1 " ¼ prob ln

i¼o L1 Y

2

#

9zi 9 4 lnt=H0 ; i ¼ 0,. . .,L1

i¼o Z 1

2 1 2 eðððxLmÞ Þ=ð2Ls ÞÞ @x ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pLs2 lnt     1 lntLm 1 lnt0:116L ¼ Uerfc pffiffiffiffiffiffiffiffiffiffiffi C erfc pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2U1:484L 2Ls2

ðE6Þ

Regarding the PD, if we take into account [31], under H1 " E ln

#  2   2  2 9zi 9 s s Ei =H ¼ ln 1 2a2 2a2 2a2

ðE7Þ

where Ei(  ) is the exponential integral function defined as: Z 1 1 x e @x EiðvÞ ¼  ðE8Þ x v

Appendix E. PFA and PD of FM

ln

p

so the variance will be h i p 2 2 var ln9zi 9 =H0 ¼ Zm2 ¼ 2ln 2 þ ¼ s2 6

ðD8Þ And for the L channels h i 2 PDSPT ¼ prob dU9zi 9 4t=H0 ; i ¼ 0. . .L1 ! " L X L UpðLdÞ=2 Uð1pÞLððLdÞ=2Þ ¼ ðLdÞ=2 d odd ! pffiffiffiffiffiffiffi!#L L t=d s , þ UpðL þ dÞ=2 Uð1pÞLððL þ dÞ=2Þ Q ðL þ dÞ=2 a a

ðE3Þ

where the second term is the mean previously calculated and the first one may be calculated as h i Z 1 1 2 2 2 2 ln xU Ueðx=2Þ @x E ln 9zi 9 ¼ ln x=H0 ¼ 2 0 Z 1 2 ln 2vUev @v ¼ Z 1 Z0 1 2 ln 2Uev @v þ 2ln 2Uln vUev @v ¼ 0 0 Z 1 2 ln vUev @v þ 2

ðD7Þ

2

2

¼ ln 2 þ 2ln 2ðln 2gÞ þ g2 þ

d odd

LY 1

2

varfln9zi 9 ¼ lnx=H0 g ¼ Efln 9zi 9

be

h i 2 PFASPT ¼ prob dU9zi 9 4 t=H0 ; i ¼ 0,. . .,L1 " ! !# " #L L L L 1 X þ Uet=2d ¼ L ðLdÞ=2 L þd=2 2

0

where g ¼0.577 denotes the Euler’s constant. On the other side, the variance may be obtained as

ðE1Þ

i¼0

and considering L large, we can make use of the central limit theorem, i.e., x will be a Gaussian random variable

On the other hand, the variance will be ( ) 2 9zi 9 var ln =H1 2a2

ðE9Þ

A. Rodrı´guez et al. / Signal Processing 92 (2012) 2293–2307

which must be estimated experimentally as an analytical expression cannot be obtained as a function of the variable s. Finally, considering the L channels, the logarithm of the statistics of FM will be a Gaussian variable whose mean and variance will be given by (58) and (60), respectively, so the PD of the method will be 0 1    Bln t L ln s2 Ei s2 C B ð2a2 ÞL C 1 2a2 2a2 C ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi s erfcB ðE10Þ   B C 2 2 @ A 9zi 9 Lvar ln 2a2 =H1

References [1] M.G. Gustafsson, Nonlinear clutter suppression using split spectrum processing and optimal detection, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 43 (1) (1996) 109–124. [2] I. Amir, N.M. Bilgutay, V.L. Newhouse, Analysis and comparison of some frequency compounding algorithms for the reduction of ultrasonic clutter, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 33 (4) (1986) 402–411. [3] X. Li, N.M. Bilgutay, J. Saniie, Frequency diverse statistic filtering for clutter suppression, Proceedings of the ICASSP 89 (1989) 1349–1352. [4] R. Drai, F. Sellidj, M. Khelil, A. Benchaala, Elaboration of some signal processing algorithms in ultrasonic techniques: application to materials ndt, Ultrasonics 38 (2000) 503–507. [5] L. Ericsson, T. Stepinski, Algorithms for suppressing ultrasonic backscattering from material structure, Ultrasonics 40 (2002) 733–734. [6] J. Saniie, D.T. Nagle, K.D. Donohue, Analysis of order statistic filters applied to ultrasonic flaw detection using split-spectrum processing, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 38 (2) (1991) 133–140. [7] P.M. Shankar, P. Karpur, V.L. Newhouse, J.L. Rose, Split spectrum processing: analysis of polarity thresholding algorithm for improvement of signal-to-noise ratio and detectability in ultrasonic signals, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 36 (1) (1989) 101–108. [8] N.M. Bilgutay, J. Sniie, The effect of grain size on flaw visibility enhancement using split spectrum processing, Materials Evaluation 42 (6) (1984) 808–814. [9] K.D. Donohue, Maximum likelihood estimation of a-scan amplitudes for coherent targets in media of unresolvable scatterers, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 39 (3) (1992) 422–431. [10] P. Rubbers, C.J. Pritchard, An overview of split spectrum processing, NDT.net 8 (8) (2003). [11] X. Li, M. Bilgutay, Wiener filter realization for target detection using group delay statistics, IEEE Transactions on Signal Proceedings 41 (6) (1993) 2067–2074.

2307

[12] Q. Tian, N.M. Bilgutay, Statistical analysis of split spectrum processing for multiple target detection, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 45 (1) (1998) 251–256. [13] J. Bilgutay, N.M. Sniie, E.S. Furgason, V.L. Newhouse, Flaw-to-grain echo enhancement, Proceedings of the Conference (Ultrasonics International 79), Graz, Austria, 1979, pp. 152–157. [14] V.L. Newhouse, N.M. Bilgutay, J. Sniie, E.S. Furgason, Flaw-to-grain echo enhancement by split-spectrum processing, Ultrasonics 20 (2) (1982) 59–68. [15] L. Vergara, J. Moragues, J. Gosalbez, A. Salazar, Detection of signals of unknown duration by multiple energy detectors, Signal Processing 90 (2) (2010) 716–719. [16] M. Karaoguz, N. Bbilgutay, T. Akgul, S. Popovics, Defect detection in concrete using split spectrum processing, Proceedings of the IEEE Ultrasonics Symposium 1998 1 (1998) 843–846. [17] I. Bosch, L. Vergara, Normalized split-spectrum: a detection approach, Ultrasonics 48 (2008) 56–65. [18] T.Q. Nguyen, S. Jayasimha, Polarity-coincidence filter banks and nondestructive evaluation, Proceedings of the ISCAS 1994 2 (1994) 497–500. [19] A. Rodriguez, L. Vergara, A new split-spectrum algorithm for dispersive materials using variable bandwidth filters, Proceedings of the IEEE Ultrasonics Symposium 2009 1 (2009) 710–713. [20] L.L. Scharf, Statistical Signal Processing, Addison Wesley, 1991. [21] J. Saniie, T. Wang, M. Bilgutay, Statistical evaluation of backscattered ultrasonic grain signals, Journal of Acoustics Society of America 84 (1) (1988) 400–408. [22] R.F. Wagner, S.W. Smith, J.M. Sandrik, H. Lopez, Statistics of speckle in ultrasound b-scans, IEEE Transactions on Sonics and Ultrasonics 30 (3) (1983) 156–163. [23] P. Karpur, O.J. Canelones, Split spectrum processing: a new filtering approach for improved signal to noise ratio enhancement of ultrasonic signals, Ultrasonics 30 (6) (1992) 351–357. [24] M.G. Gustafsson, T. Stepinski, Studies of split spectrum processing, optimal detection, and maximum likelihood amplitude estimation using simple clutter model, Ultrasonics 35 (1997) 31–52. [25] G.S. Kino, Acoustic Waves: Devices, Imaging, and Analog Signal Processing, Prentice-Hall Inc., 1987. [26] G. Cincotti, G. Cardone, P. Gori, M. Pappalardo, Efficient transmit beamforming in pulse-echo ultrasonic imaging, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 46 (6) (1999) 1450–1458. [27] R.B. Kuc, Application of kalman filtering techniques to diagnostic ultrasound, Ultrasonics Imaging 1 (2) (1979) 105–120. [28] A. Rodrı´guez, R. Miralles, I. Bosch, L. Vergara, New analysis and extensions of split-spectrum processing algorithms, NDT and E International 45 (1) (2012) 141–147. [29] S.M. Kay, Fundamentals of Statistical Signal Processing: Detection Theory, Prentice-Hall, 1998. [30] I. Pitas, A.N. Venetsanopoulos, Nonlinear Digital Filters. Principles and Applications, 1st edn., Springer, 1990. [31] S.M. Moser, The expected logarithm of a noncentral chi-square random variable, http://moser.cm.nctu.edu.tw/explog.html, November 2009. [32] A. Salazar, L. Vergara, ICA mixtures applied to ultrasonic nondestructive classification of archaeological ceramics, EURASIP Journal on Advances in Signal Processing 2010, (Article ID. 125201), doi:10.1155/2010/125201. [33] A. Salazar, L. Vergara, R. Llinares, Learning material defect patterns by separating mixtures of independent component analyzers from NDT sonic signals, Mechanical Systems and Signal Processing 24 (6) (2010) 1870–1886.