Realization of Analog Wavelet Filter Using Hybrid Genetic Algorithm for On-Line Epileptic Event Detection

As the evolution of traditional electroencephalogram (EEG) monitoring unit for epilepsy diagnosis, wearable ambulatory EEG (WAEEG) system transmits EEG data wirelessly, and can be made miniaturized, discrete and social acceptable. To prolong the battery lifetime, analog wavelet filter is used for epileptic event detection in WAEEG system to achieve on-line data reduction. For mapping continuous wavelet transform to analog filter implementation with low-power consumption and high approximation accuracy, this paper proposes a novel approximation method to construct the wavelet base in analog domain, in which the approximation process in frequency domain is considered as an optimization problem by building a mathematical model with only one term in the numerator. The hybrid genetic algorithm consisting of genetic algorithm and quasi-Newton method is employed to find the globally optimum solution, taking required stability into account. Experiment results show that the proposed method can give a stable analog wavelet base with simple structure and higher approximation accuracy compared with existing method, leading to a better spike detection accuracy. The fourth-order Marr wavelet filter is designed as an example using Gm-C filter structure based on LC ladder simulation, whose power consumption is only 33.4 pW at 2.1Hz. Simulation results show that the design method can be used to facilitate low power and small volume implementation of on-line epileptic event detector.


I. INTRODUCTION
Epilepsy is a group of neurological disorders characterized by unprovoked and recurrent seizures, which affects approximately 50 million people worldwide. Electroencephalograph (EEG) is the key tool for clinical epilepsy diagnosis by means of detecting the epileptic events, e.g. interictal spikes and spike-and-waves [1]. Traditionally, EEG monitoring is performed by a short-duration test in hospital with a low diagnostic yield. Alternatively, ambulatory EEG (AEEG) is widely used for long-term monitoring which allows to be performed in patients' home environment. Nonetheless, current AEEG systems require long cables, and hence, have The associate editor coordinating the review of this manuscript and approving it for publication was Venkata Rajesh Pamula . difficulty in reducing the size, weight and the amount of EEG data needed to be analyzed by neurologists [2]- [4]. To alleviate above difficulties, wearable AEEG (WAEEG) has been proposed to transmit EEG data wirelessly [5], [6]. However, long-term monitoring can generate huge amounts of EEG data, and make wireless transmission very powerhungry, that is unsuitable for the battery powered WAEEG which has stringent power budget. As a solution to the aforementioned problem, online data reduction for WAEEG system has been presented, in which epileptic event detection algorithm (EEDA) is used to reduce the amount of wirelessly transmitted EEG data, and thus the power dissipation [7], [8]. Fig. 1 illustrates the signal processing blocks for WAEEG system. As a precondition for saving power by data reduction strategy, EEDA should be implemented at a very low power VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ level. To achieve above goal, [9] has presented an approach to realize EEDA by continuous wavelet transform (CWT). Continuous wavelet transform can perform multiscale analysis that is well suited for detecting the nonstationary epileptic events [10], [11]. Meanwhile, CWT can be implemented by analog bandpass filter (i.e. wavelet filter) which can achieve the better performance in real-time and ultra-low power operation compared with digital counterpart [12]- [19]. A 60 pW analog wavelet filter for epileptic event detection in WAEEG system has been reported in [16]. To be synthesized by analog filter, wavelet base should be approximated by rational fraction (i.e. analog wavelet base). As a first step, the construction of analog wavelet base plays an important role in designing qualified wavelet filter to meet the requirement of WAEEG, since high approximation precision can generate low-order analog wavelet filter with high spike detection accuracy which will facilitate future low power implementation. The Maclaurin series approximation (MSA) presented in [13] is proven to be well suited to the design of wavelet filter for on-line epileptic event detection [9]. Compared with other approximation methods, MSA can generate a simple analog wavelet base with only one term in the numerator, which will greatly benefit the following high-performance filter synthesis while keeping low power dissipation. However, MSA method has low approximation precision, which will generate inaccurate wavelet coefficients and thus low spike detection accuracy. Furthermore, MSA cannot guarantee the obtained analog wavelet base stable for arbitrary filter order and time delay. To overcome the problems with MSA, the approximation method using genetic algorithm has been proposed in [20]. Although successful in generating stable wavelet filter, this approach still has trouble in enhancing approximation accuracy. This paper proposes a novel approximation method suitable for implementation of high-quality on-line epileptic event detector. Instead of direct Maclaurin series expansion, the proposed method builds a mathematical approximation model, and converts the construction of analog wavelet base to an optimization problem. Then, the hybrid genetic algorithm (GA) is employed to find the optimum approximation, that is, genetic algorithm and quasi-Newton method are combined together to obtain the near-globally-optimal solution and globally optimal solution, respectively. Also, to guarantee the stability of constructed analog wavelet base, penalty function is introduced in the fitness function of optimization model to avoid selecting unstable solutions. Finally, the optimized analog wavelet base is synthesized by Gm-C filter structure derived from LC ladder simulation. As an example, the construction procedure of analog Marr wavelet base is illustrated. In addition, the comparison of epileptic event detection obtained by proposed method and MSA is carried out. Experiment results show that the performance of proposed method is better than that of MSA. Particularly, the constructed fourth-order wavelet base achieves a better detection result compared with the seventh-order wavelet base generated by MSA, which is significant for further reduction of power consumption and chip size. Based on SMIC 1V 0.18 µm CMOS technology, the wavelet filter is designed to cover the frequency range of EEG. Simulation results show that the power dissipation of designed wavelet filter at scale a = 0.1 (i.e. 2.1 Hz) is only 33.4 pW with a 46 dB dynamic range.

II. APPROXIMATION OF WAVELET BASE IN ANALOG DOMAIN
Wavelet transform decomposes a signal f (t) into components at different scales, which can be expressed as [10] where ψ(t) is wavelet base, a and b are the scale and the time-shift parameter, respectively. Based on the definition of convolution, (1) can also be written as Apparently, the CWT at scale a can be realized by analog filter whose impulse response is 1 √ a ψ( −t a ) [13]. Fig. 2 illustrates the CWT circuit realized by analog wavelet filter bank.
For bioelectrical signal processing, Gaussian-family wavelet bases are extensively used and normally considered 33138 VOLUME 8, 2020 as the best choices. For example, the Marr wavelet transform has been applied in on-line epileptic event detection for WAEEG system [9], [16].
Generally, the Marr wavelet base at scale a in time-domain can be expressed as [13] Applying Fourier transform to (3) and denoting s = jω, the Marr wavelet base at scale a in frequency-domain can be given as Obviously, the Marr wavelet is noncausal and cannot be synthesized by analog filter directly. Therefore, a time delay t 0 should be introduced to make wavelet causal. It should be noted that the time-shift only delays the filter output, without change in amplitude. Thus, the delayed Marr wavelet base in frequency-domain at scale a can be given as [13] H (s) = −π According to filter theory, only rational fraction with strictly Hurwitz polynomial in denominator can be synthesized by analog filter. Thereby, the denominator of (5) should be approximated by a polynomial in terms of s, i.e.
H a (s) = −π 1 4 8 3 a 5 s 2 B n s n + B n−1 s n−1 + · · · + B 1 s + 1 (6) where coefficient B i is positive real number, and all of the poles locate at the left hand plane.
To construct the analog Marr wavelet base with rational form as (6), MSA method is proposed in [13] to approximate the exponential term by Maclaurin series expansion, i.e.
Different from other approximation methods producing complex numerator [12], [17], [18], the simple rational form with only one term in the numerator can yield a circuit with low complexity, the feature that is desired in ultra-low power implementation of epileptic event detection for WAEEG system [16]. More importantly, the characteristic of purely even function of s in the numerator is desired in high-quality wavelet filter design. For example, LC ladder filter structure is well suited for synthesizing analog wavelet base due to the low sensitivity to component variation, which is beneficial to future fabrication in the IC process. However, LC ladder structure can only realize the transfer function whose numerator must be a purely even or odd function of s, which means analog wavelet base with complex or arbitrary numerator cannot be synthesized directly [21].
Another merit for MSA is the automatic satisfaction with the important property of wavelets, i.e.admissibility criteria, without which the bias in computing wavelet coefficients will be introduced. Observed from (6), the analog wavelet base H a (s) equals zero at s = 0, corresponding to the admissibility property expressed in time domain as where h a (t) is the impulse response of H a (s). Although successful in many aspects, MSA method has shortcomings as below: 1) Theoretically, MSA method cannot realize high accurate approximation since excellent polynomial fitting of the denominator does not always bring out the excellent approximation of rational fraction. For instance, the seventh order analog wavelet base used for WAEEG system in [13] does not give a good approximation in high frequency region within the concerned range 0-50 dB, which will generate inaccurate CWT coefficients and may result in false spike detection and high data transmission rate. Meanwhile, the third and fifth order analog wavelet base derived from MSA both have noticeable deviation from ideal Marr wavelet, and cannot be employed in epileptic event detection [16]. And that means it is unachievable for MSA method to further reduce filter order for the purpose of minimizing power dissipation and chip size required by wearable device. 2) Practically, MSA method cannot guarantee analog wavelet base stable. The coefficients B i (i = 1, 2, . . . , n) in (6) are determined values calculated by Maclaurin series expansion shown as (7). And so, the stability of approximated transfer function depends strongly upon the inter-relationship between wavelet scale, time delay and filter order. In other words, stable analog wavelet base may not exist at several time-delay values and filter orders selected according to application requirement. For example, the third order VOLUME 8, 2020 approximation at scale a = 0.1 cannot generate stable poles at some time-shift values 13].
As an improvement of MSA, [20] presented a method converting rational approximation to an optimization problem, which can remove the possibility of instability. However, since the optimization objective is defined by polynomial fitting of the denominator in (6), the proposed method is still struggling in realizing an accurate approximation. Moreover, like MSA and other existing approximation methods, the time delay t 0 in [20] should be selected manually. Small t 0 will result in large truncation error and thus low similarity between ideal and approximated wavelet base, while large t 0 will increase the difficulty in approximation due to extra range remaining close to zero and thus require high-order wavelet filter with large chip size and high power consumption [18]. A delicate balance should be made.
To overcome above limitations, this paper proposes a novel method based on (6) to find the high accurate and stable approximation of Marr wavelet base at arbitrary wavelet scale and filter order, while generating time delay automatically.

III. PROPOSED METHOD FOR CONSTRUCTING ANALOG WAVELET BASE A. MATHEMATICAL APPROXIMATION MODEL
Generally, the approximation of wavelet base can be conducted in time or frequency domain. A majority of existing methods concentrate on approximation in time domain, however, having the limitation of manually selecting time delay t 0 first to make wavelet casual. Instead, the proposed approximation method is conducted in frequency domain.
Theoretically, the optimal approximation of magnitudefrequency and phase-frequency responses between analog and ideal wavelet base should be made. Since (4) has linear phase, the performance in approximating magnitudefrequency response and phase-frequency response reflects the waveform similarity and time shift between ideal and analog wavelet base in time domain, respectively. Hence, the proposed method only focuses on the optimal fitting of the magnitude along with frequency variation between (4) and (6) without taking phase-frequency characteristic into account. The bias in approximating phase-frequency response will lead to time shift in time domain, which can be considered as the time-delay t 0 clarified in Section 2.
Equation (6) is used as the general form of constructed analog Marr wavelet base, in which parameters B i (i = 1, 2, . . . , n) need to be determined. To evaluate the performance of approximation method, L 2 -norm is used to quantify the fitting error between ideal Marr wavelet and constructed analog wavelet base. Noted that, as an improvement, the fitting error in this paper is measured by the magnitude of overall transfer function in frequency domain, not the denominator as used in MSA and [20]. Thus, the rational approximation of wavelet base can be regarded as an optimization problem that minimizes the L 2 -norm of the error function, i.e. (9) in which the operator abs means absolute value or complex modulus, x represents the parameters B i (i = 1, 2, . . ., n) in (6) that need to be optimized. Herein, the L 2 -norm error of (9) is calculated by numerical approach with discretization, in which the discretized resolution is setting as frequency interval ω = 0.01rad/s, and the number of frequency points N is selected to cover the dominant range in frequency response.
Meanwhile, the denominator coefficient B i in (6) should be positive for filter design and the poles of H a (s) should be located at the left half of s-plane to ensure the stability of wavelet filter.
Taking all the considerations into account, the rational approximation of Marr wavelet base is converted to be a constrained optimization problem with coefficient B i in (6) to be optimized, and can be expressed as where n is the order of wavelet filter, and p i represents the pole of H a (s).
As an extra bonus, the proposed method gives a new strategy to construct analog wavelet base with time delay t 0 being selected automatically. Different from existing methods, the optimum t 0 is generated adaptively in proposed method, the reason for which can be explained as: First, the causality of constructed analog wavelet base is guaranteed by introducing constraint during optimization, that is say, the impulse response of obtained wavelet base only appears at t > 0.
Second, accurate approximation of magnitude-frequency response will lead to high similarity of impulse response between ideal and analog wavelet base.
Therefore, the causal impulse response of constructed analog wavelet base can be regarded as the time-delayed ideal wavelet base. And optimum approximation of magnitudefrequency response will generate the optimum time delay.

B. HYBRID GENETIC ALGORITHM
To resolve the optimization problem shown as (10), the hybrid genetic algorithm consisting of genetic algorithm and quasi-Newton method is employed in this paper.
Genetic algorithm is a stochastic optimization technique which mimics the process of natural evolution. As a global search method, GA can rapidly locate the region where global optimum exists, but has poor capability in locating the exact local optimum within possible region. Local search method (e.g. quasi-Newton method) can locate the local optimal solution with high speed, but strongly depends on the selection of initial solution. Thereby, the combination of GA and local search method can enhance the possibility of locating the exact global optimal solution.

1) GENETIC ALGORITHM
Based on the genetic process of biological organisms, GA searches the potential global optimum by means of bio-inspired operators, mainly including [22], [23]:

a: FITNESS FUNCTION
Fitness function is used to calculate the fitness score of each individual, which plays a critical role during optimization. The individuals with high fitness score have the chance to produce new individuals as the offspring. As for this case, the approximation result of obtained analog wavelet base is measured by the error between approximated rational fraction and ideal wavelet base. Therefore, the fitness function for constructing analog Marr wavelet base is defined by the fitting error as illustrated by (9).

b: SELECTION
Genetic algorithm selects individual genomes for later breeding according to fitness score. The individuals with high fitness score have greater probability to be selected to the next iteration, i.e. the principle of ''survival of the fittest''. There have been several selection methods, such as tournament, roulette, rank and etc.

c: CROSSOVER
A pair of parents are selected from the candidate individuals, whose chromosomes are recombined to produce new offsprings. Normally, only part of the selected individuals are performed crossover according to the preset crossover probability. By using crossover operation, the off-spring can share the characteristic of its parents, and produce high-quality chromosomes. The commonly used crossover operators are uniform, arithmetic, and etc.

d: MUTATION
To maintain genetic diversity from one generation to the next, mutation operator is usually used to introduce small changes in the chromosome, which can prevent premature and avoid local minima during the process of optimization. The popular mutation operators are bit flip, uniform, adaptive and etc.

2) QUASI-NEWTON METHOD
Quasi-Newton method is an alternative to Newton's method, which has been applied widely in optimization. Traditionally, Newton's method is frequently used in optimization problems, whose iterative formula can be given as Unfortunately, Newton's method is time-consuming in computing inverse Hessian matrix at every iteration. To reduce the computation complexity encountered by Newton's method, quasi-Newton method is proposed which constructs the approximation matrix of original Hessian matrix and its inverse matrix instead of calculating second-order partial derivative directly. Till now, several update formulas have been presented for realizing quasi-Newton method, among which Broyden-Fletcher-Goldforb-Shanno (BFGS) algorithm is the popular one [24]. In BFGS algorithm, the inverse Hessian matrix H −1 is approximated by where D is the approximate inverse Hessian matrix, s k = x k+1 − x k , and y k = g k+1 − g k .

3) CONSTRAINT CONDITION
To guarantee constructed analog wavelet base having stable poles when searching for the optimal solution of (10), a penalty term is introduced to the fitness function of the individuals who break the constraint condition. By reducing fitness score, the candidate solutions with unstable poles are penalized and have lower probability of being selected to next generation in the iteration process. The fitness function with penalty term can be generally expressed as where X defines the feasible range of optimal solution, i.e. real(p i ) < 0 in (10); r is the penalty coefficient; P(x) is the penalty function.
To optimize with random initial solution including infeasible population as required by GA, the exterior penalty function is used in this paper, which realizes r and P(x) in (13) by where gen is an increasing positive real number and equals to the number of generation in each iteration. Simultaneously, the constraint for B i in (10) can be realized by setting the lower bound on the design variables.
Then, the mathematical approximation model with penalty term can be derived from (10) and (14)

4) TERMINATION AND JUMP CONDITION
The proposed hybrid GA method can improve optimization precision and speed by integrating the advantages of genetic algorithm and quasi-Newton method. To meet above target, the optimization process should be terminated when GA converges to the vicinity of global optimum (i.e. premature convergence), and then changes to perform VOLUME 8, 2020 quasi-Newton searching method. The jump condition from GA to quasi-Newton search is defined as where c is the preset tolerance, x * represents the best individual in every generation, E i (x * ) means the highest fitness score in the i th generation, and m is a predefined number of generations to calculate the cumulative change in best fitness score.
To find the optimal approximation of Marr wavelet base modeled by (15), this paper employs GA as a pre-processor to conduct the initial search and locate the near-globally-optimal solution. When convergence is reached, the optimization process will jump out of GA and perform quasi-Newton searching based on BFGS algorithm to find the globally optimal solution.
The implementation step of proposed hybrid GA method is outlined as below.
Step 1: randomly generate the initial population; Step 2: perform selection, crossover and mutation operation to find the individual with higher fitness score; Step 3: if (16) is satisfied or the predefined maximum number of iterations is exceeded, select the best individuals in current generation as the initial solution x 0 , and jump to Step 4, otherwise jump to Step 2; Step 4: initialize the matrix D 0 (normally unit matrix I ) and the allowed maximum error ε; Step 5: if g k ≤ ε, then stop and select x k as the optimal solution, otherwise jump to Step 6; Step 6: compute the direction d k by d k = −D k g k ; Step 7: conduct line search in the direction of d k , find a stepsize α k , and compute the next solution by x k+1 = x k + α k d k ; Step 8: compute the approximated inverse Hessian matrix D k+1 by (12); Step 9: let k = k + 1, and jump to Step 5.

C. PERFORMANCE COMPARISON OF HYBRID GENETIC ALGORITHM BY TEST FUNCTIONS
The proposed hybrid GA method can greatly increase the possibility of finding optimal solution. Herein, three widely used testing functions are employed to evaluate the performance of GA, quasi-Newton method and hybrid GA. Heuristic algorithms have the characteristic of randomness. Thus, to give a fair comparison, the average performance is used for evaluation which is measured by 60 independent runs performed for each test function. The population size and maximum iterations are set to be 100 and 200, respectively.

1) SPHERE FUNCTION
Sphere function is a continuous and unimodal function, which has no local minimum but a global one. Fig. 3 plots its two-dimensional form, whose formula, search domain and   theoretical minimum point are defined in (17).
Formula : f (x, y) = x 2 + y 2 Searchdomain : −100 ≤ x, y ≤ 100 Global min imum : f (0, 0) = 0 (17) Table 1 shows the optimization results after 60 test runs for each algorithm. Obviously, quasi-Newton method can find the global minimum correctly in each test run. This can be predictable since Sphere function is convex, and local search methods are good at convex optimization. The same as quasi-Newton method, hybrid GA also gives an excellent performance. In contrast, GA cannot locate the global optimum, but mostly converges to the vicinity of minimum.

2) ROSENBROCK FUNCTION
Rosenbrock function is a unimodal and nonconvex function, whose global minimum is inside a long, narrow parabolic shaped flat valley. It is easy to find the valley, but difficult to locate the global minimum. Fig. 4 plots its two-dimensional form, whose formula, search domain and theoretical minimum point are defined in (18).  Table 2 shows the optimization results of Rosenbrock function using three algorithms. Normally, gradient-based method is inefficient for Rosenbrock function with long narrow valley. Quasi-Newton method can follow the shape of the valley   and find the global minimum using finite difference gradients. As seen from Table 2, quasi-Newton method can locate the global minimum correctly in each test run. Hybrid GA also gives an excellent performance. In contrast, GA cannot locate the global optimum, but converges to the vicinity of minimum.

3) ACKLEY FUNCTION
Ackley function is a multimodal function, which has large amount of local minima and only one global minimum. Fig. 5 plots its two-dimensional form, whose formula, search domain and theoretical minimum point are defined in (19). Global minimum : f (0, 0) = 0 As the gradient-based search method, quasi-Newton method may get trapped in local minima easily when dealing with the optimization containing many local minima such as Ackley function. Observed from Table 3, quasi-Newton method converges to local minima 59 times out of 60 independent runs when initial solutions are generated randomly. Meanwhile, GA can locate the vicinity of global minimum, but fails to converge to the accurate location, which can be seen from the best and average of optimal value in Table 3. Experiment result shows that Hybrid GA can find the global minimum correctly in each independent test run.
Apparently, hybrid GA can combine the merits of GA and quasi-Newton method, and thus enhance the performance of optimization effectively.

IV. DESIGN EXAMPLE OF ANALOG MARR WAVELET BASE USING HYBRID GA
The proposed approximation method can be used to construct analog Marr wavelet base at arbitrary filter order. To give a detailed comparison with MSA and [20], the mathematical approximation model for constructing seventh-order Marr wavelet filter at a = 1 is taken as the example. Also, the number of sampled frequency points N is set to be 701 with dominant frequency range 0-7 rad/s covered.
Then, the mathematical approximation model can be expressed by setting the parameters in (15) accordingly, i.e.
After 100 iterations, the near-globally-optimal solution obtained by GA is The L 2 -norm approximation error by GA is 1.5508. Then, quasi-Newton search using BFGS algorithm is performed with (21) as the initial solution. The globally optimal solution of (20) can be obtained as  (22) Calculation shows that the L 2 -norm approximation error by hybrid GA is decreased to 0.2695.
Substituting the parameters in (22) into (6), the seventhorder analog Marr wavelet base can be written as (23) shown at the bottom of next page. Fig. 6 shows the approximation result of proposed method compared with MSA and [20] at n = 7. It can be seen that the approximation method based on hybrid GA can improve the performance of GA greatly, and has much higher approximation accuracy in frequency domain than MSA and [20]. Table 4 shows the comparison of L 2 -norm approximation error between the methods using hybrid GA and MSA at different filter orders. Obviously, the approximation error VOLUME 8, 2020   of proposed method gets lower along with the increase of wavelet filter order, while MSA does not have this characteristic. Hence, the proposed method achieves much better performance than MSA at each filter order.
As clarified above, only focusing on the approximation in magnitude-frequency response can generate an optimum time-delay t 0 in time domain automatically. As an illustration, Fig. 7(a) gives the impulse response of constructed wavelet filter expressed by (23), as shown at the bottom of this page, compared with ideal Marr wavelet delayed by 3.3s. Fig. 7(b) shows the impulse response of constructed Marr wavelet filter by MSA and [20] compared with ideal Marr wavelet delayed by 4s. Apparently, the proposed method has much higher approximation accuracy in time domain, and yields a better fit in the domain of support, while vanishing faster with lower oscillation outside of the support interval, the feature that well suited for the characteristic of wavelet. In addition, the generated time-delay in proposed method is smaller than that used in MSA and [20] (i.e. t 0 = 4s), which gives an extra benefit for lowering the complexity of delay circuit used in WAEEG.
Simultaneously, the proposed method can guarantee constructed wavelet base stable thanks to the introduction of penalty function during optimization. Fig. 8 illustrates the poles location of seventh-order Marr wavelet base constructed by proposed method and MSA at different timedelays. Obviously, all of the poles derived from proposed method are stable (i.e. located at the left half of s-plane). However, the stability of generated poles by MSA is sensitive to the selection of time-delay.
In Table 4, it's worth noting that the approximation accuracy of proposed method at n = 4 outperforms that of MSA at n = 7. Equation (24) gives the transfer function of fourth-order wavelet base obtained by proposed method, generating 2.4s time delay in time domain. Fig. 9 shows frequency response and impulse response compared with seventh-order wavelet base by MSA. Particularly, Fig. 9(b) gives a straight comparison by aligning the impulse responses with ideal wavelet delayed by 4s. It can be seen that the fourth-order wavelet base by proposed method has better H sev (s) = −2.1741s 2 0.11s 7 + 0.45s 6 + 1.77s 5 + 3.58s 4 + 5.86s 3 + 5.68s 2 + 3.64s + 1 (23) 33144 VOLUME 8, 2020 approximation than seventh-order wavelet base by MSA both in frequency and time domain, which indicates that the chip volume and power consumption can be reduced significantly further while getting a better result in epileptic event detection.

V. PERFORMANCE ANALYSIS FOR EPILEPTIC EVENT DETECTION A. EPILEPTIC EVENT DETECTION ALGORITHM
To evaluate the performance of constructed analog wavelet base in epileptic event detection, the algorithm illustrated in Fig. 10 is employed. P a is the normalized wavelet power (NWP) at scale a denoted as (25), which is used to normalize the amplitude changes of input EEG signal recorded from different patients.
where W a represents the wavelet coefficient at scale a, and σ 2 is the variance of the signal in specific time interval. The operation procedure of detection algorithm mainly includes four stages. 1) Two Marr wavelet filters at a = 0.1 and a = 0.025 should be constructed to calculate the NWPs of input EEG signal. Herein, to facilitate low-power and low-volume implementation, the fourth-order wavelet filter defined by (24) is employed. Then, the abovementioned two wavelet filters can be designed by denormalizing (24) to desired center frequency, i.e. (26) and (27) as shown at the bottom of this page. where two time-delays will be generated automatically, i.e. t 0 = 0.24 for a = 0.1 and t 0 = 0.06 for a = 0.025, respectively. 2) Time delay blocks for τ 1 = 0.06s and τ 2 = 0.18s should be used to compensate for different time-delay values caused by two wavelet filters so as to calculate and compare the NWPs at different scales at the correct time. Compared with 0.1s and 0.3s time delay blocks used in [9], the proposed method can greatly reduce the required delay time and thus the complexity of delay circuits. 3) Candidate spikes are detected by performing the comparison where P 0.025 represents the NWP at scale a = 0.025, and β is a user-set detection threshold. 4) The artefacts and incorrect detections are rejected from candidate spikes when satisfying where P 0.1 represents the NWP at scale a = 0.1.

B. EXPERIMENT RESULTS
The epileptic EEG dataset used in detection experiment is from the Temple University Hospital (TUH) EEG Corpus, which is one of the world's largest publicly available databases of clinical EEG data [25]. Table 5 gives the summary of the information about testing EEG data, including more than 6 hours of interictal EEG data recorded from 20 patients and 60 interictal events marked by experts. The raw EEG data are sampled at a minimum of 250 Hz using a 16-bit A/D converter. The recordings have multichannel that vary between 20 and 128 channels. Wavelet analysis has been proven to be an effective tool for epileptic event detection according to three stages clarified above [9]. However, imperfect approximation to ideal wavelet base may deteriorate the detection results. Fig. 11 gives an example to illustrate this possibility.
One example of epileptic EEG data with 15s duration is shown in Fig. 11(a), and the events are marked by epileptologist. Fig. 11(b) and Fig. 11(c) depict the NWP responses at scale a = 0.025 calculated by proposed fourth-order wavelet base and the seventh-order wavelet base from MSA, respectively. The NWPs are highlighted by red arrows for real spikes from event region and black arrow for a false spike from normal EEG.
Theoretically, the spikes are supposed to have higher NWP than normal EEG at scale a = 0.025. Thus, according to rule (28), setting threshold β between the values of spikes' and normal EEG's NWP can differentiate spikes from normal background. Observed from Fig. 11(b), the NWPs of real spikes marked by red arrow are 0.068 (left) and 0.061 (right) respectively, which are higher than the false spike marked by black arrow (0.057). Therefore, based on rule (28), setting threshold β between 0.057 and 0.061 can reject the false spike and collect two real spikes successfully. However, due to the low approximation accuracy by MSA, NWPs of the real spikes and the false spike in Fig. 11(c) almost have the same value 0.07 that means false detection will be inevitable in stage 2. For example, lower threshold β will collect real spikes and also false spike as candidate spikes. Meanwhile, calculation shows that condition (29) is not satisfied at the highlighted locations with red arrow in Fig. 11(b) and the location with black arrow in Fig. 11(c). Therefore, the false spike will finally be detected by MSA incorrectly, while the proposed method can detect real spike without false detection.
The spike detection experiment is carried out by Matlab, and the Marr wavelet filters are simulated by the transfer functions. For on-line data reduction, the sensitivity and the percentage of data transmitted are employed to evaluate the performance of proposed method in epileptic event detection.
Sensitivity is defined as the percentage of the events which are marked by expert and detected correctly. Herein, to give a fair evaluation of the average performance, the total sensitivity for all tests at each user set threshold is used, i.e. total sensitivity= where M is the number of records, N i and D i are the number of events marked by expert and the number of events detected correctly in the i th record, respectively. High sensitivity is required to increase the amount of useful information so as to enhance the diagnostic yield.
The percentage of data transmitted shows the performance in data reduction, which can be expressed by percent of data transmitted = duration of data transmitted total duration of EEG data × 100% Low percentage is required to increase the amount of discarded useless information.   Normally, the total sensitivity and the percentage of data transmitted are affected by user set threshold. A small threshold will increase the number of correctly detected events and thus the sensitivity, but also the number of false detected events and thus the percentage of data transmitted. Meanwhile, a big threshold will give the opposite result. Therefore, a trade-off should be made by tuning threshold. In this paper, the threshold is set to be within the range of 0-2 in 0.05 steps. Fig. 12 shows the receiver operating characteristic (ROC) curves to evaluate the performance in epileptic event detection, in which x-axis and y-axis represent the percentage of data transmitted and total sensitivity, respectively. As illustrated in Table 6, the proposed method can correctly detect 88% of the expert-marked events while only sending 50% of the entire EEG data with 2.5s recording window, giving a 6.6% improvement compared with MSA. Table 7 gives the values of area under curve (AUC) of Fig. 12 at different recording windows. Obviously, the proposed method using fourth-order wavelet filters has the better average performance than MSA using seventh-order wavelet filters at each recording window, which is beneficial for future low-power and miniature implementation of WAEEG.

VI. DESIGN OF GM-C WAVELET FILTER BASED ON LC LADDER SIMULATION
To give a fair comparison with MSA, the Gm-C filter structure based on LC ladder simulation used in [16] is employed to synthesize the constructed fourth-order Marr wavelet base.
The transfer functions (26) and (27) have the purely even numerator and can be synthesized by LC ladder structure, which is considered as one of the most popular topologies for low power and low voltage application due to the characteristic of minimal sensitivity to inexact component value [21]. Fig. 13(a) shows the doubly terminated LC ladder prototype for the synthesis of constructed fourth-order wavelet filter.
Gm-C filter is well suited for epileptic EEG signal processing since very low characteristic frequency can be realized by Gm cell with very low transconductance. Fig. 13(b) gives the Gm-C filter structure derived from Fig. 13(a) by simulating inductors and resistor using Gm cells, whose transfer function can be expressed as (32) [21].
In this paper, the wavelet filter at a = 0.1 defined by (26) is selected as an example to elaborate the design procedure.  By using coefficient matching between (26) and (32), one can deduce the relationship as below: Numerator coefficient c in (32) can be regarded as the filter gain. The mismatch of c only introduces changes in absolute amplitude, not the relative amplitude of wavelet transform coefficient which plays an important role in threshold-dependent EEDA. Therefore, the numerator coefficient matching is not considered in (34).
The transconductance of Gm cells in Fig. 13 is selected to be identical value 100 pS to improve transconductors matching and keep the capacitance realistic for chip fabrication in low frequency application (e.g. 2.1Hz for a = 0.1). Then, the capacitances in Fig. 13 can be calculated by (33), i.e. C 1 = 6.8pF, C L1 = 10.4pF, C L2 = 22.52pF, C 2 = 6.3pF (35) Fig. 14 shows the Gm circuit realized by simple differential pair containing a cascaded load to maintain suitable bandwidth and dc operating points [26]. In a single wavelet filter, all the Gm cells share the same bias current realized by sharing transistor M 8 .
To operate in low voltage and low power application, the Gm circuit is designed in deep weak inversion, and thus the relationship between I out and V in can be deduced as [20] Then, by using Maclaurin expansion, the transconductance can be given as [16] g m = I bias 2nU t For the used CMOS process model, n and U t are 1.28 and 26mV, respectively. Therefore, to realize the transconductance of 100pS, I bias is calculated to be 6.66 pA.

VII. SIMULATION RESULTS
The fourth-order Marr wavelet filter is designed using standard SMIC 0.18µm, MIM-cap, 1 poly 6 metal CMOS process with 1V power supply. Fig. 15 shows the impulse and frequency responses of designed Marr wavelet filter at a = 0.1, respectively. Obviously, the designed wavelet filter has the satisfied responses compared with (26). The wavelet scale is realized by setting I bias = 6.67 pA to achieve the centre frequency of 2.1 Hz, corresponding to the ultra-low power consumption of 33.4 pW. The simulated filter performance is summarized in Table 8, in which the dynamic range (DR) and the signal-to-noise ratio (SNR) are defined by DR = 20 log where V im is the maximum input signal voltage, V floor is the noise floor, V irms and V noise represent the input rootmean-square and the input-referred noise voltage, respectively. Simulation results have shown that the dynamic range and SNR are 46 dB and 43dB respectively, which meet the requirements by EEG analysis [13]. The figure of merit (FOM) is also calculated by where P is the filter power consumption, V dd is the power supply, p and f c represent the number of poles and the centre frequency, respectively. As seen from Table 8, an ultra-low FOM of 0.86×10 −13 has been achieved. The scale of wavelet filter can be adjusted by changing the transconductance value of Gm cells. Fig. 16 illustrates  the impulse and frequency responses of designed wavelet filter at dyadic scales covering EEG bandwidth, i.e. 1-64 Hz. Apparently, the proposed design method can realize wavelet transform coefficients at different scales conveniently with higher accuracy.

VIII. CONCLUSION
Stable analog wavelet base with high approximation accuracy is required in low-power implementation of CWT used for on-line epileptic event detection in WAEEG system. To achieve above goal, this paper has proposed a novel method for designing stable wavelet base in analog domain. Taking the Marr wavelet as an example, the mathematical approximation model in frequency domain is constructed as an optimization problem with penalty function for stability introduced, and the hybrid genetic algorithm is employed to find the optimal approximation. Simulation results show that the proposed approximation method can yield better performance both in time and frequency domain, while generating optimum time delay automatically which is an advantage over existing approaches. The epileptic event detection experiment is conducted by using constructed fourth-order Marr wavelet base, and the result shows that an 88.0% sensitivity with 50% data reduction is achieved which is better than the seventh-order wavelet base used in existing WAEEG technique.
Besides, the fourth-order Marr wavelet filter at a = 0.1 is designed by LC ladder structure using 0.18µm CMOS process. A 2.1 Hz center frequency can be realized with 46dB dynamic range and 33.4 pW power consumption. By adjusting the transconductance of Gm cells, the wavelet filter covers the range of 1-64 Hz EEG bandwidth. Theoretical analysis and experiment results prove that the presented method can be used to reduce the complexity of existing on-line epileptic event detector used in WAEEG while keeping an acceptable detection result and data reduction rate. Leading from this paper, the future work will focus on the low power hardware implementation of other blocks in EEDA, such as delay circuit, rectifier and comparator.