Modeling and Designing of an All-Digital Resonate-and-Fire Neuron Circuit

Integrate-and-fire (IAF) and leaky integrate-and-fire (LIF) models are the popular models for spiking neurons and spiking neuron networks (SNN). They lack the dynamic properties of the ordinary differential equations (ODEs) model but are simpler. More than half of all neuron implementations were digital circuits. The majority of them were IAF and LIF models. The resonate-and-fire (RAF) model proposed by Izhikevich has both simplicity and dynamic properties. However, the RAF model does not have general structures for hardware implementation. Most realizations of the RAF neuron models were analog designs. In this study, we developed and simulated the model’s equations with digital pulses. Based on this analysis, we proposed block models of digital and all-digital RAF neurons to turn the models into digital hardware. We verified the properties of the proposed models by implementing an all-digital RAF neuron on an Intel 28 nm Cyclone V Field-programmable gate array (FPGA) device. The register-transfer level (RTL) structure of the implementation consumed fewer resources, with just 23 D flip-flops, without floating-point operations or multipliers.


I. INTRODUCTION
Integrate-and-Fire (IAF) and Leaky Integrate-and-Fire (LIF) neuron models have been mentioned in thousands of research papers and studies. Although they were invented more than 100 years ago, their simplicity allowed scientists to simulate and design AI hardware easily. A conductance-based model was invented by Hodgkin-Huxley (1952), which was considered the most complex and biophysically plausible model. Other models were proposed from 1960 to 2003, such as FitzHugh-Nagumo (1961,1962), Morris-Lecar (1981), Hindmarsh-Rose (1984), and Izhikevich (2003). All of them tried to reduce the complexity of the Hodgkin-Huxley model for software simulation and hardware implementation. In this study, we investigated over 100 papers from 2000 to 2022 that made spiking neuron models into integrated circuits or used FPGAs to emulate them: The associate editor coordinating the review of this manuscript and approving it for publication was Ludovico Minati . 1) IAF model: [1] to [27]. 2) LIF model: [28] to [64].
3) FitzHugh-Nagumo model: [65] to [71]. 4) Hindmarsh-Rose model: [72] to [79]. 5) Morris-Lecar model: [80] to [81]. 6) Spike Response Model (SRM): [82] to [83]. 7) Izhikevich model: [84] to [93]. 8) Others: [94] to [99]. As a result, a percentage diagram of popular models that were turned into integrated circuits or emulated on FPGAs is shown in Fig. 1. Fig. 1 depicts that more than fifty percent of research in the first two decades of the 21 st century used IAF and LIF models, which are 100-year-old models. We also classified hardware types in papers from [1] to [99] into three groups as in Table 1 below: In Table 1, the ''Others'' column groups publications that did not use conventional CMOS processes. These papers were implemented based on the properties of special materials  like RRAMs [27], memristors [62], biristors [61], certain FET devices [63], and even photonic components [96] to make spiking neuron cells. Even so, these special materials still needed analog CMOS circuits like amplifiers and comparators to work. The two simple models, IAF and LIF, had a ratio of digital implementation greater than fifty percent. On the other hand, most analog circuits were used to implement the FitzHugh-Nagumo, Hindmarsh-Rose, and Morris-Lecar models, which employed ordinary differential equations (ODEs). This research was only focused on realizing the neuron model, so we counted as mixed-signal designs, which used analog circuits for neuron cells, as a part of the analog CMOS group. Based on the classification in Table 1, the overall percentages of hardware types are shown in Fig. 2. Fig. 2 shows that digital implementation is the most common method to deploy spiking neuron hardware. Although many researchers stated that analog circuits were efficient in terms of power and chip area for spiking neuron cells [6], [31], the performance of this kind of hardware might be affected by device mismatches [97], leakage power [38], scaling technology [31], noise, and Process-Voltage-Temperature (PVT) variations [93]. In addition, developing analog circuits for spiking neuron cells usually takes longer than developing digital circuits. Therefore, digital implementation and FPGA emulation have been preferred in recent years. However, most of them implemented the IAF and LIF models. Because the IAF and LIF models were very simple and were invented at the dawn of neuroscientific theory, they could not mimic many dynamic properties like neurons in the human brain, for example, subthreshold oscillation and selective communication. So, spiking neuron networks based on these types of models could not take advantage of the dynamic benefits of biological neuron networks. There is a need for a new model which is simple and biophysically plausible.
In 2001, Izhikevich suggested a ''resonate-and-fire'' neuron model (RAF) [100], which added three dynamic properties to traditional IAF and LIF models. He also mentioned that the RAF model is ''an analogue'' of the IAF model. That means the RAF neuron did not only have powerful computation like the Hodgkin-Huxley model but also kept the simplicity of the IAF/LIF model. From 2001 to 2022, there were only three published integrated circuits that used this model: [97], [98], and [99], two of them are analog circuits [97] and [98]. However, it did not mean that this model was not good enough. This model was still being analyzed and simulated in many recent publications such as [101], [102], and [103]. The research [101] found that the resonating property of a neuron can optimize the learning progress of SNNs. Aside from that, recent papers have just focused on emulating the Izhikevich model (IZK) [104] on FPGAs such as [99], [90], [91], [92], and [93]. This is because the Izhikevich model was the latest neuron model, which could efficiently simulate many types of neuron behavior on a computer, including resonator behavior. Designs in these papers were typically realized using clock-driven algorithms. In order to turn the ODEs of this model into hardware, the papers simplified multiplying operations in these ODEs by using shifters [99], Euler method [89], [90], and CORDIC algorithm [91], [93]. To store the model parameters, which were real numbers, fixed-point numbers were used to replace floating-point numbers [92]. Furthermore, some dedicated finite state machines (FSMs) were also synthesized to control the model operations [88]. All of these techniques not only required many hardware resources for one neuron cell but also consumed much power due to the dynamic power of VOLUME 11, 2023 the clock-driven method. For these reasons, the original RAF model was more hardware friendly than IZK and other conductance-based models. It was a simple model similar to the IAF/LIF model but still mimicked important dynamic features of analog neurons. In this work, our study contributes to 1) Mathematical model of a digital RAF neuron with the pulse train as the damped subthreshold oscillation in a complex domain over time. 2) Block models of digital RAF and all digital RAF neurons. 3) A sample-and-hold method to realize the digital RAF neuron model into digital hardware without using multipliers, adders, and real numbers. 4) A detailed digital circuit design of the all-digital RAF neuron, which has a damped subthreshold oscillator and zero dynamic power consumption at its resting state.
To present these contributions, the contents are organized as the structure below: • Section I presents a short survey on the hardware implementation of neuron models and the existing gap is introduced as in the above content.
• Section II describes the principle operation and equations of the RAF model, then introduce related studies that turned this model into hardware.
• Section III shows a mathematical model of a digital RAF model in a complex domain. The model has a clear definition of threshold level and all-or-none spikes.
Based on the mathematical model, we propose block models of digital RAF neurons, which can be turned into not only an analog circuit model but also digital hardware.
• Section IV presents a proposed model of all-digital RAF neurons based on a true digital pulse train. Simulation results in R language are presented to clarify the operations of the proposed model.
• Section V depicts the proposed model on digital circuits and an RTL structure to turn the all-digital RAF model into digital hardware using basic logic gates and D flipflops (DFF). The concept of using a sample-and-hold circuit to emulate the resonance mechanism is discussed in this section.
• Section VI analyzes the simulation results of the Verilogbased design. This section also shows measurement results on a DE10-nano board, which is an FPGA development kit for an Intel 28 nm Low Power (LP) Cyclone V FPGA 5CSEBA6U23I7. At the end of this section, we compare the resource utilization of our proposed design to other publications in recent years.
• Section VII gives conclusions and our future plan based on this study.

II. RAF MODEL AND RELATED WORKS
According to [100], an RAF neuron was a type of neuron whose membrane potential vibrated as a damped subthresh- On the RAF neuron, the action potential was based on the relationship between the stimulus interval and the damped subthreshold oscillation of the membrane potential.
old oscillation at the first impulse and generated an action potential if the second impulse came at the first or second high state of this oscillation. This behavior is visualized in Fig. 3. The oscillatory activity of an RAF neuron over time can be described by (1): with z ∈ C is the state of the neuron, it has current part x and voltage part y. When the neuron receives the first impulse, membrane voltage starts vibrating in a damped subthreshold oscillation with frequency ω and damping constant b, which is caused by membrane impedance. Because of using complex variables and differential equations, (1) does not tell clearly how to make circuits from it. Simulations in [97] implemented Lotka-Volterra oscillator circuit, which employed subthreshold MOSFETs. This oscillator played the role of a damped subthreshold oscillator for the design. In the design, the oscillator circuit depended on thermal voltage, and its threshold depended on the supply voltage. Additionally, the mismatch between MOSFET devices had an impact on the performance of this design. Therefore, the analog RAF neuron might be sensitive to PVT variations. The device mismatch problem might be fixed by using large-size transistors, but that was not suitable for VLSI systems. Another concept in [98] simply allowed integration to occur at the high state of subthreshold oscillation and reset the neuron state at the low state of this oscillation. The concept might be similar to the RAF model, but it did not care about the interval between two impulses. Then it could not have the selective communication of an RAF neuron.
In 2003, Izhikevich suggested universal ordinary differential equations to describe the action potential of neurons: 62320 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
where v is membrane potential, u is membrane recovery variable. The equation will simulate a resonator if a = 0.1, b = 0.26, c = −60, and d = −1. With only two variables and no complex variables, the equation became famous. This model was labeled ''the simple model'' by Izhikevich, and it was also termed the ''Izhikevich model'' by other researchers.
A group of Japanese researchers came up with the idea of using two parallel shifters and a discrete state method to emulate RAF behavior in the Izhikevich model [99]. They named it a ''rotate-and-fire digital spiking neuron'' (RDN). In their design, an N-cell circular shifter represented states of the membrane potential, while an N-cell circular shifter worked as the membrane recovery variable that was described in the IZK model. Each shifter always had one cell whose value is '1' and the other cells whose value was '0'. The value '1' was shifted around all cells in each shifter at every internal clock cycle. A dedicated five-state controller was required to manage the shifting directions of the two shifters. Some reconfigurable wires were implemented to construct the reset function between the M-cell shifter and N-cell shifter. This combination divided the membrane potential oscillation and membrane recovery wave into multiple time steps. As a result, the FPGA implementation of the design used many types of both combinational and sequential components, including subtractors, counters, many adders, and lots of comparators. The advantage of the design was that it could emulate not only resonator but also up to ''14 types of neuron-like responses''. Thus, the design could be considered a realization of the IZK model, not the RAF model only. Aside from that, the design used a lot of dynamic power because the bits changed state periodically. Four flip-flops did this every internal clock cycle. This property differs from that of the IAF, LIF, and RAF models. At resting state, neurons consumed almost no dynamic power due to the ion concentration equilibrium between the intracellular and extracellular medium. In the next section, we are going to analyze the original equation of the RAF model in a complex domain over time for digital systems and propose our concept to turn the model into digital hardware.

III. DIGITAL RAF MODEL A. MATHEMATICAL MODEL
In Section I and II, the IZK model is the latest neuron model that has RAF behaviors. Most of the recent publications tried to reduce multiplying operators and floating-point numbers in implementation but it was still not efficient in power and hardware resources. Therefore, our work avoids using the IZK model and focuses on the original equation of the RAF model, which is shown in (1). Membrane potential oscillation is represented as the variable y. This oscillation is a damped oscillation. Typically, the damped oscillation is described with a sine wave function. In the view of digital design, a damped pulse train is used for this oscillation. Using Fourier Series, the oscillation is described as in (3): where A is the amplitude of membrane potential, the value of the amplitude depends on system characteristics. In analog hardware and computer simulation, the range of membrane potential is usually the same as a biology action potential, which is from −65 mV to +30 mV according to Nernst's equation [105]. In digital hardware, this amplitude should be based on logic voltages such as +5 V, +3.3 V, +1.8 V, etc. d is the duty cycle, ω is subthreshold frequency, b is damping constant that should be less than 1, m is the number of pulses, t is time, t n is the moment of the last impulse and If ω = 0 and t ≫ t n , the y(t) = 0. It means the neuron's membrane potential is zero at the resting state. It is not the same as a biological neuron's membrane potential, which is normally negative voltage. To let the y(t) have a constant negative voltage at resting level but still be simple with fewer parameters, we add an offset constant (−A/2) to (3), as follows In (5), if ω = 0 and t ≫ t n , the y(t) = −A/2. Based on conductance models, the current part through the cell membrane is an equation of membrane voltage as (6) where I is the total current of inward and outward of ion flows, C is cell capacitance and V is membrane voltage. Therefore, the x part of (1) can be written as (7) Based on [100] and [108], we assume that if there is any injected current while the membrane potential is oscillating, the state of the neuron is described by the following equatioṅ  where J injected ∈ C is the injected state. To get the equation of the total membrane potential over time, (8) can be solved as From (9), the solution of J injected dt must contain both the current part and the voltage part. If it just has a current part for the injected current, the component y(t) of z(t) is totally not affected by the injected current. Hence, we assume where I injected is a constant DC value of injected current, V (t − t n ) is the voltage integrated into the neuron's body at the moment t n . So, the equation of neuron state over time is rewritten as In (11), we see that the neuron is having two parallel operations: • Damped subthreshold oscillation of membrane potential, which is caused by Na + , K + , and leaky ion channels.
• Integration of current impulses, which is caused by Na + ion channels. The action potential, which is the neuron's response to current impulses, is the sum of these operations. Based on the operations, we build up a block model of the RAF neuron model, as shown in Fig. 4.
In Fig. 4, the damped oscillator plays the role of the neuron's dendrite membrane, and the integrator is the block of (10), which plays the role of the neuron's body. The sum function inside the model represents an axon hillock, which adds up all the electrical signals and makes action potentials. To find the potential function of the integrator, which is V (t − t n ), we employ the superposition method: assume ω = 0 and t ≫ t n , (11) in this case will be If we also use a complex variable with current and voltage parts to describe the neuron state of the IAF/LIF model, (12) is showing the IAF/LIF neuron model in a complex domain. As a result, the relationship between the real part and the imaginary part can be described using (6). So, when ion gates of the Na + channel (the inward flows) open, the membrane potential is stored in the neuron's body (the integrator), and it will leak based on the capacitor discharge equation (13) where V C (t) is the integrator's potential, and b is the neuron's time constant, which is membrane impedance. We consider the time constant b as the same as the damping constant of the neuron's membrane potential oscillation. V n is the initial voltage, which is the potential charged to the cell's capacitor at the moment t n . At the time t n , we assume that the body of the neuron adds the oscillation potential to the potential of the integrator by using (11). Therefore, the charged potential of the integrator is where τ is the interval of the injected current I injected . Thus, the action potential of the neuron is (15) includes both damped pulses and the exponent component of an integrator. It means the V A has both subthreshold oscillations and spikes. Based on the principle of dynamic polarization [106], the potentials are propagated along the neuron's axon to synaptic terminals. After this progress, the subthreshold potentials are removed from the action potentials. As a result, the action potentials become all-or-none spikes [107], which can be considered digital pulses. According to [109], a resonator does not have allor-none action potentials like an integrator. Depending on its eigenfrequency, the RAF model can be a resonator or an integrator. If its eigenfrequency is zero, an RAF neuron becomes an integrator. So, it may have all-or-none spikes in the case of an integrator but does not have that property in the case of a resonator. Subthreshold oscillation on the action potentials of a resonator is not useful for digital communication. Therefore, we suggest that a digital RAF model (DRAF) will have all-or-none spikes in both cases. To achieve this feature, we add a logistic function to convert the total membrane potentials to digital pulses, as follows where V DA are all-or-none spikes, k is the logistic growth rate. V DA keeps the subthreshold potentials if k is small such as k = 1 and k = 2. V DA is digital-like pulse if k is very large. V TH is threshold level. From Fig. 3, we assume that the threshold level V TH is above the sum between the maximum peak of the damped subthreshold oscillation and a part of the integrator's 62322 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  potential. According to the assumption, we obtain the V TH in this model by using (14) and (15) From (5), (7), (11), (13), (14), (16) and (17), we propose a block model of a digital RAF neuron (DRAF), as shown in Fig. 5. The equation system of the model is depicted in (18).
(18) has two types of parameters: fixed parameters and configurable parameters. The values of the fixed parameters (amplitude of membrane potential A, damping constant b, number of pulses m) depend on hardware characteristics including voltage level and how quickly a neuron returns to its resting state. The generating action potential depends on the configurable parameters (duty cycle and subthreshold  frequency). Therefore, these parameters work the same as the weight of the perceptron neuron model. They are used to teach a neuron depending on data and applications.
Although (18) seems to be more complicated than the existing ODE-based models, it has a clear definition of threshold level. In addition, the block model in Fig. 5 shows a general structure for the DRAF neuron that can provide concepts for hardware implementation.

B. SIMULATION 1) EVOKING DAMPED OSCILLATION
We simulate (18) by using R language with fixed parameters A = 1, m = 3, b = 0.02, k = 5000, and configurable parameters d = 0.3 (30%), ω ≈ 125.7 rad/s (frequency is 20 Hz) in 300 ms, the time resolution is τ = dt = 0.1 ms. The threshold value is V TH ≈ 1.37. An impulse, which has 0.2 amplitude and duration 0.1 ms, is injected into the cell at the time t n = 110 ms. States of the neuron in a complex domain are shown in Fig. 6. Potentials of the neuron's components (damped oscillator and integrator) are shown in Fig. 7.
In Fig. 6, the small circle marks the starting point of the trajectory. The progression of membrane potential toward the resting level can be seen along the path between the small circle and the small square. An injected current, which occurs at t n = 110 ms with interval τ = dt = 0.1 ms, pushes the membrane potential to a high level but below 1.37, which is the threshold value. After that, membrane potential evokes its damped oscillation, which is shown clearly as the curve ''Im(subOsc)'' in Fig. 7. The trajectory in Fig. 6 depicts the  boundary of the neuron state's damping and starting regions. The white region covers all returning paths of the damped oscillation. The gray area indicates whether the injected current causes bifurcations on the neuron's trajectory to cross the second-highest curve, which is the second peak of the damped wave. If it is crossed, a new damped oscillation is evoked. In Fig. 7, the resting level of the membrane potential is −0.5. With the value of the damping constant b = 0.02, and the eigenperiod of the oscillation 50 ms, the damped wave has two significant pulses at 110 ms and 160 ms. The peak value of the second pulse is half of the first pulse. Along with the subthreshold oscillation, the integrator component integrates the membrane potential into its capacitor. The curve of the ''Vc'' shows its potential. At the moment of the first injection, the integrator samples the potential of the neuron's membrane and leaks to a zero level. Fig. 8 shows the changes of the membrane potential, which is the total of the ''Vc'' and ''Im(subOsc)''.
In Fig. 8, the total membrane potential of the neuron is the image part of (11), which is labeled ''Im(cell_State)''. The potential has a damped oscillation due to the damped oscillator. Before the time t n = 110 ms, the potential is moving to the resting state, this moving corresponds to the path between the circle maker and the square maker on the neuron's trajectory in Fig. 6. During the injection moment, the potential changes its phase and rises up to a high level. The highest peak of this potential is below the threshold level V TH , then there is no spike in Fig. 9, which shows the result of the logistic function applied to the total membrane potential of the neuron. The parameter k is very large, then the ''Action'' line in Fig. 9 does not show any subthreshold wave.

2) RESONATING EVENT
In this simulation, we apply the second impulse into the neuron, the interval between the first and the second impulses is 50 ms as same as the eigenperiod of the subthreshold oscillation. The change of the neuron's state trajectory is shown in Fig. 10.
The first impulse pulls the neuron trajectory to the gray region, then evokes the damped oscillator. During the progress of returning to the resting level, the second impulse comes to the cell and once again pulls the trajectory from the square maker to the triangle maker in Fig. 10. At the moment of the second impulse, the neuron's trajectory is still high and near the boundary between the white and gray areas. Consequently, the trajectory easily rises over the threshold level and generates a spike. To explain this event more clearly, we review the potentials of the neuron's components in Fig. 11. In Fig. 11, the first impulse causes the oscillation of the damped oscillator. The integrator stores the existing potential of the damped oscillator. At the moment of the second impulse, the integrator once again samples the second peak of the damped wave while the integrator's potential is still high. The integrator's potential rises to a new peak, which is higher than the peak at the first impulse. As a result, the total membrane potential is pulled to a higher level than V TH = 1.37, as shown in Fig. 12, and evokes a spike, as shown in Fig. 13.
The spike in Fig. 13 has digital-like property due to the logistic function with a large value of the growth rate k = 5000. With the results in Fig. 9 and Fig. 13, the action potential of the DRAF neuron is having the all-or-none feature. 62324 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.    Furthermore, because the impulse interval is equal to the eigenperiod of the neuron's oscillation, we say the neuron resonates with the impulses and then fires.

3) NON-RESONATING CASES
To double-check the neuron's resonance feature, we simulate the case the impulse interval is different from the eigenperiod as in Fig. 14.
In Fig. 14, the first impulse evokes the subthreshold oscillation of the membrane potential, and the integrator catches the membrane potential as a part of its initial value for discharging. The damped oscillation is almost at the resting level with a negative value right before the moment of the second impulse. The integrator's potential, on the other hand, is at a positive value. The total membrane potential cannot cross over the threshold V TH = 1.37 to generate a spike. Positions of the neuron state in a complex domain are shown in Fig.15.
As shown in Fig. 15, at the moments right before the second and the third impulses, the neuron's trajectory is very   near to the resting state. The changing in positions due to the impulses can not rise the trajectory high enough to have spikes. Therefore, the total membrane potential just has the subthreshold oscillation, as shown in Fig. 16.
In Fig. 17, the action potential does not have any spikes. The time interval between the two impulses of this simulation is half the eigenperiod of the damped oscillation. Thus, the result verifies the frequency preference feature of the proposed DRAF model. This feature is also referred to as selective communication.

4) COINCIDENCE EVENT
This is the last simulation to confirm the operation of the proposed DRAF model. We set up two impulses having very short time intervals (5 ms) compared to the eigenperiod of the neuron (50 ms), as shown in Fig. 18.
In Fig. 18, the duty cycle of the subthreshold oscillation is d = 0.3, or 30%, and the interval of the high state is   15 ms. With a 5 ms time interval between the two impulses, the integrator component of the neuron does not leak much energy and has a chance to rise to a much higher level in a short time. Fig. 19 shows the traveling path of the neuron state in this case.
At the right before the moment of the second impulse, the neuron's trajectory in Fig. 19 is out of the white region. The boundary between the white and gray areas does not have any effect to prevent the movement of the neuron state. When the second impulse comes to the cell, the trajectory easily rises above the threshold level V TH = 1.37. Fig. 20 shows a spike of the total membrane potential in this case.
In Fig. 20, the neuron shows its capability to generate action potential when receiving coincident impulses. This feature, which is named coincidence detection, is as same as the original RAF model in [100]. The digital spike of the digital RAF neuron is shown in Fig. 21.

IV. ALL DIGITAL RAF MODEL A. MODEL
The simulations in Section III-B confirm the dynamic properties of the proposed DRAF neuron model: • Damped subthreshold oscillation of neuron's membrane potential.
• Coincidence detection. The all-or-none spikes of the proposed DRAF model allow the compatibility of the model with digital systems. Although the model uses a pulse train for the subthreshold oscillation, the damped oscillator of the block model has an output signal in the form of an analog wave, as in Fig. 7. To archive an all-digital design, we add a logistic function after the damped oscillator to convert its wave to digital pulses. As a result, we propose an all-digital block model as shown in Fig. 23 In Fig. 23, the first left logistic function (L 1 ) converts the pulses from the damped oscillator to true digital pulses. The threshold value of the function depends on the number of effective pulses of the damped oscillation.

B. SIMULATION
In this study, we simulate the model of the all-digital RAF neuron by using R language with parameter values the same as those for the DRAF neuron in Section III-B1. As a result, the damping constant b = 0.02 creates two effective pulses of the damped oscillation, as in Fig. 7. Hence, we choose the half-peak of the second pulse of the damped wave as the threshold value of the logistic function L 1 . The simulation results are shown in Fig. 22.
In Fig. 22, for the first injection at t n = 110 ms, the damped oscillator evokes two digital pulses with 50 ms of eigenperiod and 15 ms of duty cycle, as shown on the line ''Im(subOsc)''. After that, the oscillation is turned off. A long silence occurs from 200 ms to 300 ms. It denotes that the analog damped wave has been converted to a digital pulse train. We do not focus on converting the waveform of the integrator at the model level because the integrator has the 62326 VOLUME 11,2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.   function of memory. It is easy to implement the integrator using digital circuits. It is also the reason the IAF and LIF models are the most popular models used in neuromorphic hardware. To continue with the simulation, a long doublet injection is applied into the cell between 310 ms and 400 ms. The time interval of the doublet is longer than the eigenperiod of the neuron, there is no spike on the line ''Action'', which is the output action potential of the all-digital RAF neuron. The second doublet is injected between 500 ms and 550 ms. In this case, the doublet has the same 50 ms interval as the subthreshold oscillation. The total membrane potential of the cell, which is the line ''Im(cell_State)'', crosses over the threshold level ''Vth''. The result of the event is the appearance of an all-or-none spike on the line ''Action'' at the time of 550 ms.
From 650 ms to 670 ms, the third doublet enters the cell. Its time interval is 20 ms, which is shorter than the eigenperiod of the neuron. However, the second impulse of the doublet drops into the low phase of the subthreshold oscillation because the duty cycle of the oscillation is 15 ms. Consequently, there is no spike on the line ''Action'' between 650 ms and 800 ms. Finally, a very short doublet (coincident impulses) is applied at 800 ms and 805 ms. This 5 ms interval, which is within a 15 ms duty cycle of the neuron's oscillation, easily pulls the neuron's membrane potential up above the threshold level. This creates a digital spike on the ''Action'' line.
The results of the simulation show that a true digital pulse train can be used to digitally implement the RAF model. Furthermore, the proposed all-digital RAF neuron model correctly performs the RAF model's important dynamic properties. For these reasons, the proposed digital and all-digital RAF models are the most possible replacements for the IAF and LIF models in hardware implementation. In the next section, we present our work on turning the all-digital RAF model into real digital circuits.

V. HARDWARE IMPLEMENTATION A. DIGITAL CIRCUIT DESIGN
In Section IV-A, we proposed a block model of an all-digital RAF neuron, which has a digital oscillation and digital output VOLUME 11, 2023 T.-K. Le et al.: Modeling and Designing of an All-Digital Resonate-and-Fire Neuron Circuit spikes. In the model, the damped oscillator and the logistic function L 1 can be grouped into a single bock as a digitally controlled oscillator (DCO) or a numerically controlled oscillator (NCO). This DCO or NCO could be enabled or disabled to perform the function of damped oscillation. The integrator is the component that is not totally digitalized in the model.
According to (13) and (14), the integrator stores the whole neuron's state at the injection moment. From a visual perspective, Fig. 3 shows that the positions of the membrane potential and the injected current in the neuron's operation are equivalent to the positions of an analog signal and a sampling clock in an analog-to-digital process. Therefore, to implement the neuron's operation, we use a sample-andhold method, as shown in Fig. 24.
The D flip-flop in Fig. 24 uses current impulses as its main clock. Its data comes from the neuron's oscillator, which is either DCO or NCO. Based on this circuit, the phase of the neuron's oscillation is taken into the DFF at the moment of injection. The captured neuron's state is later sent into a buffer, which plays the function of a storage area like an integrator. To catch resonance and coincidence events, we just need to check the two nearest impulses. Therefore, we use 2-bit shifters as phase accumulators. From these concepts, we present in this section a schematic design of the all-digital RAF neuron in Fig. 25.
In Fig. 25, the neuron cell has an n-bit NCO whose duty cycle, and period can be programmed via the signals ''Duty[n-1:0]'', and ''Weight[n-1:0]''. The NCO uses a high-speed clock for its operation. If the NCO is designed using a ring oscillator, the high-speed clock can be removed. When an impulse is applied to the cell, it activates the signal ''Osc Enable'' via the DFF Q4. The signal releases the NCO from its reset state to generate the neuron's subthreshold oscillation. To have the function of a damped oscillation as in Fig. 22, we design a 3-bit shifter to work as an oscillation damper. The DFF Q5 and Q6 count how many pulses are generated. At the third pulse of the oscillation, the DFF Q7 is flipped, and then the NCO is reset. This reset state of the NCO is held until a new impulse arrives at the neuron. Because of the DFF Q4, an impulse cannot reactivate the NCO if the NCO is already oscillating. So, the NCO goes out of the reset state if the time interval between a new impulse and the last enabling event is longer than three times of the oscillation's period.
While the NCO is sending out pulses, the two sample-andhold circuits, Q0 and Q8, catch the phases of the subthreshold oscillation at every injection. Their mechanisms of action are not the same. The DFF Q0 stores the phases of pulses captured by excitatory injection during the high phase of the subthreshold oscillation. In contrast, the DFF Q8 just works at the low phase of the neuron's oscillation. At this phase, the Q8 samples the inverted pulses, whose levels are high. If there is any injection at the moment, we infer that the time interval of the last injected doublet is shorter than the eigenperiod of the subthreshold oscillation but longer than its duty cycle. In this case, a reset signal is sent to clear the last data in the phase accumulator of the resonance detector. After the first pulse of the neuron's oscillation, the sampled data is orderly shifted into DFFs Q1 and Q2 by the falling edges of the subthreshold oscillation. Q1 receives the latest data (the latest bit), while Q2 holds the previous data (the oldest bit). If both DFFs have a high value, the time interval between the last two impulses is approximately equal to the eigenperiod of the neuron. A logic gate AND checks that case and generates an action if both Q1 and Q2 have a high value. The action is delayed for a very short time and then fed back to the Q2 to clear the oldest data in the phase accumulator of the resonance detector. We keep the latest data to allow continuous action if the time distance from the next impulse to the last impulse is the same as the oscillation period. Logic states of the DFFs Q0, Q1, and Q2 in the mechanism of resonance detection are shown in Table 2.
In Table 2, the logic state of the subthreshold clock at the reset state is high. This state is useful for the integrator component of the neuron if its eigenfrequency is zero, as explained in (12). The first impulse in step 1 evokes the NCO, and the second impulse in step 3 samples the second pulse of the NCO. The falling edge of the subthreshold clock on the transition from step 3 to step 4 pulls both Q1 and Q2 into a high state. Consequently, the line ''RS_action'' changes to a high state and then resets Q2 to a low state. The process adds an all-or-none spike to the ''action potential'' of the neuron. Now consider the coincidence event, which occurs when two or more impulses arrive at the neuron during its high phase of oscillation. During this high phase, we can assume that the damped oscillator of the all-digital RAF model is being stopped and only the integrator is working. Thus, the neuron is working like a traditional IAF neuron. In this case, we cannot use the resonance detector's buffer, which is constructed from the DFFs Q1 and Q2, due to the missing shifting clock. A DFF Q9 is added after the Q0 to form 62328 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  a new buffer. The injected pulse is its working clock. This new buffer accumulates impulses by shifting the logic value '1' into its flip-flops. This mechanism is similar to how a digital IAF neuron works. We can assume that the neuron generates an action when either all flip-flops or a certain number of the flip-flops have a high value. Then, the buffer can be extended to a higher number of flip-flops depending on specific applications. In this study, we would like to detect the coincidence event quickly, which means just two impulses at the same time. Therefore, only Q9 is needed together with Q0 to construct a 2-bit buffer. Like the buffer of the resonance detector, a logic gate AND is also used to help the neuron determine the action potential based on the values of Q0 and Q9. The action is also delayed for a short time and then used to clear the previous data, which is stored in the DFF Q9.
The output action potential of the design in Fig. 25 is the sum of the actions of the resonance detector and the coincidence detector. The inputs of the design have two types of injections: • Excitatory spikes: These spikes are summed via an OR gate to become a spike train.
• Inhibitory spikes: These spikes use an XOR gate to cancel excitatory spikes. The two inverters on this input are used to reduce the racing condition between this input and the output of the OR gate on the excitatory input.
Due to the use of the XOR gate to sum all impulses, not only the excitatory spikes but also the inhibitory spikes can evoke the pulse train of the NCO. The function emulates the post-inhibitory spikes feature of the RAF neuron, which is a notable property of the RAF model [100]. The overall design just uses basic logic gates (NOT, AND, OR, and XOR) and D flip-flops. The detailed schematic shows the ability to turn the all-digital RAF neuron into real VLSI hardware using standard CMOS technology processes. To verify the functions of the proposed design, we implement the schematic on an FPGA kit using Verilog in the following section.

B. IMPLEMENTATION ON FPGA 1) DESIGNING OF AN NCO
As mentioned in Section V-A, the NCO can be implemented by using digital self-oscillators such as ring oscillators. VOLUME 11, 2023   To prove the functions of the all-digital RAF design, we use an FPGA, which is a clock-driven device. Therefore, the design of the NCO is implemented by using an 8-bit binary up counter. So, a high-speed clock is required to create pulses from the NCO. The schematic of the NCO is shown in Fig. 26.
The NCO in Fig. 26 is reset at a high level. At the reset state, its 8-bit up counter stops working, and both DFFs Q0 and Q1 are also held in the reset state. As a result, the output signal ''OSC'' has a high level. The dynamic power of the design is defined by where α is a constant, V is voltage, and f is switching frequency. Now, all flip-flops inside the NCO are held in a steady state. The switching frequency of the NCO is zero. So, the dynamic power of the design is zero. The NCO starts working when all of the ''Reset'' signals are pulled low. At each rising edge of the ''High speed clock'', it counts and compares the counting value to the values of ''Weight'' and ''Duty'' using simple comparators made of XNOR and AND gates. The duty cycle and period of the ''OSC'' depend on the results of these comparators.

2) DESIGNING OF DELAY CELL
In a CMOS VLSI system, a delay cell can be easily constructed by using two inverters in series. It is not easy to do that on an FPGA system. Therefore, we use the design in Fig. 27 to emulate a delay cell.

3) RTL STRUCTURE
Applying the above NCO and Delay cell to the design of the all-digital RAF neuron in Fig. 25, we synthesize one all-digital RAF neuron cell on an FPGA Cyclone V 5CSEBA6U23I7 using Quartus Prime Lite. The RTL structure of the cell is shown in Fig. 28. According to Table 3, we need 20 ALM units and 23 registers to implement one all-digital RAF neuron. Based on the total resource of the Cyclone V 5CSEBA6U23I7, we estimate a spiking neuron network (SNN), which has 2000 neuron cells, can be realized on the FPGA. In the case of using ALMs only, the proposed RTL uses 54 ALM units for one neuron cell (0.13% of the total ALMs in the Cyclone V 5CSEBA6U23I7).

A. SIMULATION
We simulate the RTL in Fig. 28 in Verilog using the below values: • Frequency of the high speed clock: f iClk = 5MHz.  Fig. 29. Fig. 29a shows the simulation results in the case only excitatory spikes arrive at the cell. The two excitatory inputs are ''siExcit[0]'' and ''siExcit [1]''. The total signal of all inputs is ''sExcit''. When the time interval between two impulses is long enough, the enable signal of the NCO ''enOsc'' is held in a high state. Then, the neuron's oscillation output ''oClkSub'' stops working. The neuron rests with zero dynamic power during this period.
In the cases of long doublet and short doublet, the neuron does not generate any spike on the output action ''oAct''. Around the time of 1 ms, the time interval of the doublets is about 50 µs. Therefore, the doublets resonate with the neuron's oscillation. Action spikes on the line ''oAct'' are the results of the events. Fig. 29 also shows clearly a coincidence event near the time of 2 ms. The doublet, whose time interval is shorter than 20 ms, causes one spike.
In Fig. 29b, inhibitory spikes are applied to the neuron on the line ''siInhib''. On the line ''siExcit [1]'', some inhibitory spikes arrive at the cell as the same time as excitatory spikes before the time of 1 ms. These spikes cancel each other out, and the neuron is kept in a resting state. When there are only inhibitory spikes, they can evoke the damped oscillation 62330 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.   of the NCO, and the neuron can send out action spikes at resonance events.

B. EXPERIMENTAL MEASUREMENT
We verify the important functions of the proposed all-digital RAF neuron model by setting up an experimental system on a DE10-nano development kit, which uses the Cyclone V 5CSEBA6U23I7. The Cyclone V family is based on Intel's 28 nm CMOS technology. The core's supply voltage is 1.1 V. The I/O pins of the kit were configured to work with 3.3 V TTL logic. A high-speed, multi-channel oscilloscope is used to capture the all-important signals of the cell. The setup of the board and measuring machine is shown in Fig. 30. The captured waveform is shown in Fig. 31.
In the Fig. 31, we capture the neuron's subthreshold oscillation, the output action potential, the total input spike train, and the output logic of the DFF flip-flop, which plays the role of the sample-and-hold circuit at the high phase of the oscillation (the Q0 in Fig. 25). The results on the waveform confirm the functions of the design, including the coincidence detector, the resonance detector, and the damped oscillation of the neuron (the rest state). Thanks to the resonance effect,   the response time of the neuron is within one neuron's eigenperiod.

C. COMPARISON RESULTS
A comparison of resource utilization among recent publications is listed in Table 4. The proposed all-digital RAF neuron requires the least resources compared to the other models. The most related work [99] utilizes many types of resources, including counters, comparators, adders, and even subtractors. The number of DFFs in [87] and [92], which use the IZK model, is more than twenty times our architecture. That means the proposed design of the ADRAF neuron saves energy and chip area by 70% to 90%, according to Table 4. Moreover, their implementations have to use multipliers and DSP blocks. Realizing these components in a VLSI system will increase design cost, chip area, and power consumption.
The proposed all-digital RAF neuron just has two simple configurable parameters: eigenperiod and duty cycle. The design in Section V-B uses two 8-bit buses for period and duty. Thus, we can estimate that an SNN with 2000 ADRAF cells requires 4 KB to store the neuron configurations. Also, because the pulse train is slowed down, the neuron's dynamic power goes down to zero when it is at rest. This feature cannot occur in ODE-based models, which continuously calculate their neurons' states on every system's clock.
In a comparison with the analog implementation of the RAF neuron model [97], our all-digital RAF neuron model is free from PVT variations. Table 5 shows the comparison results. The analog RAF neuron in [97] was simulated using an AMI 0.35 µm CMOS process, while the ADRAF was tested successfully with an Intel 28 nm CMOS process. It means the ADRAF uses 12 times smaller transistors without any issues. Due to its digital components and fast response time, the neuron is sufficient for high-speed applications in speech and visual processing.

VII. CONCLUSION
The resonate-and-fire neuron model has not only the important dynamic features of a biological neuron but also the simplicity of a mathematical model like the IAF neuron model. However, it is still not a common model in hardware implementation, especially in digital integrated circuits. By using math and simulation, we constructed the first block models of both digital and all-digital RAF neuron models. We also presented a digital circuit design to turn the all-digital RAF model into real hardware. The design is simulated using Verilog and measured on an FPGA Cyclone V. All results confirmed that the proposed models and design have the most of salient properties of the original RAF model: • Damped subthreshold oscillation of the neuron's membrane potential.
Our work also showed a detailed RTL structure of the all-digital RAF neuron. By doing so, our study solved the weakness of the spiking neuron model in both analog and digital designs. Using digital components, the ADRAF neuron solves the device mismatch and low-speed issues in analog design. Due to the damped oscillation, the ADRAF neuron solves the dynamic power issue in digital design. The implementation of the ADRAF utilizes the fewest hardware resources compared to recent publications. It allows us to create a large-scale neuron network on a medium-resource FPGA and make tiny VLSI systems.