Comparison of Time Domain and Frequency-Wavenumber Domain Ultrasonic Array Imaging Algorithms for Non-Destructive Evaluation

Ultrasonic array imaging algorithms have been widely developed and used for non-destructive evaluation (NDE) in the last two decades. In this paper two widely used time domain algorithms are compared with two emerging frequency domain algorithms in terms of imaging performance and computational speed. The time domain algorithms explored here are the total focusing method (TFM) and plane wave imaging (PWI) and the frequency domain algorithms are the wavenumber algorithm and Lu’s frequency-wavenumber domain implementation of PWI. In order to make a fair comparison, each algorithm was first investigated to choose imaging parameters leading to overall good imaging resolution and signal-to-noise-ratio. To reflect the diversity of samples encountered in NDE, the comparison is made using both a low noise material (aluminium) and a high noise material (copper). It is shown that whilst wavenumber and frequency domain PWI imaging algorithms can lead to fast imaging, they require careful selection of imaging parameters.


Introduction
Ultrasonic arrays have been widely used for non-destructive evaluation (NDE) in recent years [1]. Ultrasonic arrays have advantages over single-element devices such as the ability of one array to produce beams at a range of angles or focal depths, and the intuitive imaging of the interior of a structure that they produce as an output. In array imaging the overall approach, independent of imaging algorithm is to: (1) fire the array elements in some sequence to generate an ultrasound beam in the structure, (2) the reflected ultrasound signals are captured by the same array elements and (3) these received signals are processed to form the image. This paper discusses two categories of imaging algorithm, time domain and frequency-wavenumber (f-k) domain methods. In the time domain, or delay-and-sum (DAS) algorithms [2], the amplitude of an image pixel is the summation of the amplitudes of the received signals after the application of time delays, typically set to represent a specific wave propagation path. These imaging algorithms are relatively simple to implement and can be configured for a wide variety of inspection cases. Among the possible time domain imaging algorithms, the total focusing method (TFM) [3] and plane wave imaging (PWI) [4] are the focus of this study. The TFM imaging algorithm requires full matrix capture (FMC) in which all transmit-receive signals are captured and hence an N-element array requires N transmission cycles. All the transmitter-receiver pairs are used in the summation to produce a synthetic focus at each pixel. Variants of the approach have been used to inspect structures with multiple layers [5], through nonplanar surfaces [6] and considering wave mode conversions at a defect or an interface [7]. Because the time domain approaches are based on a simple summation, the algorithms are readily parallelisable using graphics processing units (GPU) and field-programmable gate arrays (FPGA). Such approaches have been used for fast TFM imaging [8,9] and can achieve imaging speeds of 150 frames/s for an image with 200 × 200 pixels and a 32 element array [9].
The time domain PWI algorithm uses a discrete number of plane waves with various steering angles to inspect a structure [4]. Relative to the FMC, PWI typically results in a smaller number of transmit cycles and, as the energy associated with each transmitted plane wave is greater than that of an individual element emission, larger penetration depths can be achieved. By optimising the steering angles of the transmitted plane waves, this method can achieve a high frame rate with adequate image quality. The PWI algorithm has also been developed for multimodal imaging cases [10] and inspecting the structures through an irregular surface [11].
The development of f-k domain imaging algorithms was motivated by fast image formation, capitalising on the efficiency of the Fourier transform algorithm and a reduced number of computations. In general, these algorithms transform the array data into the frequency-wavenumber (f-k) domain where a mapping/migration operation is performed, followed by an inverse Fourier transform to the image domain. Stolt [12] first derived what is now known as Stolt's f-k migration algorithm for applications in seismology. After that, the f-k domain imaging algorithm has been developed in both the medical imaging and NDE fields. In medical ultrasound imaging field, Lu et al. [13] proposed a method that uses a single transmitted plane wave to achieve very high imaging frame rate and demonstrated 3750 frames/s in 200 mm depth biological tissues. Chen et al. [14] improved image resolution by extending this to use multiple plane waves at different angles. Garcia et al. [15] applied Stolt's method to PWI by modifying the exploding reflector model and achieved good contrast-to-noise ratio and lateral image resolution. Inspired by medical ultrasound imaging, Stepinski et al. [16] implemented the f-k imaging algorithm in the NDE field and showed good results could be achieved on typical test structures. Skjelvareid et al. [17] used the f-k domain imaging algorithm to inspect multi-layered structures. Hunter et al. [18] proposed the wavenumber imaging algorithm which uses an FMC array data set to achieve a high image resolution and increase the image speed by a factor of the number of array elements compared to the TFM imaging algorithm. More recently, Moghimirad et al. [19] developed a frequency domain algorithm that generates focused virtual sources and this achieved large penetration depth and 20 times faster imaging speed than the time domain equivalent.
Velichko et al. [20] theoretically analysed the connections between the TFM and wavenumber imaging algorithms by expressing them as a linear superposition of phase-delayed transmitter-receiver signals in the frequency domain with different focusing coefficients. They demonstrated that the wavenumber method can provide lower side lobes in the image than the TFM as the transformation process acts as a frequency-independent spatial filter. In terms of computational performance, although the wavenumber algorithm is superior, it requires a large amount of computational memory. Zhang et al. [21] and Anton et al. [22] proposed methodologies to assess the suitability of ultrasound imaging algorithms for NDE, especially for inspecting materials with large grains which can cause strong ultrasound backscattering that results in noise in ultrasound images. A general conclusion is that the TFM imaging algorithm is the most robust algorithm with largest flexibility for inspecting a range of different types of defects [21,22]. Merabet et al. [23] theoretically compared the time domain and f-k domain PWI and concluded that, for computational performance implemented using Matlab, the f-k method is faster than time domain PWI both in 2D and 3D imaging. In terms of image quality, the f-k method improves the lateral resolution because of the spatial-frequency filter effect mentioned above. It is noted that all above comparisons were based on defects at a specific location (generally direct below the array centre). Hence, whilst these papers give insights to imaging algorithm performance in specific NDE inspection examples, they do not present a general picture of the algorithm performance on the diverse range of inspection scenarios encountered.
The motivation of this paper is to comprehensively compare the most commonly used ultrasound imaging algorithms in NDE field with respect to the performance on defects at different locations and in materials with different backscattering noise. Hence, this paper represents a useful guide to imaging performance and how to choose and optimise these imaging algorithms in practice. It is hoped that this represents a further step towards the industrial uptake of these imaging techniques. The algorithms are first overviewed, the methodology of how to choose imaging parameters for time domain and f-k domain PWI algorithms is then proposed, the chosen imaging parameters are finally used to form images for the comparison. Note that given their widespread industrial application, this paper is limited to the application of 1D linear arrays for 2-D imaging on a 2-D structure.

Overview of Ultrasound Imaging Algorithms Compared
In this section, the four chosen imaging algorithms are briefly overviewed. The TFM imaging algorithm and the time domain PWI algorithm are implemented by mapping the delayed and summed time domain signals to the image domain directly. In the wavenumber imaging algorithm and the f-k domain PWI algorithm, the time domain signals are first transformed to the f-k domain using a Fourier transform (FT), the different mapping methods are then applied in the f-k domain and the final image is generated using an inverse Fourier transform (IFT). Figure 1 shows the 1-D array and 2-D test structure geometry as well as the coordinate system used in both the experimentally measured and simulated images. In an FMC data set captured by a 1-D linear array with N e elements from a 2-D structure, e(u, v, t) represents the time domain signal transmitted by the element at (u, 0) and received by an array element at (v, 0). In a plane wave data set with N p transmitted plane waves, f (θ, v, t) represents the time domain signal transmitted with a plane wave steering angle of θ and received by an array element at (v, 0). In the simulation presented in this paper, we synthesize the plane waves from the FMC data set by (1), where, c is the longitudinal wave speed in the test structure material.

TFM Imaging Algorithm
The intensity of the TFM image at pixel located at (x, z) is summation of the appropriately delayed time domain data across all transmitter-receiver combinations and can be expressed as [3], where, τ is the wave propagation time from an array element to a pixel at (x, z) and, where, d is the distance from an array element to a point at (x, z).

Wavenumber Imaging Algorithm
In forming a wavenumber image, the f-k domain data set, E(k u , k v , ω), is first generated from an FMC array data set using a 3D FT by [18], where, ω is the angular frequency, k u and k v are the wavenumbers in the f-k domain described in more detail in Section 4.1. E(k u , k v , ω) is then mapped into the image Fourier domain to generate the Fourier spectrum at each k u using an inverse Stolt mapping. The Fourier spectrum of the final image is the summation of the Fourier spectra at all k u , expressed as [12], where, S −1 [·] denotes the inverse Stolt mapping, k x and k z are the wavenumbers in the image Fourier and are described in more detail in Section 4.1. The relationship between the wavenumbers in the f-k domain and the image Fourier domain is [12], and this leads to, The intensity of the final image is the 2D IFT of S(k x , k z ) expressed as,

Time Domain PWI Algorithm
In the time domain PWI algorithm, the image intensity at a pixel located at (x, z) is the summation of the images from a number of plane waves with different steering angles [4], where, τ is the wave propagation time from an array element to a pixel at (x, z) and,

Lu's f-k Domain PWI Algorithm
In Lu's f-k domain PWI algorithm [13], the spectrum of f (θ, v, t) is first calculated using a 2D FT with respect to v and t, as, F(θ, k v , ω) is then remapped in the image Fourier domain to generate the spectrum of the image using the inverse Stolt mapping, where, The intensity of the final image is the 2D IFT of the summation of all S(θ, k x , k z ) over k x and k z ,

Experimental Setup
A 5-MHz linear array with 64 elements (manufactured by Imasonic, Besancon, France), of parameters shown in Table 1, was used in both the experiments and simulations. The array probe was coupled (using Sonagel-W250, Sonatest Ltd., Milton Keynes, UK) to the top surface of the specimens as shown in Figure 1, its x-location adjusted depending on the defect under test, and an FMC captured. A commercial array controller (Micropulse MP5PA, Peak NDT Ltd., Derby, UK) was used to capture the complete set of time domain signals from every transmitter-receiver pair of the ultrasonic array (i.e., the FMC data set), as shown in Figure 2a. The excitation signal from the array controller is a negative square wave with a pulse width of 80 ns, which leads to the frequency ranging from 0 to 12.4 MHz. The captured data were then exported and processed using MATLAB (The MathWorks Inc., Natick, MA, USA). The experimental data used in this paper are openly available in Supplementary Materials.  Figure 2b shows the details of specimens' geometry and defect locations. An aluminium specimen of 50 mm thick in z direction was used for the low backscattering noise case and a copper specimen of 32 mm thick was used for the high backscattering noise case. The aluminium specimen had three well-spaced 1-mm-diameter side-drill-holes (SDHs) at z = 13, 25 and 38 mm and the copper specimen had three well-spaced 2-mm-diameter SDHs at z = 10, 15 and 22 mm. Note that the data sets used for generating plane wave images were synthesized from the FMC data using (1). Figure 2c shows an example of the experimentally measured pulse-echo signal from an array element, in which the signal reflected from the back-wall is labelled. The back-propagated version of this signal is shown in Figure 2d. In the image from a defect (either simulated or experimentally measured), the array performance indicator (API) is a measure of the size of the point spread function and hence relates to image resolution [3], (15) where, N pi is the number of pixels with an amplitude above some threshold (e.g., −6 dB relative to the peak amplitude at the defect), ∆x and ∆z are pixel sizes and λ c is the wavelength at the centre frequency. The signal to noise ratio (SNR) governs the detectability of a defect and was calculated from, where, I d is the peak amplitude from the image with a defect present and rms(I f ) is the root-mean-squared (RMS) of a region of the image (e.g., a 5 × 5 mm square box) centred on the defect location, but with the array placed over a region without a defect.

Simulation Model
A 2-D linear systems modelling approach [21,24] is used to simulate the FMC data set from a defect and hence produce noise-free images for comparison purposes. In the frequency domain, for a wave path from transmitter element at (u, 0) to a defect at (x d , z d ) and back to a receiver element at (v, 0), the received signal can be expressed as [21,24], where, A 0 is the signal spectrum of the transmitted signal, B is the directivity function of an array element [25], α is the measured material attenuation coefficient as shown in Figure 3 is the scattering matrix [27] of the defect at the wave incident angle of β i and scattering angle of β s . The simulated FMC data sets, e(u, v, t), can be obtained using 1-D IFT on G(u, v, ω) with respect to ω. In the simulation, the defects were modelled as the SDHs with 1 mm and 2 mm diameters and their scattering matrices were obtained from [28]. Note that the longitudinal wave speed in aluminium and copper were measured at 6300 m/s and 4600 m/s respectively. In order to make the simulation close to the experimental measurements, the transmitted signal, A 0 , is taken as the back-propagated signal from the back-wall reflection [29], as shown in Figure 2d. This signal has a bandwidth of 0.6-12 MHz at the amplitude of −40 dB of the maximum. Furthermore, note that before forming all images FMC data sets were filtered using a digital filter with a central frequency of 5 MHz (100% fractional bandwidth at 40 dB).

The Selection of Imaging Parameters
In order to make a fair comparison of the chosen imaging algorithms, reasonable imaging parameters for each imaging algorithm should be found first to ensure good imaging performance is achieved. The TFM imaging algorithm has no free parameters so no optimisation is required. Hence, this section focuses on how to find good imaging parameters for the wavenumber imaging algorithm as well as the time domain and f-k domain PWI algorithms. For the wavenumber imaging algorithm and f-k domain PWI algorithm, the first requirement is to properly define the grid for the image domain. The maximum steering angle and angle increment of plane waves in the time domain and f-k domain PWI algorithms are then discussed and optimised.

Image Grid Definition in the Frequency Domain Algorithms
In the time domain algorithms, the image region and pixel size can be defined arbitrarily, although a smaller pixel sizes lead to smoother images. However, in the frequency domain algorithms, image artefacts produced by Fourier domain wrapping effects [30] must be avoided. Here we assume a Cartesian coordinate system and a rectangular image domain with equally distributed pixels in both x and z directions. The wavenumbers in the f-k frequency domain, i.e., k, k u and k v , and those in the image Fourier domain, i.e., k x and k z , are then distributed as, where, f s is the sampling frequency of the time domain data, p is the array element pitch, N k , N u , N v , N x and N z are the sizes of the k, k u , k v , k x and k z domains respectively. In practice, N k is the number of points used in the Fast Fourier transform (FFT) of the time domain signals, N u and N v are the numbers of array elements (although zero padding can be used to extend this), N x and N z are the numbers of pixels defined in the reconstructed final image along x and z directions. From (18), k varies from −(π f s )/c to (π f s )/c, k u and k v from −π/p to π/p determined by the array element pitch, k x and k z span from −π/∆x to π/∆x and −π/∆z to (π/∆z, respectively, determined by the image pixel size. From (6), k x can by up to 2k u and k z to 2k. In order to efficiently use the information from an FMC array data set, the size of the wavenumber space in the image Fourier domain should be greater than that in the f-k domain. This means k x ≥ 2k u and k z ≥ 2k leading to the requirement for an image pixel size of, ∆x ≤ p/2, ∆z ≤ c/(2 f s ). (19) In practice, most of the transmitted ultrasonic energy is within a limit bandwidth, hence, to reduce computational cost, ∆z in (19) can be taken as the largest frequency in the transmitted bandwidth (i.e., f bw max ≤ f s ).

Angle Parameters of the Transmitted Plane Waves in the PWI Algorithms
Here, the choice of maximum steering angle, θ m , and angle increment, ∆θ, of the transmitted plane waves in the PWI algorithms is investigated for both the aluminium and copper samples. Note that uniform angular spacing of the plane waves is used and so the angular increment, ∆θ, can be calculated from number of plane waves, N p , by ∆θ = (θ Optimal m )/(N p − 1). Plane wave images with θ m varying from 0 • to 80 • and at a fixed the angle increment of ∆θ = 1 • were first generated, and the API (15) and SNR (16) measured from the images.
It is known that when the steering angle of a plane wave is out of range of the geometrical angles between the end array elements and a defect (defined as ϕ 1 and ϕ 2 in Figure 1), the emitted wave does not propagate to the defect as a plane wave, hence it generates unwanted artefacts in the image [31]. Hence, it is common for the plane waves with the steering angle greater than the subtended angles to be discarded. This logic leads to a simple route to choosing θ m which can simply be set separately in the positive and negative direction as, ϕ 1 , ϕ 2 or to a common value of ϕ c = max{ϕ 1 , ϕ 2 }. The θ m that leads to the best API, defined as θ API m , can also be found. Then, using, θ Optimal m = max{θ API m , ϕ c }, the angle increment was varied from ∆θ = 1 • to ∆θ = 20 • with ∆θ Optimal chosen to be as large as possible without significantly compromising image quality. Figure 4a compares the API of the time domain plane wave images extracted from the experimental measurements and simulations on a 1 mm diameter SDH at (0, 13) and (0, 25) in the aluminium specimen. The API decreases with increased θ m (for ∆θ = 1 • ) until a plateau in API performance is reached, i.e., θ API m = 56 • and 40 • for the SDHs at (0, 13) and (0, 25), respectively. Note that these values are close to the geometrical angles as discussed above and shown in Figure 1, i.e., ϕ 1 = ϕ 2 = tan −1 ((N e − 1)p/2z d )) = 56.8 • and 38.4 • for the SDH at (0, 13) and (0, 25), respectively. This confirms that these geometrical angles can be used as a simple way to set the angular range [31]. Figure 4a also shows that the simulations are in good agreement with the experiments (RMS difference between the model and experimental results less than 10% in all cases), meaning that imaging performance and optimisation can be predicted through the simulation.

Case for the PWI Algorithms on the Aluminium Specimen
With regard to the imaging performance of time domain PWI imaging algorithm, Figure 4b shows the experimentally measured peak amplitude from the SDHs as a function of θ m . Peak amplitude is seen to increase with θ m , until it reaches a plateau value at a similar angle to that at which the API plateaued. Figure 4b shows that the experimentally measured noise varies only by a small amount with θ m and Figure 4d shows the experimental SNR which follows the peak amplitude variation seen in Figure 4b. The plateau point in Figure 4d was measured as θ SNR m = 37 • and 16 • for the SDHs at (0, 13) and (0, 25), respectively. As the SNR in Figure 4d is uniformly high (≥35 dB), minimisation of the API was used here to select the maximum steering angle of the transmitted plane waves, i.e., θ   It can be seen that the API remains almost constant with N p (or ∆θ), whereas the SNR increases with increased N p (or decreased ∆θ) until a plateau. Hence, there is no optimal value but the choice of N p (or ∆θ) is a compromise between detection performance (maximise SNR) and inspection time (minimise the number of firings). The approach taken here is to arbitrarily select the number of plane waves, e.g., N p = 21.  , is always close to the maximum geometrical angle, ϕ c = max{ϕ 1 , ϕ 2 } [31]. This means that the optimal imaging parameters depend on the location of the defect, which is typically unknown. The approach suggested here is to consider defects in a range of locations and select image parameters based on a defect in the worse possible location, which from Table 2   Note that ∆θ c is taken at the case of N p = 21. Figure 4a-d also show that the performance of the f-k domain PWI algorithm with angular range follows the same trend as those from the time domain PWI with SNR differing only by ±2 dB. As for the time domain PWI algorithm, the SNR in Figure 4d is high, the imaging performance is optimised by consideration of the API and/or the subtended angles, leading to, θ Optimal m = max{θ API m , ϕ c } = θ API m = 57 • and 39 • for the SDHs at (0, 13) and (0, 25), respectively. These optimal parameters are in good agreement with those for the time domain PWI algorithm. Figure 5a,b also show that the imaging performance of the f-k domain PWI algorithm with angular increment is very similar to the time domain imaging version, with the API having almost no variation with N p (or ∆θ) and the SNR increasing with increased N p . Repeating the same procedure used to generate  Figure 6a,c shows results of time domain PWI and compares the API extracted from the experimental measurements and simulations from a 2 mm diameter SDH in the copper specimen as a function of θ m and N p . As seen in Figure 6a the API varies only slightly with θ m (for ∆θ = 1 • ) which is a distinctly different trend from the aluminium specimen (i.e., Figures 4 and 5). This difference is due to the combined effect of the larger defect geometry used in this case to ensure detection (i.e., a 2 mm SDH) and the multiple scattering from the grain structure which distorts the reflected signals and adds backscattered noise [32]. The effect of defect geometry is explored further in Figure 6e which compares the simulated API as a function of θ m for SDHs with various diameters. As shown, there is a clear difference between small SDHs (e.g., 1 mm) which are good approximations to point scatterers and larger SDHs (e.g., 2 mm) which result in an extended image and a different trend of API as a function of θ m . Despite these differences, it can also be seen that the θ m at which the performance plateaus is similar (54 • ) across the defects and we note that this angle is close to the geometrical angle (53 • ).
The API was also measured from the images after the data sets had been processed using a digital filter with a centre frequency of 1.5 MHz (100% fractional bandwidth at −40 dB) and the results are shown in Figure 6b,d. In this case, with a lower centre frequency, the 2 mm SDHs are now closer to point sctterers and, in line with this the simulation results show a good agreement with the experimentally measured ones. From Figure 6a,c, it is shown that the variation of API with θ m and N p (or ∆θ) is indistinct and hence, unlike the low-noise Aluminium, API is not a good metric to use to select optimal parameters for this noisy material. Figure 7a,b show the SNR of the time domain PWI images from a 2 mm diameter SDH at different depths in the copper specimen, as a function of θ m and N p , respectively. As for the aluminium case, the SNR can be seen to increase with θ m (for ∆θ = 1 • ) to a plateau at around 30 • , which is smaller than the geometric angles of 53 • and 42 • from the defect at (0, 15) and (0, 22), respectively. However, to make choosing imaging parameters simpler, the maximum steering angle was selected as the geometrical angle, θ Optimal m = ϕ c , without compromising the reduction of SNR. Again, following a similar trend to the aluminium, the SNR can be seen to increase as N p increases (or ∆θ decreases). As before, there is no optimal value but the choice of N p (or ∆θ) is a compromise between detection performance (maximise SNR) and inspection time (minimise the number of firings). For comparison purposes, the approach taken here is to select the number of plane waves as the same as the case for the aluminium specimen, i.e., N p = 21. Table 3 shows the experimentally extracted time domain PWI parameters for 2 mm SDH defect in the copper specimen at different locations when N p = 21. It can be seen from Table 3 that the worst optimal case considered was at (25,22) with a SNR of 22 dB for which the PWI parameters were, θ  Note that ∆θ c is taken at the case of N p = 21.

Imaging Performance Comparison
In this section, the computational time from all chosen imaging algorithms is compared. The calculated API and SNR are then used to compare these algorithms in terms of imaging performance. For the scenario where the defect location is unknown and a large inspection area is needed, we fix θ m to vary from −80 • to 80 • and use various N p .

Computational Time
The computational costs of the chosen imaging algorithms are made up of the sum of costs from 1-D and 2-D linear interpolations, 2-D and 3-D FFT, as shown in Table 4 [18,23]. One obvious difference seen in Table 4 between the time domain and frequency domain algorithms is the use of the FFT in the imaging reconstruction. In the wavenumber imaging algorithm, a 3-D FFT is first used to process an FMC data set. The 2-D slices of the Fourier transformed data matrix are then taken at each wavenumber sample to calculate the contributions for each image pixel through linear interpolation. Finally, the contributions at each image pixel are summed and the final image is reconstructed through the 2-D inverse Fast Fourier transform (IFFT). In the f-k domain PWI algorithm, the 2-D interpolation and the 2-D FFT are used to speed up imaging.
The computation performance of the imaging algorithms inevitably depends to some extent on their implementation. To compare their computation efficiency here, the algorithms were implemented in MATLAB (MathWorks, Natick, MA, USA) on a standard computer with Quad-Core Intel i3-8100 CPU and 8 GB of RAM. The images were generated from the aluminium sample using the chosen imaging algorithms within various image areas (N x × N z ). The pixel size was fixed as 0.1 mm, the number of plane waves used for the PWI algorithms is 21 (N p = 21). Figure 8 compares the recorded computational time for the chosen imaging algorithms. Also shown are scaled predictions using the formula shown in Table 4 and the least square fit method, made using the experimentally obtained scale factors α 1 = 1.1 × 10 −8 , α 2 = 1.6 × 10 −8 , β = 0.65 × 10 −8 , γ = 0.5 × 10 −8 .
As shown in Figure 8, the computational costs of the TFM and the time domain PWI algorithms increase with the number of image pixels. The wavenumber and the f-k domain imaging algorithms have almost no change in computation time when the number of image pixels is small and then the computation time increases for a large number of image pixels. For larger imaging areas (or larger numbers of pixels), the optimised PWI algorithms have lower computational costs than the TFM and the wavenumber imaging algorithms due to a smaller number of firings (N p ≤ N e ). The efficiency of the frequency domain algorithms can be seen to be superior to the time domain imaging algorithms for large image dimensions. It is also noted that the frequency domain algorithms require a heavy memory load compared with the time domain imaging algorithms. For example, based on the computer resource used here, for imaging an area with N x × N z = 400 × 300 pixels, the minimum required memory is 26, 1100, 26 and 43 MB for the TFM, the wavenumber, the time domain PWI and f-k domain imaging algorithm, respectively. As is apparent, the largest memory requirement is from the wavenumber imaging algorithm and it is due to a 3D FFT operation in the imaging process [31]. The time domain imaging algorithms can be implemented in a memory efficient fashion by storing only a single time trace (and the final image) in memory at one time. Conversely, the frequency domain imaging algorithms must operate on the entire 3D matrix of echo data simultaneously. It should be noted that the computation efficiency can be increased using GPU parallel processors for both time domain and frequency domain image algorithms.  Figure 9 compares the experimentally extracted API and SNR using the chosen imaging algorithms from the 1 mm diameter SDHs in the aluminium specimen for a range of defect locations. Note that ∆θ = 1 • (N p = 161) is used here as a reference to provide for the best plane wave imaging performance. Due to its large subtended angle, ϕ 1 − ϕ 2 , the defect nearest to the array centre and the top surface, i.e., at (0, 13), has the best API for each imaging algorithm. The TFM and the time domain plane wave images show similar APIs which are between 5% and 55% lower than the other two types of images depending on defect location, with the difference being largest at the extreme, x d = 30 mm. For each defect located at x d ≤ 20 mm, the SNR results from all images are similar with variation within 5 dB and this indicates approximately equal performance from all four algorithms. For each imaging algorithm, the SNR at x d = 30 mm is lower than those at the other locations due to lower directivity of array element and longer ray path.  Figure 10 compares the PWI performance for the cases of different numbers of waves which selected according to; (i) ∆θ = 1 • for best performance (N p = 161); (ii) N p = 21; (iii) the minimum number of plane waves required to cover all defects (N p = 5). As shown, for the same number of plane waves, the SNRs from time domain plane wave images and f-k domain plane wave images are similar with a variation less than 5 dB. The SNR reduction due to fewer plane waves depends on the defect location. For a fixed defect under a fixed imaging algorithm, the difference of the reduction between the cases of N p = 5 and N p = 21 and the cases of N p = 21 and N p = 161 varies from 0 dB to 8 dB. This means that the contribution from an extra plane wave to SNR is higher when number of plane waves is low. Furthermore, it should be noted that for low noise materials, even using 5 plane waves can still reach an SNR above 25 dB.  Figure 11 compares the API and SNR extracted from the chosen imaging algorithms from a 2 mm diameter SDH in the copper specimen. As discussed in Section 4.2.2, the signal can be expected to be distorted due to the multiple scattering from the material grain structure and hence it is difficult to use API to analyse imaging performance. For all defects, the SNR results from the TFM images, the time domain plane wave images and the f-k domain plane wave images show similar SNRs (i.e., within 5 dB) which is at least 5 dB higher than the SNRs from the wavenumber images. This may be due to the number of wavenumbers being 64 (corresponding to k u in Equation (5)) in the wavenumber imaging algorithm, which is less than the number of plane waves (corresponding to θ in Equation (12)), i.e., 161, in the f-k domain plane wave imaging algorithm. Hence noise is better suppressed from more averages in the f-k domain plane wave imaging algorithm than in the wavenumber imaging algorithm.  Figure 12 compares the PWI performance for the cases of different number of waves which selected according to; (i)N p = 161; (ii) N p = 21; (iii) N p = 5. Again, the SNR reduction due to the use of fewer plane waves depends on the defect location. For the copper specimen, more plane waves are needed when compared to the aluminium case to achieve an SNR above 20 dB. Again, the contribution from an extra plane wave to SNR is higher when number of plane waves is low. Note that when N p = 5, the SNRs from the defects located at either x d = 25 mm or z d = 22 mm are less than 20 dB. This is because of the high attenuation meaning that the defects farthest away from the probe have lower SNRs.

Conclutions
The imaging performance of the TFM algorithm, the wavenumber imaging algorithm, the time domain and f-k domain PWI imaging algorithms were investigated and compared.

•
In order to reduce image artefacts, in the wavenumber and f-k domain PWI imaging algorithms, the pixel size in the array lateral direction is required to be less than a half of the pitch of an array element while that in the depth direction less than a half of the ratio between the wave speed and the highest frequency of the transmitted signals.

•
The API performance in the PWI algorithms depends on the subtended angle between an image point and the ends of an array and can be predicted using the proposed simulation model in the single scattering regime. However, when the multiple scattering occurs, the image of the defect is distorted and the SNR is reduced, often making the API unsuitable used for selecting imaging parameters.

•
There is no optimal value for the number of plane waves but the choice of number of plane waves is a compromise between detection performance (maximise SNR) and inspection time (minimise the number of firings). When number of plane waves is high, e.g., Np = 161, for low noise material, all chosen imaging algorithms have similar SNR performance, i.e., all SNRs within 5 dB. However, for high noise material, the TFM imaging algorithm, the time domain PWI algorithm and the f-k domain PWI algorithm have similar performance with SNR at least 5 dB higher than that obtained using the wavenumber imaging algorithm. • 5 plane waves can be used for imaging low noise materials, e.g., aluminium specimens with SNR above 25 dB for a 1 mm SDH defect. However, for imaging materials with high backscattering, e.g., copper specimens, the multiple scattering distorted the API and 21 plane waves were required to achieve SNR greater than 25 dB for a 2 mm SDH defect.