The work describes an algorithm for 3D radio holographic images synthesis in a subsurface sounding system designed by the authors. The algorithm was parallelized and implemented for multicore processors and NVIDIA GPU. Computational complexity of the algorithm was estimated. Numerical experiments were performed for estimating the algorithm performance and e ciency. Experiments were carried out on multicore processors with vector instructions set and NVIDIA GPU. The developed system can locate inhomogeneity (size of which can be a few centimeters) in doublelayered reinforced concrete at the depth of a few tens of centimeters.
A number of problems, which are solved using the subsurface sounding [
Image synthesis algorithm
The algorithm comprises the following steps [
1. There is a set of beat-frequency waveforms fsk(t)g; k = 1; K, for di erent positions f(xk; yk; zk)g of antenna.
2. Compute the complex spectrums of these waveforms using the Fourier
transform
practical application of the system, it is essential to use e ective software and
hardware for fast processing of radio signals and showing results in the real time.
The authors see two ways of using hardware to implement fast signal
processing. The rst one is to use specialized computers with signal processors and
programmable logic devices. The second one is to use general-purpose personal
computers with high-performance coprocessors. Combination of these two ways
can be used [
3. Perform phasing for the whole set
Sk0( ) = Sk( ) e j 0( ); I(x; y; z) = k X Sk0( )k;
K where 0( ) = 2 (f0 f =2) . Complexity of this step is O(N K).
4. For all coordinates of the required volume (x; y; z), calculate the total intensity where k is the total delay of the echo from the (x; y; z) point considering possible refractions at the media interfaces. The complexity of this step is O(M K), where M = X Y Z is the total number of voxels in the synthesized image, X; Y; Z are the voxel resolutions for each coordinate.
The total algorithmic complexity of the algorithm is O(K(N log(N ) + M )). To achieve the desired resolution of synthesized images, we have to impose certain restrictions onto all components of this expression. Application of the system with the serial sounding allows one to use parallelization for each single signal processing (on each step of the Fourier transform computation and accumulation) and to nd values independently for di erent points using T threads. Thus, the complexity can be reduced to O( KT (N log(N ) + M )).
Parallelization and numerical experiments results Application of the multithreaded system allows one to process the data for each thread independently. So, application of T threads reduces time complexity by T times. However, this also increases space complexity by requirement to store copy of the data for each thread. Combination of these methods can be used. In the experiment, the synthesized image has the following parameters: the resolution in voxels is 100 100 300, the beat-frequency waveforms set K has the length of 2048, and the number of quantization steps is 216. Both the serial and parallel methods were applied to the same image. Table 1 shows the computation times for the serial method and parallel one with four threads. It is seen that the computation time increases linearly with the number of signals.
Another way to increase the implementation performance is to use optimization for the certain processor architecture. For example, contemporary Intel processors support the vector instructions sets. This allows one to increase the performance of oating-point operations signi cantly.
The results of the automatically vectorized program for the Intel Core i7 processor using the Intel C++ Compiler 14.0 are shown in Table 1. The thread parallelization was performed using the OpenMP technology. The number of signals is K = 2048.
The algorithm was also implemented for GPU using the NVIDIA CUDA technology. The experiments were carried out using the NVIDIA GTX 780Ti (2880 cores, 875 MHz). Table 3 shows that this implementation allows one to produce an image with 15 frames per second.
The experiments were also performed using a PC with the Intel Core i7 3610M (2.4 GHz, 9 logical cores) processor with 16 GB of RAM. The algorithm parameters were: M = 106 and N = 262144. The time shown in tables is the minimal time of 1000 experiments with the same conditions. The FFT algorithm was taken from the FFTW and NVIDIA cuFFT libraries. The problem of the low-sized subsurface object visualization motivates creating devices with characteristics extremely close to their theoretical physical limits. Some of the fundamental bases of those limits are a wavelength of the probing signal, aperture sizes, and signal losses in the probed mediums. At the same time, methods of additional image ltering facilitate reaching these limits, which is useful in solving various applied problems.
For example, consider a problem of probing di erent building structures. The purpose of linear ltering consists of highlighting the low-sized objects on a background of high-intensity re ections from reinforcements and wall surfaces. Would the ltering method based on the minimal prior information (e.g. that the wall is considered to be at, the reinforcement bars to be long, etc.) be still e ective?
The principle of the lter operation is the following. Figure 1 (top) shows the simulated two-dimensional image of the reinforcement lattice and the absolute value of its spatial frequency spectrum.
It can be shown that the absolute value does not depend on the shift of the lattice with respect to the picture frame. The shift occurs only in the phasefrequency spectrum. The uctuations of the absolute values of frequency components depend on the lattice's step size. The signi cant volume of spatial harmonic curves for a rectangular reinforcement lattice is concentrated along the coordinate axes. With rotation of the lattice, the spectrum rotates about the vertical axis.
Horizontally oriented ribs correspond to the harmonics located alongside one of the axes, and the vertically oriented ribs correspond to the harmonics located alongside the other one. Extending this conclusion to the third dimension, we can be assured that any orthogonal plain in the 3D image corresponds to the spectrum concentrated alongside a single plain as well.
This fact suggests using a linear 3D lter with the spatial-frequency spectrum that equals 0 for any coordinates that have at least one component at zero and equals 1 in any other case. The lter is supposed to be e ective in clearing any plains and lines from the images that are orthogonal or parallel to any of the coordinate axes. However, the rate of distortion of intensity level and the form of low-sized objects remain vital.
Figure 1 (bottom) shows the simulated image of a low-sized object and its spectrum. As we can see, the signi cant part of spectrum of the object's image is spread along all the frequencies. This fact proves that our lter can be e ective.
Figure 2 shows summary images and their spectra before (top) and after (bottom) ltering. The object's intensity is slightly reduces and the cross-like beams appeared, but these are admissible distortions.
The algorithm for the lter is implemented in the frequency domain by using forward and backward Fourier transform with transitional multiplication of the spectrum by the spatial-frequency function (see Fig. 3). The complexity of the Fourier transform is O(M 3 log(M )), the complexity of multiplication is O(M 3). Thus, the algorithmic complexity of ltration is O(M 3). Analysis of algorithmic complexity of the algorithm suggested for synthesis of 3D radio holographic images was carried out. Application of the serial and parallel processing was considered. The experiments were carried out for performance increasing by application of the parallel computing and vector instructions set. It was shown that by using the GPUs, it is possible to achieve performance su cient for the real-time synthesis of the images. We have shown possibility of using the secondary linear ltering for the radio-frequency 3D images to reduce the masking e ect of wall surfaces and reinforcement lattice as well as highlight the low-sized masked objects. The developed system can locate inhomogeneity (size of which can be a few centimeters) in double-layered reinforced concrete at the depth of a few tens of centimeters.
Acknowledgment This work was partially supported by Russian Foundation for Basic Research (Project No. 15-01-00629).