A Time-Domain Analog Spatial Compressed Sensing Encoder for Multi-Channel Neural Recording

A time-domain analog spatial compressed sensing encoder for neural recording applications is proposed. Owing to the advantage of MEMS technologies, the number of channels on a silicon neural probe array has doubled in 7.4 years, and therefore, a greater number of recording channels and higher density of front-end circuitry is required. Since neural signals such as action potential (AP) have wider signal bandwidth than that of an image sensor, a data compression technique is essentially required for arrayed neural recording systems. In this paper, compressed sensing (CS) is employed for data reduction, and a novel time-domain analog CS encoder is proposed. A simpler and lower power circuit than conventional analog or digital CS encoders can be realized by using the proposed CS encoder. A prototype of the proposed encoder was fabricated in a 180 nm 1P6M CMOS process, and it achieved an active area of 0.0342 mm2/ch. and an energy efficiency of 25.0 pJ/ch.·conv.

Since AP has a bandwidth of 100 Hz to 10 kHz [25], a high-speed data transmission is needed. For example, a 10 bit, 20 ksps/ch., and 100 channel neural recording system requires 10 bits × 20 ksps × 100 ch. = 20 Mbps bandwidth data transmission. If the system requires 1000 simultaneously recording channels, the data bandwidth becomes 200 Mbps, which is unrealistic for implantable applications. Hence, a data compression technique is inescapably required for multiple-channel neural recording systems [25].
Compressed sensing (CS) [26,27] is a data reduction technique that can be realized by using a simple operation. In a measurement system based on Nyquist-Shannon sampling theorem, superfluous sampling is required in spite of the sparse information in the signal. CS is a mathematical framework that can ensure accurate reconstruction from fewer measured data, which is observed using simple matrix-vector multiplication. Since CS encoder does not require any additional circuits such as feature extractor to compress the data by other compression methods (e.g., spike detector), which causes an increase in chip area, CS-based measurement systems can be expected to reduce the chip area and power consumption [28].

Theoretical Background of CS
In this section, a theoretical background of signal compression and reconstruction based on CS is introduced. CS is a mathematical framework that ensures accurate data reconstruction from fewer measured data than that is required for the conventional Nyquist-Shannon-based signal acquisition, and has been established by Donoho [26], Candes [27], and Tao [27]. Figure 1a shows the neural signal measurement (encode) process based on CS. In CS theory, an input signal vector v ∈ R N×1 can be represented as v = Bs, where B ∈ R N×N is a basis for representing v, and s ∈ R N×1 is coefficient vector [36]. When the number of non-zero elements of s is K N, vector v is K-sparse on the basis B. If v has sparsity on the arbitrary basis B, v can be represented by using vector c ∈ R M×1 (M ≤ N), which has fewer dimensions, as where A ∈ R M×N is a sensing matrix for an incoherent sampling [37] . It is known that Bernoulli matrix in which all the entries are either +1 or −1 can be used as the sensing matrix [28]. A compression ratio (CR) can be defined as CR = N/M, and c becomes uncompressed data when N = M.  The reconstruction (decode) process is shown in Figure 1b. As effective signal reconstruction methods, l p -norm minimization [38] and block sparse Bayesian learning (bSBL) [39] are widely known. l p -norm minimization derives a sparse vector by minimizing l p -norm a p = ∑ N n=1 |a n | p (0 ≤ p < 1). bSBL can improve reconstruction performance by applying Bayesian learning. In this study, the classical l 1 norm minimization [26,27] is used to reduce the amount of calculation. The input vector v can be reconstructed by solving the convex optimization problem about l 1 -norm as argmin s∈R N×1 ŝ 1 subject to c = ABŝ, where l 1 norm of vector a ∈ R N×1 is defined as a 1 = ∑ N n=1 |a n | [28]. Since the basis B is not required for the signal encode process, the CS encoder does not require any feature extraction for signal compression. Therefore, the CS-based measurement system can reduce the system complexity. Figure 2a shows the concept of the proposed time-domain analog spatial CS encoder. The CS encoder executes the following product-sum operation per a clock cycle between an input vector v and a row of matrix A as

Time Domain Analog Signal Processing for CS Encoder
where c i is an element of result vector c, a ij is an element of the matrix A, v j is an element of v, and 1 ≤ i ≤ M is the number of elements in the vector c. In spatial CS architecture, frame rate should be larger than twice of signal bandwidth, because CS is applied only for spatial domain and each frame is encoded. Therefore, if a input signal bandwidth including AP is 10 kHz [25], a frame rate corresponding to sampling frequency, f s , should be greater than 20 kHz to guarantee the Nyquist-Shannon sampling theorem. The product-sum operation is executed by using cascaded voltage-to-delay-time converters (VTCs), which convert the control voltage into a time-domain signal. A product, a ij × v j , is set to each control terminal of the VTCs, and a conversion cycle starts at the rising edge of CLK. The timing diagram of the encoder is shown in Figure 2b. When CLK rises, the first VTC starts a voltage-to-delay-time conversion. After a delay of t d1 , which corresponds to the product a ij × v j , the first VTC raises its output D 1 . Then D 1 rises, and the second VTC starts the conversion in the same manner. Thus, the delay time generated by the VTCs are accumulated, and the total delay time, t d,ci , corresponding to c i appears between CLK and D OUT . Finally, t d,ci is converted into a digital code c i by the time-to-digital converter (TDC).  Conventional analog CS encoders, including the voltage domain [33] and digital implementation shown in Figure 3, tend to be power hungry because of the following reason: In the analog implementation, the number of samples is reduced in front of the ADC by using a product-sum circuit composed of a voltage-domain mixer and a switched-capacitor adder. Similar to the time-domain implementation, the analog CS encoder can execute the product-sum operation per a clock cycle. However, the switched-capacitor adder requires an operational amplifier, which satisfies the fast settling condition, resulting in increased power consumption. The digital implementation shown in Figure 3b requires full-sampled data before compression. It is similar to the other digital signal compressors which require large-scale memory to store the data before compression. When the product-sum circuit is realized by using a single accumulator to reduce hardware cost, a system clock frequency of M × N × f s is required for the 1-frame encode. For example, if f s = 20 kHz, M = N = 20 (uncompressed), and the required system clock frequency is 8 MHz. Since the dynamic power consumption of the clock synchronization circuit is proportional to the system clock frequency, a higher system clock frequency is undesirable. On the other hand, the proposed CS encoder is an analog circuit, which can be mainly composed of logic elements, and therefore, the proposed encoder can essentially reduce its power consumption. Moreover, since the operation of the proposed encoder is based only on delay propagation, the total number of transition cycles in the proposed encoder can be lowered than that in the conventional digital implementation. Therefore, the power consumption of the proposed CS encoder can be significantly reduced.  Figure 4 shows the block diagram of the proposed CS encoder. It comprises 5 measurement units with 20 electrodes, and consequently, the system can simultaneously measure 100 channels of a neural signal. In this design, each input signal of the channel is represented as a pseudo-differential signal, which is converted by VTC j+ and VTC j− (1 ≤ j ≤ 20), and the control voltages of VTC j+ and VTC j− are set as input voltage v j and reference voltage V CM for v j , respectively. The control voltages are kept constant by the sample and hold (S/H) circuits during 1-frame conversion. The S/H circuits sample their inputs (V CM or v j ) when the sampling clock φ is high. Multiplication with ±1 is executed by using choppers as a time-domain multiplier. The resulting encoder output is represented by the time difference between the rising edges of D OUT+ and D OUT− , and is converted into a digital code by the TDC.

Overview of the Entire Operation
The timing diagram of the CS encoder is shown in Figure 5. This encoder executes a product-sum operation between a row of the sensing matrix and the input vector per a single clock cycle. VTC j+ and VTC j− of each channel are assigned to a positive or negative delay line, respectively, by the chopper. Delay accumulation of the positive and negative delay lines are asynchronously executed, and the product-sum output appears in the relationship between the delay time t d and the control voltage V CTL of the VTC, which can be represented by a linear function as where α and β are constants. When the delay times of the VTC j+ and VTC j− are defined as respectively, the time difference between the rising edges of D OUT+ and D OUT− is directly proportional to the product-sum output as  The upper limit of each VTC's delay is determined by the stage number of the VTCs N and the frame rate f s . Since the period of the frame is T s = 1/ f s , M-times product-sum operation must be executed in T s . Therefore, the upper limit of the time for the product-sum operation is Since T conv.,max is the total delay time of the VTCs, the upper limit of each VTC's delay can be expressed by using Equation (1) as Delay propagation in the CS encoder The required number of bits for TDC N bit must be determined so as not to degrade the reconstructed data. To determine N bit , a system level simulation by using MATLAB was performed as shown in Figure 6, which plots the reconstructed SNR vs. the N bit on each compression ratio CR (= N/M). The reconstructed SNR is defined as where v i is the original data, andv i is the reconstructed data [28]. Note that noise-less TDC and VTC are assumed in the simulation. A reconstructed SNR saturation on the higher N bit is caused by CR, and thus, dropping the reconstructed SNR on a lower N bit denotes the insufficiency of N bit . According to the result, N bit was decided as 10 bit in this design. The time resolution of the TDC T LSB can also be derived from Equation (1) and N bit as In this design, T LSB was decided as 1.63 ns. Note that its input-referred noise is 4.86 µ V RMS , and is sufficiently lower than amplitude of AP (50 to 500 µ V PP [25]). In time-domain analog circuits, jitter corresponds to the noise in voltage-domain analog circuits, and it is defined as a timing deviation from the true operation timing. The total value of jitter is mainly determined by thermal noise represented by voltage or current source [40,41]. Since thermal noise follows Gaussian distribution, jitter value also follows Gaussian distribution. Therefore, jitter value can be discussed in statistics.  total comprises of the output-referred jitter of VTC array σ 2 VTC,total and the TDC input-referred jitter σ 2 TDC , where σ 2 VTC,total is the total accumulated jitter of each VTC in the VTC array, and σ 2 TDC is determined by the variations in the operation timings of the flip-flops and a ring oscillator in the TDC. Note that σ 2 TDC corresponds to a sampling noise of ADC. Assuming that all the jitter deviations follow the Gaussian distribution, the overall jitter can be derived by the sum of squares as In this design, the condition for satisfying the required accuracy of the CS encoder is defined as at least. As mentioned below, σ 2 TDC is sufficiently smaller than σ 2 VTC,total in this design. Thus, σ 2 TDC is considered to be constant. The following subsections describe the architecture and design consideration of TDC and VTC.

TDC
TDC can convert the time difference between two input clock edges into a digital code with high time resolution by using a delay time of logic elements. Figure 8 shows the block diagram of a delay-line-based and a ring-oscillator-based TDC. The delay-line-based TDC, which is shown in Figure 8a, is composed of a delay line, a D-FF array, and a decoder. The TDC converts the time difference between the rising edges of the START and STOP signals. When the START signal rises, the rising edge propagates through the delay line. When the STOP signal rises, the D-FF array captures the delay propagation state, and the captured state is converted into a digital code by the decoder. Although this TDC has the simplest architecture, 2 N bit stage delay lines are required to realize N bit bit resolution. If the desired resolution is 10 bit, a 1024-stage delay line is required. On the other hand, a ring-oscillator-based TDC, which is shown in Figure 8b, is composed of a ring oscillator, a binary counter, a D-FF array, and a decoder. In the ring-oscillator type TDC, the D-FF array captures the phase of the ring oscillator for fine conversion, and the counter measures the ring oscillator output for coarse conversion. Therefore, the ring-oscillator-based TDC requires fewer delay line stages than the delay-line-based TDC. In this design, the ring-oscillator-based TDC is employed to reduce the number of delay line stages. The following paragraphs describe the TDC jitter and the number of ring oscillator stages N. In ring-oscillator-based TDC, the delay time of the ring oscillator's delay cell t d,Ring becomes the time resolution of the TDC T LSB . The frequency of the ring oscillator f o can be expressed as where N Ring is the number of ring oscillator stages. The TDC jitter σ 2 TDC can be expressed as where σ 2 Ring is the ring oscillator jitter and σ 2 D-FF is the timing variation when the D-FF array captures the ring oscillator phase. Since σ 2 D-FF is a sufficiently smaller constant value than σ 2 Ring , σ 2 TDC ≈ σ 2 Ring . Also, σ 2 Ring can be expressed as [41] where κ is a proportionality constant that is determined by the circuit parameters, and ∆t is the measurement time. According to [41], κ is determined by the size of the transistor for the ring oscillator (W and L), the number of ring oscillator stages N Ring , and the current noise power spectral density in which are input to ring oscillator nodesī 2 n /∆ f . Since the time resolution of TDC T LSB is already determined by Equation 3, the transistor size for the ring oscillator cannot be modified to satisfy T LSB . In addition,ī 2 n /∆ f is also not changeable in this design; only N Ring can be modified. However, the TDC jitter value hardly changes regardless of changing N Ring [41]. Therefore, σ 2 TDC ≈ σ 2 Ring is regarded as a constant.
In this design, since TDC consumes most of the power in the entire CS encoder, its design should be optimized to achieve lower power consumption. The power consumption of TDC P TDC can be expressed as P TDC = P Ring + P D-FF + P Cnt , where P Ring , P D-FF and P Cnt are the power consumptions of the ring oscillator, the phase capturing D-FF, and the counter for coarse conversion, respectively. Figure 9 plots the simulated power consumption vs. N Ring during measurement. P D-FF is sufficiently smaller than the others. P Ring is almost constant and occupies the majority of P TDC . When N Ring is small, P Cnt becomes larger than P Ring . According to the above results, increasing N Ring reduces the power consumption of the TDC, and it approaches the value of P Ring . Finally, the number of ring oscillator stages is determined as N = 64 in this design. To estimate the value of σ 2 Ring , 100-times transient-noise simulation is performed. According to the result, the TDC jitter with N Ring = 64 becomes sufficiently small as 3 σ Ring = 221.6 ps.

VTC
To realize the proposed CS encoder, a transfer function with a high linearity is required between the delay time and the control voltage. Therefore, an integrator-based architecture is composed of a capacitor and current source, as shown in Figure 10, is employed in this design. A conversion trigger signal A controls the state of VTC. When the node voltage of A is logical low, VTC is set as a reset state and the capacitor's terminal voltage becomes V CTL − V CM as shown in Figure 11a, where V CM = V DD /2 is a reference voltage. The voltage-to-delay-time conversion starts at the rising edge of A.
Finally, when v c− crosses the threshold voltage of the comparator, V CM , the conversion is completed and the output of the comparator becomes high (Figure 11c). To ensure high linearity over a wide range of V CTL , the delay time is controlled only by the initial voltage of the capacitor, and the integrating current and threshold voltage of the comparator are constant. Note that the comparator is composed of a simple logic inverter to reduce its operating power.   The jitter requirement for VTC can be derived by substituting T LSB and σ TDC into Equation 4 as σ VTC ≤ 1.61 ns. Figure 12 shows the equivalent circuit model for the integration state, where I t is an integrating current, i n is the noise current of the current source, C t is the integrating capacitance, and r o is the output resistance of the current source. The relationship between the delay time t d and the control voltage V CTL can be expressed as In this design, the measurement unit shown in Figure 4 has 20 measurement channels, and each measurement channel includes two VTCs (VTC j+ and VTC j− ). Thus, the total accumulated jitter for the product-sum operation σ 2 VTC,total in Equation 4 can be represented by using a square-sum of all the VTCs' jitter σ 2 td as Therefore, the ratio between C t and I t is derived by introducing the control voltage range and the maximum delay time into Equation 6, and C t can be determined as satisfying σ VTC ≤ 1.61 ns.  Figure 13 plots the 100-times transient-noise simulation result of 3 σ VTC,total as a function of C t . In this design, the capacitance of the integrator was determined as C t = 521 fF and I t = 1.0 µA to satisfy the jitter requirement.

Design Constraints in Proposed Architecture
In this subsection, design constraints and trade-off of the proposed architecture are discussed. The chip area and the total number of channels are design constraints of the proposed CS encoder. These constraints provide maximum capacitance in the VTCs, and hence it decides realizable minimum jitter of the VTCs. On the other hand, if desired SNR and sampling rate are given, specifications of time resolution of the TDC and the VTC jitter are also determined. Since the specification of the VTC jitter cannot be less than the realizable one, the rest of design constraint is how balancing between the number of channels per measurement unit N unit and power consumption. A smaller N unit relaxes jitter requirement for VTC and TDC, and thus low-power implementation could be realized, while higher CR cannot be achieved because realizable maximum CR is same with N unit . In contrast, a larger N unit requires low jitter for VTC and TDC, increasing power consumption. Especially, since a VTC jitter is limited by its capacitor size, the low jitter requirement for the VTC cannot be realized for much larger N unit . Therefore, considering the above discussion, in this design, a moderate N unit is set as 20, and resulting maximum power consumption is 6 µW/ch. which can be expected as lower value than previous studies [33][34][35].

Measurement Results and Discussion
The proposed 100-ch. time domain analog CS encoder was fabricated in a 180 nm 1P6M CMOS process, as shown in Figure 14. The active area of the prototype encoder is 1.85 mm × 1.82 mm. The frontend for each measurement channel comprises the electrode, the low-noise amplifier (LNA), the sample and hold circuit (S/H), and the two VTCs. The active area with TDC is 0.0331 mm 2 /ch., and without TDC is 0.0272 mm 2 /ch.  Figure 15 shows the evaluation environment for the proposed CS encoder. The measurement system, which is shown in Figure 15a, comprises of the prototype CS encoder, a prototype evaluation board, an FPGA board, a power supply, and a PC. A logic analyzer is used for the development and debugging of the environment. The prototype is controlled by the control logic and the micro controller (MCU), which are embedded in the FPGA. The test signal of the encoder is provided for each measurement unit as a time-interleaved voltage signal from DAC. The test signal which simulates spontaneous neuronal activity is generated on MATLAB. The details are described in Appendix A.

Prototype evaluation board
Prototype CS encoder PC (for system control) The evaluation procedure for the CS encoder is shown in Figure 16. In this measurement, spatial test signal (input vector) were prepared for each frame and directly stored in the S/H circuits of the prototype at the beginning of the conversion frame. The test signals are generated for simulating the AP waveforms from a neural probe array, and each signal represents a 2-dimensional input voltage distribution. The acquired data from the prototype CS encoder is transferred to the PC and then reconstructed by a MATLAB-based program. Note that the reconstruction by solving the convex optimization was realized using CVX [42]. As the basis for the reconstruction, a discrete cosine transformation (DCT) matrix was selected. Data compression method based on spatial DCT has been proposed in Ref. [43] with 1/69 times data reduction at 6% root mean square error. Using DCT for data compression can imply that multi-channel APs is inherently sparse on 2D frequency domain, and thus DCT matrix can be potentially used as basis of spatial CS reconstruction. Finally, the reconstructed signal quality is evaluated by calculating the reconstructed SNR defined in Equation (2). Figure 17 shows the 100-ch. reconstructed temporal waveforms from the compressed data encoded by the prototype CS encoder at CR = 4, and Figure 18 plots the temporal change in the reconstructed SNR. Note that amplitude of waveforms shown in Figure 17 indicates input voltage for VTCs. The reconstructed SNR at t = 11.1 ms (when the input signal becomes peak amplitude) was 15.3 dB. Since the reconstructed SNR depends on the input signal, non-sparse signal degrades the SNR. If a improved SNR is required, CR should be relaxed which means an increase in the number of sampling for each frame, M. In this design, the frame rate of the encoder is defined to support uncompressed condition (M = N). Therefore, the proposed encoder can control CR from 1.0 to 20 without any hardware changes. The reconstructed SNR vs. compression ratio (CR) at t = 11.1 ms is plotted in Figure 19. Note that the measurement results of Figures 17-19 include all noises induced by the prototype chip and the measurement environment. In practical applications, influence for a spike sorting [44] must be discussed. Not to affect spike sorting performance, enough reconstructed SNR before and after spike is required for spike detection and spike classification [44]. As shown in Figures 17 and 18, spike amplitude and waveform are successfully reconstructed at CR = 4, and reconstructed SNR before and after spikes (indicated with allows in Figure 18) are around 10 dB. Note that the measurement results of Figures 17 and 18 include noise of the measurement environment. SNR for a spike detection and a classification requires more than 10 dB [45], recovered data do not affect spike sorting performance. The SNR was saturated at 20.0 dB with a CR lower than 3. The saturated SNR is lower than the expected MATLAB simulation result as shown in Figure 6. In addition, the SNR at CR = 4 is 15.3 dB, and it dropped over 3 dB compared to the SNR simulated using MATLAB (19.5 dB).
From the simulation result shown in Figure 6, it is equivalent to a degradation of 0.7 effective number of bits (ENOB), and it is considered a systematic variation in the gain of the VTC's, which degrades the dynamic range of the CS encoder. Since the gain variation of VTCs is not compensated in this design, it could be affected by the process, voltage, and temperature (PVT) variation. Therefore, a gain compensating technique for VTC is desired for improving the SNR. To achieve further improvement of CR and reconstructed SNR in practical in vivo measurement, an optimized basis could be used, which is obtained from uncompressed (CR = 1) recorded data by dictionary learning algorithm such as K-SVD [46]. As other solution for improving reconstructed SNR, an optimization techniques for sensing matrix have been proposed in Ref. [47]. Indeed, optimization techniques for sensing matrix can improve reconstruction performance. However, extra registers, which almost consumes 20% active area of the prototype, to contain the optimized sensing matrix is required. Thus, the technique have not been applied in this prototype, and the sensing matrix is generated on the chip by using simple linear-feedback shift register (LFSR).      A performance comparison of the proposed CS encoder with those developed in previous works for neural recording applications is summarized in Table 1. The prototype CS encoder achieved the lowest power consumption and the smallest area compared to the encoders from previous works. Especially, the power efficiency of the CS encoder improved by about 10-times compared to the digital CS encoder.

Conclusions
In this paper, a low-power energy-efficient neural signal acquisition system, which uses the novel time-domain analog spatial CS encoder, is proposed. In this technique, the product-sum operation for the CS encoder can be executed by accumulating the delay time information. Since a major part of the proposed CS encoder can be realized by using logic elements, it can reduce power consumption and chip area compared to conventional analog or digital CS encoders. Some design parameters for the proposed encoders were considered and optimized by a trade-off between noise and power consumption.
The 100-ch. neural signal acquisition system employing the proposed time-domain CS encoder was fabricated in a 180 nm 1P6M CMOS process, and its active area is 0.0331 mm 2 /ch. A 100-ch. CS encode experiment was performed using the prototype CS encoder, and it achieved a reconstructed SNR of 15.3 dB and conversion energy efficiency of 25.0 pJ/ch.·conv. at f s = 20 ksps and CR = 4. The prototype CS encoder achieved the lowest power consumption and the smallest area compared to the encoders in other previous works for neural recording applications. Therefore, the proposed time-domain spatial CS encoder is suitable for exponentially increasing multi-channel neural recording applications.
Future works are to generate a basis which is optimized for measured spatial APs and thoroughly to evaluate the performance of spike sorting with reconstructed spatial information.
Acknowledgments: This work has been supported in part by MEXT Grants-in-Aid for Scientific Research (15H05525) and by VDEC, the University of Tokyo in collaboration with Cadence Design Systems, Inc., Synopsys, Inc., and Mentor Graphics Corp.
Author Contributions: T.O. proposed the time-domain CS architecture, performed implementation and evaluation, and wrote the paper; I.A. designed and led the project, and wrote the paper.

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

Appendix A. Test Signal Generation
A real neural dataset which meets the conditions of our prototype chip including electrode density, electrode arrangement and the number of electrodes could not be found. Hence, in this study, the test signal which simulates spontaneous neuronal activity is generated by using MATLAB.
Upon pseudo AP signal generation, at first, membrane voltage signal of the neuron V m is generated by using algorithm based on Ref. [48]. Then, V m is converted to membrane current I m . Finally, extracellular potential is obtained from I m . An equivalent circuit of nerve cell's membrane is shown in Figure A1a, where C m is a membrane capacity, R Na and R K are sodium and potassium resistance, respectively, E Na and E K are potassium and sodium potentials, respectively, R L and E L are leakage resistance and potential, respectively [49]. Practical parameters of the equivalent circuit shown in Figure A1a are referred to Ref. [49]. N neurons are randomly placed around electrodes, which are arranged at equal intervals, as shown in Figure A1b. Note that neuron density is determined by referring to Ref. [21] as 150 neurons/mm 2 , averaged firing rate of neurons is 0.3 Hz, and the average number of firing neuron in the sensing area per 25 ms is 6 times. A membrane current of i-th (1 ≤ i ≤ N) neuron can be expressed as where I m,i is an ionic current, C m is a membrane capacity. When I m,i is assumed to be a point sink current source, extracellular potential of j'th(1 ≤ j ≤ 100) electrode can be obtained as where σ is constant conductivity, and r ij is distance between current source and electrode [50].