Wave-based sensor, actuator and optimizer

Universidade do Minho Escola de Engenharia Paulo Jorge da Silva Santos Wave-based sensor, actuator and optimizer outubro de 2023 UMinho | 2023 Paulo Jorge da silva santos Wave-based sensor, actuator and optimizer Paulo Jorge da Silva Santos Wave-based Sensor, Actuator and Optimizer Tese de Doutoramento Programa Doutoral em Sistemas Avançados de Engenharia para a Indústria (AESI) Trabalho efetuado sob a orientação do Universidade do Minho Escola de Engenharia outubro de 2023 Professor Doutor Jorge Miguel Nunes dos Santos Cabral Doutor Filipe Manuel Serra Alves DIREITOS DE AUTOR E CONDIÇÕES DE UTILIZAÇÃO DO TRABALHO POR TERCEIROS Este é um trabalho académico que pode ser utilizado por terceiros desde que respeitadas as regras e boas práticas internacionalmente aceites, no que concerne aos direitos de autor e direitos conexos. Assim, o presente trabalho pode ser utilizado nos termos previstos na licença abaixo indicada. Caso o utilizador necessite de permissão para poder fazer um uso do trabalho em condições não previstas no licenciamento indicado, deverá contactar o autor, através do RepositóriUM da Universidade do Minho. Licença concedida aos utilizadores deste trabalho Atribuição-NãoComercial-SemDerivações CC BY-NC-ND https://creativecommons.org/licenses/by-nc-nd/4.0/ ii Acknowledgements I would like to express my sincere and special thanks to my advisor and friend, Dr. Professor Jorge Cabral. His exceptional guidance, dedication, and wisdom have profoundly influenced my academic and personal development throughout this journey. His constant support and encouragement to challenge myself were fundamental to the successful completion of this thesis. I am deeply grateful. I would also like to express my deep gratitude to my other two supervisors, Dr. Filipe Serra alves and Dr. José Azevedo Gonçalves, for their support and guidance throughout my research journey. I would like to thank the Fundação para a Ciência e Tecnologia (FCT) and Bosch Car Multimedia for funding my PhD (grant PD/BDE/142901/2018). I would like to extend a special and heartfelt thanks to João Carvalho, who has been a friend and a constant source of support, encouragement, and inspiration throughout my thesis. His unwavering presence, insightful advice, and genuine interest in my academic and personal development were truly invaluable. I am deeply grateful. Also, I would like to extend my sincere thanks to my dear friends Miguel Esteves, Francine Oliveira, Rodrigo Marinho, Eduardo Mendes, and Luís Novais, for their unwavering support and friendship during this challenging journey. Your encouragement, understanding, and ongoing support were crucial in helping me navigate the ups and downs of this research endeavor. To my two brothers, Pedro Santos and Tiago Santos, thank you for always being by my side and providing me with the motivation and inspiration to keep going. Finally, I would like to thank my parents, Maria Silva and Álvaro Santos, for their unconditional love, unwavering support, and endless encouragement. Their sacrifices and belief in me were fundamental in making this thesis a reality. Paulo Jorge da Silva Santos October 16th, 2023 iii STATEMENT OF INTEGRITY I hereby declare having conducted this academic work with integrity. I confirm that I have not used plagiarism or any form of undue use of information or falsification of results along the process leading to its elaboration. I further declare that I have fully acknowledged the Code of Ethical Conduct of the University of Minho. iv Resumo A presente tese explora a utilização de ondas para abordar dois desafios significativos na indústria automóvel. O primeiro desafio consiste no desenvolvimento de um sistema de cancelamento ativo de ruído (ANC) que possa reduzir os ruídos não estacionários no compartimento de passageiros de um veículo. O segundo desafio é criar uma metodologia de conceção ótima para sensores de posição indutivos capazes de medir deslocamentos lineares, rotacionais e angulares. Para abordar o primeiro desafio, foi desenvolvido de um sistema ANC onde wavelets foram combinadas com um banco de filtros adaptativos. O sistema foi implementado em uma FPGA, e testes demonstraram que o sistema pode reduzir o ruído não estacionário em um ambiente acústico aberto e não controlado em 9 dB. O segundo desafio foi abordado através de uma metodologia que combina um algoritmo genético com um método numérico rápido para otimizar um sensor de posição indutivo. O método numérico foi usado para simular o campo eletromagnético associado à geometria do sensor, permitindo a maximização da corrente induzida nas bobinas recetoras e a minimização da não-linearidade no sensor. A minimização da não-linearidade foi conseguida através do desenho (layout) das bobinas que compõem o sensor. Sendo este otimizado no espaço de Fourier através da adição de harmónicos apropriados na geometria. As melhores geometrias otimizadas apresentaram uma não-linearidade inferior a 0,01% e a 0,25% da escala total para os sensores de posição angular e linear, respetivamente, sem calibração por software. O sistema ANC proposto tem o potencial de melhorar o conforto dos ocupantes do veículo, reduzindo o ruído indesejado dentro do compartimento de passageiros. Isso poderia reduzir o uso de materiais de isolamento acústico no veículo, levando a um veículo mais leve e, em última análise, a uma redução no consumo de energia. A metodologia desenvolvida para sensores de posição indutivos contribui para o estado da arte de sensores de posição eficientes e económicos, o que é crucial para os requisitos complexos da indústria automóvel. Essas contribuições têm implicações para o desenho de sistemas automotivos, com requisitos de desempenho e considerações ambientais e económicas. Palavras-chave: Cancelamento ativo de ruído, otimização, sensor de posição indutivo v Abstract This thesis explores the use of waves to tackle two major engineering challenges in the automotive industry. The first challenge is the development of an Active Noise Cancelling (ANC) system that can effectively reduce non-stationary noise inside a vehicle’s passenger compartment. The second challenge is the optimization of an inductive position sensor design methodology capable of measuring linear, rotational, and angular displacements. To address the first challenge, this work designs an ANC system that employs wavelets combined with a bank of adaptive filters. The system was implemented in an FPGA, and field tests demonstrate its ability to reduce non-stationary noise in an open and uncontrolled acoustic environment by 9 dB. The second challenge was tackled by proposing a new approach that combines a genetic algorithm with a fast and lightweight numerical method to optimize the geometry of an inductive position sensor. The numerical method is used to simulate the sensor’s electromagnetic field, allowing for the maximization of induced current on the receiver coils while minimizing the sensor’s non-linearity. The non-linearity minimization was achieved through its unique sensor’s coils design optimized in the Fourier space by adding the appropriate harmonics to the coils’ geometry. The best optimized geometries exhibited a non-linearity of less than 0.01% and 0.25% of the full scale for the angular and linear position sensors, respectively. Both results were achieved without the need for signal calibration or post-processing manipulation. The proposed ANC system has the potential to enhance the comfort of vehicle occupants by reducing unwanted noise inside the passenger compartment. Moreover, it has the potential to reduce the use of acoustic insulation materials in the vehicle, leading to a lighter vehicle and ultimately reducing energy consumption. The developed methodology for inductive position sensors represents a state-of-the-art contribution to efficient and cost-effective position sensor design, which is crucial for meeting the complex requirements of the automotive industry. Keywords: Active Noise Cancelling (ANC), inductive position sensor, optimization vi Table of Contents Resumo v Abstract vi List of Figures x List of Tables xv List of Algorithms xvii List of Acronyms xviii 1 Introduction 1 1.1 Motivation ..................................... 3 1.2 Scope ....................................... 4 1.3 Researchquestions................................. 5 1.4 Researchmethodology ............................... 5 1.5 Thesisstructure................................... 7 2 Theoretical context 9 2.1 Signalprocessing.................................. 9 2.1.1 SignalProperties.............................. 10 2.1.2 FourierTransform ............................. 12 2.1.3 WaveletTransform ............................. 16 2.2 Computational Electromagnetism . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.1 Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.2 Methodofmoments ............................ 29 vii 4.70 FMM vs MoM - Excitation coil current distribution at 1 GHz . . . . . . . . . . . . . . 162 5.1 Hardware accelerated ANC system board . . . . . . . . . . . . . . . . . . . . . . 166 5.2 Manufactured inductive position sensors’ PCB featuring the optimized geometries . . . 167 5.3 LIPSsensortestsetup ............................... 168 xiv List of Tables 3.1 Computing a very rough approximation to the inverse of a number by reversing the order ofthebits...................................... 76 3.2 FPGA resources utilization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.3 Results - average attenuation for each noise type . . . . . . . . . . . . . . . . . . . 86 4.1 Comparison of different inductive position sensors . . . . . . . . . . . . . . . . . . 105 4.2 Comparison of different inductive position sensors - results . . . . . . . . . . . . . . 105 4.3 Example of fixed and optimizable parameters for induced current maximization . . . . 111 4.4 Example of fixed and optimizable parameters for non-linearity minimization . . . . . . 112 4.5 APS optimization using MoM - Limits and initial parameters for induced current maximization ...................................... 124 4.6 APS optimization using MoM - Limits, initial and optimized parameters for non-linearity optimization (amplitudes) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 4.7 APS non-linearity optimization (amplitudes and phases) - Genetic Algorithm (GA) parameters........................................ 129 4.8 APS optimization using MoM - Limits, initial and optimized parameters for non-linearity optimization (amplitudes and phases) . . . . . . . . . . . . . . . . . . . . . . . . 130 4.9 APS optimization using MoM - Limits, initial and optimized parameters for harmonic contentoptimization ................................ 131 4.10 APS optimization using FHM - Limits and initial parameters for induced current maximization ...................................... 137 4.11 APS optimization using FHM - Limits, initial and optimized parameters for non-linearity optimization (amplitudes and phases) . . . . . . . . . . . . . . . . . . . . . . . . 139 4.12 LIPS optimization using FHM - Limits and initial parameters for induced current maximization ...................................... 143 xv 4.13 LIPS optimization using FHM - Limits, initial and optimizable parameters for non-linearity optimization (amplitudes and phases) . . . . . . . . . . . . . . . . . . . . . . . . 145 4.14 LIPS non-linearity optimization (amplitudes and phases) - GA parameters . . . . . . . 147 4.15 LIPS harmonic content optimization - GA parameters . . . . . . . . . . . . . . . . . 148 4.16 LIPS optimization using FHM - Limits, initial and optimizable parameters for harmonic contentoptimization ................................ 149 A.1 APS Geometric optimization - Error measurement minimization . . . . . . . . . . . . 184 A.2 APS mainresults.................................. 186 A.3 Signal components’ amplitudes of the APS sensors’ receiver coils . . . . . . . . . . 188 A.4 LIPS geometric optimization - Error measurement minimization . . . . . . . . . . . . 189 A.5 LIPS geometric optimization - Harmonic content minimization . . . . . . . . . . . . 192 A.6 LIPS mainresults.................................. 194 A.7 Signal components’ amplitudes of the LIPS sensors’ receiver coils . . . . . . . . . . 196 xvi List of Algorithms 1 MLFMMAlgorithm ................................. 47 2 Simulated Annealing Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3 Particle Swarm Optimization Algorithm . . . . . . . . . . . . . . . . . . . . . . . 53 4 GeneticAlgorithm.................................. 55 5 Tournament Selection Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6 CrossoverAlgorithm................................. 56 7 Gaussian Mutation Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 8 ElitismAlgorithm .................................. 57 xvii List of Acronyms ABS Antilock Braking System. AI Artificial Intelligence. AM Amplitude Modulation. AMR Anisotropic Magnetoresistance. ANC Active Noise Cancellation. APS Angular Position Sensor. BCLK Bit Clock. BRAM Block of Random-Access Memory. CEM Computational Electromagnetism. CIC Cascaded Integrator-Comb. DFT Discrete Fourier Transform. DSP Digital Signal Processor. DWT Dyadic Wavelet Transform. EA Evolutionary Algorithm. EMF Electromotive Force. ESP Electronic Stability Program. FDTD Finite Difference Time Domain. xviii FEM Finite Element Method. FFT Fast Fourier Transform. FHM Fast Harmonic Method. FIR Finite Impulse Response. FLANN Functional Link Artificial Neural Network. FMM Fast Multipole Method. FPGA Field Programmable Gate Array. FxLMS Filtered-x Least Mean Square. FxNLMS Filtered-x Normalized Least Mean Square. GA Genetic Algorithm. GMR Giant Magnetoresistance. GMRES Generalized Minimum Residue. I2SInter-IC Sound. IDFT Inverse Discrete Fourier Transform. IIR Infinite Impulse Response. LIPS Linear Position Sensor. LMS Least Mean Square. LRCLK Left/Right Clock. LU Lower–Upper. MEMS Microelectromechanical System. MIMO Multiple-Input Multiple-Output. MLFMM Multilevel Fast Multipole Method. xix MoM Method of Moments. NLMS Normalized Least Mean Square. PCB Printed Circuit Board. PCM Pulse Code Modulation. PDM Pulse Density Modulation. PL Programmable Logic. PS Processing System. PSO Particle Swarm Optimization. RAM Random-Access Memory. RCGA Real-Coded Genetic Algorithm. SA Simulated Annealing. SD Serial Data. SISO Single-Input Single-Output. SNR Signal-to-Noise Ratio. TCS Traction Control System. TDM Time-Division Multiplexing. TMR Tunnel Magnetoresistance. UWPT Undecimated Wavelet Packet Transform. WPT Wavelet Packet Transform. WS Word Select. xx Chapter 1 Introduction Waves are a fundamental aspect of the natural world and play a crucial role in our ability to sense our surroundings. We see and hear because of the existence of waves, electromagnetic and mechanical. Electromagnetic waves do not require a medium to propagate as they can travel through a vacuum, such as space. They range in frequency from radio waves through infrared and visible light to X-rays and gamma rays. Unlike electromagnetic waves, mechanical waves, such as sound waves, require a medium to propagate. They are created by the vibration of an object, such as a vocal cord, guitar string, or drumskins, which creates a pressure change or disturbance in the medium through which they travel. This pressure change travels through the medium as a wave, and we perceive these waves as sound. We can hear sound waves propagating through the air (or water, if we are in a dive) from around 20 Hz to 20 kHz. A sensor is a device designed to detect or measure physical properties or changes in the environment, converting this information into a signal that can be analyzed or monitored [1]. Essentially, sensors serve as an interface between the physical world and a system that processes the data they gather. In particular, sensors can detect or measure several types of waves, such as acoustic or radio waves. They play a critical role in numerous fields, including engineering and physics, as they enable the measurement of diverse physical quantities. For instance, a microphone, a specific type of sensor, transforms acoustic waves into electrical signals, allowing for the measurement of acoustic pressure—even those outside the audible range, such as ultrasound and infrasound. Moreover, electromagnetic or mechanical waves can be indirectly used to measure other physical quantities, such as position, distance, and velocity. As demonstrated in this thesis, linear or angular positions can be indirectly measured through electromagnetic waves. A coil generates an electromagnetic wave, a conductive target modulates it based on its 1 2 Chapter 1. Introduction position, and another coil subsequently detects it. An actuator is a component or device that converts input energy, such as electrical, hydraulic, or pneumatic, into mechanical motion or action, thereby facilitating control and movement in systems or processes [1]. Actuators perform the opposite function of sensors by typically converting a signal, often electrical, into a physical quantity like generating an electromagnetic or mechanical wave. A common example is a speaker, which transforms an electrical signal, typically an audio signal, into an acoustic wave that we perceive as sound. Both sensors and actuators are widely used in machinery, robotics, automation, and automotive systems to execute specific actions or tasks, and can be considered transducers since they generally convert energy from one form to another. Sensors and actuators are essential components in the rapidly-evolving automotive industry. As new challenges arise, such as electric and autonomous vehicles, increased environmental concerns, and the need for more efficient and reliable technology, there is a growing demand for high-performance, lowcost sensors and actuators. An emerging technology in the automotive industry that uses sensors and actuators is Active Noise Cancellation (ANC). ANC systems use sensors (microphones) and actuators (loudspeakers) to replace bulky and heavy noise-reduction materials, thus reducing consumption and weight, and ultimately resulting in improved passenger comfort. This way, passengers can enjoy a quieter and more comfortable ride in their electric and autonomous vehicles. Another example can be seen in the evolving positioning sensor technologies. Traditional permanent magnet-based sensors, composed of rare raw materials, are replaced by more economical and environmentally friendly inductive technologies immune to constant magnetic fields. This is particularly important for electric cars due to the strong magnetic fields caused by high currents. Overall, the automotive industry is constantly evolving, and sensors and actuators play a vital role in meeting new challenges and demands. They are essential components that enable the development of electric and autonomous vehicles, as well as improve the efficiency and performance of traditional vehicles. As technology advances, we can expect to see even more innovative uses of sensors and actuators in the automotive industry, resulting in improved safety, performance, and environmental sustainability. Chapter 1. Introduction 3 1.1 Motivation The University of Minho and Bosch Car Multimedia Braga established a research and technology development collaboration back in 2013. During the first phase (2012-2015), HMIExcel (Human Machine Interface Excellence) was conducted, focusing on developing multimedia solutions for automobiles in the man-machine interface. Building on this success, in 2015, the University of Minho and Bosch Car Multimedia Portugal launched a new initiative called INNOVATIVE Car HMI , which aimed to generate knowledge and technologies that transfer into global advancements in future car solutions. In July 2018, the partnership moved into its third phase with the Sensible Car Program, further concentrating on developing smart critical sensors to meet the capabilities required of autonomous vehicles. Driven by their impressive track record in innovative automotive technology and research, I opted to work in this partnership with Bosch Car Multimedia and the University of Minho. This collaboration allowed me to contribute to cutting-edge advancements in the field while benefiting from their extensive knowledge and resources, and provided the opportunity to work alongside industry experts and researchers, fostering a stimulating environment for learning and professional growth. As a result of my involvement in the partnership, I explored two emerging technologies in the automotive industry. The first was the ANC system, which enhances vehicle occupant comfort by reducing unwanted noise inside the passenger compartment. This technology also has the potential to reduce the use of acoustic insulation materials, leading to lighter vehicles and ultimately reducing energy consumption. However, current ANC systems in the automotive industry mainly focus on reducing stationary noise, leaving non-stationary noise unaddressed. This gap in the technology motivated research into finding possible solutions using spectral methods, making the development of an ANC system for reducing nonstationary noise one of this thesis’s main goals. The other main goal was the inductive position sensor technology based on planar coils printed directly on Printed Circuit Boards (PCBs), which has been gaining attention in the research community in recent years. These sensors offer several advantages, most notably their robustness in the presence of an external continuous magnetic field, crucial for applications in electric vehicles, among others. Furthermore, their production is relatively low-cost, as the coils can be printed directly onto the PCB, eliminating the need for magnets. Despite these strengths, the performance of these sensors depends on the coil geometry, and the current development process often involves time-consuming trial-and-error methods. To address 10 Chapter 2. Theoretical context in analyzing signals that are not stationary or that change over time, with applications such as image compression, noise reduction, and feature extraction [4]. Fourier analysis and wavelets are widely-used and essential mathematical tools in signal processing. They are employed in various fields to extract useful information from signals. Both have advantages and disadvantages and are often used together for comprehensive signal analysis. It is important to note that signal processing is a vast and complex field with many different tools that can be used to extract useful information from signals. Fourier analysis and wavelets are just two examples of the many tools used in signal processing. Understanding the basics of these two tools provides a foundation for understanding decisions and implementations throughout this thesis. 2.1.1 Signal Properties In this section, some of the properties of signals are introduced. These properties are an essential base for the discussion the Fourier integrals further ahead in the chapter. 2.1.1.1 Linearity An operator, L, is said to be linear if the superposition theorem applies; that is, if the input, f(t), is a weighted sum of multiple signals, the output, Lf (t), will also be defined by a weighted sum of responses to each of the inputs [2, 3]. Analytically, if Lf1(t)is the output of f1(t), and Lf2(t)is the output of f2(t), then: af1(t) + bf2(t)⇒aLf1(t) + bLf2(t)(2.1) where in the above equation, aand bare constants that can be complex numbers. 2.1.1.2 Time invariance Linear time-invariant operators are used to implement signal processing operations like signal transmission, stationary denoising, and predictive coding. An operator is said to be time-invariant if a time shift in the input f(t)results in the same exact time shift in the output signal [2,3]. Chapter 2. Theoretical context 11 g(t) = Lf (t)⇒g(t−τ) = Lf (t−τ)(2.2) 2.1.1.3 Impulse response Linear time-invariant systems can be characterized entirely by their response to a Dirac impulse. If f(t) is continuous, its value at tcan be given by Equation 2.3, where δ(t−u)is the Dirac impulse located at time u[2,3]. f(t) = ˆ+∞ −∞ f(u)δ(t−u)du (2.3) Because of L’s continuity and linearity, the system’s output, Lf (t), is the convolution of the input of the system, f(t), with the system’s impulse response, Equation 2.4, where his the impulse response of L (h(t) = Lδ (t)). Lf (t) = ˆ+∞ −∞ f(u)h(t−u)du (2.4) The time-invariance demonstrates that Lδ (t−u) = h(t−u), resulting in: Lf (t) = ˆ+∞ −∞ f(u)h(t−u)du =ˆ+∞ −∞ h(u)f(t−u)du =h(t)∗f(t)(2.5) A linear time-invariant filter is comparable to a convolution with the impulse response h(t), not requiring that f(t)be continuous. This formula is true for every signal f(t)for which the convolution integral converges [4]. 12 Chapter 2. Theoretical context 2.1.1.4 Causality A filter is said to be causal if its output depends only on the current and past input values. If the filter depends on future values, it is non-causal [2, 3]. Equation 2.5, represents a convolution between an input signal, f(t), with an impulse response, h(t). In this example if Lf (t)does not depend on the values f(u)for u>t, meaning that h(u)=0for u < 0. This impulse response, h(t), is said to be causal [4]. 2.1.1.5 Stability A stable operator, L, is one in which the output, Lf (t), does not diverge for small inputs of f(t), or more specifically, one whose output’s amplitude is limited for limited input’s amplitude. Basically, the stability property ensures that Lf (t)is bounded if f(t)is bounded [4]. Since: |Lf (t)| ≤ ˆ+∞ −∞ |h(u)||f(t−u)|du ≤sup u∈R|f(u)|ˆ+∞ −∞ |h(u)|du (2.6) It is sufficient that ´+∞ −∞ |h(u)|du < +∞, even if h(u)is an unbound function. Thus, h(u)is said to be stable if it is absolutely integrable [4]. 2.1.2 Fourier Transform In this subsection, to prevent problems with convergence, the Fourier integral is first defined in space L1(R). Then, it is expanded to the space L2(R). 2.1.2.1 Transfer functions Eigenvectors of convolution operators are complex exponentials ejωt. Lejωt =ˆ+∞ −∞ h(u)ejω(t−u)du =ejωt ˆ+∞ −∞ h(u)e−jωudu =H(ω)ejωt (2.7) Chapter 2. Theoretical context 13 Where the eigenvalue: H(ω) = ˆ+∞ −∞ h(u)e−jωudu (2.8) is the Fourier transform of h(u)[4]. Because complex sinusoidal waves, ejωt, are eigenvectors of timeinvariant linear operators, it is tempting to decompose any function f(t)as a sum of these eigenvectors. The eigenvalues H(ω)can then be used to express Lf (t)[4]. 2.1.2.2 Fourier transform in L1(R) The Fourier integral, Equation 2.9, is an integral that measures how much of a particular frequency ωis present in a given function f(t). If f(t)∈L1(R), then the integral converges, and f(t)is an absolutely integrable signal, Equation 2.11 [4]. Similarly, it may be possible to reconstruct the original signal f(t) by applying the Inverse Fourier Transform, Equation 2.10. If f(t)∈L1(R)and F(ω)∈L1(R) then one can apply the Inverse Fourier Transform to recover the original signal f(t)[4]. For example, if f(t) = e−|t|is absolutely integrable (Equation 2.11), then F(ω) = 2/ (1 + ω2)is also absolutely integrable, so the Inverse Fourier Transform (Equation 2.10) works. However, if f(t) = δ(t), is absolutely integrable as well, then F(ω)is equal to 1 in the entire domain, so it is not absolutely integrable, and the Inverse Fourier Transform fails to recover the original signal. F(ω) = ˆ+∞ −∞ f(t)e−jωtdt (2.9) f(t) = 1 2πˆ+∞ −∞ F(ω)ejωtdω (2.10) |F(ω)| ≤ ˆ+∞ −∞ |f(t)|dt < +∞(2.11) The Inverse Fourier Transform, Equation 2.10, recovers f(t)by integrating an infinite set of sinusoids ejωt against the Fourier Transform of the function, F(ω). By doing so, one can prove that the hypothesis F(ω)∈L1(R)implies that f(t)is continuous [4]. Therefore, the reconstruction (Equation 2.10) is not 14 Chapter 2. Theoretical context proved for discontinuous functions. To extend the Fourier Transform to hold any function f(t)∈L2(R), it is necessary to extend the definition of the Fourier transform from L1(R)to L2(R)[4]. One of the most important properties of the Fourier Transform is that for any two functions f(t)and h(t) that are elements of L1(R), the function g(t)resulting from the convolution f(t)∗h(t)is also an element of L1(R)[4]. This can also be expressed as complex sinusoidal waves, ejωt, being eigenvalues of convolution operators. Furthermore, the response Lf (t) = g(t) = f(t)∗h(t)of a linear timeinvariant system can be calculated from its Fourier Transform G(ω) = F(ω)H(ω). By converting G(ω)into the original domain using the Inverse Fourier Transform, Equation 2.10, one can obtain the Lf (t). This implies that each frequency component ejωt of amplitude and phase F(ω)is amplified or attenuated, and phase shifted, by H(ω). As a result, such a convolution is known as frequency filtering , and the filter’s transfer function is H(ω)[4]. 2.1.2.3 Fourier transform in L2(R) The Fourier transform of the indicator function f= 1[−T,T ]that has discontinuities at t=±T, therefore, its Fourier Transform is not integrable: F(ω) = ˆT −T e−jωtdt = 2sin (Tω) ω(2.12) Thus, the Inverse Fourier Transform, Equation 2.10, does not apply in this case. This prompts the expansion of the Fourier Transform to the space L2(R)of functions f(t)with finite energy, as shown in Equation 2.13. By doing so, one can access all of the resources offered by the presence of an inner product in the Hilbert space L2(R)[4]. The inner product of f(t)∈L2(R)and g(t)∈L2(R)is given by Equation 2.14, where ∗denotes complex conjugate, and the norm of f(t)∈L2(R)is given by Equation 2.15. |F(ω)|2≤ˆ+∞ −∞ |f(t)|2dt < +∞(2.13) ⟨f, g⟩=ˆ+∞ −∞ f(t)g∗(t)dt (2.14) Chapter 2. Theoretical context 15 ∥f∥2=⟨f, f⟩=ˆ+∞ −∞ |f(t)|2dt (2.15) If f(t)and h(t)are in L1(R)∩L2(R), then: ˆ+∞ −∞ f(t)h∗(t)dt =1 2πˆ+∞ −∞ F(ω)H∗(ω)dω (2.16) Also, the resulting norm of f(t)∈L1(R)∩L2(R)is: ˆ+∞ −∞ |f(t)|2dt =1 2πˆ+∞ −∞ |F(ω)|2dω (2.17) Equations 2.16 and 2.17 are called the Parseval and Plancherel formulas, respectively [4]. 2.1.2.4 Convergence In the last two sections, the fundamental concepts have been described, enabling the determination of the class of signals for which the Fourier integrals, Equations 2.9 and 2.10, converge. Additionally, this ensures that ˆ f(t), the result of applying the Inverse Fourier Transform to F(ω), is a valid representation of f(t). For ˆ f(t)to be a valid representation of f(t),f(t)must be either square integrable (Equation 2.13) or absolutely integrable (Equation 2.11). In the latter case, f(t)must also be continuous. 2.1.2.5 Regularity and decay |F(ω)|decay is determined by f(t)’s worst single behaviour. For instance, because f(t) = 1[−T,T ]is discontinuous at t=±T,F(ω)decays like |ω|−1. In this case, it is also worth mentioning that f(t) is a normal function for t=±T. The decay of |F(ω)|does not provide this information. Therefore, to determine the local regularity of a signal f(t), it is necessary to decompose it using time-localized waveforms rather than sinusoidal waves ejωt [4]. 2.1.2.6 Uncertainty principle Due to the Heisenberg uncertainty principle, it is well known that it is impossible to know a particle’s position and momentum simultaneously with arbitrary precision [4]. One can only measure one of them. 16 Chapter 2. Theoretical context Similarly, the uncertainty principle does not permit the existence of a function f(t)with highly localized energy in time that has a Fourier Transform F(ω)with energy concentrated in a low-frequency interval. For example, the Dirac δ(t−u)has a support limited to t=u, but its Fourier Transform e−jωu has an energy distribution that is uniform across all frequencies. It is known that only if f(t)has regular temporal variations does |F(ω)|decay fast at high frequencies. As a result, the energy of f(t)must be dispersed over a relatively broad domain. In order to reduce the time spread of f(t), it can be scaled by s < 1while maintaining the overall energy unchanged [4]: fs(t) = 1 √sft s(2.18) The Fourier transform Fs(ω) = √sF (sω)is dilated by 1/s, meaning that what was gained in time localization was lost in frequency localization. It can be seen as a trade-off between time and frequency localization [4]. The Heisenberg uncertainty principle restricts time and frequency energy concentrations [4]. 2.1.2.7 Compact support If a signal, f(t), which is not identically null has a compact support (is zero outside a finite interval), its Fourier Transform, F(ω), cannot be zero on an entire interval. Again, if F(ω), which is not identically null, has a compact support, then f(t)cannot be zero on a whole interval, which shows that a function of compact support, whose Fourier Transform has a compact support, does not exist [4]. It is impossible to know an event’s time and frequency simultaneously with arbitrary precision. 2.1.3 Wavelet Transform The wavelet transform is a signal processing technique that computes the signal spectrum over time, like the windowed Fourier transform. However, the wavelet transform uses a compact support basis function called wavelet. This compact support makes them particularly suitable for handling non-stationary signals. There are several wavelet-based functions, such as Haar (see Equation 2.19) and Daubechies [5]. The Haar wavelet is also known as a particular case of the Daubechies wavelet and is the simplest possible wavelet. Chapter 2. Theoretical context 17 Figure 2.1: Translation-scale boxes in Gabor space Ψ(t) =              1if 0≤t < 1/2 −1if 1/2 ≤t < 1 0otherwise (2.19) The wavelet function (Ψ) is scaled and translated (represented as s and u in Equation 2.20) before being convoluted with the signal of interest f(t)(see Equation 2.21). This translation and scaling in wavelet transform plays similar roles to time and frequency in the windowed Fourier transform. Unlike the windowed Fourier transform that has the same resolution for all time-frequency boxes, the wavelet transform has different time-frequency (translation-scale) resolutions (Figure 2.1), making it more suitable to analyze the frequency-temporal dynamics of the signal of interest. Ψu,s(t) = 1 √sΨt−u s(2.20) Wf(u, s) = ˆ+∞ −∞ f(t)Ψ∗ u,s(t)dt (2.21) The wavelet transform can be computed using a bank of filters. The filter banks’ technique [6–8], also known as the fast wavelet transform algorithm, decomposes a discrete signal (e.g., f[n] ) into two half-sized signals using a filtering and subsampling procedure. This procedure is cascaded as shown in Figure 2.2, 18 Chapter 2. Theoretical context Figure 2.2: Dyadic wavelet transform Figure 2.3: Packet wavelet transform where his a low-pass filter and gis a high-pass filter, d0is the wavelet approximation coefficient, and d1 is the wavelet detail coefficient. Figure 2.2 illustrates the Dyadic Wavelet Transform (DWT) in which the wavelet approximation coefficient is cascaded as opposed to the wavelet detail coefficient. This wavelet transform results in translation-scale boxes illustrated in Figure 2.1 [4] (higher temporal resolution at high frequencies , but with low spectral resolution; and low temporal resolution at low frequencies , but with high spectral resolution). Figure 2.3 depicts a Wavelet Packet Transform (WPT) where both the approximation coefficient and the detail coefficient of the wavelet are decomposed, resulting in translation-scale boxes with the same resolution. There is a variant of the WPT called Undecimated Wavelet Packet Transform (UWPT) [9]. As the name suggests, this approach does not perform decimation after filtering; instead, the wavelet’s filter coefficients are upsampled. Using this technique, the wavelet coefficients are doubled at each level of decomposition, in contrast to the WPT, where the number of wavelet coefficients is equal to the number of samples of the input signal. Figure 2.4 illustrates the UWPT where the dilated filters, hand g, are upsampled. Figure 2.4: Undecimated wavelet packet transform Chapter 2. Theoretical context 19 Equations 2.22 and 2.23 represents the filter upsampling process by inserting zeros between each pair of coefficients of h(l)and g(l), where j is the decomposition level, and L is the total number of coefficients of the original filters h(l)and g(l). hj[l] =  h[0],0, ..., 0 |{z} 2j−1 , h[1],0, ..., 0 |{z} 2j−1 , ..., h[L−1],0, ..., 0 |{z} 2j−1  (2.22) gj[l] =  g[0],0, ..., 0 |{z} 2j−1 , g[1],0, ..., 0 |{z} 2j−1 , ..., g[L−1],0, ..., 0 |{z} 2j−1  (2.23) 2.2 Computational Electromagnetism In this section, the exploration of fundamental principles of Computational electromagnetics (CEM) is undertaken, with a focus on designing optimal inductive position sensors. The goal is to develop a customized simulation tool for streamlining the optimization process. The insights provided are essential for creating an efficient tool that can perform thousands of electromagnetic simulations within a reasonable timeframe. CEM is a technology that applies numerical methods in the field of electromagnetics, allowing for the modeling of electromagnetic field interactions with matter using Maxwell’s equations [10,11]. Numerous numerical methods have been widely adopted in CEM, such as the Finite Difference Time Domain (FDTD) method [12–14], the Finite Element Method (FEM), the Method of Moments (MoM), and the Fast Multipole Method (FMM) [15]. The FDTD and FEM methods are based on differential equations, while MoM and FMM rely on integral equations. To introduce CEM, one first explores Maxwell’s equations, which consist of four fundamental equations describing the behavior of electric and magnetic fields. These equations serve as the foundation of classical electromagnetism and offer a mathematical framework for comprehending various electromagnetic phenomena. Subsequently, the MoM and FMM methods are presented among the four mentioned methods, as they offer greater advantages in terms of speed and computational resources for the needs of the inductive position sensors optimization process. 26 Chapter 2. Theoretical context ∇2 E+k2 E=−jωµ  J(2.48) ∇2 B+k2 B=−µ∇×  J(2.49) This is done by applying the vector calculus identity 2.50. In the case of Equation 2.48, it was assumed that there are no charges (∇·  E= 0) [16–18]. ∇×∇×  F=∇∇·  F−∇2 F(2.50) 2.2.1.6 Scalar and vector potentials Maxwell’s equations provide all the knowledge required for fully identifying an electromagnetic phenomenon, but they can offer difficulties when the solutions are more complex. To make it easier to find a solution, auxiliary functions known as potentials (electric scaler potential and magnetic vector potential) are defined, simplifying the search of a solution in some situations [16–18]. ∇·∇×  F= 0 (2.51) ∇×(∇ϕ) = 0 (2.52) The magnetic vector potential,  A, can be introduced through Maxwell’s Equations 2.37. By considering the vector calculus identity 2.51, the following is obtained:  B=∇×  A(2.53) Similarly, via Maxwell’s Equation 2.39 and applying the vector calculus identity 2.52, the function of the electric scaler potential, ϕ, is presented as:  E=−∇ϕ−jω  A(2.54) Chapter 2. Theoretical context 27 Magnetic vector potential (  A) and electric scaler potential (ϕ) are not well-defined functions [16–18]. This is demonstrated by introducing an arbitrary function, ψ, to create a new pair of potentials,  A′, ϕ′, where:  A′= A−∇ψ(2.55) ϕ′=ϕ+jωψ (2.56) By applying Equations 2.53 and 2.54 to the new potentials defined by Equations 2.55 and 2.56, and considering the vector calculus identity 2.52, the following is obtained:  B=∇×  A′=∇×  A−∇×∇ψ=∇×  A(2.57)  E=−jω  A′−∇ϕ′=−jω  A+jω∇ψ−∇ϕ−jω∇ψ=−jω  A−∇ϕ(2.58) The differential equations that must be satisfied for the potentials in isotropic media are direct consequences of Maxwell Equations 2.36 and 2.40. By applying Equations 2.53 and 2.54, along with the constitutive relations, Equations 2.43 and 2.44, the following is obtained: ∇×1 µ∇×  A=ω2ϵ A−jωϵ∇ϕ+ J(2.59) jω∇·ϵ A+∇·(ϵ∇ϕ) = −ρ(2.60) For isotropic and homogeneous media, it can be further simplified by applying the vector calculus identity 2.50 [16–18]: ∇2 A+k2 A−∇∇·  A+jωµϵϕ=−µ J(2.61) ∇2ϕ+jω∇·  A=−ρ ϵ(2.62) 28 Chapter 2. Theoretical context Gauge invariance is a concept that applies to many areas of physics, and several authors have proposed restrictions on this degree of freedom [16]. In this work, the Lorentz gauge condition has been utilized: ∇·  A+jωµϵϕ = 0 (2.63) Imposing the Lorentz gauge condition (Equation 2.63), the differential equations for the potentials, Equations 2.61 and 2.62, are simple: ∇2 A+k2 A=−µ J(2.64) ∇2ϕ+k2ϕ=−ρ ϵ(2.65) 2.2.1.7 Boundary conditions For a given problem, there are many solutions that satisfy Maxwell’s equations. However, just one of them is the solution to the real-world problem. To obtain it, in addition to knowledge of differential equations, it is required to know the boundary conditions associated with the domain. Those conditions can be deduced from the integral form of Maxwell’s equations (Equations 2.27 to 2.30). Where in its general form at the interface between two media, given by: ˆn× E1− E2= 0 (2.66) ˆn· D1− D2=ρs(2.67) ˆn× H1− H2= Js(2.68) ˆn· B1− B2= 0 (2.69) where ˆnis the unit vector normal to the interface, as depicted in Figure 2.7. Also, only two equations, one from Equations 2.66 and 2.69, and the other from Equations 2.67 and 2.68, are independent. The boundary conditions can be simplified in a particular case where one of the media is a perfect conductor [17]. Since an ideal conductor cannot sustain internal fields, Equations 2.66 and 2.69 becomes: Chapter 2. Theoretical context 29 Medium 2 n Medium 1 μ1,ε1 μ2,ε2 Figure 2.7: Interface between two media ˆn× E= 0 (2.70) ˆn· B= 0 (2.71) where  Eand  Bare the outside fields to the conductor, and ˆnpoints away from the conductor [17]. 2.2.2 Method of moments MoM [19] was designed to solve the integral field equations, and it is frequently employed to analyze antenna devices because the antennas are typically thin wires and thin sheets. The MoM discretizes the 3D devices by decomposing them into a 2D surface mesh and, similarly, decomposes the thin wire devices into a straight 1D segment mesh. This discretization and the open domain nature of the MoM (meaning that the specification of boundary conditions is not required) are the main reasons for being applied in the antenna domain. The basis and test functions are then applied to each element [20,21]. The choice of these functions is central, as errors can occur due to this decision [22]. Pulse, triangular, and sinusoidal functions are the most commonly used basis and test functions for thin wires. Green’s functions, the core of MoM, are the impulse response of an inhomogeneous linear differential operator, which aids in analyzing the current on each element and the strength of each moment. Typically, there are many issues associated with MoM, such as low-frequency breakdown and singularities [23]. However, by employing singularity extraction, the accuracy of calculating Green’s function integral is enhanced. Also, MoM discretization results in a large and dense matrix, making it time-consuming and unsuitable for significant problems. Nevertheless, this could be overcome by employing preconditioners and iterative approaches [24]. In this thesis MoM was used to model thin-wire structures with speed, precision, stability, and versatility by using triangular functions as basis and test functions [10]. 30 Chapter 2. Theoretical context The MoM has several strengths over other methods, namely FDTD and FEM, which are methods based on differential equations. The main strengths are: i) the open domain, meaning that it does not require boundary conditions in the domain edge in contrast with the FDTD and FEM, which requires it to be closed domain methods; ii) the discretization of thin wires devices with unidimensional segments, which in contrast with FDTD and FEM that require the discretization of the entire volume where the device is in; iii) the memory usage that is much less than in FDTD and FEM; and iv) the processing time that is also shortened compared to the processing time needed by FDTD and FEM. The derivation of the terms used in the MoM to form the matrix equation (manipulated to solve the thinwire-based inductive sensors) starts with the imposition of the Lorenz gauge in the potential formulation of the electromagnetic field equations. In which the imposed Lorenz gauge, Equation 2.72, yields Equation 2.73. Where, ∇is the vector differential operator, j2=−1,ωis the angular frequency, µis the magnetic permeability, ϵis the electric permittivity, ϕis the electric scalar potential,  Ais the magnetic vector potential,  Jis the current density, and kis the wavenumber. ∇·  A=−jωµϵϕ (2.72) ∇2 A+k2 A=−µ J(2.73) In the context of Green’s functions, one can solve the non-homogeneous linear equation, Equation 2.73, which can be represented generically by Equation 2.74. In this equation, Lis a linear operator, Φis an unknown function, and gis a force function. As a result, for Equation 2.74, a solution can be expected in the form of Equation 2.75, where the core of this integral operator, G, is the Green’s function. In this case, Green’s function has the form of Equation 2.76 due to the Loperator being ∇2+k2and because it is a three-dimensional problem. Moreover, the integral operator in Equation 2.75 can be seen as the sum of the influences produced by all sources (at r′positions) in the Ωdomain at the rposition. In the case of Equation 2.73, where the unknown function Φis the magnetic vector potential  Aand the known force function gis µ J, making the respective substitutions in Equation 2.75 leads to Equation 2.77 [19]. Furthermore, because µ, in most cases, defines the permeability of a linear, isotropic, and homogeneous medium that a scalar can represent, µcan be extracted from the integral, as shown in Equation 2.77. Chapter 2. Theoretical context 31 (a) Curved thin wire (b) Discretization Figure 2.8: Curved thin wire discretization LΦ = g(2.74) Φ(r) = ˆΩ G(r, r′)g(r′)dr′(2.75) G(r, r′) = e−jk|r−r′| 4π|r −r′|(2.76)  A(r) = −µˆΩ G(r, r′) J(r′)dr′(2.77) Using the potential formulation of the electromagnetic field equations, the electric field  Ecan be obtained through Equation 2.78a. By rearranging this equation, Equation 2.78a, the term −jωµϵϕ can be replaced with ∇·  Afrom Equation 2.72, which leads to Equation 2.78b. Equation 2.78c is obtained by replacing the magnetic vector potential  Awith the right-hand side of Equation 2.77.  E=−jω  A−∇ϕ(2.78a) =−j ωµϵ k2−∇∇· A(2.78b) =j ωϵ k2−∇∇·ˆΩ G(r, r′) J(r′)dr′(2.78c) 32 Chapter 2. Theoretical context Figure 2.9: Triangle functions with half triangle on first and last segment Figure 2.8a depicts a thin curved wire, and Figure 2.8b shows this thin curved wire subdivided into N segments with N+ 1 endpoints. Although the segments, in this example, are the same length, they can be of different lengths. However, because the tangent vector is a piecewise continuous function, the segments must be small enough not to disturb (distort) the curvature of the wire. Since the wire is very thin, it can be assumed that the current  Jis a filamentary current  I. This current  Iis expanded using base functions (un), Equation 2.79, where Inare the unknown coefficients and unis an arbitrary basis, but most frequently is a pulse, triangular or sinusoidal function. Among the three, the triangular function, Equation 2.80, is the most used, being depicted in Figure 2.9 [19]. As represented in Figure 2.9, the first and last segments are assigned a ”half triangle” , allowing the solution to take any value at the ends. Otherwise, if a ”full triangle” is assigned, the solution will be zero at the ends, which is the case of an antenna.  I(r) = N X n=1 Inun(r)(2.79) un(r) =                r −rn |rn+1 −rn|, r ∈rn,rn+1 rn+2 −r |rn+2 −rn+1|,r ∈rn+1, rn+2 (2.80) In the inductive sensor design, the excitation coil is sourced by a delta-gap source, as shown in Figure 2.10a [19]. The delta-gap source treats the source as if the field exists only between the coil terminals, having zero value outside. Whenever it is assumed that the field exists on a single wire segment and is zero on the others, it means that the resulting excitation vector will only have non-zero elements for basis functions supported on that segment. This sourcing method typically yields less accurate results for input impedance simulation, but it does very well at simulating radiation patterns, which is one of the reasons Chapter 2. Theoretical context 33 (a) Delta-Gap (b) Frill Figure 2.10: MoM’s thin wire feeding techniques for its selection. The other reasons are its implementation simplicity and reduced computational cost, when compared to other techniques, namely frill (Figure 2.10b [19]). Equation 2.81a is obtained by applying thin-wire discretization to Equation 2.78c. This is accomplished by substituting the current density  J, in Equation 2.78c, with the discretized filamentary current  Iin Equation 2.79. Also, since Inis independent of the integration variable r′and the order of operations can be changed, PN n=1 Incan be moved outside the integral, leading to Equation 2.81b (or Equation 2.81c, where gmis given by Equation 2.82).  E(rm)≈j ωϵ k2−∇∇·ˆ∆r′ n G(rm, r′) N X n=1 Inun(r′)dr′(2.81a) = N X n=1 In j ωϵ k2−∇∇·ˆ∆r′ n G(rm, r′)un(r′)dr′(2.81b) = N X n=1 Ingm(2.81c) gm=j ωϵ k2−∇∇·ˆ∆r′ n G(rm, r′)un(r′)dr′(2.82) To solve the unknown current amplitudes In,Nequations must be obtained from Equation 2.81c by multiplying Equation 2.81c with weight (or test) functions vn, and integrating them along the wire’s length. 34 Chapter 2. Theoretical context As a result, the inner product of each weighting function and gmis formed, which leads to Equation 2.83b. In this way, a set of Nsimultaneous equations written in matrix form are obtained, as shown in Equation 2.84 or 2.85, where Zmn, the impedance matrix, is given by Equation 2.86, and Vm, the source vector, is given by Equation 2.87. The unknown current’s value is then calculated by solving the matrix equation, Equation 2.84 or 2.85 [19]. vn, N X n=1 Ingm=vn, E⇔(2.83a) N X n=1 Invn,gm=vn, E(2.83b)      ⟨v1,g1⟩ ··· ⟨v1,gN⟩ . . ..... . . ⟨vN,g1⟩ ··· ⟨vN,gN⟩           I1 . . . IN      =     Dv1, E1E . . . DvN, ENE      (2.84) [Z] [I] = [V](2.85) Zmn =jωµ ˆ∆rm vm(r)·ˆ∆r′ n un(r′)G(r,r′)dr′dr −j ωϵ ˆ∆rm∇·vm(r)ˆ∆r′ n∇′·un(r′)G(r,r′)dr′dr (2.86) Vm=ˆ∆rm vm(r)· E(r)dr (2.87) In the matrix elements, Equation 2.86, the calculations are carried out using an M -point numerical GaussLegendre quadrature formula, resulting in Equation 2.89 [19], where wpand wqare the quadrature weights. Additionally, Equation 2.89 includes ∆pand ∆q, which are related to the divergence operation. This is because the divergence of a triangular function is 1/∆lor −1/∆lon a segment of length ∆l when used as a basis and test function, depending on whether it is the first or second segment of the Chapter 2. Theoretical context 35 element, as demonstrated in Equation 2.88. It is important to mention that Equation 2.89 only applies to non-overlapping elements because a division by zero occurs when applying it to overlapping elements. For this reason, singularity extraction must be applied for these overlapping elements, leading to Equation 2.92, where S1and S2are given by Equations 2.90 and 2.91, respectively, with abeing the wire radius. Equations 2.92 and 2.89 provide the elements for the impedance matrix, Z, for the overlapping and non-overlapping segments, respectively. ∇·un=                ∇· r −rn |rn+1 −rn|, r ∈rn,rn+1 ∇· rn+2 −r |rn+2 −rn+1|,r ∈rn+1, rn+2 =             ∇· r −rn ∆l ∇· rn+2 −r ∆l =             1 ∆l∇·(r −rn) 1 ∆l∇·(rn+2 −r) =                  1 ∆l∂(rx−rxn) ∂x +∂(ry−ryn) ∂y +∂(rz−rzn) ∂z  1 ∆l"∂rxn+2 −rx ∂x +∂ryn+2 −ry ∂y +∂rzn+2 −rz ∂z # =                1 ∆l∂rx ∂x +∂ry ∂y +∂rz ∂z  1 ∆l−∂rx ∂x −∂ry ∂y −∂rz ∂z  =       1 ∆l −1 ∆l (2.88) Zmn =1 4π M X p=1 M X q=1 wp(rp)wq(r′ q) jωµvm(rp)·un(r′ q)±j ωϵ∆p∆qe−jkrp−r′ q| |rp−r′ q|(2.89) 42 Chapter 2. Theoretical context (a) Octree level 1 (b) Octree level 2 (c) Octree level 3 Figure 2.14: Octree levels 1 to 3 Equation 2.107, T(un,ˆ k),H(k, ˆ k,rab), and R(vm,ˆ k)are, respectively, the radiation , transfer , and receive functions, depicted in Equations 2.108, 2.102, and 2.109. In the radiation function, Equation 2.108, the term [1 −ˆ kˆ k]is included, which is applied to the vector resulting from the integral term. Suppose the resulting vector is called x, and applying the term [1 −ˆ kˆ k]on it results in x −ˆ k(ˆ k·x). Analyzing this reveals that this term operates by removing the ˆ kcomponent from the vector x, leaving only the ˆ θ and ˆ ϕcomponents. Consequently, when pre-computing the radiation and receive functions, only the ˆ θ and ˆ ϕcomponents are stored [19]. L(vm, un) = ˛1 R(vm,ˆ k)·H(k, ˆ k,rab)T(un,ˆ k)dS (2.107) T(un,ˆ k) = h1−ˆ kˆ ki·ˆ∆r′ n un(r′)ejkˆ k·rr′adr′(2.108) R(vm,ˆ k) = jωµ ˆ∆rm vm(r)ejkˆ k·rbr dr (2.109) The integration on the sphere in Equation 2.107 ( [19]) was computed using the Gauss-Legendre quadrature rule in θand Simpson’s rule in ϕ, where the quadrature weights ws(ˆ k)were pre-multiplied with the transfer functions. 2.2.3.1 Multilevel Fast Multipole Method The complexity of the one-level FMM is O(N3 2), but it can be extended and applied recursively, leading to Multilevel Fast Multipole Method (MLFMM) [27], which has a complexity of O(NlogN). In a nutshell, Chapter 2. Theoretical context 43 (a) 3D view (b) Top view Figure 2.15: Linear Position Sensor (LIPS) on octree level 3 the MLFMM is an octal tree structure of Llevels that house the device. This octal tree structure’s first level is a single cubic box recursively divided into eight identical child boxes at subsequent levels. Figure 2.14 (adapted from [19]) shows it on a three-level octal tree, with the first level’s box divided into eight identical child boxes on the second level. Each second-level box is subdivided into eight smaller child boxes on the third level. The near and far regions are delimited similarly to one-level FMM. Like in the one-level FMM, in MLFMM, all neighbour boxes at level lthat share at least one vertice with the local group box are considered to be part of the near group, and there are a maximum of 33= 27 near group boxes, including the local group. The far group boxes at level lare all the child boxes of all parent boxes at level l−1, which are neighbours (that share at least one vertice) of the local group parent box and do not belong to the near group zone at level l. Meaning there are a maximum of 63−33= 189 far group boxes. Furthermore, because the transfer function (Equation 2.104) is a function of only k,ˆ k, and rab, where kand ˆ khave the same values across the entire octal tree level, it is clear that there are only 73−33= 316 different possible values of rab in a particular level of the octal tree. In other words, on each level of the octal tree, it is only, required to compute a maximum of 316 unique transfer functions. On the MLFMM, the Znear matrix is assembled according to the thin-wire elements inside each near region at the highest level only. In other words, this octal tree’s highest level dictates how sparse the Znear will be, meaning that the highest number of levels the octal tree has, the more sparse the Znear matrix will be. However, a higher sub-division level comes at a cost. The closer the centre boxes of a local group and a far group are, the higher multipoles are required to compute the far interactions accurately. Nevertheless, a higher number of multipoles could potentially present numerical issues. Due to the spherical Hankel function, h(2) l(k|rab|)on the transfer function of Equation 2.104, becomes highly oscillatory for fixed k|rab|and increasing l, [19]. Due to that, using MLFMM to simulate significantly small devices’ in comparison to their operation wavelength can be challenging. 44 Chapter 2. Theoretical context (a) 3D view (b) Top view Figure 2.16: LIPS on octree level 4 (a) 3D view (b) Top view Figure 2.17: LIPS on octree level 5 Figures 2.15, 2.16 and 2.17 show an inductive sensor inserted on an octal tree cubic structure. The octal tree in this example has five levels (L= 5), with levels 1 and 2 omitted and levels 3, 4, and 5 depicted in Figures 2.15, 2.16 and 2.17, respectively. In light-green is a local group, and in light-blue and light-red are the near and far regions, respectively, in the correspondent level. The MLFMM algorithm has three stages: aggregation, translation, and disaggregation (see Figure 2.19, adapted from [28] and Algorithm 1). The aggregation phase is the first to be executed and starts at the highest level on the octal tree. It begins by combining the radiation patterns of all basis functions inside each box at the highest level, following steps 5-8 in Algorithm 1. Where tθand tϕare the aggregated field vectors on the ˆ θand ˆ ϕcomponents vector in the unit sphere, Tθand Tϕare the radiation functions (see Equation 2.108), Iis the current vector where Zfar matrix is multiply by using this algorithm, and Wl,L is an interpolation matrix. This interpolation matrix upsamples the result of multiplication of the transfer function with the current vector to the rate needed at the highest level on the octal tree. The use of interpolation is explained further in the current subsection. As the algorithm progress into the lowers levels of the octal tree (upward pass), the radiation fields of the child boxes are aggregated and passed up to their parents using interpolation and phase shift, steps 14-17 in Algorithm 1. Where Wn+1,n is the interpolation matrix that upsamples from the rate at level n+ 1 to the rate at level n, and ejkˆ kn·(rc−rp) is the phase shift term that shifts the phase of the aggregated child field vectors, tθ c,n+1 and tϕ c,n+1, that irradiates at child’s box centre’s position, rc, as it radiates at parent’s box centre’s position, rp. The Chapter 2. Theoretical context 45 Figure 2.18: LIPS on octree Znear mask (white square means 1 and black square means 0) operation (⊙) in lines 15 and 16 denotes term-wise multiplication. Note that levels one and two do not have far groups, so the aggregation phase ends at level three. In the translation phase, the aggregated fields are translated into incoming fields in the downward pass, starting at level 3 and ending at the highest level of the octal tree. This is accomplished by performing all transfers between all far groups at each level of the octal tree, as shown in steps 26, 27, 35, and 36 of Algorithm 1. Where sθand sϕare the local field vectors, and Hrepresents the transfer function depicted in Equation 2.104. The local fields are then passed down to the child boxes in the so-called disaggregation phase, steps 32 and 33 in Algorithm 1, which is done simultaneously with the translation phase. In this phase, the parents’ local fields are phase-shifted, e−jkˆ kn·(rc−rp), and downsampled, WT n,n−1, where T denotes matrix transpose. Finally, steps 40-49 of Algorithm 1 disaggregate the local fields sθand sϕ at the highest level of the octal tree using the receive functions Rθand Rϕfrom Equation 2.109. This disaggregation is done by downsampling the local fields (steps 41 and 42) and multiplying them with the appropriate receive function (step 46). As previously stated, the unit sphere surface integral is computed using the Gauss-Legendre quadrature in the θdimension and Simpson’s rule in the ϕdimension, where the quadrature weights ws(ˆ k)were pre-multiplied with the transfer function. As a result, the integration is performed by simply adding all discretized points in the unit sphere, as shown in step 46 of Algorithm 1. 46 Chapter 2. Theoretical context (a) Aggregation (b) Translation (c) Disaggregation Figure 2.19: MLFMM algorithm The interpolation matrices, Wl,L in steps 6 and 7 and Wn+1,n in steps 15 and 16, are used to save memory by storing the transfer functions and the aggregated field vectors at a lower sampling rate, meeting the Nyquist rate, and upsampling them to the higher rate as the process progresses into the lower levels on the octal tree. This bandwidth increase is due to the bandwidth of the product being the sum of the individual bandwidths. In the disaggregation phase, these interpolation matrices Ware transposed to act as downsampling. The interpolation matrices are formed using Lagrange Polynomials, where the weights are calculated using Equation 2.110, with xrepresenting the position of the points to upsample and x′ representing the position of the upsampled points. Wmn = N Y k=1 k=n x′ m−xk xn−xk (2.110) Chapter 2. Theoretical context 47 Algorithm 1 MLFMM Algorithm 1: {Aggregation Phase} 2: for all non-empty cubes c in level L do 3: Initialize all elements in vector tθ c,L to zero 4: Initialize all elements in vector tϕ c,L to zero 5: for all thin-wire elements e in cube c do 6: tθ c,L =tθ c,L +Wl,L Tθ e,cIe 7: tϕ c,L =tϕ c,L +Wl,L Tϕ e,cIe 8: end for 9: end for 10: for n = L -1 to 3 do 11: for all non-empty cubes p in level n do 12: Initialize all elements in vector tθ p,n to zero 13: Initialize all elements in vector tϕ p,n to zero 14: for all child cubes c in parent cube p do 15: tθ p,n =tθ p,n +ejkˆ kn·(rc−rp)⊙Wn+1,ntθ c,n+1 16: tϕ p,n =tϕ p,n +ejkˆ kn·(rc−rp)⊙Wn+1,ntϕ c,n+1 17: end for 18: end for 19: end for 48 Chapter 2. Theoretical context 20: {Disaggregation Phase} 21: for all non-empty cubes c in level 3 do 22: Initialize all elements in vector sθ c,3to zero 23: Initialize all elements in vector sϕ c,3to zero 24: for all non-empty cubes f far from c do 25: sθ c,3=sθ c,3+H(k, ˆ k3,rfc)⊙tθ f,3 26: sϕ c,3=sϕ c,3+H(k, ˆ k3,rfc)⊙tϕ f,3 27: end for 28: end for 29: for n = 4 to L do 30: for all non-empty cubes c in level n do 31: sθ c,n =WT n,n−1e−jkˆ kn·(rc−rp)⊙sθ c,n−1 32: sϕ c,n =WT n,n−1e−jkˆ kn·(rc−rp)⊙sϕ c,n−1 33: for all non-empty cubes f far from c do 34: sθ c,n =sθ c,n +H(k, ˆ kn,rfc)⊙tθ f,n 35: sϕ c,n =sϕ c,n +H(k, ˆ kn,rfc)⊙tϕ f,n 36: end for 37: end for 38: end for 39: for all non-empty cubes c in level L do 40: sθ c,f =WT f,Lsθ c,L 41: sϕ c,f =WT f,Lsϕ c,L 42: for all thin-wires elements e in cube c do 43: Ve= 0 44: for all unit sphere points p do 45: Ve=Ve+sθ c,f Rθ p,e,c +sϕ c,f Rϕ p,e,c 46: end for 47: end for 48: end for Chapter 2. Theoretical context 49 2.3 Global optimization Global optimization is a technique used to obtain the best or one of the best solutions from a set of possible solutions. It is particularly beneficial when dealing with complex, non-linear problems that have multiple local optima, such as the optimization of inductive position sensors. A variety of global optimization algorithms are available, each with its unique strengths and weaknesses. Examples of popular global optimization algorithms include Simulated Annealing (SA), Particle Swarm Optimization (PSO), and Genetic Algorithm (GA). SA is a global optimization algorithm that draws inspiration from the physical annealing process employed in metallurgy. This algorithm simulates the gradual cooling of a material to enhance the quality of its structure. Similarly, the SA algorithm gradually reduces the temperature of the search space to find the global optimum [29]. PSO is a global optimization algorithm that imitates the collective behaviour of a swarm of birds or a school of fish. Each particle in PSO represents a candidate solution that moves through the search space and updates its position based on the best solution found so far [30]. GA is a global optimization algorithm based on the principles of natural selection and genetics. This algorithm simulates the process of evolution by generating a population of candidate solutions and iteratively improving them using a set of genetic operations, including selection, crossover, and mutation [31]. These algorithms are widely used in various engineering and scientific fields. This section discusses the aforementioned global optimization algorithms: SA, PSO, and GA, keeping in mind the optimization of inductive position sensors. Ultimately, only one of these algorithms is selected and employed in the optimization tool. 2.3.1 Simulated Annealing algorithm SA [32] is a simple and general algorithm for finding a global minimum inspired by an analogy to physical systems. It is a Monte Carlo search method named from the heating-cooling methodology of metal annealing. SA operates by modeling a random walk on the domain of interest, Ωspace. The random walk steps represent the states of the system and its elements ω∈Ω. This algorithm searches the state space 50 Chapter 2. Theoretical context Algorithm 2 Simulated Annealing Algorithm 1: {Initialize state} 2: ω=ω0 3: {Initialize temperature} 4: T=T0 5: for ntemp = 1 to Ntemps do 6: for nstep = 1 to Nsteps do 7: Select at random a neighbour ω′∈η(ω) 8: △E=E(ω′)−E(ω) 9: {Evaluate the change in energy, △E} 10: if △E≤0then 11: {Downhill move} 12: ω=ω′ 13: else 14: Generate a random number Runiform on [0,1] 15: if R < e−△E/Tthen 16: {Uphill move} 17: ω=ω′ 18: end if 19: end if 20: end for 21: {Next lowest temperature} 22: T=Temp [ntemp] 23: end for for low-energy states. At each moment during the search, there is only one current state ω, unlike PSO and GA described further, which search multiple states simultaneously. Based on this current state, the algorithm randomly selects a neighbor ω′∈η(ω), where η(ω)are the neighbors of ω. Then, a bias is incorporated in favor of motions that reduce energy. This is achieved by calculating the energy of ω′, E(ω′), and comparing it to the energy of the current state ω,E(ω). Let △E=E(ω′)−E(ω); when △E > 0, the walk is considered as going uphill, and when △E≤0, the walk is considered as going downhill. To influence the bias, downhill moves are always accepted, while uphill moves are accepted only occasionally. This approach is essential to prevent the system from becoming trapped in local minima. Uphill moves of size △Eare allowed with x△Eprobability, where x∈[0,1] is a control parameter. It is important to note that when x= 1, all moves are approved; however, when x= 0, only downhill moves are accepted. For intermediate xvalues, the likelihood of accepting an uphill move decreases as xlowers. A random walk is used in SA to decrease xfrom an initial value around one to a final value near zero. As the algorithm progresses, less time is spent traveling uphill. In more conventional physical notation, the parameter xis specified in terms of the temperature Tby Chapter 2. Theoretical context 51 x=e−1/T, mimicking a state with variable temperatures where the temperature affects the decisionmaking probability at each step [29]. Algorithm 2 resumes all the above described. 2.3.2 Particle Swarm Optimization PSO was initially proposed by Eberhard and Kennedy [33–35] as an alternative to traditional optimization techniques, such as gradient descent and SA. PSO draws inspiration from the social behaviour of bird flocks, fish schools, and insect swarms, and aims to mimic the collective intelligence of these animals. The PSO algorithm simulates the movement of particles (individuals) in an n-dimensional search space, where each particle (individual) represents a candidate solution to a particular problem. In the initial iteration, the particle positions and velocities are randomly initialized in the search space. Subsequently, in the following iterations, the velocity of each particle is determined by the best position achieved by any neighbouring member of the swarm, the best position achieved by the particle, and the particle’s previous position. The best position found by any member of the swarm is referred to as the global best position. The algorithm terminates when the quality of the current solution is satisfactory or when the maximum number of iterations is reached. Otherwise, the velocities and positions of the particles are updated, and the next iteration begins. The updated formulas for the velocity and position of the particle are given below: Vn(t+ 1) = w·Vn(t) + c1φ1Xbest n−Xn+c2φ2Xbest −Xn(2.111) Xn(t+ 1) = Xn(t) + Vn(t+ 1) (2.112) Where c1and c2are positive constants, and φ1and φ2are two random variables with a uniform distribution between zero and unity. The Equation 2.111 incorporates the inertia weight w, which represents the impact of the previous velocity vector, Vn(t), on the new velocity vector, Vn(t+ 1). To prevent particles from moving too quickly from one area in the search space to another, the velocity is restricted by an upper bound limit Vmax. At each iteration, the best solution value, Pbest n, and its corresponding position, Xbest n, are updated for each particle when a better solution is found. Similarly, the best swarm solution, 58 Chapter 2. Theoretical context inductive position sensors. Finally, global optimization techniques were discussed, focusing on SA, PSO, and GA. These algorithms were examined having in mind the optimization of inductive position sensors, and only one of them would be selected for use in the optimization tool developed in this thesis. In conclusion, this chapter has provided the necessary theoretical background for the development of an Active Noise Cancellation (ANC) based on wavelets and a customized simulation tool aimed at optimizing inductive position sensors. The concepts and methods introduced in this chapter will be applied and further explored in the subsequent chapters of this thesis. Chapter 3 Mechanical waves application – ANC An Active Noise Cancellation (ANC) system is used to reduce unwanted noises, either periodic or random [46]. There are many applications for these systems, ranging from headphones for aviation or general usage [47] to the automotive sector. The latter uses ANC systems to eliminate stationary noises originating from the engine, road, or the external environment [48]. The working principle consists of a transducer, usually a speaker that reproduces a cancellation signal, which is a sound wave (Mechanical wave) with the same amplitude as the signal to be cancelled, but with an inverted phase. The two sound waves (the signal to be cancelled and the cancellation signal) combine to form a new wave, in a process called interference, resulting in destructive interference as the two signals cancel each other out [46]. 3.1 State-of-the-art Over recent years, ANC systems have been extensively researched and several systems were developed for the automotive industry to improve the acoustic comfort of vehicle occupants by reducing engine, road, and wind noise [49–55]. These systems work by analyzing and identifying the frequency components of unwanted ambient noise and generating an opposing anti-noise signal to cancel it out. The most common type of ANC system used in automotive applications is the feedback ANC system, which uses microphones to measure the ambient noise and speakers to generate an opposing anti-noise signal to cancel it out. These systems are typically based on adaptive filtering algorithms such as the Least Mean Square (LMS) [56] and the Normalized Least Mean Square (NLMS) [57,58], which continuously adapt the filter coefficients in real-time to reduce the noise. However, these systems estimate the secondary path, 59 60 Chapter 3. Mechanical waves application – ANC the acoustic path from the microphones to the speakers, only during the calibration phase. If the acoustic properties change for some reason, such as a change in vehicle interior layout or the introduction of new materials that affect the sound transmission, the performance of the ANC system can drop significantly, and in some cases, it can even become unstable and generate additional noise. To overcome this problem, some researchers, such as Padhi et al., propose a technique for evaluating the secondary path during the ANC activity [59]. Recent advancement in ANC systems uses frequency-domain techniques, such as wavelets, to decompose the noise signal into its frequency components, allowing for the developing of more targeted and effective control algorithms. Wavelet-based ANC systems are able to reduce noise in specific frequency ranges more effectively and provide the ability to separate noise sources and treat them individually [52,60–66]. This is an alternative to the popular frequency-domain filtering method based on the Fast Fourier Transform (FFT) algorithm. Another trend in recent ANC systems is using artificial neural networks to improve their performance. This results in a more effective noise reduction than traditional ANC systems. Studies have shown that this approach can be more effective than methods such as LMS, NLMS, and others [60, 67–72]. Researchers proposed using hybrid systems that combine different ANC methods to achieve better performance [59, 73]. These hybrid systems usually consist of a sinusoidal noise cancellers subsystem, narrowband ANC subsystem and broadband ANC subsystem. The conventional hybrid ANC system is capable of suppressing mixed noise. However, in non-linear environments, the attenuation performance of the conventional hybrid ANC system can be significantly weakened due to its lack of non-linear components. An improved hybrid ANC system, the hybrid Functional Link Artificial Neural Network (FLANN) system, has been proposed to overcome this limitation [70,72]. This system adds the FLANN structure to the conventional hybrid ANC system. The use of smart materials, such as piezoelectric materials, is another area of research in ANC systems. These materials can actively change their properties in response to external stimuli, creating active noise cancellation systems that can adapt to changing noise levels in real-time, providing a more effective and efficient noise reduction [53]. Another aspect of ANC systems is the number of microphones and loudspeakers. A Single-Input SingleOutput (SISO) ANC system works well in attenuating noise at a specific point in space, such as headphones, but is less effective in large areas like the passenger compartment of a vehicle. In such cases, Chapter 3. Mechanical waves application – ANC 61 (a) An ANC region defined by its spatial boundaries, represented by a black area, comprises a circular array of microphones with a radius R1 and a circular array of loudspeakers with a radius R2 (b) Primary noise field in free-field (200 Hz) (c) Noise field after ANC activation (free-field) (d) Primary noise field in the reverberant environment (200 Hz) (e) Noise field after ANC activation (reverberant environment) Figure 3.1: Noise cancellation performance of a Multiple-Input Multiple-Output (MIMO) ANC system, studied by Zhang et al. 62 Chapter 3. Mechanical waves application – ANC a MIMO ANC system is needed. This system uses a constellation of microphones and loudspeakers to capture and cancel out unwanted noise. In a study by Zhang et al., an ANC system composed of eight microphones and eight loudspeakers arranged in a circular configuration was proposed [74]. The researchers evaluated the system’s performance in both a free-field and reverberant environment, with and without the ANC system activated. The results (Figure 3.1 [74]) showed that the proposed MIMO ANC system effectively reduced noise in both environments, with a noticeable noise reduction when the ANC system was activated. Recent studies have proposed using global optimization algorithms to design and identify the optimal configuration of an ANC system. For example, Long et al. used a Genetic Algorithm (GA) to find an optimal solution for a sub-band ANC system for broadband noise cancellation [75], and Porghoveh et al. developed an ANC system for reducing engine noise within smart cubic vehicle enclosures using a global optimization framework based on Particle Swarm Optimization (PSO) [49]. Overall, ANC technology is constantly evolving to meet new challenges and demands in the automotive industry to provide a more comfortable and pleasant ride for the passengers while also reducing the need for bulky and heavy acoustic insulation materials, which contributes to a more lightweight and fuel-efficient vehicle consequently reducing the environmental impact. 3.2 Classical feedback ANC system In a typical ANC system, an adaptive filter called FxLMS [56] is used. The core of this algorithm is the LMS method [76–79]. Equation 3.1 is used to calculate the coefficients ( ˆw) for the LMS Finite Impulse Response (FIR) filter, which approximates the response of the primary path (from the noise source to the error microphone, see Figure 3.2) to minimize the error signal ( e[n] ). Where ˆw [n] is the current LMS FIR filter coefficient vector, µis the step size, xs[n] is the delayed signal reference, e[n] is the error signal, T denotes matrix transpose, and ˆw [n+1] is the LMS FIR filter coefficients vector for the new iteration. The error signal results from the sum of the noise signal ( d[n] ) with the anti-phase signal ( y’[n] ) (see Equation 3.2). The FIR filter output signal, y[n] , is the outcome of transforming the reference signal ( x[n] ) accordingly to the LMS’s calculated coefficients, as described by Equation 3.3. Then, the anti-phase signal is obtained by multiplying the filter output signal ( y[n] ) by -1. The error signal will converge to zero as the FIR filter response converges to that of the primary path (acoustic environment). Chapter 3. Mechanical waves application – ANC 63 FxNLMS Algorithm Secondary Path Primary Path Acoustic Environment Noise Power Amplifier Second Path Estimation x[n] d[n] e[n] -y[n] Reference Input Signal Coeff. Update xs[n] y[n] X-y[n] -1 FIR Filter x[n] Error Microphone Secondary Loudspeaker Noise Source Figure 3.2: Filtered-x Least Mean Square (FxLMS) block diagram The LMS algorithm has the limitation of requiring virtually no delay between the emitted anti-phase signal and the input of the coefficients update algorithm. This path is commonly referred to in the literature as a secondary path. For that reason, Widrow et al. [56] proposed the FxLMS to account for the secondary path’s delays. By considering the effect of the secondary path on the system, the algorithm can correctly calculate the coefficients that minimize the error. For this reason, the same delay must be created in the reference signal that feeds the coefficients update algorithm. ˆw[n+ 1] = ˆw[n] + µxs[n]eT[n](3.1) e[n] = d[n]−y′[n](3.2) y[n] = ˆwTx[n](3.3) The adaptive filter LMS is susceptible to the input power changes of the reference signal. For example, a sudden volume change may turn the algorithm unstable [80]. To solve this problem, one common solution is the normalization of the reference signal before using it to calculate the coefficients. Using this method, called NLMS, the estimation of the learning rate is much more effective, which allows the error to converge to zero. Equation 3.4 is the formula for updating the NLMS coefficients. 64 Chapter 3. Mechanical waves application – ANC ˆw[n+ 1] = ˆw[n] + µ ∥xs[n]∥2xs[n]eT[n](3.4) 3.3 Proposed hardware accelerated ANC system The present section describes the system overview of the wavelet-based ANC system proposed in this thesis, and the arrangement of the various constituent elements (secondary loudspeakers, noise and error microphones). The positioning of the elements was crucial for the overall performance of the proposed ANC system. 3.3.1 System overview The block diagram depicted in Figure 3.3 represents a MIMO approach, with two error microphones, two secondary loudspeakers, and one reference microphone. The Microelectromechanical System (MEMS) based microphones used have a Pulse Density Modulation (PDM) output signal. This justifies the need for the PDM to Pulse Code Modulation (PCM) converter right after the microphones. The noise source is entirely independent of the system. The error ( e[n] ) and reference ( x[n] ) signals are processed and divided into several sub-band channels using Undecimated Wavelet Packet Transform (UWPT). Each channel is an input for a dedicated Filtered-x Normalized Least Mean Square (FxNLMS) algorithm. Four FIR filters apply the secondary path estimation to the reference signal, one for each secondary loudspeaker and error microphone combination. In addition, there are two other FIR filters to apply the third path estimation (the path from the secondary loudspeakers to the reference microphone) to the anti-phase signal ( y[n] ). The goal is to eliminate the presence of the anti-phase signal on the reference signal (as the reference microphone shares the same space as the secondary loudspeakers) and obtain a signal that contains almost only information from the noise source. Finally, the outputs of all FxNLMS algorithms are added and then shifted to obtain the anti-phase signal. To further elaborate, each output channel of the UWPT is directly connected to the input of a distinct FxNLMS module, which processes the sub-bands of the input signal derived from the UWPT outputs in parallel. The number of FxNLMS modules running in parallel influences the system’s efficiency. The higher the number of modules running, the better the cancellation of the noise will be. However, since each Chapter 3. Mechanical waves application – ANC 65 NLMS NLMS Coeff. Update FIR Filter FIR Filter FIR Filter FIR Filter FIR Filter FIR Filter Wavelet Decomposition Wavelet Decomposition Wavelet Decomposition Wavelet Decomposition Wavelet Decomposition Wavelet Decomposition Wavelet Decomposition Secondary path estimation (right speaker; right mic) Secondary path estimation (right speaker; left mic) Secondary path estimation (left speaker; right mic) Secondary path estimation (left speaker; left mic) ∑ ∑ Tertiary path estimation (right speaker; noise mic) Tertiary path estimation (left speaker; noise mic) Multiply by -1 Multiply by -1 Left secondary loudspeaker Right secondary loudspeaker Left error microphone Right error microphone Reference microphone NLMS ∑ -y1[n] -y2[n] y2[n] y1[n] e1[n] e2[n] x[n] xs22[n] xs21[n] xs12[n] xs11[n] PDM to PCM PDM to PCM PDM to PCM Figure 3.3: ANC block diagram with wavelets reference signal channel has its unique power, it is necessary to normalize each channel independently. In this way, all error signals from all channels converge to zero at almost the same rate. It was decided to implement the UWPT rather than the Dyadic Wavelet Transform (DWT) or Wavelet Packet Transform (WPT). These techniques have drawbacks that limit their usage. As an example, the DWT has two limitations: i) the output rates on each output channel are different and; ii) there is an output channel that contains half of the spectrum of the input signal (i.e., higher frequencies). For that reason, the FxLMS algorithm will operate with a high bandwidth input, and consequently, this would result in lower performance. In the case of the WPT, the drawback is the low output sampling rate. The FxLMS algorithm requires an adequate sampling rate of the input signal (which is related to the update of the coefficients). This sampling rate is particularly important to actuate on the higher spectral components of the input signal. A lower sampling rate means that the FxLMS will not be able to cope with such fast-transitioning signals. For those reasons, the UWPT was the selected wavelet implementation. It has all outputs with the same sample rate and with the same rate as the input signal. Also, the input signal’s bandwidth is equally divided by all outputs. In this way, the FxNLMS algorithm is constantly fed with an adequate sampling rate, thus improving its overall performance. 3.3.2 Arrangement for random noise cancellation One way to interpret the FxNLMS algorithm is to bear in mind that its objective is to constantly search for the coefficients (of the FIR filter) that estimate the response of the primary path. Nevertheless, it 66 Chapter 3. Mechanical waves application – ANC Noise Speaker Reference Microphone Secondary Loudspeaker Secondary Loudspeaker Error Microphone Error Microphone 75cm 40cm 4.3cm 50cm Figure 3.4: Example of ANC components arrangement is not always possible to find a physically realistic solution. For example, the delay of the primary path could be smaller than the delay of the secondary path. Such a particular solution would not represent a typical geometry of an ANC system. To improve the stability of the system, the constituent components were arranged, as illustrated in Figure 3.4. The delay of the primary path (the distance of the reference microphone to the error microphone) is larger than the secondary path (distance from the secondary loudspeaker to the error microphone). 3.3.3 Field Programmable Gate Array (FPGA) implementation The implementation of the ANC system comprises several distinct modules, with the primary focus on five key components: the PDM to PCM converter, the FIR filters, the UWPT, the FxNLMS, and the Inter-IC Sound (I2S) output module. The PDM to PCM converter module transforms the output signal from the MEMS microphones, in PDM format, into a PCM format that the implemented ANC system can process. The FIR filter modules serve multiple purposes (see block diagram of Figure 3.3), such as applying second path estimations to the reference signal ( Second path estimation blocks in the block diagram) (x[n]), third path estimations ( Third path estimation blocks in the block diagram) to the anti-phase signals (y1[n]and y2[n]), and, within the FxNLMS module, primary path estimations to each sub-band of the x[n]signal. The UWPT module (identified as Wavelet Decomposition in the block diagram) decomposes error signals (e1[n]and e2[n]), and reference signals (x[n],xs11 [n],xs12 [n],xs21 [n], and xs22 [n]) into several sub-bands. The FxNLMS module continuously estimates the primary path, with one module dedicated to each sub-band decomposition. Since the overall system performance primarily depends on the number of parallel FxNLMS algorithms, the modules have been designed to optimize the balance between FPGA logic Chapter 3. Mechanical waves application – ANC 67 Figure 3.5: FPGA implementation block diagram resource usage and performance. The I2S output module supplies an I2S amplifier with the anti-phase signals (−y1[n]and −y2[n]). Given the requirement for a considerable number of modules running in parallel, an FPGA is a suitable choice, as it enables hardware acceleration and promotes hardware parallelization. It is also crucial to select an FPGA with enough resources to implement the ANC system. Furthermore, these resources must accommodate all the FxNLMS modules necessary for achieving effective random noise reduction, as previously mentioned in Section 3.3.1. To optimize random noise cancellation through a larger number of parallel FxNLMS modules, the XILINX FPGA XC7Z045 was chosen. This FPGA offers one of the highest numbers of Digital Signal Processors (DSPs) available on the market within its family. Figure 3.5 depicts the FPGA ANC implementation block diagram. In this block diagram, the PDM to PCM conversion and the I2S output modules are not present for simplicity. Here, the diagram inputs, mu , left_mic , right_mic , and noise_mic are respectively the step size µ, the left error microphone signal 74 Chapter 3. Mechanical waves application – ANC port of the BRAM_A , where each bit represents the value of a given coefficient for the associated filter. The final result of the filter is truncated (the 18 most significant bits). In the BRAM_A the coefficients of the 64 filters are stored, having 1-bit wide word, being its content depicted in Figure 3.7. As already mentioned, the output of the BRAM_A is directly connected to the SEL bit of the arithmetic unit. The access to the BRAM’s content is done at the system clock, being the read address port ( ADDR_RD ) managed by the six least significative bits of the Counter_A , required to address the 64 coefficients. The acquisition of the input signal (18-bit) happens when the Counter_A (incrementing at the system clock) overflow to 0 value. When the overflow occurs, the input signal is fetched by a flip-flop in order to be stored on the BRAM_B , implemented as a circular buffer of 64 elements. The write address port ( ADDR_WR ) is updated at the input sampling frequency and its value corresponds to CounterB−1, which is the address where the incoming samples will be stored. Note that this counter is also working at the input signal’s frequency, 47.97 kHz. Whenever the end of the BRAM_B is reached, the incoming sample overwrites the oldest stored sample. Reading samples from the BRAM_B is done through the read address port ( ADDR_RD ). The Counter_B value is used as a reference for the first value to be read (write address value plus one). To the reference, an offset is added, which is related to the sequential nature of the filtering process. Since the reading is done at the system clock, all the samples stored in memory are read from the oldest to the newest (in this order) for the filtering process. On the filters accumulator flip-flop, there is a reset signal. This is needed for the convolution of each new input sample. Moreover, the BRAMs elements work as a pipeline system, but two clock cycles are required to output the memory requested content. For this reason, the reset signal of the accumulator flip-flop is set when Counter_A is equal to 2. This delay and the internal delay of 1 clock cycle of the accumulator flip-flop result in a final clock count of 67. It takes 67 system clocks to obtain a processed sample. 3.3.3.4 FxNLMS module The FxNLMS module interface is depicted in Figure 3.10 and features 11 inputs and 2 outputs. The anc_on input enables or disables the output’s anti-phase signals, y1 and y2 . The mu input governs the convergence speed of the FxNLMS algorithm, while inputs E1 and E2 correspond to the error microphone signals. The reference signal is represented by X , and the inputs XS11 , XS12 , XS21 , and XS22 are the reference signal filtered by the second path estimation coefficients. Chapter 3. Mechanical waves application – ANC 75 clk rst anc_on mu[5:0] E1[1151:0] y1[17:0] E2[1151:0] X[1151:0] XS11[1151:0] XS12[1151:0] XS21[1151:0] XS22[1151:0] y2[17:0] FxNLMS Figure 3.10: FxNLMS module interface The module has two main components, as illustrated in Figure 3.3: the FIR filters and the coefficient update. In the first stage at the coefficients update block, there is a normalization operation. This normalization was crafted to avoid mathematical divisions, thus saving DSP resources (details can be found later in this section). Each UWPT outputs 64 channels and each channel is processed by a single FxNLMS module. As such, there are 64 FxNLMS per UWPT, and since a MIMO configuration was used, two FIR filters are needed for each FxNLMS algorithm (one for each cancelling speaker). In total are required the instantiation of 128 FIR blocks to filter the reference signal, x[n], using the on-line primary path estimation coefficients, ˆw[n]. Each FIR filter has 1000 coefficients of 48-bit. The number of coefficients was limited to the available resources of the FPGA. The implementation of FIR filters inside the FxNLMS module is similar to the one described in the FIR filter module section. The coefficients are stored in BRAMs as well as the last 1000 samples of the reference signal x[n]. The storage of the input signal is done as described in the FIR filter module section. As the FPGA’s DSPs have 18+25 bit inputs and 48-bit outputs each and the coefficients are 48-bit, three DSPs per filter are required. As previously mentioned, the total number of coefficients, 1000, was determined by the available resources of the FPGA, which has a total of 545 BRAMs with a capacity of 36 kbit each. These BRAMs can be configured as memory of various sizes, such as 32K of 1-bit word, 16K of 2-bit word, 8K of 4-bit word, 4K of 9-bit (or 8-bit) word, 2K of 18-bit (or 16-bit) word, 1K of 36-bit (or 32-bit) word, or 512 of 72-bit (or 64-bit) word. Additionally, each BRAM can be divided into two completely independent 18 Kb BRAMs. Given these restrictions and the 48-bit word size of the coefficients, one complete BRAM was configured as 1k of 32-bit word, and half a BRAM as 1k of 16-bit word, necessitating the use of 1.5 BRAMs. Moreover, 76 Chapter 3. Mechanical waves application – ANC Table 3.1: Computing a very rough approximation to the inverse of a number by reversing the order of the bits Number Inverse Decimal Binary Binary Decimal 1 00000001.0000000 00000001.0000000 1.0000000 2 00000010.0000000 00000000.1000000 0.5000000 4 00000100.0000000 00000000.0100000 0.2500000 8 00001000.0000000 00000000.0010000 0.1250000 16 00010000.0000000 00000000.0001000 0.0625000 32 00100000.0000000 00000000.0000100 0.0312500 64 01000000.0000000 00000000.0000010 0.0156250 100 01100100.0000000 00000000.0000010 0.0156250 storing the last 1000 samples of the reference signals, x[n],xs11 [n],xs12 [n],xs21 [n], and xs22 [n], which have an 18-bit word size, required 2.5 BRAMs, configured as five halves of 1k of 18-bit word. Altogether, this amounted to 1.5 BRAMs per FIR filter per FxNLMS module, plus 2.5 BRAMs, resulting in 5.5 BRAMs per FxNLMS module. When accounting for the 64 FxNLMS modules running in parallel, a total of 352 BRAMs out of the available 545 were needed. The NLMS coefficients update block implements Equation 3.4. The algorithm starts by calculating the noise signal energy, ∥xs[n]∥2, of the last 1000 samples. This is done by adding all the 1000 elements after being multiplied by themselves (see Equation 3.6). For this operation, the same strategy as described in the FIR filter module section was used, the multiplications are performed serially (pipeline) to minimize the DSPs number. ∥xs[n]∥2=Xxs[n]2(3.6) After calculating the energy of the xs[n] signal, it is necessary to invert the result. To avoid the division operation, the order of the bits was reversed around the decimal point, followed by a shift of one bit to the left, as shown in Table 3.1. Then, all bits to the left of the first 1, counting from right to left, are set to be zero. In this way, a very rough approximation of the inverse of a number is obtained. Once the inverse calculation of the xs[n] signal energy is done, the next step is to multiply it by the constant µ(see Equation 3.4), which is always a power of two. This multiplication is done by shifting log2(µ)bits to the left. The rest of the calculations are done serially to use as few DSPs as possible. Chapter 3. Mechanical waves application – ANC 77 Table 3.2: FPGA resources utilization Resource Utilization Available Utilization LUT 82061 218600 37.54 % LUTRAM 4855 70400 6.90 % FF 108629 437200 24.85 % BRAM 365 545 66.97 % DSP 900 900 100.00 % IO 26 250 10.40 % BUFG 3 32 9.38 % To conclude, Table 3.2 shows the total resources used by the proposed ANC in a XILINX FPGA XC7Z045 . 3.3.3.5 I2S output module The output module represents the final stage of the ANC system and is responsible for transmitting the processed PCM signal to an external I2S stereo amplifier. The I2S interface is a serial communication protocol explicitly designed for transmitting digital audio data between devices. It has gained widespread use in audio systems due to its simplicity, synchronization capabilities, and low jitter characteristics. The I2S protocol necessitates three signals for operation: •Bit Clock (BCLK): Provides the clock signal for synchronizing data transmission between the output module and the I2S stereo amplifier. •Word Select (WS): Also known as the Left/Right Clock (LRCLK), this signal indicates whether the current data sample belongs to the left or right audio channel. •Serial Data (SD): Transmits the actual audio data in a Time-Division Multiplexing (TDM) format. The output module obtains the processed PCM signal from the previous stage of the ANC system (FxNLMS module) and formats the PCM signal in accordance with the I2S protocol requirements. This process involves organizing the data into left and right audio channels and ensuring the correct bit depth and sample rate. In this particular ANC system, the external I2S stereo amplifier is designed to accept a bit depth of 20 bits and supports various sample rates, such as 44.1 kHz, 48 kHz, 96 kHz, and 192 kHz. However, for this specific implementation, an 18-bit bit depth is employed, accompanied by a sample rate of 47.97 kHz, which is close to the accepted 48 kHz. The amplifier can accommodate this sample rate due to the synchronization capabilities inherent in the I2S protocol. 78 Chapter 3. Mechanical waves application – ANC 3.3.3.6 Second and third path impulse responses The impulse responses of the second and third paths (see Figure 3.2 and Subsection 3.3.1, respectively) are calculated in the Processing System (PS). First, the PS plays a wave file on the left secondary loudspeaker that contains a pre-recorded swept-frequency sine. At the same time, the signals captured by all three microphones, the reference microphone and the two error microphones, are recorded. Next, the same procedure is repeated for the right secondary loudspeaker, and three more captured microphone signals are stored. Then, the PS calculates all the six impulse responses using the reproduced signal and the six recorded signals. Finally, it sends all six impulse responses to the respective FIR filter modules on the Programmable Logic (PL) side. The procedure described here, in this sub-section, is the first operation that the ANC system does when it is turned on. 3.4 System evaluation This section starts with the presentation of the methodology used to assess the ANC system developed during this research, followed by the results and a brief discussion. Chapter 3. Mechanical waves application – ANC 79 (a) FPGA system board (b) System hardware Figure 3.11: Hardware accelerated ANC system 3.4.1 Test methodology The tests of this ANC system were done in a spacious and uncontrolled acoustic environment, more concretely an academic auditorium (Figure 3.12), where a loudspeaker was used as the noise source. This loudspeaker, including its amplifier and controller, are entirely independent of the ANC system. The ANC system controls two loudspeakers to generate the anti-phase signal and uses the reference microphone and the two error microphones to close feedback path. These elements were disposed as illustrated in Figure 3.4. As the source of ”noise”, four signals were selected to evaluate the system’s performance: O Fortuna , a piece of classical music composed by Carl Orff and directed by André Rieu; an urban street environment sound captured in a busy city; an ambient noise of a crowded restaurant; and pink noise. O Fortuna was chosen because it has a rich spectrum and a high dynamic range, suitable for showing the performance of this ANC system. Urban street and crowded restaurant ambient noises were selected because they are unpredictable and are more realistic use cases, and pink noise because it is random and has a full spectrum that fades along with the frequency, suitable for evaluating the system’s dynamic range. 80 Chapter 3. Mechanical waves application – ANC Figure 3.12: ANC in the academic auditorium 3.4.2 Results The results are organized as follows. For each source of ”noise” ( O Fortuna , urban street, crowded restaurant, and pink noise) there are four associated graphics (Figures 3.13, 3.14, 3.15 and 3.16). The first graphic (top) depicts the noise (time representation) captured by one of the two error microphones: the red curve was captured when the ANC system was off and the blue curve when the ANC system was turned on . Note that the temporal signals were time-shifted so that they became overlaid. The calculation of the gain was done using these signals values. The second and third graphics (central) are spectrograms of the time domain signals. To draw the spectrograms (exclusively), a WPT up to level 9 with the Daubechies wavelet having two vanishing moments ( db2 ) was used. This allowed a good translatescale resolution, being the scale ordered by its frequency dominance. The fourth and last graphic shows the noise spectrum (FFT) when the ANC system is off (blue curve) and when the ANC system is on (red curve). Figure 3.13 shows the results obtained when the ”noise” source is the classical piece, O Fortuna . Using the data of the time domain curves, the ANC system was able to reduce an average of 9.45 dB of ”noise”. The central spectrograms depict the spectral noise evolution over time, when the system is on and when the system is off , respectively. As can be seen, the ANC system was capable of reducing an entire band and working with the same performance in low and high amplitude noise. The bottom graphic shows the FFT of the time domain curves for their entire duration (15s). There was a decrease in intensity for almost all frequencies below 7 kHz. The system does not show any cancellation effect above this frequency. Chapter 3. Mechanical waves application – ANC 81 Figure 3.14 depicts the result obtained when urban environment was used to generate the noise signal. In this case, the attenuation gain was, on average, 9.09 dB. Considering the time-domain signals, the noise’s attenuation is noticeable when the ANC system is working. Regarding the spectrograms, an attenuation in the whole spectrum is visible along with the entire signal FFT, as expected. Between the 8th to 10th second, an engine of a car sound was present, and it is visible that the system attenuated this high-intensity sound. Looking at the FFT, the system reduced almost all frequencies up to 7 kHz, as was already observed in the previous result. Figure 3.15 illustrates the ANC system response to the noise of a quotidian sound of a crowded restaurant. The trend from the previous results remains present. There is a noticeable reduction of the noise (temporal signals), more precisely 9.06 dB on average. The spectrograms show an attenuation for the entire lower part of the spectrum, and the FFT support this. The last result (Figure 3.16) examines the ANC system’s response to pink noise. By analyzing the time domain of the spectrograms, it is evident that the noise was attenuated, though to a lesser degree compared to the previously discussed results. In this instance, the average attenuation reached 6.34 dB. By inspecting both the spectrogram and the FFT graphics in Figure 3.16, it becomes clear that the system was particularly effective in reducing frequencies below 7 kHz. Consistent with the responses observed for the other three noise sources, nearly all frequencies up to 7 kHz experienced a reduction. 82 Chapter 3. Mechanical waves application – ANC (a) Time domain and spectrogram of the error microphone’s signal (b) Frequency domain of the error microphone’s signal Figure 3.13: Experimental results of the proposed ANC system - classical music as a ”noise” source Chapter 3. Mechanical waves application – ANC 83 (a) Time domain and spectrogram of the error microphone’s signal (b) Frequency domain of the error microphone’s signal Figure 3.14: Experimental results of the proposed ANC system - city ambient sound as a noise source 90 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology the system can infer the target’s position, distance, or orientation. This technology offers advantages such as robustness in harsh environments, resistance to external magnetic fields, contactless operation, low cost, and most importantly, the absence of rare raw elements in its material composition. These benefits make it suitable for various applications, including automotive systems, industrial automation, and medical devices. Despite its advantages, the performance of inductive position sensors is closely related to the shape of the sensor’s coil, emphasizing the need for a dedicated optimization tool for the coil’s design. To the best of the authors’ knowledge, there are currently no such tools available in the market. 4.1 State-of-the-art A position sensor measures a relative, linear or angular position concerning a reference. In most cases, this type of sensor is a transduction mechanism that converts physical quantities, for example, distance travelled, into electrical quantities. The market is filled with several position sensors designed to address specific requirements. They are widely used in robotics, mechatronics, automation, machine tools, automotive, avionics, aerospace, medical devices, and consumer electronics. The increasing demand for these sensors results in a market requirement for more reliability, higher performance, and lower power sensors. For applications in harsh environments, such as automotive or industrial applications, may limit technology choices. Usually, they require sensor robustness to various influential factors, such as temperature, mechanical vibrations, mechanical tolerances, and external magnetic fields. In addition to the above factors, the cost of manufacturing the position sensors must be kept competitive. Furthermore, the position sensors must meet strict specifications regarding resolution, accuracy, hysteresis, robustness and stability. The following sections describe several physical principles used in developing position sensors. These are classified into magnetic, optical, capacitive and inductive based. The magnetic category is further subdivided into individual technologies. 4.1.1 Magnetic-based position sensors Magnetic-based position sensors are used primarily in the automotive industry due to their robustness, contactless, and non-destructive sensing. The two main operating principles of magnetic-based position sensors are Hall-effect, where a voltage is produced by magnetic deflection of current-carrying electrons, Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 91 and magnetoresistance, where a change in resistance occurs in a material exposed to an external magnetic field. Magnetic-based position sensors are commonly found in applications with permanent-magnet machines (especially, Hall-effect sensors) like, Antilock Braking System (ABS), Traction Control System (TCS) or Electronic Stability Program (ESP) [81–84]. Among the most prominent manufacturers of this type of sensors are Bosch, Honeywell International, Inc., STMicroelectronics N.V., Infineon Technologies AG, and TDK Corporation. A Hall-effect sensor is a transducer that measures a magnetic field based on the Hall-effect principle— discovered in 1879 by American physicist Edwin Herbert Hall [85]. It results from Lorentz’s force, which deflects moving charge carriers in the presence of a magnetic field applied perpendicularly to the current flow. When a magnetic field Bis applied in the direction parallel to the current flow, Lorentz force deflects the moving charge carrier towards the surface. Consequently, charges of opposite sign accumulate at two surfaces/edges orthogonal to the direction of current flow, resulting in a differential voltage between them. A critical parameter determining the sensitivity of Hall-effect sensors is their building materials’ carrier mobility [86]. Hall-effect sensors present several advantages: simple device architecture, easy manufacturing, low cost, the possibility of scaling down, integration with Complementary Metal Oxide Semiconductor (CMOS) circuits, and can provide contactless measurement [87]. Nevertheless, Hall-effect sensors suffer from some drawbacks: i) the output signal is very weak compared to other technologies, namely magnetoresistive sensors, which require electronics to amplify the output signal; ii) the output signal drifts with temperature due to the thermal instability of the devices, which requires complex electronics to compensate the drift; iii) the sensitivity to external magnetic disturbances can cause undesired effects such as erroneous data measurement [87]. Magnetoresistance-based sensors, similar to the previously mentioned Hall-effect-based sensors, can also provide contactless measurements. Magnetoresistance-based sensors are sub-categorized into Tunnel Magnetoresistance (TMR), Anisotropic Magnetoresistance (AMR) and Giant Magnetoresistance (GMR). They all share a similar operational property, dependence of the electrical resistance on the angle between the electric current’s direction and the direction of the magnetic field, but these have different structures [87]. Certain materials’ resistance is influenced by the angle formed between the magnetization of the material and the applied current’s direction [87]. This phenomenon, called AMR, was initially observed by William Thomson (Lord Kelvin) in 1856. At the time, he noted that the resistance is highest when the magnetization 92 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology is parallel to the current and lowest when it is perpendicular [88]. AMR sensors have the disadvantage of being less sensitive than the other magnetoresistance technologies, such as GMR and TMR sensors. Nonetheless, they are much simpler to make, have more design flexibility regarding device shape and resistance, and perform better at low frequencies in terms of Signal-to-Noise Ratio (SNR) [89, 90]. Because of their robustness, AMR sensors are used in linear positioning systems [91], and an Angular Position Sensor (APS) is described in [92] using a quadruple permalloy layer AMR sensor. The resistance of a thin-film structure composed of two ferromagnetic layers separated by a non-ferromagnetic spacer layer depends on the relative alignments of the magnetizations of the two ferromagnetic layers. This effect is known as the giant magnetoresistive effect, discovered independently by two scientists who were awarded the 2007 Nobel Prize in Physics for its discovery, Albert Fert and Peter Grunberg [93,94]. One of the two ferromagnetic layers of a typical GMR sensor is manufactured to have a fixed direction of magnetization. This layer is known as the pinned (or reference) layer. The other one is the free layer, which is made of a soft magnetic substance that, under the effect of the applied magnetic field, has its magnetization direction rotates accordingly [87]. The GMR effect is a product of the spin-dependent scattering of electrons. For example, suppose the magnetization of the pinned layer is parallel to that of the free layer. In that case, the device’s resistance is relatively small because electrons with a spin direction opposite to the magnetization are more weakly scattered. On the other hand, if the magnetization directions of the pinned and free layers are antiparallel to each other, the device’s resistance increases substantially due to electrons with both up and down spin polarizations being strongly scattered [87]. Since the free layer rotates at an angle under the influence of the applied field, the resistance of the GMR device is proportional to the sine of that angle. The electrical resistance can be determined by running a continuous current through the device and measuring its voltage. Due to their excellent sensitivity, resolution, and low sensor resistance, GMR sensors are highly appealing [95]. This lower resistance also leads to higher operational bandwidth, lower power consumption for a given sensor current, and reduced noise levels [87]. Similar to the GMR sensor, a tunnelling magnetoresistance sensor comprises two ferromagnetic layers. However, instead of being separated by a conductor layer, it is separated by a tunnelling barrier (an extremely thin insulator, typically a few nanometers thick). If the insulating layer is thin enough, electrons can tunnel from one ferromagnet into the other. Since this process is forbidden in classical physics, the TMR is a strictly quantum mechanical phenomenon. Like the GMR sensor, the resistance relies Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 93 on the magnetic vectors’ relative alignment in the free and fixed layers, affecting the flow of conduction electrons between the layers. The parallel configuration results in minimum resistance, and the antiparallel configuration results in maximum resistance. The resistance is measured by passing a current through the layer structure and measuring its voltage. The concept behind the TMR effect was initially demonstrated in 1975 by Julliere. Julliere demonstrated that a junction of two ferromagnetic materials separated by a thin insulating film exhibits the TMR effect [96]. Similar to the GMR sensors, one of the ferromagnetic layers of a typical TMR sensor is free to rotate with the applied field. In contrast, the magnetization of the other layer is pinned. Like GMR sensors, the resistance of the TMR stack is proportional to the cosine of the angle between the free and pinned layers. Consequently, to create a TMR sensor with a linear response to an applied magnetic field, the magnetization of the reference layer must be pinned perpendicular to the free layer. TMR sensors have a higher sensitivity and resistance and consume less power than GMR sensors. However, they suffer from higher noise. TMR sensors are also more expensive and difficult to manufacture since they require a high-quality, pinhole-free, extremely thin tunnelling barrier. 4.1.2 Optical-based position sensor Optical sensors have better resolution than Magnetoresistance-based sensors [97]. Although many structures may be applied, optical-based position sensor present immunity to external magnetic fields [97]. Conventionally, optic sensors are developed using either a shadowing effect or linear gratings. Conventional optic sensors are often discarded from most applications, despite high-resolution requirements, due to their high manufacturing cost from the production of lithographic patterns for encoders [97]. Diffractive optical sensors allow for a chosen distribution of light, and with Fourier diffractive optics, a similar resolution is obtained, but more resistant against vibrations and assembly tolerances [98]. The optical sensors require a full-proof enclosure. This is due to their sensibility to environmental pollution, and thus they are prone to be discarded in most industrial, and automotive applications [97]. 4.1.3 Capacitive-based position sensors Capacitive sensors are an alternative type of sensor that can be used to measure position. They can provide high-precision measurements in environments where magnetic fields may cause interference. In addition, capacitive position sensors have the advantages of simple design, small size, low cost, low power 94 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology consumption [99–102], and superior accuracy compared to magnetoresistance-based sensors [103]. However, they are susceptible to environmental pollution and require special care when operating in dusty environments with moisture and grease [99]. As a general operational concept, these sensors use capacitance to detect displacement. They operate based on a principle similar to that of a parallel-plate capacitor. As one plate rotates/moves, the dielectric is modulated by the rotational movement causing a change in the capacity. Typically, capacitive position sensors employ an array of capacitors that transmit signals based on electric field coupling between fixed and moving plates [104–109]. Their basic structure allows contactless readings but is susceptible to condensation, and electrostatic build-up [103, 109]. In addition to their ability to work without a direct power connection, capacitive sensors typically exhibit lower power consumption than alternative sensors. Although optical encoders are comparatively more accurate, position sensors made with capacitors can gather high-precision and long-range displacement measurements due to the absence of any inherent limitation in capacitor array length. Also, using the electrode planes as the sensing unit significantly reduces the accuracy requirements of the measurement mechanism [104]. Nevertheless, they are vulnerable to environmental conditions, such as moisture and temperature, and present more significant non-linearity due to the electric field bending effect and stray capacitance, affecting accuracy [103,110–112]. Zhang et al. [113] proposed a new capacitive array absolute positioning system that provides two coarse and two fine orthogonal signals, allowing absolute positioning measurement. This intricate device, however, needs several electrical connections. As a result, many of the benefits of capacitive position sensors are lost in the design, which is also unsuitable for portability and downsizing, reducing its applicability in real-world situations. 4.1.4 Conclusions This section reviewed several position sensor technologies, which are vital components in numerous industries such as automotive, robotics, mechatronics, automation, machine tools, avionics, aerospace, medical devices, and consumer electronics. The market demands sensors with higher reliability, performance, and lower power consumption while maintaining competitive manufacturing costs. Furthermore, position sensors must meet strict specifications regarding resolution, accuracy, hysteresis, robustness, and stability. Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 95 Magnetic-based position sensors, including Hall-effect and magnetoresistance-based sensors, are prevalent in the automotive industry due to their robustness, contactless, and non-destructive sensing. Opticalbased position sensors offer high resolution and immunity to external magnetic fields, but their susceptibility to environmental pollution and higher manufacturing costs limit their applications in industrial and automotive sectors. Capacitive-based position sensors provide an alternative solution for environments where magnetic fields may cause interference, featuring simple design, small size, low cost, low power consumption, and superior accuracy compared to magnetoresistance-based sensors. However, they are sensitive to environmental conditions, such as moisture and temperature, and exhibit significant nonlinearity due to the electric field bending effect and stray capacitance. In conclusion, each position sensor technology has its advantages and disadvantages, making them suitable for specific applications and environments. 4.2 State-of-the-art in inductive-based position sensors Inductive sensors are adaptable to a wide range of applications due to their extensive design options. Nevertheless, inductive sensors designed with planar coils stands out above the rest in terms of its application. This design opens the possibility of impressing the coils directly on PCBs or flexible materials, offering solutions with a lower weight, better mechanical stability, and a compact design when compared to solenoidbased [114–126]. Furthermore, since planar coils can be printed on PCBs, they are highly repeatable, predictable, and economically efficient, facilitating assembly and integration processes. Combining the production method with the high reliability of inductive-based technologies delivers robust solutions offering good thermal behaviour, flexible design, and high replication capability for various applications–including automotive, healthcare, robotics and electronics. In addition, PCB-based planar coils are well-suited for various sensing applications due to their device size (smaller for low-power applications), cost-effectiveness (more affordable), and harsh operational conditions (e.g. dust and oil) constraints. There are already devices on the market that utilize this technology in their transduction mechanisms, such as the IPS2200 magnet-free inductive position sensors from Renesas [127]. However, New designs continue to emerge. The remainder of this subsection is devoted to presenting some of these. Anandan et al. [128] proposed a shaft angle sensor, Figure 4.2, which uses four identical flexible planar coils as excitation coils. These four coils, arranged in cylindrical quadrants, induce a voltage on a fifth coil 96 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology Figure 4.2: The flexible planar coil-based sensor for thru-axis angle detection proposed by Anandan et al. Figure 4.3: The schematization of the non-contact angle sensor based on the Foucault current proposed by Kumar et al. placed next to them. This fifth coil, coupled to the shaft, is free to move and its induced voltage, related to its azimuthal angle, is measured and processed to determine the angle. In addition, in the processing phase, an algorithm is needed to identify the quadrant of the measured angle. According to the authors, the sensor has a maximum full-scale error of 1% (equivalent to 3.6°) and a resolution of 0.15°. The excitation coils are fed with a 10 kHz sinusoidal signal. In [129], Kumar et al. proposed a different shaft angle sensor. The sensor consists of a conductive rotating shaft with a groove and a stationary part with four identical flexible coils (see Figure 4.3). The coils are fed with a 60 kHz sinusoidal signal, and there is a 3 mm air gap between the shaft and the coils. As a brief description of its operating principle, the groove modifies the inductance value of the stationary coils according to their position. Thus, these stationary coils, connected to a signal conditioning circuit, determine the circuit’s output voltage, which is proportional to the coils’ inductance change. The authors claim that the sensor has a resolution of 0.08°and a precision of 0.02%. The maximum non-linearity of the sensor is 0.25% for the entire circle range (0.9°). Kumar proposes another shaft angle sensor in [130] (Figure 4.4). The working principle is the same as in the previous work. A surface groove introduced on the shaft changes the inductance of the stationary coils according to their angular position. The difference is the stationary part of the sensor that is composed of two layers of six flexible square-planar coils each. The idea is to use a successive approximation algorithm to determine the quadrant in which the surface groove is present. Then a fine measurement is done Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 97 Figure 4.4: The Foucault current-based angle sensor proposed by Kumar et al. Figure 4.5: The absolute magnetic rotary sensor proposed by Zhang et al. using the coils from the two layers that belong to that particular quadrant. The coils are fed with a 200 kHz sinusoidal signal, and there is a 0.5 mm air gap between the shaft and the coils. According to the authors, the sensor has a resolution of 0.1°and a non-linearity of 0.9%. Zhang et al. [131, 132], proposed an angle sensor composed of two inductors and a code disc (Figures 4.5 and 4.6). Both inductors are embedded with a big planar spiral copper coil and four smaller copper coils using Microelectromechanical System (MEMS) technology (Figure 4.6). Two circles of regular copper sheets are listed on the surface of the code disc. The outer has 177 strips and the inner 176 strips. The excitation signal of the sensor is 4 MHz. According to the authors, the sensor’s resolution is 0.5°. In [133], Kumar proposed another angle sensor with an eccentric disc-shaped rotor and a stator in a PCB with four circular coils (see Figure 4.7). The rotor rotates off-centre such that the inductances of 98 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology Figure 4.6: The working principle of the absolute magnetic rotary sensor proposed by Zhang et al. Figure 4.7: Angle sensor with an eccentric disc-shaped rotor proposed by Kumar et al. the coils change as the angle changes. The air gap between the rotor and the stator is 0.5 mm. In this configuration, the authors claim that the sensor’s resolution and worst-case non-linearity are 0.06°and 0.7%, respectively. The design of a multiperiod bipolar inductive absolute angle sensor is proposed in [134]. The sensor structure comprises a stator and a rotor, where the last is a fan-shaped target with four outer blades and one inner blade, Figure 4.8. The receiving coils are on the stator, and there are two groups: the outer group and the inner group. The outer group comprises two orthogonal sinusoidal-shaped with four cycles that match the fan-shaped target at the outer edge. The inner group comprises two orthogonal sinusoidal-shaped with one cycle corresponding to the centre’s fan-shaped target. The coil with more sinusoidal cycles provides high measurement accuracy, and the coil with a small number of sinusoidal cycles identifies the quadrant of measurement. Additionally, on the stator, there is an excitation coil driven by a 4 MHz sinusoidal signal that induces the Foucault currents on the target sensed by the receiving coils. (As a side note, a more detailed description of its operation is provided in the following subchapter, as this sensor is similar to the ones described in the present work.) Zhang et al. state that the measurement error, the accuracy and the resolution of the sensor are 0.04°, 0.1°, and 0.005°, respectively, in the full range of 360°. Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 99 Figure 4.8: Multiperiod bipolar inductive absolute angle sensor proposed by Zhang et al. Figure 4.9: Planar inductive-based sensor proposed by Gao et al. Goa et al. propose another planar inductive-based sensor in [135] similar to the previous work. In this work, the PCB stator consists of three circular excitation coils and two sets of two orthogonal sinusoidalshaped receiver coils (see Figure 4.9). The rotor, printed on PCB, has two segmented copper circles that match the receiver’s coils’ position. There are 64 segments on the outer circle and 63 on the inner circle, the exact number of sinusoidal cycles in the respective receiver’s coils. In this way, there are two fine sense channels that, together, make an absolute angle sensor. This absolute capability is achieved by representing the two sensing signals (the inner and the outer) with complex numbers where the real and imaginary parts are the values of the respective orthogonal receiver’s coils. Multiplying these two complex numbers together results in a signal with one cycle per 360°, the same as if the sensor had just one sinusoidal cycle shape in the entire 360°, but with the precision of the 64/63 cycles. Exciting the excitation’s coils at 125 kHz, the authors claim that the sensor has 0.000278°of full-range error. In [136], Babu et al. proposed a LIPS using a U-shaped moving magnetic core as the target and stationary 106 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 4.3 Designing inductive-based position sensors Several works have been carried out in recent years to obtain a better performance in the main characteristics of the position sensors, such as safety, reliability, sensitivity, resolution, stability, thermal compensation, energy consumption and size [110,131,144–151]. The position sensors are typically divided into five main groups: optical [152], capacitive [144,152], magnetic [144,152], and inductive [144,152]. A wide variety of position sensors are currently available on the market. Depending on the application, they can indirectly measure other physical quantities converted into motion [153]. The application range of position sensors goes from the automotive industry, robotics, control applications, medical equipment, instrumentation, and military applications to the aerospace industry [131,144,145]. This section only covers the inductive sensors based on Foucault currents since linear and angular optimized position sensors were designed in this thesis. Foucault currents are closed circuits of induced electric current that circulate in a conductive target, usually perpendicular to the magnetic flux (time-varying). In their simplest form, Foucault currents concentrate near the adjacent surface of the excitation coil, and their density decreases with depth. Being Lenz’s law, the principle of operation of Foucault’s current sensors [154]. There are two main versions of Foucault’s current-based sensors for position measurement, frequencybased and amplitude-based (see Figures 4.18a and 4.18b, respectively). The frequency-based version has an oscillator circuit connected to the sensor coils (Figure 4.19), where the frequency of this oscillator is used for the target’s position calculation. Figure 4.20 depict the block diagram of the amplitude-based version. This version includes an excitation coil that generates the magnetic field and receiver coils connected to an Amplitude Modulation (AM) demodulator circuit. The target’s position dictates the amplitude Coil 1 Coil 2 Coil 3 Target (a) Frequency-based Excitation coil Receiver coil 1 Receiver coil 2 Target (b) Amplitude-based Figure 4.18: Angular position sensor Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 107 Figure 4.19: Frequency-based APS block diagram Figure 4.20: Amplitude-based APS block diagram of these AM signals. Both versions of Foucault’s current-based sensors have a coil energized by an alternating current, facing a conductive target. This alternating current generates an alternating magnetic field perpendicular to the direction of the electric current, which in turn induces currents in the conducting metal target, known as Foucault currents. These Foucault currents, also alternating currents, generate a new magnetic field that adds to the original one. The difference between the two types of sensors starts here. In the case of the Foucault current-based sensor version with the oscillator circuit (frequency-based version), the interaction between the two magnetic fields will change the coil’s impedance value [153,155]. In turn, this coil is a passive component of an electronic oscillator circuit, in which the oscillation frequency is related to the inductance value. Thus, by knowing the oscillating frequency, it is possible to calculate the inductance value of the coil. Moreover, since this inductance value is related to the relative position between the coil and the metal target [145, 147, 153, 154, 156], it is possible, as well, to calculate the displacement of the sensor. In the case of the amplitude-based version, which has dedicated receiver coils, the interaction between the two magnetic fields (the field generated by the excitation coil and the 108 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology field generated by the Foucault currents on the target) dictates the induced signal amplitude on the receiver coils. Consequently, by knowing the amplitude of the induced signal related to the relative position between the coil and the metal target, it is possible to calculate the displacement of the sensor. The sensors based on Foucault currents provide a non-contact measurement system, are immune to external magnetic fields, have high performance even in hostile environments, and are robust [147,153, 157,158]. With the increasing number of sensors in a single application, compact, low-cost, and low-power solutions are gaining relevance [110,145,153]. However, the skin effect limits high-frequency Foucault’s current-based sensor’s resolution and stability [150], their usual compact shape and insensitivity to the environment have become very attractive for industrial applications [159–161], leading to their deployment in several automotive systems. 4.4 Proposed methodology for inductive position sensor design This section proposes and describes an automatic geometry generator for inductive position sensors that optimizes the sensor’s coils for high induced current while simultaneously reducing the overall sensor’s non-linearity. Finding the best coils’ geometry configuration is required to accomplish these goals. The automatic geometry generator tool optimizes the sensor’s coils’ geometry for those two goals (more can be added) by performing two optimizations independently. In the first optimization iteration, the tool searches for one possible design that maximizes the induced current in the receiver’s coils and fulfils the user constraints. Typically, those constraints are related to the selected PCB (e.g., stack size and the number of layers), the sensor’s maximum allowed size and operating frequency. Those constraints could be imposed (e.g., the operating frequency), or an allowed range could be given (e.g., the sensor’s coils’ minimum and maximum permitted dimensions, including the target coil). In the second optimization iteration, the automatic geometry generator picks the optimized geometry from the first optimization iteration and harmonically deforms the two receiver coils to minimize the sensor’s non-linearity. This is done by searching the amplitudes and phases of added harmonics components to the receiver’s coils’ forms (this is one of the main contributions of this thesis). The automatic geometry generator tool utilizes a global search algorithm, the Real-Coded Genetic Algorithm (RCGA), and an electromagnetic field simulator to optimize the geometry of the sensor’s coils. In Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 109 -20 -10 0 10 20 mm -20 -15 -10 -5 0 5 10 15 20 mm Excitation coil Receiver coil 1 Receiver coil 2 Thin wire target (a) Angular position sensor -30 -20 -10 0 10 20 30 mm -5 0 5 mm Excitation coil Receiver coil 1 Receiver coil 2 Thin wire target (b) Linear position sensor Figure 4.21: Examples of position sensor geometries the implemented RCGA, a population of individuals, each containing a single chromosome, evolves over a predetermined number of generations. Each chromosome consists of multiple genes, where the entire chromosome represents a sensor coil’s geometry and each gene corresponds to one optimizable parameter. For the electromagnetic field simulator, three methods were implemented: the Method of Moments (MoM), a modified version of the MoM named Fast Harmonic Method (FHM) (detailed in Section 4.5), and the Multilevel Fast Multipole Method (MLFMM). Among these three methods, only one is employed throughout the entire optimization process, with the selecting being made during the tool’s configuration phase. The tool was primarily designed for linear and angular inductive position sensors, Figures 4.21a and 4.21b, respectively. It optimizes the sensor’s geometry (excitation, receivers, and target coils) to achieve higher induced current on the receiver’s coils and lower the sensor’s non-linearity. However, it is not intended to limit the scope of the tool’s applicability. It could be easily extended to other types of inductive sensors and other optimization goals like flux maximization and mechanic tolerances. The flowchart of Figure 4.22 overviews how the automatic geometry generator tool operates. The process begins with collecting the user’s input settings that define the non-optimizable parameters and their respective fixed values. Correspondingly, the tool prompts the user to specify the optimizable ones and their limits. As an option, the user can provide an initial geometry. If provided, it is incorporated into the Genetic Algorithm (GA)’s first generation as an individual. This option allows for improving an existing 110 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology Start Creates the GA initial population with valid random geometry values Last generation? Selection (Tournament Selection) Crossover Mutation (Gaussian Mutation) Elitism (Select the best individuals from the latest and current generation to pass to the next generation) Next generation Select the best individual of the last generation Shows the optimized geometric parameters End Set GA parameters: Number of individuals; Number of generations; Mutation probability; Mutation standard deviation; Elitism ratio. Yes No Fitness (Field Simulation) Field simulator initialization Fitness (Field Simulation) Set the sensor's fixed parameters Set the limits of the sensor's optimizable geometric parameters Set the initial geometry Figure 4.22: Optimization algorithm flowchart geometry. Lastly, the last settings before the optimization process begins are the GA’s settings. These are the number of individuals (geometries) that compose the population, the maximum number of generations that the population will be able to evolve, how often a mutation occurs (mutation probability) and its non-normalized standard deviation, and the elitism ratio that defines the ratio of individuals from the previews generation that will be part of the next one. With this configuration as input, the automatic geometry generator tool creates the first generation of -25 -20 -15 -10 -5 0 5 10 15 20 25 mm -5 0 5 mm Excitation coil Receiver coil 1 Receiver coil 2 Thin wire target Figure 4.23: Geometry generated in the first generation Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 111 -25 -20 -15 -10 -5 0 5 10 15 20 25 mm -5 0 5 mm Excitation coil Receiver coil 1 Receiver coil 2 Thin wire target Figure 4.24: Geometry generated by the optimization tool Table 4.3: Example of fixed and optimizable parameters for induced current maximization Fixed parameters Description Value PCB trace width 165.000 µm PCB copper thickness 18.000 µm Thickness between PCB layers 1 and 2 70.000 µm Thickness between PCB layers 2 and 3 320.000 µm Thickness between PCB layers 3 and 4 70.000 µm Receiver coils’ lenght 50.000 mm Receiver coils’ width 10.000 mm Target airgap 2.000 mm Frequency of the excitation coil’s power source 10.700 MHz Voltage of the excitation coil’s power source 707.107 mV RMS Maximum current of the excitation coil’s power source 7.071 mA RMS Optimizable parameters Description Minimum Maximum Turns in the excitation coil 1 3 Layers in the excitation coil 1 4 Gap between turns of excitation coil 165.000 µm 1.000 mm Excitation coil inner lenght 50.660 mm 101.320 mm Excitation coil inner width 10.660 mm 21.320 mm Target lenght 1.000 mm 50.000 mm Target width 1.000 mm 21.320 mm induced position sensors’ geometries (individuals). If the user supplies an initial geometry as a starting point, it is inserted as the first individual in the first generation, and the remaining individuals are generated randomly. If not, the first individual is also generated randomly. All randomly generated individuals are valid designs (e.g., all geometries where the receiver’s coils invade the excitation coil area and geometries outside the limits of the user-defined optimizable parameters are not allowed). As a curiosity, Figure 4.23 depicts an example of a random geometry generated in the first generation in the second iteration (low non-linearity goal). As an example of a valid use case, the geometry illustrated in Figure 4.24 was achieved through a series of parameters set. In this example, the fixed parameters in the first optimizing iteration (high induced 112 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology Table 4.4: Example of fixed and optimizable parameters for non-linearity minimization Receiver coil Harmonic Minimum Maximum Amplitudes 1 0 -5.000 mm 5.000 mm 2 -5.000 mm 5.000 mm 3 -5.000 mm 5.000 mm 4 -5.000 mm 5.000 mm 5 -5.000 mm 5.000 mm 2 0 -5.000 mm 5.000 mm 2 -5.000 mm 5.000 mm 3 -5.000 mm 5.000 mm 4 -5.000 mm 5.000 mm 5 -5.000 mm 5.000 mm Phases 1 2 -90° 90° 3 -90° 90° 4 -90° 90° 5 -90° 90° 2 2 -90° 90° 3 -90° 90° 4 -90° 90° 5 -90° 90° current goal) were: the PCB trace width, the PCB copper and layers’ thickness, the receiver’s coils’ length and width, the target airgap (airgap space between the receiver’s coils plane and target coil plane), the frequency and voltage of the excitation coil’s power source, and the maximum current allowed in the excitation coil. The optimizable parameters were: the total number of turns in the excitation coil in each layer, the total number of layers that compose the excitation coil, the gap between the excitation coil’s turns’ traces, the excitation coil’s inner length and width (the length and width of the innermost turn), and the length and width of the target’s coil. All of them and their respective values are summarised in Table 4.3. Also, in this example, the optimization tool was fed with an initial geometry, the geometry depicted in Figure 4.21b. In the second optimization iteration (low non-linearity goal), the tool is fed with the geometry generated in the first iteration as initial geometry. This time the optimizable parameters are the amplitudes and phases of the harmonics components to be physically added to the receiver’s coils. These components are the continuous component and all harmonics from 2 through 5. Table 4.4 summarizes the optimizable parameters and depicts their respective limits. As a side note, the fundamental harmonic amplitude is adjusted accordingly to avoid superposition over the excitation coil. As stated before, the global search algorithm implemented is the RCGA to search through all possible solutions for optimal design. In essence, the RCGA algorithm evolves a population of candidate designs, Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 113 where every individual has a chromosome that encodes the optimizable parameters in several genes. These parameters are evaluated against the desired objective of higher induced current or lower nonlinearity, where a fitness function gives a score. Any potential offspring from a pair of parents would inherit their parent’s genes and any new gene they acquired during the evolution process. If the proposed solution did not fit the goal, it will not produce offspring and eventually will die out. The individuals that survived would pass on their genes to further generations. In addition, a mutation may occur. The algorithm continues until it comes across a solution that has met the criteria. Essentially, the RCGA is composed of genetic operators (selection, crossover and mutation) that generate a series of populations (sensor’s coils’ geometries) whose individuals will have evolved to satisfy a particular propose (higher induced current or lower non-linearity). A more detailed description of the GA is presented in Section 2.3.3. The first evolutionary operator applied to the population is the Selection , as shown in the flowchart in Figure 4.22. The primary objective of this operator is to probabilistically advance the most suitable solutions to the next generation, while discarding the less optimal ones. There are various strategies available for selection, one of which is tournament selection. In the implemented tournament selection algorithm, the chromosomes, representing the sensor’s geometries, are randomly reordered into two distinct arrays, each containing the entire population in a different sequence. Chromosomes in corresponding positions within the arrays compete with each other, and the winner is the one with the highest fitness value. This process is repeated for all positions within the arrays. The second evolutionary operator to be applied is the crossover . This operator is the primary search tool for GAs since it combines chromosomes with relevant genetic information, creating a new population. The implemented crossover algorithm randomly selects the genes (optimizable parameters) from parent one (sensor’s geometry) to be inherited by offspring one (crossed sensor’s geometry). At the same time, the remains come from parent two. In turn, offspring two inherits both parents’ rejected genes. Parent one and two are randomly selected from the population and removed to ensure they are not selected again. This process repeats until all chromosomes (sensor’s geometries) have been crossed. The subsequent evolutionary operator is mutation . As described in Section 2.3.3, it is introduced to prevent early convergence and to ensure diversity. The implemented version of this operator modifies the offspring genes (optimizable parameters) randomly by altering their value in compliance with a unidimensional gaussian curve, resulting in a new individual (sensor’s geometry). The user-specified mutation 114 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology standard deviation value determines the gaussian curve, and the mutation probability, which is also userspecified, dictates how frequently this mutation can occur. The elitism operator guarantees that the best preview individuals (sensor geometries) are included in the following generation. This operator is necessary for the most refined individual of the previous generation to be present in the next. The ratio of the individuals of the prior generation and the current generation that will be part of the next generation is given by the elitism ratio (user-specified), and the best individuals from both generations (previous and current) are selected. With the next generation formed, the process repeats until a limited number of generations is achieved (see Figure 4.22). When this occurs, the best individual from the last generation is selected and presented as the optimized sensor geometry. 4.4.1 Fitness function As mentioned previously, the optimization tool has three methods to solve the field equations, the MoM, a modified version of the MoM named FHM (detailed in Section 4.5), and the MLFMM. Each of them is called in the fitness function. The flowchart in Figure 4.25 depicts the implementation of the fitness function that employs the MoM to solve the field equations. This function begins by generating the geometry encoded in the individual’s chromosomal genes, followed by the discretization using perfectly conducting thin wires. After that, it generates the system of linear equations to be solved, as described in Section 2.2.2. The system of linear equations is then solved, using the Generalized Minimum Residue (GMRES) method solver, and stored. Depending on the optimization objective, increasing the induced current in the receiver’s coils or reducing the sensor’s non-linearity, the sensor’s coils’ geometry must be simulated with different target positions. To minimize non-linearity, the sensor must be simulated with the target at several positions to know the sensor’s non-linearity overall the target excursion (linear and angular). Thus, after each simulation, the target is repositioned, and the simulation process is repeated, ending when the target reaches the final position. Finally, the fitness score is calculated according to the optimization objective through the knowledge of the current [i]distribution for different target positions. The flowchart of the fitness function that implements the FHM is depicted in Figure 4.26. Similar to the MoM’s fitness function, this function starts by generating the geometry encoded in the individual’s Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 115 Fitness (MoM Simulation) Generate the geometry Discretizes the geometry using thin wires Solves the matrix equation Calculates the fitness score End [z][i]=[v] Generate the matrices [z][i]=[v] Target at the final position? Store [i] Move the target to the next position Yes No Figure 4.25: MoM simulation algorithm flowchart chromosomal genes, followed by the discretization using perfectly conducting thin wires. After, it generates the system of linear equations in the same way as in the MoM, the only difference being that the impedance matrix [z]and the source matrix [v]are converted to the frequency domain. In this new domain, the higher frequencies are removed to reduce the size of the matrices and, therefore, the processing time needed to solve the system of linear equations. This reduction is only possible because inductive position sensors are typically much smaller than the operating wavelength. The system of linear equations in the frequency domain is then solved using the GMRES solver, and the current [i]is converted back to the original domain and stored. A similar way as in the MoM’s fitness function, the target is also moved to the next position, and the simulation process is repeated until it reaches the last position. Finally, the fitness score is calculated. The last fitness function employs the MLFMM algorithm to solve the field equations. Figures 4.27 and 4.28 show the algorithm’s initialization and the fitness function flowcharts, respectively. The MLFMM initialization algorithm initializes the octal tree structure and related matrices. As the structure of the octal tree and the matrices related to it are independent of the sensor geometry, they only need to be computed once. In this way, all matrices that meet these criteria are initialized in the Figure 4.27 flowchart. These are the radially oriented unit vectors on the sphere (ˆ k), the interpolation matrices (W), and the transfer functions (H) of all possible combinations of distinct values of ˆ kand rab. All transfer 122 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology Subsequently, in the second iteration, the optimized geometry from the first iteration was used as the starting point. Furthermore, the second iteration was performed independently twice, each time with different objectives. The first objective aimed to minimize the sensor’s non-linearity, while the second one focused on reducing the harmonic content. Again, both objectives used the optimized geometry from the first iteration. Further details are discussed in the following sections. Since three electromagnetic simulators were developed, the remaining of this section describes the results obtained with each one. 4.6.2 Optimization using the MoM This section summarizes the main results obtained with the optimization tool implemented with the MoM to solve the electromagnetic field equations. The optimization tool was used to optimize the geometry of an APS, consisting of an excitation coil, two receiving coils and one target coil, with a measurement range of 360°. The excitation coil and both receiver coils were printed on the same PCB, defined as the stator, while the target coil was printed on a separate PCB identified as the rotor. The goals were to maximize the induced current in both receivers’ coils and minimize the APS non-linearity. To accomplish this, the optimization procedure was divided into two iterations. The first iteration used the optimization tool to optimize the sensor’s geometry to obtain a maximization of the induced current on both receivers’ coils. The second iteration takes the optimized geometry obtained in the first iteration, and uses it as the starting point to minimize the sensor’s non-linearity, using the same optimization tool. In addition, the second iteration was repeated 14 times, with different fitness functions and GA parameters. The goal was to determine which fitness functions and GA parameters best fit the objective in terms of convergence and geometric solution. A diagram of the optimization procedure can be seen in Figure 4.30. Geometry 1 , depicted in Figure 4.31, is the initial human-made geometry most suitable for producing acceptable results. This geometry was the starting point to maximize the induced current, and Geometry 2 results from geometry optimization. Geometry 3 is the optimized version of Geometry 2 when maximum nonlinearity is used as a fitness parameter to be optimized. Geometry 4.x ( Geometry 4.x refers to 12 geometry optimizations, from Geometry 4.1 through Geometry 4.12 , that employ the same fitness function but different GA parameters) and 5 begin with the same geometry as Geometry 3 but with extra optimizable parameters. As in Geometry 3 , the maximum sensor non-linearity is the fitness parameter for Geometry 4.x . However, rather than the maximum non-linearity of Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 123 Figure 4.30: Optimization flowchart 0 2 mm 20 10 20 10 mm 0 mm 0 -10 -10 -20 -20 (a) Perspective view -20 -10 0 10 20 mm -20 -15 -10 -5 0 5 10 15 20 mm Excitation coil Receiver coil 1 Receiver coil 2 Thin wire target (b) Top view Figure 4.31: APS initial geometry - Geometry 1 eight equally spaced angular positions, like in Geometry 3 , Geometry 4.x has 16. Finally, Geometry 5 has the same optimizable parameters as Geometry 4.x . The only difference is the fitness parameter, which is the harmonic content in the two receivers’ coils signals. The goal is to reduce the harmonic content and, consequently, the non-linearity. A more detailed description of each optimization and respective results are given in the following subsections, one for each optimization. 4.6.2.1 APS peak-to-peak induced current optimization This subsection describes the use of the optimization tool to generate a geometry that maximizes the induced current in the sensor receiver’s coils. As previously stated, this process begins with Geometry 1 as a starting point. It is necessary to select which geometric parameters will be fixed and which will 124 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology Table 4.5: APS optimization using MoM - Limits and initial parameters for induced current maximization Fixed parameters Description Value PCB trace width 165.000 µm PCB copper thickness 18.000 µm Thickness between PCB layers 1 70.000 µm and 2 Thickness between PCB layers 2 320.000 µm and 3 Thickness between PCB layers 3 70.000 µm and 4 Receiver coils’ outermost radius 20.250 mm Target airgap 2.000 mm Frequency of the excitation coil’s 10.700 MHz power source Voltage of the excitation coil’s 707.107 mV RMS power source Maximum current of the excitation 7.071 mA RMS coil’s power source Optimizable parameters Description Initial Minimum Maximum Optimized Geometry 1 Geometry 2 Turns in the excitation coil 1 1 3 1 Layers in the excitation coil 4 1 4 4 Gap between turns of excitation coil 335.000 µm 165.000 µm 1.000 mm 732.212 µm Excitation coil inner radius 21.500 mm 1.000 mm 50.000 mm 20.428 mm Receiver coils’ innermost radius 14.500 mm 1.000 mm 50.000 mm 1.000 mm Target inner radius 14.500 mm 1.000 mm 50.000 mm 1.000 mm Target outer radius 21.750 mm 1.000 mm 50.000 mm 20.320 mm Target angular wing size 180.000° 0.000° 360.000° 205.363° Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 125 be optimizable, as well as their upper and lower limits in the case of the optimizable parameters. As fixed geometric parameters, the PCB trace width, the PCB copper thickness, the PCB layers stack, the receiver’s coils’ outermost radius and the conductive target airgap (airgap between the stator and rotor PCBs) were imported to the tool. Some were selected for obvious reasons, such as copper thickness and PCB layers’ stack, resulting from PCB fabrication constraints. Others are due to mechanical constraints, such as the receiver’s coils’ outermost radius and the target airgap. For optimizable parameters, eight parameters were selected: the total number of turns in the excitation coil in each layer, the total number of layers that compose the excitation coil, the gap between the excitation coil’s turns’ traces, the excitation coil’s inner radius (the radius of the innermost turn), the receivers coils’ innermost radius (the shortest distance that receiver’s coils’ traces has to the sensor’s centre), the inner and outer radius of the target’s coil, and the angular size of the target. Table 4.5 summarizes these parameters and shows the upper and lower limits, as well as the initial values of the optimizable parameters used in Geometry 1 (initial geometry). After selecting and configuring the fixed and optimizable parameters, the fitness function and GA parameters must also be coded and configured, respectively. The fitness function is made up of MoM simulations, in which the APS’s coils (the excitation coil together with the two receiver coils and the target coil) are all immersed in a medium with a relative electrical permittivity and a relative magnetic permeability of one. Also, the excitation coil is powered by a differential potential of 0.7071 V RMS at 10.7 MHz, being the maximum induced current on the receiver’s coils used as a fitness value. For this, two MoM simulations were performed, being the target positioned 180° apart in the two areas where the induced current is maximum. The fitness value is calculated by adding those two maximums values (peak-to-peak). Also, the fitness function sets the fitness value to zero if the optimizable geometric parameters yield a physically impossible geometry or if the excitation current exceeds a threshold value, which was set to 7.071 mA RMS . Finally, the GA parameters are configured to produce a population of 100 individuals that will evolve over 100 generations with a mutation probability of 1/8. This implyes that, on average, one gene per chromosome will be mutated, as there are a total of eight genes per chromosome (eight configurable parameters per individual or geometry). The non-dimensionalized standard deviation mutation is 1/30, and the elitism ratio is 1/2. This value indicates that half of the individuals from the previous generation, the best ones, will be part of the next generation. Being the optimization tool fully configured as described, the optimization procedure took approximately 126 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 10 20 30 40 50 60 70 80 90 100 Generation 0 0.5 1 1.5 2 2.5 3 Max. peak-to-peak induced current (mA) Figure 4.32: APS optimized for maximum induced current - Convergence curve 0 2 mm 20 10 20 10 mm 0 mm 0 -10 -10 -20 -20 (a) Perspective view -20 -10 0 10 20 mm -20 -15 -10 -5 0 5 10 15 20 mm Excitation coil Receiver coil 1 Receiver coil 2 Thin wire target (b) Top view Figure 4.33: APS optimized for maximum induced current - Geometry 2 14 wall time hours to optimize Geometry 1 and get Geometry 2 on a high-end 2014 workstation with 64 GB of Random-Access Memory (RAM). Figure 4.32 shows the convergence curve, which depicts the evolution of the best individual’s fitness value (highest peak-to-peak induced current) over 100 generations, and Figure 4.33 illustrates the optimized geometry, Geometry 2 . Geometry 2 reached a fitness value of 2.908 mA (peak-to-peak), whereas Geometry 1 has 0.941 mA (peak-to-peak), indicating that Geometry 2 has approximately a three times stronger induction current than Geometry 1 . 4.6.2.2 APS non-linearity optimization (amplitudes) The next step is to reduce Geometry 2 ’s maximum non-linearity after optimizing the induced current in the receiver’s coils. Geometry 2 , with a maximum non-linearity of 0.663°, is used as the starting point. The goal is to see if the optimization tool can discover a solution that considerably reduces sensor nonlinearity by introducing some harmonic distortion into the receiver’s coils design. Because the receiver’s coil’s signals had a significant harmonic amplitude in the continuous component and the third and fifth harmonics, those components were employed to generate the geometric harmonic distortion. In other Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 127 Table 4.6: APS optimization using MoM - Limits, initial and optimized parameters for non-linearity optimization (amplitudes) Receiver Harmonic Initial Minimum Maximum Optimized value coil Geometry 2 (Geometry 3) Amplitudes 1 0 0.000 mm -19.250 mm 19.250 mm 0.000 mm 1 19.250 mm - - 16.576 mm 3 0.000 mm -19.250 mm 19.250 mm 1.140 mm 5 0.000 mm -19.250 mm 19.250 mm -1.534 mm 2 0 0.000 mm -19.250 mm 19.250 mm -0.007 mm 1 19.250 mm - - 17.762 mm 3 0.000 mm -19.250 mm 19.250 mm 0.115 mm 5 0.000 mm -19.250 mm 19.250 mm -1.366 mm words, a geometric harmonic distortion was created in the continuous component and in third and fifth harmonics to counterbalance the presence of those on the receiver’s coil signal. Table 4.6 depicts the optimization parameters, including their initial values as well as their upper and lower limits. It is worth noting that the fundamental harmonic amplitude is lowered proportionally to accommodate the existence of the additional harmonics and avoid geometric constraints. The maximum non-linearity that the sensor can have in a single complete turn was used as fitness value. For this, eight MoM simulations were executed in the fitness function, being the target position in each of them angularly equally separated. The angle detected by the sensor is calculated using Equation 4.2, where αis the measured angle, and xand yare the amplitude signals from receiver coils 1 and 2, respectively. The non-linearity is defined as the difference between the true angular target position and the measured angle, with the fitness value (the maximum non-linearity) being the highest of the eight nonlinearities. It is also important to note that the sensor’s zero location was chosen such that the maximum non-linearity was as minimal as feasible. This was accomplished by subtracting the average non-linearity value from all non-linearity values. α=                            tan−1y x, if x > 0 tan−1y x+π , if x < 0and y ≥0 tan−1y x−π , if x < 0and y < 0 +π 2, if x = 0 and y > 0 −π 2, if x = 0 and y < 0 (4.2) 128 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 10 20 30 40 50 60 70 Generation 0.3 0.4 0.5 0.6 Maximum non-linearity (º) Figure 4.34: APS optimized for minimal non-linearity - Convergence curve 0 2 mm 20 10 20 10 mm 0 mm 0 -10 -10 -20 -20 (a) Perspective view -20 -10 0 10 20 mm -20 -15 -10 -5 0 5 10 15 20 mm Excitation coil Receiver coil 1 Receiver coil 2 Thin wire target (b) Top view Figure 4.35: APS optimized for minimal non-linearity - Geometry 3 For the settings chosen for each GA parameter, the GA population, like the receiver’s coils induced current optimization, consists of 100 individuals who will evolve over 70 generations. The non-dimensionalized standard deviation mutation and elitism ratio are again adjusted at 1/30 and 1/2, respectively. The only difference is that the mutation probability is set to 1/6. This ensures that on average, one gene per chromosome gets modified since there are six genes per chromosome (six optimizable parameters per individual or geometry). The optimization tool took approximately 142 wall time hours on the same workstation used to optimize Geometry 2 , and Geometry 3 with the specifications mentioned in this subsection was obtained. Figure 4.34 depicts the convergence curve, which shows the evolution of the best individual’s fitness value (highest non-linearity) over 70 generations, and Figure 4.35 depicts the optimized geometry, Geometry 3 . Figure 4.35 shows that an unexpected geometry, a more ”squared” geometry, was achieved. Despite being more ”squared” , Geometry 3 has less than half the non-linearity of Geometry 2 . Geometry 2 had a maximum non-linearity of 0.663°, while Geometry 3 had a maximum non-linearity of 0.283°. Note that the convergence curve employs the maximum value among the eight MoM simulations performed within Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 129 Table 4.7: APS non-linearity optimization (amplitudes and phases) - GA parameters Geometry Population Mutation Probability Non-dimensionalized standard deviation mutation 4.1 50 1/10 1/30 4.2 100 1/10 1/30 4.3 150 1/10 1/30 4.4 50 1/100 1/30 4.5 100 1/100 1/30 4.6 150 1/100 1/30 4.7 50 1/1000 1/30 4.8 100 1/1000 1/30 4.9 150 1/1000 1/30 4.10 50 1/1000 1/10 4.11 100 1/1000 1/10 4.12 150 1/1000 1/10 the fitness function. The non-linearity values shown here are the true maximum values. This explains the discrepancy between the non-linearity of the convergence curve and the non-linearity presented here. The same goes for the following subsections. 4.6.2.3 APS non-linearity optimization (amplitudes and phases) This subsection discusses another optimization that takes the same geometry, Geometry 2 , as a starting point and the same fitness parameter, the sensor’s maximum non-linearity. The only changes in the fitness function are the total number of MoM simulations performed for each fitness function execution, which is 16 instead of the eight used in Geometry 3 , and the addition of four more optimizable parameters. This allows the adjustment of the initial phase of the receiver’s coils’ third and fifth harmonics. By introducing those changes, the four additional optimizable parameters and twice the MoM simulations extend the optimization processing time. However, as demonstrated later, it can discover a better solution with a substantially smaller maximum non-linearity. Along with those changes, and differently from the previous section, there are 12 specific optimizations named Geometry 4.1 to Geometry 4.12 . These optimizations have the same fitness function but different GA parameters, namely mutation probability, non-dimensionalized standard deviation mutation and population. The other two GA parameters, the elitism ratio and the total number of generations that the population will evolve, are the same across all optimizations, and they are set to 1/2 and 100, respectively. Briefly stated, Table 4.7 contains all GA parameters for all geometric optimizations. 130 Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology Table 4.8: APS optimization using MoM - Limits, initial and optimized parameters for non-linearity optimization (amplitudes and phases) Receiver Harmonic Initial Minimum Maximum Optimized coil Geometry 2 Geometry 4.2 Amplitudes 1 0 0.000 mm -19.250 mm 19.250 mm 0.003 mm 1 19.250 mm - - 18.695 mm 3 0.000 mm -19.250 mm 19.250 mm 0.232 mm 5 0.000 mm -19.250 mm 19.250 mm -0.154 mm 2 0 0.000 mm -19.250 mm 19.250 mm -0.003 mm 1 19.250 mm - - 18.379 mm 3 0.000 mm -19.250 mm 19.250 mm 0.570 mm 5 0.000 mm -19.250 mm 19.250 mm -0.133 mm Phases 13 0° -90° 90° 7.440° 5 0° -90° 90° -4.131° 23 0° -90° 90° 4.732° 5 0° -90° 90° 13.301° 10 20 30 40 50 60 70 80 90 100 Generation 0 0.1 0.2 0.3 0.4 0.5 0.6 Maximum non-linearity (º) Figure 4.36: APS optimized for minimal non-linearity - Convergence curve The optimization tool required approximately an average of 107, 213 and 321 wall time hours to work with a population of 50, 100 and 150, respectively, on a high-end 2020 workstation with 128 GB of RAM to optimize Geometry 2 , and obtain the respective optimized geometry with the specifications listed in this subsection. Figure 4.36 depicts the convergence curve of the best geometric optimization, Geometry 4.2 , among the 12, which shows the evolution of the best individual’s fitness value (highest non-linearity) over 100 generations. Figure 4.37 depicts the respective optimized geometry. As a result, Geometry 4.2 achieved 0.021° of maximum non-linearity, saying that Geometry 4.2 has non-linearity reduced by 32 times when compared to Geometry 2 , which has 0.663° of maximum non-linearity. Table A.1 in Appendix A shows the geometric results (the values of the optimizable parameters) for all 12 optimized geometries, and the non-linearity values are depicted in Table A.2. Chapter 4. Electromagnetic waves application - Inductive position sensor design methodology 131 0 2 mm 20 10 20 10 mm 0 mm 0 -10 -10 -20 -20 (a) Perspective view -20 -10 0 10 20 mm -20 -15 -10 -5 0 5 10 15 20 mm Excitation coil Receiver coil 1 Receiver coil 2 Thin wire target (b) Top view Figure 4.37: APS optimized for minimal non-linearity - Geometry 4.2 Table 4.9: APS optimization using MoM - Limits, initial and optimized parameters for harmonic content optimization Receiver Harmonic Initial Minimum Maximum Optimized coil Geometry 2 Geometry 5 Amplitudes 1 0 0.000 mm -19.250 mm 19.250 mm 0.002 mm 1 19.250 mm - - 18.629 mm 3 0.000 mm -19.250 mm 19.250 mm 0.350 mm 5 0.000 mm -19.250 mm 19.250 mm 0.105 mm 2 0 0.000 mm -19.250 mm 19.250 mm -0.003 mm 1 19.250 mm - - 18.611 mm 3 0.000 mm -19.250 mm 19.250 mm 0.333 mm 5 0.000 mm -19.250 mm 19.250 mm -0.138 mm Phases 13 0° -90° 90° 3.462° 5 0° -90° 90° -18.951° 23 0° -90° 90° 3.330° 5 0° -90° 90° 89.994° 4.6.2.4 APS harmonic content optimization Another strategy to minimize the sensor non-linearity involves reducing the amplitude of the harmonics present in the receiver coils’ signals. Essentially, this approach uses the same initial geometry ( Geometry 2 ) the same optimizable parameters, upper and lower limits values, and GA configurations as Geometry 4.2 . The only difference is the fitness function that uses the harmonic content instant of maximum nonlinearity as a fitness value. In this optimization, the fitness value is calculated by adding the continuous component and the third and fifth harmonics amplitudes of the two receiver’s coils signals together. As previously stated, the goal is to reduce the continuous component and the third and fifth harmonics amplitudes in the receiver’s coils’ signals and, consequently, the sensor’s non-linearity. Using the same workstation, the optimization procedure took around 394 wall time hours. Figure 4.38