Time-Domain Electromagnetic Identification Based on Rectangular Grooves

Received 31 August 2022, accepted 13 September 2022, date of publication 20 September 2022, date of current version 27 September 2022. Digital Object Identifier 10.1109/ACCESS.2022.3208172 Time–Domain Electromagnetic Identification Based on Rectangular Grooves PETR KADLEC , (Member, IEEE) Lerch Laboratory of EM Research, Department of Radio Electronics, Faculty of Electrical Engineering and Communication, Brno University of Technology, 616 00 Brno, Czech Republic e-mail: [email protected] This work was supported by the Czech Science Foundation under Grant 20-01090S. 1 2 3 4 5 6 7 8 9 10 11 ABSTRACT Unique identification of goods and products is an integral part of today’s automated world. Conventional optical systems cannot be operated in environments of poor visibility. The radio–frequency identification systems require expensive development, production, and installation of passive tags. In this paper, we propose a system that stores the information necessary for the distinction of individual items into a set of grooves with different geometric properties. The individual items are then identified based on observing the response to the electromagnetic wave in the time domain. The proposed identification system benefits from the cooperation of a semi–analytical computational scheme based on the Cagniard–DeHoop Method of Moments and a global optimization algorithm that solves the inverse problem of grooves characterization. The proposed system is validated on several computational examples. Also, the resilience of the proposed system to the influence of the noise added to the observed response is investigated. Finally, the influence of reflected signals on the accuracy of the system is assessed. 12 13 INDEX TERMS Radiofrequency identification, inverse method, particle swarm optimization, evolutionary computation. I. INTRODUCTION14 The identification of goods and products is addressed in15 many industries ranging from ordinary shops, including16 logistics, to high–end technology laboratories. Traditional17 optical bar/QR code systems are reliable and make it easy18 to distinguish between large numbers of individual items.19 Nevertheless, the use of optical systems is limited to (clean)20 environments through which the light can propagate without21 serious scattering effects [1].22 Radiofrequency identification (RFID) has become one of23 the most widely used and profit–generating wireless sys-24 tems [2]. A conventional RFID system consists of an active25 reader device and a passive tag attached to an item that is to26 be identified [3]. The design, manufacturing, and attachment27 of the passive tag are unavoidable steps when using RFID28 technology.29 The contemporary RFID research focuses on the follow-30 ing areas. The development of chipless RFID sensors offers31 The associate editor coordinating the review of this manuscript and approving it for publication was Sotirios Goudos . both increased reliability and cost-effectiveness [4], [5], [6]. 32 Development of efficient RFID readers arrangement algo33 rithms should enhance the capacity of the identification sys34 tems [7]. The pattern reconfigurable readers increases the 35 RFID coverage area, can be used for the on-body imple36 mentation, and indoor localization applications [8]. Authors 37 in [9] and [10] try to integrate the RFID systems to Internet of 38 Things applications. Another important topic is ensuring the 39 security of the RFID systems including used hardware [11] 40 or software protocols [12], [13]. 41 Nevertheless, the information necessary to distinguish 42 individual items can be stored in a set of grooves that are 43 scratched in an attached conducting tag or directly to a con44 ducting surface of the item. Individual items are distinguished 45 by altering the geometrical properties of the selected number 46 of grooves (e.g. their depths and mutual positions). Such a 47 system would have several advantages. The set of several 48 grooves can be scratched on the metallic surface using a very 49 cheap technology e.g. widely available Computer Numerical 50 Control (CNC) machining. Last but not least, it can be oper51 ated in low visibility conditions. 52 100104 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ VOLUME 10, 2022 P. Kadlec: Time–Domain Electromagnetic Identification Based on Rectangular Grooves Wakuami et al. [14] proposed an identification system53 based on the principle of magnetic sensing. In the current54 paper, we propose an identification system based on the55 same principle that distinguishes individual items based on56 sensing the response to an electromagnetic (EM) plane wave.57 According to our knowledge, no identification system uses58 the sensing of time-domain (TD) scattered EM fields from a59 set of grooves to distinguish between a set of items.60 An essential part of such a system is a tool calculating61 the EM fields scattered by a metallic surface with rect-62 angular grooves. Most of the computational methods are63 based on integral formulation and therefore apply to the64 time–harmonic fields only [15], [16], [17], [18], [19]. The65 Finite-Difference Time-Domain approach was used to solve66 for the EM fields scattered by grooves in [20]. The authors67 in [21] proposed a method to characterize the properties of68 a single groove on a metallic surface by the angular dis-69 tribution of the scattered light generated by a set of plane70 waves. Authors in [22] derived a method based on Geomet-71 ric/Physical Optics to analyze the terahertz–wave scattering72 characteristics of objects with multiple small–scale grooves73 having the same dimensions. A semi–analytical method to74 solve for the electromagnetic (EM) scattering by a 2D groove75 on a perfectly electric conducting (PEC) surface was fully76 derived in the time domain in [23]. This approach was then77 extended to any number of grooves in the book [24, Ch. 13.2].78 Recently, we published a feasibility study [25] in which we79 successfully solved the inverse problem of characterization of80 a single groove. More specifically, the geometrical properties81 of a rectangular groove were retrieved from the observed82 voltage response using an optimization method. In this con-83 tribution, we further extend this approach - properties of84 two grooves are determined based on the voltage response85 computed by the method introduced in [24] for a groove con-86 figuration proposed by an optimization algorithm. Through-87 out the paper, the computational method will be called the88 forward solver while the optimization algorithm will be89 referred to as the inverse solver. A comparative study of one90 local and four global optimization algorithms namely the91 quasi–Newton Broyden–Fletcher–Goldfarb–Shanno (BFGS)92 algorithm [26], and Particle Swarm Optimization (PSO) [27],93 Differential Evolution (DE) [28], Self–organizing Migrating94 Algorithm (SOMA) [29], and Genetic Algorithms (GA) [30]95 is performed to find the most appropriate inverse solver. The96 limitations of the potential identification system are criti-97 cally discussed based on the results of various examples.98 We investigate the influence of Additive White Gaussian99 Noise (AWGN) and the presence of echoed signals on the100 reliability of the identification system.101 The paper is organized as follows. Section II describes102 the inverse problem of characterization of geometrical prop-103 erties of two rectangular grooves. It is formulated as a104 single-objective optimization problem. Section III reviews105 the numerical computational tools that solve the formulated106 problem. Several instances of the grooves characterization107 problem are defined, solved, and discussed in Section IV.108 Finally, Section V concludes the paper with a critical dis109 cussion of the advantages and limitations of the proposed 110 identification system. 111 II. PROBLEM DESCRIPTION 112 The problem under consideration consisting of two rectan113 gular grooves etched in the PEC surface, denoted by GA114 and GB, is shown in Fig. 1. Here, the position is located 115 using the Cartesian reference frame with the origin Oand 116 the (standard) base {ix,iy,iz}. The centers of the pertaining 117 apertures are located at {0,0,0}and {xB,0,0}, respectively. 118 The grooves are considered to have the same width wand can 119 vary in their depths dAand dB, respectively. Grooves GAand 120 GBoccupy domains {−w/2<x<w/2,−∞ <y<∞,121 0<z<−dA}, and{xB−w/2<x<xB+w/2,−∞ <y<122 ∞, ..0<z<−dB}, respectively. 123 The grooves are surrounded by a linear, isotropic and loss124 free medium D0that is characterized by its electric permittiv125 ity and permeability µ. Both EM parameters are considered 126 to be scalar, real-valued, and positive. Then, the EM wave 127 propagates through D0with speed c0=(µ)−1/2>0 and 128 corresponding wave admittance Y0=(/µ)1/2>0. 129 The corrugated surface is illuminated by a pulsed EMplane 130 wave that is defined by its pulse shape, ei(t), and the angle of 131 incidence, θ. The superscript istands for ‘‘incident’’ and t132 denotes the time coordinate. 133 The proposed identification system is supposed to distin134 guish individual users (or shortly IDs) based on sensing the 135 scattered EM fields in a form of the voltage response mea136 sured over the grooves. Every ID then has its own unique set 137 of grooves. The degrees of freedom are geometrical proper138 ties of the grooves namely depths dA, and dB, and their mutual 139 distance xB. To recognize a correct ID (i.e. values dA,dB,140 and xB) based on the observed voltage Vopretends an inverse 141 problem. The inverse problem is solved using the global opti142 mization algorithm (inverse solver), that proposes a candidate 143 solutions whose correctness is evaluated by the CDH–MoM 144 (forward solver) until the correct ID is not determined. 145 The inverse problem of determining the geometrical prop146 erties of the grooves can be formulated as a single– objective 147 optimization problem:148 min uf(u)= Nt X k=1Vc k(u)−Vo k149 s.t. u∈0(1) 150 where k=1,2,3,...,Ntdenotes individual time samples, 151 Ntis the total number of samples considered. Furthermore, 152 symbol Vkstands for the k-th sample of the voltage induced 153 across the grooves. Superscripts cand odistinguish the volt154 age being ‘‘computed’’ (by a forward solver) and ‘‘observed’’ 155 (e.g. by measurement, or by a forward solver with an added 156 AWGN). The decision space vector uconsists of three param157 eters of the analyzed configuration {dA,dB,xB}. They are 158 highlighted in red in Fig. 1. The symbol 0denotes the 159 decision space (all feasible combinations of dA,dB, and xB). 160 VOLUME 10, 2022 100105 P. Kadlec: Time–Domain Electromagnetic Identification Based on Rectangular Grooves FIGURE 1. Two rectangular grooves in a PEC surface (the variables of the decision space are highlighted in red). Objective function (1) simply minimizes the absolute error161 between the observed voltage and the one computed based162 on parameters proposed by the forward solver (optimization163 algorithm).164 III. COMPUTATIONAL METHODS165 The whole process of solving the inverse problem (see (1))166 is governed by the inverse solver. Any single–objective con-167 strained optimization algorithm can be used as the inverse168 solver. The inverse solver proposes estimations of the deci-169 sion space variables u={dA,dB,xB}. The forward solver170 then computes the voltage Vcfor the particular uto eval-171 uate it using (1) whose value indicates the quality of the172 proposed u.173 Concerning the efficiency and robustness of the solu-174 tion, we opted for the PSO algorithm as the inverse solver.175 As shown by the results of our previous work [25] and176 in Fig. 6, it outperforms the other state–of–the–art algo-177 rithms BFGS, DE, SOMA, and GA significantly. The forward178 and inverse solvers are briefly reviewed in the following179 subsections.180 A. FORWARD SOLVER181 The full derivation of the method based on the reci-182 procity theorem of the time-convolution type [31] and the183 Cagniard–DeHoop Method of Moments (CDH–MoM) [32]184 can be found in [24, Ch. 13.2]. Pursuing this approach, the185 voltage response of the surface Vcon the EM pulse ei(t) can186 be obtained via the marching–on–in–time technique:187 Vm=Y−1 1·Hm− m−1 X k=1Ym−k+1−2Ym−k 188 +Ym−k−1·Vk(2)189 for m=1,2,...,Ntbeing the indexes of the time samples.190 Symbol Y−1 1denotes the inverse of the admittance matrix for191 the first time step t1=1t. The excitation of the system (2)192 is given by 193 Hm=−2Y0ei(tm) 0,(3) 194 where Y0is the free–space admittance. The elements of 195 2×2 matrices Ykin (2) are obtained:196 Yk=2Y0 c01t8A(tk)9(tk)/2 9(tk)/28B(tk),(4) 197 where 9(t)is given by 198 9(t)=1 2πcosh−1c0t xBHt−xB c0,(5) 199 where symbol H(t)stands for the Heaviside unit step func200 tion. Then, the main diagonal components of (4) are given by 201 8A,B(t)=h8(t)+8− A,Bi/2.(6) 202 The TD function for 8(t)reads 203 8(t)=hϒB(w,t)−2ϒB(0,t)+ϒB(−w,t)i/w2,(7) 204 where 205 ϒB(x,t)=1 2πx2 2cosh−1c0t |x|+c0tx 2 c2 0t2 x2−1!1/2 206 −c0t|x|tan−1  c2 0t2 x2−1!1/2 207 ·Ht−|x| c0+c0tx 2H(x)H(t),(8) 208 for all x∈Rand t>0. The limit for x=0 is then 209 ϒB(0,t)=c2 0t2 4πH(t).(9) 210 Finally, the second component of (6) is given by 211 8− A,B=c0t 2wH(t)212 +1 w ∞ X n=1c0t−2ndA,BHc0t−2ndA,B.(10) 213 Throughout this study, we use a TD causal bipolar pulse 214 with the triangular signature defined as [23] 215 ei(t)=2 twtH(t)−(2t−tw)H(2t−tw)216 +(2t−3tw)H(2t−3tw)217 −(t−2tw)H(t−2tw)(11) 218 to excite the grooves. In (11), symbol twdenotes the base 219 length of the triangle pulse. The pulses impinging grooves 220 GAand GBare shown in Fig. 2. The delay between the time 221 of arrivals depends on the distance of grooves xBand on the 222 incident angle θ. The pulse (11) can be easily generated on a 223 standard signal generator [33]. 224 100106 VOLUME 10, 2022 P. Kadlec: Time–Domain Electromagnetic Identification Based on Rectangular Grooves FIGURE 2. Signature of the excitation pulse eitat grooves GAand GB. The accuracy of the forward solver is validated by a225 direct comparison with the full–wave solver CST Microwave226 Studio. The voltage responses obtained by the CDH–MoM227 and CST for three different configurations of the grooves228 ID 1-3 are shown in Fig. 3. Here, voltages over grooves GA 229 and GBobserved by both computational methods are in an230 excellent agreement. However, the calculation of one config-231 uration of grooves takes only a few ms using CDH–MoM232 while one run of CST takes about 250s. This difference in233 computation times brings a significant acceleration of the234 whole identification process because the forward solver is235 called by the inverse solver several times during the identi-236 fication process.237 The responses shown in Fig. 3indicate the variation of the238 observed TD responses when changing grooves parameters239 dA,dB, and xB. The responses are completely different when240 two of three parameters are changed like in case of ID 1 and241 ID 2. The responses have a similar shape if only one of the242 parameters is slightly changed like in case of ID 1 and ID 3.243 Nevertheless, as shown by the results presented below (please244 refer to Fig. 7) the sum of differences between those two245 curves provides still enough information to the PSO inverse246 solver to distinguish between the two IDs provided by these247 two configurations.248 To further show that the proposed system will be able to249 safely distinguish between IDs, we introduce the quantity of250 the objective function margin:251 1=10log min ∀u∈0,u6=ID Nt P k=1 |V(u)−VID| Nt P k=1 |VID| ,(12)252 Quantity 1simply shows what is the difference between253 the observed voltage for the particular ID VID and the other254 ID that produces the most similar voltage response. Larger255 values of 1should mean easier identification of the partic-256 ular ID. A zero signal would have the margin 1=0dB.257 FIGURE 3. Comparison of voltage responses observed by CDH–MoM (solid curves) and CST (dashed curves) solvers for three problem instances ID 1 (dA/w=5.0, dB/w=5.0, xB/w=7.0), ID 2 (dA/w=8.0, dB/w=5.0, xB/w=5.0), and ID 3 (dA/w=5.0, dB/w=5.0, xB/w=6.0). FIGURE 4. The objective function margin 1for all possible IDs of the proposed system. The distribution of the margin quantity over all IDs of the pro258 posed system is shown in Fig. 4. The value of 1= −10dB 259 can be expected for the whole decision space 0what should 260 be sufficient for a reliable identification. 261 B. INVERSE SOLVER 262 Any single–objective optimization algorithm can be used as 263 the inverse solver of problem (1). The effectiveness of the 264 optimization process varies when using different optimizers 265 on the same problem according to the well–known ‘‘No 266 free lunch’’ theorem [34]. Therefore, we have tried different 267 state–of–the–art optimization algorithms on the single groove 268 characterization problem [25] and as well in the current study 269 where we face the problem to characterize two grooves. 270 Namely, we compare the results obtained by quasi–Newton 271 BFGS method [26], PSO [27], DE [28], SOMA [29], and 272 GA [30]. As shown by the results (please, refer to Sec. IV and 273 [25, Sec. IV]), PSO outperforms the other algorithms for this 274 VOLUME 10, 2022 100107 P. Kadlec: Time–Domain Electromagnetic Identification Based on Rectangular Grooves type of optimization problem. For the sake of brevity, we shall275 limit our description to the PSO algorithm only.276 PSO was introduced by Kenedy and Eberhart [35]. It fully277 exploits the intelligence of the swarm of particles. By a278 particle, we mean a decision space vector uthat moves in279 individual iterations of the algorithm based on information280 about the decision space obtained so far by the whole swarm.281 PSO starts with a random generation of positions u(and282 velocities v) for individual particles:283 up=umin +r(umax −umin)(13)284 Both uand vare vectors of length Dwhich is the number of285 decisions space variables (i.e. optimized parameters). Index286 pdenotes the p-th particle, umin and umax are the lower and287 upper limits of the decision space 0, and rdenotes the random288 vector of values chosen from interval h0,1iwith a uniform289 probability. Symbol has the meaning of the element–wise290 multiplication of two vectors with the same lengths.291 Every particle then changes its position in every iteration i292 according to the formula293 up(i)=up(i−1)+1tvp(i)(14)294 where 1tis the time step that equals 1 and is present for the295 sake of physical correctness, only.296 The most important equation of the PSO algorithm is the297 velocity update formula that reads:298 vp(i)=wvp(i−1)+c1r1pbp−up(i−1) 299 +c2r2gb −up(i−1).(15)300 Here, w,c1, and c2are the user–defined parameters con-301 trolling the flow of the algorithm namely the inertia weight,302 the cognitive learning factor, and the social learning factor,303 respectively. Symbol rdenotes the random number from304 interval h0,1i, again. Finally, symbols pb and gb stand for the305 personal and global best position, respectively. Every particle306 keeps its own pb where it found so–far the best value of the307 objective function. Global best is then the best position visited308 by the whole swarm.309 Equation (15) redirects the particles based on three tenden-310 cies that are balanced partially by the user’s choice (values of311 parameters w,c1, and c2) and partially by the randomness312 (values r1, and r2): 1) a particle remains to move in its313 previous direction (controlled by w), 2) a particle is attracted314 to its pb position (controlled by c1,r1), and 3) a particle is315 attracted to the global best position gb (controlled by c2,r2).316 The proper setting of w,c1, and c2balances the exploitation317 and exploration properties of the algorithm.318 After the new position for each particle is set accord-319 ing to (14), these positions are checked if they are feasible320 (i.e. inside of the decision space 0). If any of them is out321 of 0, one of the absorbing, reflecting, or invisible boundary322 conditions has to be applied to the violating particles [36].323 Then, the new positions u(i)are evaluated using the objec-324 tive function f(1). The personal best positions and the325 global best positions are updated for the particles that find326 better positions. This process repeats until the stop condition 327 is met. Usually, the maximal number of iterations Iis used 328 as one of the conditions along with the sufficient value of the 329 objective function. The pseudocode of the PSO algorithm is 330 summarized in Alg. 1. 331 Algorithm 1 Pseudocode of the PSO Algorithm Input: Set of parameters S, objective function f, decision space 0 Output: Vector u∈0with minf(u) 1: Generate random particles and velocities u,v∈0 2: Compute f(u) 3: while i≤Ido 4: Update v(i)and u(i) 5: Check u(i)∈0 6: Compute new f(u(i)) 7: Update pb and gb 8: i=i+1 9: end while IV. RESULTS 332 We performed several experiments to assess the feasibility 333 of the proposed methodology in the identification problem. 334 Accordingly, the system consisting of the forward and inverse 335 solver as described in the previous section is asked to find 336 the geometrical properties of two grooves. The quality of the 337 search is expressed through the decision space error [25]: 338 DER ="D X d=1ubest d−u∗ d2#1/2 ,(16) 339 where superscript ∗marks the true (optimal) decision space 340 vector and best indicates the decision space vector found 341 by the optimization algorithm. Symbol Dis the size of the 342 decision space (D=3 for our problem). The DER value is 343 simply the Euclidean distance between the found and correct 344 solutions. It should be noted that DER is normalized with 345 respect to the grooves’ width. We take w=1.0mm in the 346 examples that follow. 347 All the tests were performed 100–times to eliminate 348 statistical anomalies that can occur when using stochas349 tic optimization algorithms such as PSO. All the tests 350 are carried out for the decision space 0defined as 351 follows: {3.0<dA/w<10.0},{3.0<dB/w<10.0}, and 352 {1.0<xB/w<10.0}. The limits of the proposed 0are 353 selected taking into account the capabilities of today’s sig354 nal generators [33, Ch. 5] and the spatial possibilities - the 355 grooves used for the identification should not occupy an 356 exceedingly large area or/and be too deep. 357 The parameters of the excitation pulse are presented in 358 Fig. 2. The voltage responses used to compute objective 359 function (1) have Nt=101 samples. A single run of 360 the forward solver takes approximately 4.0 ms on a stan361 dard PC. The optimization algorithms implemented in FOPS 362 (an in–house MATLAB toolbox [36]) were used. All the 363 100108 VOLUME 10, 2022 P. Kadlec: Time–Domain Electromagnetic Identification Based on Rectangular Grooves tests could use 1000 objective function evaluations during364 one run (i.e. 20 agents and 50 iterations). Default values of365 controlling parameters as defined in FOPS were used for366 individual optimizers (please refer to [36] for further details).367 The objective function of a problem instance having368 d∗ A/w=5.0, d∗ B/w=7.0, and x∗ B/w=5.0 is shown in369 Fig. 5. There, three cuts are shown to better visualize the 3D370 function (xB=x∗ Bfor the dA-dBcut and so forth). It can371 be seen, that the objective function is relatively smooth and372 valley–shaped for all three cuts. The objective function value373 increases most rapidly with a change of the dAvalue. But374 the most important finding is that the objective function has a375 distinguishable difference between the fvalue for the optimal376 u∗and for uin a close distance from u∗. Let us consider the377 distance 1 ×was a close one. This distance is crucial for the378 intended identification system.379 A. COMPARATIVE STUDY380 The effectiveness of the optimization algorithm as an inverse381 solver is the most important parameter for the potential iden-382 tification system. Too slow convergence or a large variation in383 the results would disqualify the system to be used for real–life384 applications. Therefore, we assessed five algorithms to see385 which one converges most quickly and surely to the global386 optimum. Namely, we compared BFGS, PSO, DE, SOMA,387 and GA algorithms.388 Results of the DER metric for individual algorithms are389 summarized in the form of standard boxplots in Fig. 6. The390 BFGS algorithm achieved the worst results which is not391 surprising as it was started at the random position in 0for392 every run. Quasi–Newton methods tend to converge to the393 closest local optimum. PSO outperforms also the other global394 algorithms (DE, SOMA, GA). The worst results obtained395 by PSO are at the same level of DER as the best results396 of the other algorithms. The median value obtained by PSO397 is approximately two orders of magnitude better than the398 median values of DE, SOMA, and GA. The excellent results399 of PSO are caused by the valley–shaped nature of the objec-400 tive function and the low number of decision space variables401 D=3. Usually, DE and GA outperforms PSO in case of402 higher dimensions D10 [37].403 In this study, we consider an identification system with404 individual items characterized by three geometrical proper-405 ties dA,dB, and xB. Individual users are mutually differen-406 tiated by changing one of the three parameters by at least407 1×w. Therefore, the critical value of the DER metric meaning408 that the inverse solver was successful or not is DER =0.5409 (it is shown as a red dashed line in Fig. 6). Only PSO and410 DE achieved DER well below this particular value for all411 independentruns whichmeans thatthese twoalgorithmswere412 successful in all runs. Using PSO, there is still a margin of413 two orders in magnitude. Therefore, a lower number of agents414 and iterations could be used to make the identification process415 faster.416 B. RELIABILITY 417 Another watched parameter of every system is its reliability. 418 Thanks to the fact, that our identification system has only 419 three degrees of freedom, we tested the PSO inverse solver for 420 all the possible combinations that are feasible according to 0421 (see Sec. IV). Considering that all the parameters are sampled 422 with 1×wour system can distinguish 8×8×10 users. All the 423 resulting 640 problem instances were solved 100 times by the 424 PSO inverse solver. The mean value of DER metric is plotted 425 as a sliced 3D graph in Fig. 7.426 We can see the shades of dark blue colors in the whole 427 range of 3D space 0. The DER grows to larger values only 428 for the lower limit of GAdepth dA/w=3.0. These errors 429 are caused by the numerical issues of the forward solver. 430 The numerical discrepancies can be neglected by reducing 431 the width of the spatial support of the excitation pulse or 432 by decreasing the time step 1t. The first option may run 433 into the limits of available signal generators. The second 434 option will lead to a slowdown of the identification process. 435 Nevertheless, the capacity of the identification system can 436 be enhanced by adding more grooves which requires only a 437 minor change in the forward solver. 438 C. RESILIENCE TO NOISE 439 In a real scenario, the influence of noise corrupting the 440 ‘‘observed’’ signal Vohas to be taken into account. Therefore, 441 the response Vowas corrupted by the additive white Gaussian 442 noise. AWGN with SNR = {5.0, 10.0, 15.0, 20.0, 25.0, 30.0}443 (in dB) was added to Vo. The statistical results are shown 444 in Fig. 8. Here, we can observe that the proposed system can 445 identify the correct user in 100% of cases for SNR ≥15.0 dB. 446 Only a few outliers (5 of 100 trials) are misidentified for 447 SNR =10.0dB. 448 D. ECHO SIGNAL INFLUENCE 449 In this experiment, we analyze the impact of a delayed and 450 attenuated signal (e.g. a reflected signal) on the accuracy of 451 identification. For these purpose, the observed signal Vois 452 composed of two voltages: the primary response Vp(with no 453 delay or attenuation) and echoed response Ve. The echoed 454 voltage is attenuated 10-times and is delayed by a value τ.455 An example of the thus distorted voltage signal is shown 456 in Fig. 9.457 The DER values for different values of the delay τc0/w=458 {5.0,10.0,20.0,50.0}are shown in Fig. 10. The influence of 459 the echoed signal increases with the decreasing delay τwhen 460 the echo signal changes the strong initial part of the observed 461 pulse. On the contrary, the observed voltage is corrupted by 462 zeros for the most time samples when τbecomes large in 463 comparison to spatial parameters of grooves. Fig. 10 shows 464 that the identification is virtually flawless for τc0/w≥10. 465 Therefore, the proposed approach can easily handle unwanted 466 signals due to scattering by objects located in the close vicin467 ity of the device under test. 468 VOLUME 10, 2022 100109 P. Kadlec: Time–Domain Electromagnetic Identification Based on Rectangular Grooves FIGURE 5. Three cuts of the objective function for problem instance: d∗ A/w=5.0, d∗ B/w=7.0, and x∗ B/w=5.0. FIGURE 6. Achieved decision space error DER for different optimization algorithms. FIGURE 7. Decision space error DER for PSO algorithm for all possible instances of problem (1) when individual parameters dA,dB, and xBare sampled with step 1 ×w. E. SENSOR POSITION INFLUENCE469 In the current state, the TD response Vis observed just470 above the aperture of the grooves. The real-life system would471 FIGURE 8. Decision space error DER for PSO algorithm for various SNR values of noise added to the ‘‘observed’’ signal. FIGURE 9. The observed signal Vowhen it is composed of a primary signal Vpand an echoed signal Vewith normalized parameter τc0/w=2.0. need a sensor to be placed at a certain distance from the 472 screen. As the reflected field above the analyzed structure 473 can be attributed to the voltage response (that equals to the 474 100110 VOLUME 10, 2022 P. Kadlec: Time–Domain Electromagnetic Identification Based on Rectangular Grooves FIGURE 10. The observed DER values for different delay τof the echoed signal Ve. FIGURE 11. The electric field TD pulses observed at SP1 (xS/w=10.0, yS/w=20.0), and SP2 (xS/w=50.0, yS/w=50.0) for ID1 (solid curves) and ID2 (dashed curves). Data computed by CST Microwave Studio. equivalent magnetic-current surface density) it could be eval-475 uated through a space-time convolution with the pertaining476 2-D Green’s function [31]. This continuation of the voltage477 response to the upper space will then, indeed, change both478 the amplitude and shape of the TD response. Nevertheless,479 the same effect on the observed pulse can also be achieved by480 changing the properties of the initial pulse ei(t).481 The effect of the sensor lateral position xSand height482 ySabove the screen on the observed response is shown in483 Fig. 11. Here, the electric field TD pulses observed at two484 different positions of the sensor (SP1, SP2) for two different485 IDs (ID1, ID2, see Fig. 3) are compared. The TD pulses were486 obtained using the CST Microwave Studio. Especially the487 later parts of the pulses are well distinguishable for the same488 SP and different ID. The objective function margin between489 the two IDs 1is well below −10 dB for both locations490 (1= −6.95dB for SP1, and 1= −7.85 dB for SP2).491 As shown by the results presented earlier, this value of 1is492 sufficient for the optimization algorithms to securely identify 493 the correct ID. 494 V. CONCLUSION 495 This paper proposes the identification system based on a 496 solution to the grooves characterization problem. The indi497 vidual items are distinguished by a unique set of grooves 498 etched on their surface or an attached conducting tag. The 499 grooves have different depths and mutual positions. The char500 acterization problem is solved by employing two powerful 501 computational tools: a semi–analytical scheme based on the 502 Cagniard–DeHoop method of moments is used as the forward 503 solver and a stochastic optimization technique is used as the 504 inverse solver. Results of the comparative study show that 505 PSO is the most efficient optimization algorithm from the 506 chosen set of tested algorithms. Furthermore, we demon507 strated that the proposed identification approach works well 508 in the presence of noise signals. 509 The proposed system has several advantages compared 510 to traditional optical or RFID systems. There is no need to 511 design passive tags or bar codes. The grooves can be etched 512 by a widely–available CNC technology directly into the con513 ducting surface. The system can be operated even in places 514 with poor visibility. Next, the number of distinguishable 515 items grows exponentially with the number of used grooves. 516 The results show, that 640 items can be easily distinguished 517 using only two grooves with spatial step 1mm. On the other 518 hand, the number of grooves cannot grow to infinity because 519 the accuracy of the optimization algorithm would decrease 520 thanks to the curse of dimensionality. Also, the dimensions of 521 the grooves have to be selected with regard to the capabilities 522 of available signal generators. 523 ACKNOWLEDGMENT 524 The author would like to thank Dr. Martin Štumpf for provid525 ing the code serving as the forward solver, for valuable dis526 cussions of the problem, and proofreading of this manuscript. 527 REFERENCES 528 [1] M. Zhou, Q. Wang, T. Lei, Z. Wang, and K. Ren, ‘‘Enabling online robust 529 barcode-based visible light communication with realtime feedback,’’ IEEE 530 Trans. Wireless Commun., vol. 17, no. 12, pp. 8063–8076, Dec. 2018. 531 [2] A. Lozano-Nieto, RFID Design Fundamentals and Applications.532 Boca Raton, FL, USA: CRC Press, 2017. 533 [3] A. Haibi, K. Oufaska, K. E. Yassini, M. Boulmalf, and M. Bouya, ‘‘Sys534 tematic mapping study on RFID technology,’’ IEEE Access, vol. 10, 535 pp. 6363–6380, 2022. 536 [4] L. Shahid, H. Shahid, M. A. Riaz, S. I. Naqvi, M. J. Khan, M. S. Khan, 537 Y. Amin, and J. Loo, ‘‘Chipless RFID tag for touch event sensing and 538 localization,’’ IEEE Access, vol. 8, pp. 502–513, 2020. 539 [5] P. Fathi, N. C. Karmakar, M. Bhattacharya, and S. Bhattacharya, ‘‘Potential 540 chipless RFID sensors for food packaging applications: A review,’’ IEEE 541 Sensors J., vol. 20, no. 17, pp. 9618–9636, Sep. 2020. 542 [6] A. Subrahmannian and S. K. Behera, ‘‘Chipless RFID: A unique technol543 ogy for mankind,’’ IEEE J. Radio Freq. Identificat., vol. 6, pp. 151–163, 544 2022. 545 [7] P. Yan, S. Choudhury, and R. Wei, ‘‘A machine learning auxiliary approach 546 for the distributed dense RFID readers arrangement algorithm,’’ IEEE 547 Access, vol. 8, pp. 42270–42284, 2020. 548 [8] M. A. S. Tajin and K. R. Dandekar, ‘‘Pattern reconfigurable UHF RFID 549 reader antenna array,’’ IEEE Access, vol. 8, pp. 187365–187372, 2020. 550 VOLUME 10, 2022 100111 P. Kadlec: Time–Domain Electromagnetic Identification Based on Rectangular Grooves [9] A. Ghasempour, ‘‘Internet of Things in smart grid: Architecture, applica-551 tions, services, key technologies, and challenges,’’ Inventions, vol. 4, no. 1,552 p. 22, Mar. 2019.553 [10] G. Zhang, S. Tao, W. Xiao, Q. Cai, W. Gao, J. Jia, and J. Wen, ‘‘A fast and554 universal RFID tag anti-collision algorithm for the Internet of Things,’’555 IEEE Access, vol. 7, pp. 92365–92377, 2019.556 [11] A. Rivadeneyra, A. Albrecht, F. Moreno-Cruz, D. P. Morales, M. Becherer,557 and J. F. Salmeron, ‘‘Screen printed security-button for radio frequency558 identification tags,’’ IEEE Access, vol. 8, pp. 49224–49228, 2020.559 [12] L. Gao, L. Zhang, F. Lin, and M. Ma, ‘‘Secure RFID authentication560 schemes based on security analysis and improvements of the USI proto-561 col,’’ IEEE Access, vol. 7, pp. 8376–8384, 2019.562 [13] F. Zhu, ‘‘SecMAP: A secure RFID mutual authentication protocol for563 healthcare systems,’’ IEEE Access, vol. 8, pp. 192192–192205, 2020.564 [14] H. Wakaumi, T. Komaoka, and E. Hankui, ‘‘Grooved bar-code recogni-565 tion system with tape-automated-bonding head detection scanner,’’ IEEE566 Trans. Magn., vol. 36, no. 1, pp. 366–370, Jan. 2000.567 [15] K. Barkeshli and J. L. Volakis, ‘‘Scattering from narrow rectangular filled568 grooves,’’ IEEE Trans. Antennas Propag., vol. 39, no. 6, pp. 804–810,569 Jun. 1991.570 [16] D. D. Reuster and G. A. Thiele, ‘‘A field iterative method for computing571 the scattered electric fields at the apertures of large perfectly conducting572 cavities,’’ IEEE Trans. Antennas Propag., vol. 43, no. 3, pp. 286–290,573 Mar. 1995.574 [17] Y. Shifman and Y. Leviatan, ‘‘Scattering by a groove in a conducting575 plane–A PO-MoM hybrid formulation and wavelet analysis,’’ IEEE Trans.576 Antennas Propag., vol. 49, pp. 1807–1811, 2001.577 [18] E. Howe and A. Wood, ‘‘TE solutions of an integral equations method578 for electromagnetic scattering from a 2D cavity,’’ IEEE Antennas Wireless579 Propag. Lett., vol. 2, pp. 93–96, 2003.580 [19] M. Bozorgi, ‘‘A generalized method for scattering from wide cavities with581 specified wave functions,’’ IET Microw., Antennas Propag., vol. 15, no. 1,582 pp. 69–79, Jan. 2021.583 [20] A. Kumar Agrawal, A. Ninawe, and A. Dhawan, ‘‘Non-uniform narrow584 groove plasmonic nano-gratings for SPR sensing and imaging,’’ IEEE585 Access, vol. 9, pp. 10136–10152, 2021.586 [21] Z. Cao, F. Cui, F. Xian, S. Pei, J. Li, and H. Ye, ‘‘Toward characterization of587 a rectangular groove on a metallic surface by multi-angle light scattering,’’588 IEEE Access, vol. 8, pp. 60210–60217, 2020.589 [22] X. Meng, L. X. Guo, C. L. Dong, and Y. C. Jiao, ‘‘GO/PO method for590 the terahertz scattering computation of objects with multiple small-scale591 grooves,’’ IEEE Access, vol. 7, pp. 40738–40745, 2019.592 [23] M. Štumpf, ‘‘Time-domain electromagnetic scattering by a two-593 dimensional narrow groove—A solution based on the Cagniard–DeHoop594 method of moments,’’ IEEE Antennas Wireless Propag. Lett., vol. 21,595 no. 3, pp. 586–589, Mar. 2022.596 [24] M. Štumpf, Metasurface Electromagnetics: The Cagniard-DeHoop Time-597 Domain Approach (Electromagnetic Waves). Rijeka, Croatia: SciTech,598 2022.599 [25] P. Kadlec, M. Štumpf, and T. Doležal, ‘‘Towards time-domain characteriza-600 tion of a rectangular groove,’’ in Proc. Int. Symp. (ELMAR), 2022, pp. 1–4.601 [26] R. Fletcher, Practical Methods of Optimization. Hoboken, NJ, USA: Wiley, 602 2013. 603 [27] J. Robinson and Y. Rahmat-Samii, ‘‘Particle swarm optimization in elec604 tromagnetics,’’ IEEE Trans. Antennas Propag., vol. 52, no. 2, pp. 397–407, 605 Feb. 2004. 606 [28] K. Price, R. M. Storn, and J. A. Lampinen, Differential Evolution: A Prac607 tical Approach to Global Optimization. Berlin, Germany: Springer, 2006. 608 [29] P. Kadlec and Z. Raida, ‘‘Self-organizing migrating algorithm for opti609 mization with general number of objectives,’’ in Proc. 22nd Int. Conf. 610 Radioelektronika, 2012, pp. 1–5. 611 [30] O. Cosma, P. C. Pop, and I. Zelina, ‘‘An effective genetic algorithm for 612 solving the clustered shortest-path tree problem,’’ IEEE Access, vol. 9, 613 pp. 15570–15591, 2021. 614 [31] T. Adrianus and D. Hoop, Handbook of Radiation and Scattering of Waves.615 New York, NY, USA: Academic, 1995. 616 [32] M. Štumpf, Time-Domain Electromagnetic Reciprocity in Antenna Model617 ing. Hoboken, NJ, USA: Wiley, 2020. 618 [33] I. Oppermann, M. Hämäläinen, and J. Iinatti, UWB: Theory and Applica619 tions. Hoboken, NJ, USA: Wiley, 2004. 620 [34] D. H. Wolper and W. G. Macready, ‘‘No free lunch theorems for optimiza621 tion,’’ IEEE Trans. Evol. Comput., vol. 1, no. 1, pp. 67–82, Apr. 1997. 622 [35] J. Kennedy and R. Eberhart, ‘‘Particle swarm optimization,’’ in Proc. IEEE 623 ICNN, vol. 4. Nov./Dec. 1995, pp. 1942–1948. 624 [36] M. Marek, P. Kadlec, and M. Čapek, ‘‘FOPS: A new framework for the 625 optimization with variable number of dimensions,’’ Int. J. RF Microw. 626 Comput.-Aided Eng., vol. 30, no. 9, Sep. 2020, Art. no. e22335. 627 [37] M. Marek and P. Kadlec, ‘‘Another evolution of generalized differential 628 evolution: Variable number of dimensions,’’ Eng. Optim., vol. 54, no. 1, 629 pp. 61–80, Jan. 2022. 630 PETR KADLEC (Member, IEEE) received the 631 B.Sc., M.Sc., and Ph.D. degrees in electrical 632 engineering from the Brno University of Tech633 nology (BUT), Brno, Czech Republic, in 2007, 634 2009, and 2012, respectively. He is currently 635 a Researcher with the Lerch Laboratory of 636 EM Research, Department of Radio Electronics, 637 BUT. He has coauthored several computational 638 toolboxes, including the Antenna Toolbox for 639 MATLAB (AToM), and the Fast Optimization Pro640 cedureS (FOPS). His research interests include global optimization methods 641 and computational methods in electromagnetics. 642 643 100112 VOLUME 10, 2022