The Spec-Radiation Method as a Fast Alternative to the Re-Radiation Method for the Detection of Flaws in Wooden Particleboards

: For real-time evaluation of non-destructive air-coupled ultrasonic testing of wood-based materials, efﬁcient and reliable calculation methods from ultrasonic holography are essential. Presented here is the spec-radiation method as a fast alternative to the re-radiation method. The spec-radiation method offers a more accurate and up to 88% faster evaluation than the re-radiation method for the determination of ﬂaws in particleboards. Flaws of sub-wavelength sizes can be identiﬁed and their shape and location can be determined with this method. The spec-radiation method produces a better reproduction of the sound ﬁeld than the re-radiation method, especially in the area of the measuring plane.


Introduction
Panel-shaped, wood-based materials are becoming increasingly important in today's world. Whether it is in the furniture industry, where furniture with processed wood is becoming more common as an approach to sustainability, or in the construction industry, where wood is considered a cost-effective, durable material unlike expensive high-performance materials such as Carbon Fiber Reinforced Polymer (CFRP) or Glass Fiber Reinforced Polymers (GFRP). Most big name furniture manufacturers want to offer their products at lower prices and this is usually only possible through material savings, without compromising the load-withstanding capacity of the material. This means that a 100% control of the processed wooden particleboards is necessary. This must be non-destructive, fast, and reliable. This is what non-destructive ultrasonic testing is used for today, which has been continuously developed since Sokolov [1] first used it for material testing. In contact ultrasonic testing, both transmitter and receiver are in contact with the test material. The difference in impedance between the transmitter/receiver and the test material often requires a layer of contact oil. However, this can leave residues on the test material and even be absorbed by the wood material and change its properties (Schafer [2]). The contact force also has a great influence on the accuracy of the results (Gyekenyesi et al. [3]). These led to the development of non-contact testing methods, with air serving as the ambient or contact medium (Fang et al. [4]). Usually, water is used as contact or ambient medium, due to its suitable impedance properties to many materials (Zhang et al. [5], Jasiuniene et al. [6], Mitri et al. [7]). However, due to the large difference in impedance with air as the contact medium, most of the ultrasonic energy is reflected at the interface between the media (Sanabria et al. [8]). To counteract this problem, development of powerful ultrasonic transducers (Hillger et al. [9]) and suitable matching layers between the ultrasonic transducer and the air (Chimenti [10], Álvarez Arenas [11]) were developed, allowing meaningful investigations to the area of the measuring plane. We show this through two experiments with flaw imitations in two dimensions. One has a larger diameter than the wavelength λ and the other has a smaller one. Large distances (approximately 500 mm) between the transmitter through the particleboard to the receiver have to be bridged with the sound attenuation in air depending on the frequency (Jakevicius and Demcenko [39], Stößel [40]). For this, we use a low excitation frequency of 50 kHz. The presented spec-radiation method is up to 88% faster then the re-radiation method. We compare the re-radiation method with the spec-radiation method on the basis of measured data from a transmission test of a wooden particleboard with and without flaws. We show the influence of the measuring grid point distance on both methods. Using an academic piston transducer as an example, we investigate the calculation speed of both methods. This publication is structured as follows: In Section 2 the materials and methods are introduced, with the spec-radiation method in Section 2.1, the re-radiation method in Section 2.2, the experimental setup in Section 2.3, and the practical implementation of the spec-radiation method in Section 2.4. In Section 3 the results are reported. In Section 3.1 the experimental data with a measured grid point distance of 2 mm is studied and the influence of the evaluation window size is investigated. In Section 3.2 the effects of coarse measuring grids are explored. In Section 3.3 are the results from the use of a coarse measuring grid to enhance detectability using the method from Schmelt et al. [22], and in Section 3.4, the difference in the computing speed is investigated. In Section 4 is the discussion.

Materials and Methods
In this section the methods of this publication are presented as well as the experimental conditions. In Section 2.1, the spec-radiation method based on the angular spectrum method is presented. Section 2.2 presents the re-radiation method based on the Rayleigh-Sommerfeld diffraction integral. In Section 2.3, the experimental setup is described. Finally, Section 2.4 presents the practical implementation of the spec-radiation method. The differences to the re-radiation method are elaborated.

The Spec-Radiation Method
The spec-radiation method is based on the angular spectrum method. The following derivation is mainly based on Goodman [41], who described it for light, and on Liu and Waag [38], who described it for acoustics.
The Helmholtz equation describes the propagation of acoustic pressure waves in an incompressible fluid: where x, y, z are the space coordinates, t is the time, c is the speed of sound in the corresponding medium, which is assumed to be constant here, and p(x, y, z, t) describes the acoustic pressure field in space and time. If Equation (1) is transferred to the Fourier domain, the following results are obtained: with the space coordinate vector x = [x, y, z] and: where ω is the time dependent angular frequency. The acoustic sound pressure p(k, ω) can be written as follows: Here, the wave vector k = [k x , k y , k z ] and k x , k y , k z are the spatial angular frequencies and i is the complex unit. Equation (4) satisfies Equation (2) when the exponent factor also satisfies Equation (2), and we get: with k x = 2π f x and k y = 2π f y . The sound pressure distribution in the plane z: This now describes the propagation in the z-direction with the spatial frequencies f x and f y . Two cases must be distinguished: The first case describes homogeneous waves and the second case describes exponentially decaying evanescent waves. If the spectrum of a plane, e.g., the measuring plane with the coordinates x = x , y = y and z = z , is known at z = 0: with F as the Fourier transform operator and Equation (8) substituted into Equation (5), resulting in: To display the result in the spatial domain and not in the Fourier domain, an F −1 is necessary: Depicted here is that only factor e −ik z z is responsible for the propagation of the sound wave. Equation (10) becomes with propagation function H(k x , k y , z): The propagation function can now take two forms, depending on whether the sound waves are to be calculated forward in time or backward. For both forward or backward propagation, the propagation function is now defined as:

The Re-Radiation Method
The re-radiation method, which is based on the Rayleigh-Sommerfeld diffraction integral, is derived in this work according to Sanabria et al. [18], Marhenke et al. [19,20,21], Schmelt et al. [22,23], Tsysar and Sapozhnikov [24], Schmelt et al. [25], Schmerr and Song [42], Sommerfeld [43]. Here, the derivation is equal to the derivation of the spec-radiation method from Equation (1) to Equation (6). As Delen and Hooker [26] described, the inverse Fourier transform of Equation (13) is: and thus describes the derivative of free-space Green's function, with R = (x − x ) 2 + (y − y ) 2 + z 2 as the distance description. A spatial representation is thus derived here. Sommerfeld [44] found that the right term in Equation (13) describes oscillating dipoles. Transferred to acoustics, this corresponds to oscillating point sources. This makes it easier to understand the parameters in the equation. If Equation (13) is inserted into Equation (10), the result is: with x = [x , y , (z = 0)]. Here, V describes the propagation parameter, which for z ≤ 0 is V = −1 and for z > 0 is V = 1. This simplified interpretation of the re-radiation method through representation in the spatial domain is, however, at the expense of additional computational operations, including an additional 2D Fourier Transformation (FFT) that must be performed. A faster calculation is also possible using the spec-radiation method.

Experimental Setup
In the experiment, air was used as the coupling medium between the transmitter and the particleboard and the receiver. As a transmitter, an Airmar pulse ultrasonic transducer was used. The used model AT50 has an active diameter of 45 mm and operates at 50 kHz. Therefore, the wavelength in air at 50 kHz is approx. λ = 6.9 mm. With an amplitude of 80 V, the transmitter was sending 10 sinusoidal oscillations with a repeating frequency of 1 Hz. The distance between the transmitter and the particleboard was always 280 mm and, therefore, greater than the near field length, which is, according to Krautkrämer and Krautkrämer [45], approximately 72 mm, of the transmitter in air at this frequency. During the whole measurement, the transmitter was stationary. The particleboard was a commercially available wooden medium-density fiberboard (MDF). The dimensions of the particleboard measured at height = 800 mm, width = 800 mm, and depth = 20 mm. For wooden particleboards, the choice for suitable flaw imitations are pieces of paper (Marhenke et al. [19,21], Schmelt et al. [22,23]) instead of Teflon tape which is historically used for CFRP and GFRP (Fahr [46]). This is through the assumption that paper has nearly the same material properties as the particleboard. Here, the thin air film between the piece of paper and the particleboard serves as flaw imitation. The air film causes a similar impedance change in the material properties like a real air inclusion at the top of the particleboard. The pieces of paper used as flaw imitations were circular with a diameter of ø20 mm, which is larger than the wavelength, and ø4 mm, which is smaller than the wavelength. The wooden particleboard was also immovable. A Knowles microphone, SPU0410LR5H-QB, was used as a receiver, due to its recording frequency up to 80 kHz and its small dimension with the height = 1.1 mm, width = 3 mm, length = 3.76 mm, and an acoustic port of ø0.25 mm. Since both transmitter and particleboard were immobile, the receiver must move to obtain the information for a C-scan. It was attached to an XYZ traversing stage from RoboCylinder with the maximum traversing path of x = 0-400 mm, y = 0-200 mm, and z = 0-200 mm. During the whole measurement, the XYZ traversing stage did not move in the z direction. The entire system was controlled by a standard computer. Via USB, the XYZ stage unit was connected to the computer and controlled by a LabView program. With a National Instruments (NI) system, the transmitter was controlled and the signal of the microphone recorded. The NI system consisted of a chassis PXIe-1085, with an inserted card for the communication with the computer (PXIe-8398), a card to record the measurement data (PXIe-5171R with a maximum sampling frequency of 250 MHz and 8 channels), and a card (PXIe-5423) as a frequency generator to control the transmitter. For the amplification of the control signal of the transmitter, the power amplifier 1040L of the company E&I was used to generate the required amplitude of 80 V. Figure 1 depicts the experimental setup. The experimental setup is identical to that in Schmelt et al. [22], for completeness it is explained here. The only difference is that in the present study, only the experiments with a flaw imitation ø20 mm > λ and a flaw imitation ø4 mm < λ and a measurement from a flawless particleboard are performed.    The program starts in Figure 2 by loading the measured transient sound pressure field p(x , y , 0, ω) in the measuring plane at z = 0 into the working memory. Afterwards, important parameters are set, such as the sound velocity of air, which can be calculated according to Cramer [47]. In these carried out experiments, the speed of sound was approximately 343 m/s. In the next section, a decision has to be made whether an enhancement of the detectability is desired on the basis of whether a coarse measuring grid is to be used or not. If an enhancement of the detectability is desired, the data is interpolated by a spline interpolation to a finer virtual grid and so p(x , y , 0, t) is interpolated to p(x Virt , y Virt , 0, ). Afterwards, the transient sound pressure field is brought into the time frequency domain by a one-dimensional FFT, whereby the complex amplitudes in the spatial domain are obtained. These amplitudes then go into the block "Radiation Technique", depicted in Figure 3. Here, the re-raditaion and the spec-radiation techniques use different codes. For the re-radiation technique, a grid with the spatial coordinates is created first. With this grid, the term can then be calculated from Equation (13). Next, the 2D FFT of this term is calculated to determine the propagation function H(k x , k y , z). A 2D FFT of the complex amplitudes in the spatial domain is performed as well to transfer the complex amplitudes into the Fourier domain. Then, the multiplication of H(k x , k y , z) and p(k x , k y , 0, ω) is performed, which corresponds to the convolution, and the complex amplitudes in the Fourier domain at the distance z are obtained. The complex amplitudes are transferred from the Fourier domain back into the spatial domain utilizing a 2D inverse Fourier Transform (iFFT). For the calculation with the spec-radiation method, a grid of the spatial angular frequencies k x and k y is created first. Then, a two-dimensional FFT of the complex amplitudes in the spatial domain is performed to obtain the complex amplitudes in the Fourier domain. In contrast to the re-radiation method, the propagation function H(k x , k y , z) can be calculated directly through the grid of the spatial frequencies according to Equations (5) and (12), making a 2D FFT not necessary. The wavenumber step is defined with the size of the measurement field to ∆k x = 2π

Practical Implementation
x max and ∆k y = 2π y max . By multiplying the complex amplitudes in the Fourier domain with the propagation function, i.e., the convolution, the complex amplitudes in the Fourier domain at the distance z are obtained. The complex amplitudes are transferred from the Fourier domain back into the spatial domain utilizing a 2D iFFT. The results can now be displayed as the amplitude or phase in the spatial domain to identify flaws. To compute the time data, the results must be brought back into the time domain with a corresponding iFFT. It is easy to see that the spec-radiation method is faster in the numerical calculation than the re-radiation method because only one 2D FFT is necessary. The advantages of the spec-radiation method will be discussed in the following section.

Results
In this section, the results are reported. In Section 3.1, experimental data with a measured grid point distance of 2 mm is studied. The influence of the evaluation window size on the re-radiation method as well as on the spec-radiation method is investigated. In Section 3.2, the effect of coarse measuring grids are explored. Section 3.3 discusses the results from the use of a coarse measuring grid to enhance detectability using the method from Schmelt et al. [22]. In Section 3.4, the computing speed is determined using an example of an idealized piston transmitter.

Results with 2 mm Grid Point Distance
In this section, the grid point distance is 2 mm in the x and y direction. First, the experimental result without a flaw imitation is evaluated. The size of the evaluation window corresponds to the measurement window and has 101 × 101 grid points. A single frequency (50 kHz) is evaluated here to draw a straightforward comparison. In Figure 4, the computed pressure distribution on the surface of the particleboard utilizing the spec-radiation method is shown. In Figure 4a, the amplitude is depicted in dB related to the maximum value in the plane. The red area represents the position of the transmitter. In Figure 4b, the corresponding phase distribution is depicted. The position of the transmitter can be seen at the red circle area. In Figure 5, the same measurement is evaluated using the re-radiation method. The result is again given in amplitude (in Figure 5a) and angle (in Figure 5b). When comparing Figures 4 and 5, no difference can be seen at first glance. Upon closer inspection, however, a difference in both phase and amplitude is noticeable. In order to present this difference, an difference image DI of the difference function of the amplitudes was generated (see Figure 6), according to: where |p(x, ω) spec | is the pressure amplitude determined by the spec-radiation method and |p(x, ω) re | is the pressure amplitude determined by the re-radiation method. The difference image is now scaled to −40 dB to make the difference clearer. To find the reason for these differences, several horizontal planes with a distance-step of 1 mm were calculated utilizing both methods. Each method is used to calculate the forward and backward propagation, to and from the measuring plane away. The measuring plane is at z = 0 m. The result is given in Figure 7.    It can be seen in Figure 7a that the results of the spec-radiation and the re-radiation method are not identical. Due to a pole at the measurement plane in the re-radiation method, significant differences can be found, with distance to the measurement plane being similar but not equal. The maximum value of both methods changes in relation to the size of the evaluation window. To show this, the measurement window, with a size of 101 × 101 grid points, was inserted into a larger evaluation window which is produced by zero padding. The grid point distance remains dx = dy = 2 mm. In Figure 7b, the respective results for 128 × 128 and 151 × 151 grids are given. With increasing grid size, the difference become smaller. The results depicted in Figure 7c show that for a further increase in the measurement grids, the differences nearly vanish, with only the difference around the measurement plane remaining. The possible reason for this pole in the re-radiation method is the calculation. The re-radiation method starts from Equation (14), where the term for the free space Green's function for the center of the coordinate system (x = y = 0) at a distance for z = 0 with R(0, 0, 0) = 0 takes the following value: On the other hand, Equation (12) is for H(k x , k y , z = 0) = 1, which returns the measurement results again with the spec-radiation method. Since a pole location at the position of the measuring plane does not make sense, due to physical reasons, the spec-radiation method reflects the reality better than the re-radiation method in the region of the measurement plane. The calculated sound fields at the surface of the flawless particleboard with a calculation window size of 256 × 256 grid points are depicted for the spec-radiation method in Figure 8 and for the re-radiation method in Figure 9, respectively. These are optically almost identical, with only phase differences visible at the corners. However, these differences are irrelevant, because this area around the measurement window was only filled with zeros and, therefore, depending on the numerical deviations, large differences in a phase can be expected there. The difference image DI of the two amplitudes is shown in Figure 10. The scaling had to be extended to −60 dB in order to make the minute differences recognizable. In Conta et al. [48], the influence of zero-padding for the fictitious increase of the frequency resolution is described in detail.
The results are now considered identical for a calculation window size of 256 × 256 at a sufficient distance from the measuring plane. In Figure 11, the amplitude distribution of the measurement with a ø4 mm flaw on the top (the side closer to the microphone) of the particleboard is depicted.
In the middle of the red area, i.e., in the middle of the transmitter influence area, a yellow spot can be seen in a size of approximately ø4 mm. This identifies the ø4 mm flaw which is smaller than the wavelength λ = 6.9 mm in air. In Figure 12, the ø20 mm flaw is well identified by the blue area in the influence area of the transmitter. This shows that with the spec-radiation method, it is possible to identify flaws both larger and smaller than the wavelength by evaluating a single frequency reference. The resolution limit for both the re-radiation method and the spec-radiation method, without further processing, is theoretically λ/2, according to Wolf [49] and Simonetti [50]. To visualize the result from the spec-radiation method, we provide a Video Supplementary Materials S1 within the supplementary materials, where the amplitude distribution of each plane above the particleboard can be seen in 3D space. The influence of the flaw and the diffraction behind can be seen.

Influence of the Grid Point Distance
For the efficiency of non-destructive testing with air-coupled ultrasound (ACU), the measurement point distance in the experiment is also an important parameter. The more measurement points are required, the more time for the data acquisition or more equipment (e.g., more receivers) is needed. Using the measurement with the ø20 mm flaw imitate, the influences of the grid point distance is investigated for both the re-radiation method and the spec-radiation method. Here, and for further investigations, the measurement window size is constant at approximately 200 mm × 200 mm, and the evaluation window size is approximately 2.5 times larger, due to zero-padding. Only the grid point distance dx and dy will change. In order to get a good comparison with the re-radiation method, the amplitude distribution on the surface of the particleboard is first calculated and displayed for different grid point distances using the re-radiation method. Figure 13a depicts this for the 2 mm grid. As expected, this is optically identical to the result of the spec-radiation method in Figure 12 from Section 3.1, with a grid point distance of 2 mm. There are artifacts in the corners that can be misinterpreted as another sound source and thus complicate the evaluation with a 8 mm grid point distance (see Figure 13b). However, if the approximate position of the transmitter is known, the flaw is still detectable. With a grid point distance of 12 mm, these errors occur at a higher rate and are closer to the evaluation area (see Figure 14a). However, with knowledge of the approximate range in which the transmitter is located, the flaw can still be identified since a separation from the artifacts is possible. From a grid point distance of 16 mm, the artifacts are so close to the evaluation area that it is no longer possible to ignore them (see Figure 14b). The flaw can at most be identified by exact knowledge of its position, everything else would be an estimation. When calculating the same situations with the spec-radiation method, these artifacts do not occur. (a) (b) Figure 14. Amplitude at the top of the wooden particleboard with a ø20 mm piece of paper. Evaluated with the re-radiation method for the frequency 48-52 kHz, with (a) dx = dy = 12 mm and (b) dx = dy = 16 mm a grid of the size 512 mm × 512 mm. Figure 15a depicts the result for a grid point distance of 8 mm calculated by the spec-radiation method. The transmitter and the flaw can be seen clearly. There are no artifacts that could be mistaken for sound sources. Even with a grid point distance of 12 mm, both the transmitter and the flaw are still clearly visible (see Figure 15b). In the area of the flaw, however, it is noticeable that the amplitude drop is not as strong and, therefore, only a yellowish coloring and no blue coloring appears. While the re-radiation method identifies the ø20 mm flaw with a grid point distance of 16 mm (see Figure 14b) through guessing or by exact knowledge of the position of the flaw, the spec-radiation method identifies both the transmitter and the flaw unambiguously (see Figure 16a). At a grid point distance of 18 mm, the flaw can no longer be identified (see Figure 16b). The transmitter and its position can still be identified. These results show that the spec-radiation method gives clearer results than the re-radiation method using rough measurement grids. The extent to which the results with large grid point distances can be improved and thus increase the detectability of flaws with the spec-radiation method using the technique from Schmelt et al. [22], as with the re-radiation method, is examined in the following Section 3.3.

Results with Enhancement of the Detectability
For the re-radiation method, it has already been shown in [22] that from the data used here, a grid point distance of dx = dy = 18 mm is sufficient to identify the ø20 mm flaw if the time data of the measurement are interpolated to a finer grid before processing. The resulting new grid then contains the measured points and interpolated points. This results in a virtual grid point distance that is smaller than the one originally used for the measurement. In [22], the grid point distance was changed from 18 mm to 2 mm. The command "interp3" with the "spline" interpolation in the program Matlab was identified as a good compromise in effort and accuracy. The same refinement method has been applied to the data from Figure 16b. The sound pressure field is reconstructed by the spec-radiation method, depicted in Figure 17. This shows that the refinement method can be applied here as well. The evaluated data allow a clear identification of the flaw, unlike in Figure 16b. This is expected, since the re-radiation method also produces positives results for this grid point distance with a virtual grid point distance of 2 mm. This is because the calculation methods were not changed and the time data was interpolated before processing. This describes a situation similar to Section 3.1, where the input data and the calculation window are large enough and identical.

Determination of the Computing Speed
To compare the speed and efficiency of the spec-radtiation method compared to the re-radiation method in the numerical calculation obtained using the process shown in Figure 3, an idealized circular piston transducer was investigated. The amplitude distribution in the plane is depicted in Figure 18. In the red area, the amplitude is 1, and in the blue area the amplitude is 0. The signal is a pulse of two 50 kHz sinusoidal waves and the grid point distance is 2 mm. The surrounding medium was air. The propagation was calculated up to a distance of 1 m in increments of 1 mm, which requires the calculation of 1000 planes. In Figure 19, the amplitude distributions of the two calculations are given. These look almost identical. Differences only become prominent when the scaling is extended to −60 dB (see Figure 20). The larger differences become visible especially at the borders. During the re-radiation method, it appears as if sound enters the volume from the sides and the sound from the piston transducer cannot spread uninhibitedly to the sides. This is not the case with spec-radiation. Corresponding to the pole point of the re-radiation method described in Section 3.1 (see Figure 7), it does not give the default values in layer 0, whereas the spec-radiation method reproduces this exactly. For the calculation, a standard computer was used with a processor of the type: Intel(R) Core(TM) i5-6500 CPU @ 3.2GHz with 32GB RAM and a 64 bit operating system with Windows 10 Pro. The calculation of the 1000 planes with the re-radiation method needed a time of 141 s and the calculation with the spec-radiation method needed 16.6 s. The spec-radiation method took about 12% of the time compared to the re-radiation method. To visualize the result from the spec-radiation method, we produced the Video Supplementary Materials S2, where the amplitude distribution of each plane can be seen in 3D space. Therefore, the spec-radiation method is faster and more accurate than the re-radiation method in calculation.

Discussion
In this paper, a new innovative method called spec-radiation for the detection of flaws in wooden particleboards was presented and evaluated. With this method, flaws near the theoretical resolution limit of λ/2 can be identified. We showed in the analytical derivation that the re-radiation method can be easily interpreted by its representation with spatial coordinates, but requires a lot of computational effort. While the spec-radiation method is more complicated to interpret by using the spatial angular frequency coordinates, it is much more efficient to compute. A comparison of the two methods showed that the spec-radiation method can be 88% faster in computation than the re-radiation method. We showed this using an example of an idealized circular piston transmitter. We also showed that the re-radiation method is less accurate than the spec-radiation method, especially in the peripheral areas (see Figure 20). Another important aspect when it comes to time efficiency is the number of measuring points needed. We showed that with the re-radiation method, it is possible to identify the flaw up to a certain grid point distance, but that there are artifacts that make the interpretation less clear (see Figures 13 and 14). In order to get good interpretable results with the re-radiation method despite a large grid point distance, one has to rely on techniques, as in Schmelt et al. [22]. For the spec-radiation method, it was shown that a large grid point distance does not produce the same artifacts as with the re-radiation method and that a ø20 mm flaw can still be clearly identified with a 16 mm grid point distance. For even larger grid point distances, we showed that the technique from Schmelt et al. [22] can also be used to increase detectability. We also showed that the re-radiation method has a pole near the measuring plane, which makes the calculation results inaccurate close to the measurement plane. This pole position does not exist in the spec-radiation method and at the position of the measuring plane exactly the measuring plane is returned. With the help of the attached videos, we showed that we can determine and interpret the amplitude sound fields in space as well as the transient sound fields with the spec-radiation method. We are confident that this is a big step towards real-time analysis in the non-destructive air-coupled ultrasonic testing of wood-based materials.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: