Hardware architecture for real-time EEG-based functional brain connectivity parameter extraction

Functional brain connectivity has been shown as one of the main methods to understand the functions of the brain. Typically, FC creates a graph amongst different areas of the brain where each area of the brain represents a node of the graph whilst the arcs between them represent the connectivity strengths between different brain regions (nodes). The fundamental importance of such an approach is that it can characterize the information communications between different parts of the brain quantitatively in terms of graph network. Two global characteristics of the brain that are important for any neurodegenerative/neurodevelopmental disease (such as Alzheimer's, epilepsy, autism and Parkinson's [1–8]) are: (1) the ability of the brain to integrate information from a distant source and; (2) specialized local processing. The graph-theoretic parameters that could quantitatively describe such properties of brain are shown in the table 1.

The PLI is a functional connectivity measurement commonly used to analyze EEG signals [4, 6, 15] introduced by [16] as a new measure to address the problem of volume conduction in functional connectivity. The PLI is calculated using (1), where φ is the instantaneous phase of the signal in the channels i and j, across a window with a total of L samples.

$$\begin{equation} PLI(i,j) = \left| \frac{1}{L} \sum_{l = 1}^{L} \textrm{sign}(\phi_i(l) - \phi_j (l)) \right|. \end{equation} \tag{ 1 }$$

The matrix obtained using (1) is an N-by-N matrix where the N corresponds to the number of nodes in the network (total number of electrodes in the EEG). Every 'i, j' element in the matrix is known as an edge, and it represents the synchronization between these two i and j nodes.

The node degree is a measurement used to represent the interconnection of one node with the other nodes in the network. This provides information about how closely interconnected the nodes in the network are. The node degree can be estimated using (2), counting the number of edges (a_ij ) greater than a certain threshold (i.e. zero) for a particular node i. Furthermore, a weighted degree ($k_{i}^w$) could also be calculated but in this case, the weights (w) of the edges are accumulated for a particular node. The weighted degree not only establishes that there is a connection between two nodes but also provides information about how strong this connection is.

$$\begin{equation} k_i = \sum_{j = 1}^{N} a_{ij}\ ;\ k_{i}^w = \sum_{j = 1}^{N}w_{ij}. \end{equation} \tag{ 2 }$$

The density can be obtained using (3), which counts the number of edges (w) greater than a threshold in the upper triangular elements of the matrix and divided by the total possible number of connections in the network. This ratio provides information on how dense the connection between the different brain regions is.

$$\begin{equation} D = 2 \sum_{j,i\in N \ j\gt\,i} upper(w_{ij} \gt threshold) / (N^2 - N). \end{equation} \tag{ 3 }$$

The clustering coefficient is used to measure the cliquishness between a group of nodes in the network, which in turn determines how clustered the network is [17]. The calculation of the clustering coefficient is done as in (4), calculating the geometric mean (t_i ) of the weights (w) for all the possible triangles around a specific node i.

$$\begin{equation} t_i = \frac{1}{2} \sum_{j,h\in N} \sqrt[3]{w_{ij}w_{jh}w_{ih}} \ ;\ CC_i = \sum_{i\in N} \frac{2t_i}{k_i(k_i-1)}. \end{equation} \tag{ 4 }$$

Also, it is possible to obtain the mean clustering coefficient as in (5) which gives a more general idea of how clustered the whole network is.

$$\begin{equation} CC_{mean} = \frac{1}{N}\sum_{i\in N}CC_i. \end{equation} \tag{ 5 }$$

Transitivity is a variation of the mean clustering coefficient measurement used to understand the segregation of the whole network. The transitivity is calculated using (6), similar to the clustering coefficient but in this case, the values of the numerator and denominator of every node are accumulated before performing the division. This normalization helps the transitivity to be less susceptible to nodes with a lower degree [13] and therefore its advantage over the mean clustering coefficient.

$$\begin{equation} T = \frac{\sum_{i\in N}2t_i}{\sum_{i\in N}k_i (k_i -1)}. \end{equation} \tag{ 6 }$$

The characteristic path length (CPL) is a measurement of integration that describes the interconnection among the different nodes and how efficient the communication between them is. The CPL is defined in (8). The shortest distance (d_ij ) between two nodes can be obtained using (7), where f(w) is a function to convert the edge values into distance values. A commonly used function is the inverse value of the weight.

$$\begin{equation} d_{ij} = ShortestPath_{ij} (f(w)) \end{equation} \tag{ 7 }$$

$$\begin{equation} CPL = \sum_{i \in N} \frac{\sum_{j\in N \ j\ne i}d_{ij}}{N(N-1)}. \end{equation} \tag{ 8 }$$

Once the edge values are converted into distances, the shortest distance between the two nodes is found. The most common algorithm used to find the shortest distance between two nodes is Dijkstra's algorithm (DA) presented in [18]. The steps to perform the DA are the following:

• (a)  
  Create a set of unvisited nodes, which includes all nodes except node i.
• (b)  
  Set the tentative distance d between the node i and every other node j to the edge length between them, i.e. d_i,j = l_i,j . Throughout the algorithm, d_i,j will hold the shortest path length from i to j through visited nodes. Currently, only node i is visited, so the initial value of d_i,j is l_i,j .
• (c)  
  In the set of unvisited nodes, find the node k with the shortest tentative distance.
• (d)  
  For every node j, update the tentative distance as follows: $d_{i,j}(new) = minimum( d_{i,j}$, $d_{i,k} + l_{k,j})$. The existing shortest path from i to j through visited nodes has a path length of d_i,j . If a path through node k (which has length $d_{i,k} + l_{k,j}$) is shorter, this would become the new shortest path between the nodes.
• (e)  
  Remove node k from the set of unvisited nodes. Repeat steps 3–5 until the set of unvisited nodes is empty.
• (f)  
  After visiting all nodes, d_i,j now holds the shortest path length from i to j through visited nodes.

The eccentricity is the maximum distance between two nodes in the network. In addition to this measurement, there is also the radius and diameter of a network, which corresponds to the minimum and maximum distance respectively for the whole graph. The distance between the nodes is calculated using (7) as in the CPL.

Since this is the first attempt, to our knowledge, for implementing such a system in hardware, there is no existing design specification for such a system. However, given the visionary application described in section 1, we considered the following conservative logical design constraints. The overall sensing, parameter extraction and the stimulation process (see figure 1) need to be completed within 1 ms considering the neuronal firing rate is O(ms). This leaves only a fraction of a ms to complete the FC computation and corresponding graph-theoretic parameters extraction part, which is the main focus of this paper. Therefore, we conservatively estimated our real-time definition in this scenario as $\lt$0.5 ms. In terms of duration of operation over a day, we do not envisage that the system will need to operate continuously for more than 18 h in the worst case and its battery could be recharged during the sleeping time of a subject. Considering the current battery capacity of a wearable system is 4200 mAh at 3.82 V [19], this means that the ceiling of the current consumption of the system ${\leqslant}$233 mA h⁻¹. Considering this figure, the energy consumption of the system should be ${\leqslant}$890 mW.

The most computationally expensive operation needed in the clustering coefficient calculation is the computation of the geometric mean of all the triangles t_i . This is an iterative process that involves the multiplication of the weights for every possible triangle permutation as in (4) followed by cubic root evaluation. As can be seen in [20, 21] such computation needs significant arithmetic resources. To avoid that, we opt to take a different approach by completely avoiding evaluation of cubic root and simply using the multiplication between the weights as the numerator in (9) as in essence the cubic root will only contribute in terms of a scaling factor and not in terms of the actual nature of the t_i and hence for the value of clustering coefficients.

$$\begin{equation} t_i = \frac{1}{2} \sum_{j,h\in N} w_{ij}w_{jh}w_{ih} \ ;\ CC_i = \sum_{i\in N} \frac{2t_i}{k_i(k_i-1)}. \end{equation} \tag{ 9 }$$

To verify the impact of such modification in the definition of t_i , we carried out a correlation analysis of the clustering coefficient calculated using (4) and (9) withdrawing five-minute real EEG signals from ten different subjects using data obtained from [22]. The data was recorded by 23-channels at 256 Hz with 16-bit resolution. The results are shown in table 2. The average correlation between these signals is 0.923, which indicates that even though the scale of the signal is different, the morphology of the signals remains very close to the original and therefore could be used in the characterization of FC matrices.

To reduce the complexity even further, we utilize the symmetry properties of the PLI matrix. This leads to the computation involving only the elements of the upper triangle of the PLI matrix in the geometric mean calculation of t_i resulting in a reduction of computational complexity by half. This optimized geometric mean of t_i was then used to calculate clustering coefficients as well as the degree k_i (which was also used in the calculation of the clustering coefficient).

Since the transitivity calculation is similar to the clustering coefficient, we reuse the geometric mean and the degree for deriving transitivity. To validate that the modifications done to the clustering coefficient do not affect the calculations of the transitivity, the same correlation analysis has been done and its results are shown in table 2. The average correlation is 0.9405, which demonstrates that not using the cubic root on the triangle calculation does not have a substantial impact on the transitivity either.

Exploiting the symmetric properties of PLI, the summation of the binary degree k_i is equal to two times the summation of the non-zero values in the upper triangular elements of the matrix. Thus, it is possible to calculate the density as in (10) with the already calculated degree, reducing the arithmetic resources required.

$$\begin{equation} D = \frac{\sum_{i\in N}k_i}{N(N-1)}. \end{equation} \tag{ 10 }$$

The architecture used for the PLI calculation is the same as the one presented in [10] and it is shown in figure 3. The main components of this architecture are the Analytical signal, Phase calculation and PLI calculation modules. The Analytical signal module performs a fast Fourier transform (FFT) to the EEG signal followed by a convolution with a mask function that doubles the power for frequencies higher than zero, keep the same power when the frequency is equal to zero, and makes the power equal to zero for the rest of the frequencies. After this, an inverse FFT is performed to obtain the analytical signal. The analytical signal is a complex signal in which the real component is the original signal (x(t)) and the imaginary component ($\tilde{x}(t)$) corresponds to the Hilbert transform. Then, the instantaneous phase is obtained by calculating the arctan function of the analytical signal. This is done using a COordinate Rotation DIgital Computer, which is a low complexity computing technique commonly used to perform arithmetic operations in hardware. Finally, the PLI calculation module uses these instantaneous phases to calculate the PLI value as in (1). This architecture was extended to 19-channels in this work instead of 16-channels as originally done in [10].

Zoom In Zoom Out Reset image size

As explained in section 3, the clustering coefficient module is considered as a core module that calculates not only the clustering coefficient but also the geometric mean of the triangles and the degree (binary and weighted). The block diagram of this module is presented in figure 4, consisting of three basic processing blocks: geometric mean, degree and weighted degree. The geometric mean t_i was calculated first as in (9), multiplying and accumulating the weights of all the possible permutations for every node i in the matrix. At the same time, the degree (binary and weighted) are calculated in the respective blocks (explained in a later subsection). Once the geometric mean of all the permutations in node i have been done, the values are divided by $k_i (k_i-1)$. The same process is repeated for every node and the results are accumulated until the final clustering coefficient value is obtained.

Zoom In Zoom Out Reset image size

The degree and the weighted degree are calculated as in (2), accumulating the number of edges greater than zero and the weights of the edges respectively for every node in the matrix. Both degrees are calculated in the clustering coefficient module. By doing all triangles permutations for every node, we guarantee that all the edges of this node will be used at least once, therefore it is possible to reuse these weight values in the calculation of the degree. By doing this it is not necessary to include an extra copy of the PLI values and the extra logic to perform the degree calculation in a separate process. The only disadvantage of this is that the degree calculation would take longer than it would have done using a separate process, but this does not have a critical impact on the performance since the degree result is still available before the other calculations, e.g. the clustering coefficient.

As mentioned before, the density module reuses the previously computed values of the degree in the clustering coefficient module, and thus consists only of an accumulator and divider as shown in figure 5. The density is calculated as in (10), accumulating the degree values of every node followed by a division by the total number of elements outside the main diagonal. Since the total number of elements is constant, we precompute the inverse of it and multiply this inverted value with the output of the accumulator to perform the division operation. This approach reduces the total arithmetic complexity as opposed to using an arithmetically complex full divider.

Zoom In Zoom Out Reset image size

The transitivity calculation is almost identical to the clustering coefficient with the exception that the numerator and denominator are accumulated before performing the division. Thus, the triangle and the degree values from the clustering coefficient calculation are reused here and the only extra hardware needed is an accumulator for the numerator, a multiplier and accumulator for the denominator and the final divider to perform the calculation. The architecture for the transitivity is shown in figure 6.

Zoom In Zoom Out Reset image size

To perform the CPL calculation it is necessary to map the weights of every edge in the matrix to an equivalent distance. Here, we decided to perform the mapping as in (11). With this mapping, the weight values closer to one which represents a strong connection on the brain connectivity matrix will have a small distance value. On the other hand, values close to zero will have a greater distance value representing a weak connection between the nodes.

$$\begin{equation} f_{ij}(w_{ij}) = 1 - w_{ij}. \end{equation} \tag{ 11 }$$

Once the mapping is done, the shortest path is found using DA discussed in section 2. The shortest path finding unit in figure 7(a), performs the distance mapping shown in (11) and compares the actual minimum distance from the node against the particular node distance plus the current node distance, as it is done in step 4 in section 2. The Minimum finding unit finds the node with the minimum distance that will be used as the new node (step 5) and the current node distance is updated with this node distance. The process is repeated until the distance from the initially selected node and the rest of the nodes is found. Then, the whole process is repeated with a different initial node and the minimum distance from that node with respect of the other is found. When all the minimum distances from all the nodes with respect to the other nodes is completed, these values are accumulated and later divided by 'N²−N' as in (8), obtaining the final CPL value.

Zoom In Zoom Out Reset image size

To calculate these three measurements, it is necessary to calculate first the shortest distance between nodes. This is already done during the calculation of the CPL, then these values are stored and reused here. Then the calculation is just as simple as finding the maximum values for every node (eccentricity) and finding the maximum and minimum of the whole network (diameter and radius respectively). This is done by checking node distance and using a comparator to determine if the new value becomes a maximum/minimum against the previous value examined.

The validation was done using a five-minute, 19-channel EEG signal recorded at 256 Hz, from ten different subjects. These records were obtained from the database of [22]. These EEG signals are the same used in the parameter optimization for section 3.

The proposed system was simulated using the ModelSim software and the output of this system was compared with the output obtained using the toolbox from [13]. The window size used in the PLI calculation was 128 samples. Since the EEG used was sampled at 256 Hz, it is possible to calculate two windows for every second of the EEG. Therefore, a five-minute EEG signal will generate a total of 600 PLI matrices for every subject. From every PLI matrix, the graph-theoretic parameters mentioned in section 2 are calculated, followed by calculation of the relative error between the hardware architecture and the MATLAB toolbox. The same process was repeated for every one of the ten subjects.

5.2.1. Error calculation

The results of the relative error are shown in figure 8 as a histogram for every parameter calculated across all the windows and all the subjects. Based on the distribution of the histogram fit in this figure, it can be observed that the center of the curve is located close to zero, which indicates that for most of the values obtained the relative error is concentrated in this area. The degree is the one with the maximum relative error range that goes from −7.97% to 7.82% (mean = −0.076%, 3σ = 7.9%). Despite this wide range, most of the errors obtained were concentrated in a range of about −2% to 2%. In addition, the RMSE value was calculated for all the parameters. The maximum RMSE found among all the parameters is 0.123 for the degree, a small value that also proves the correct behavior of the system.

Zoom In Zoom Out Reset image size

5.2.2. Time calculation

The number of clock cycles required to perform the calculation of every connectivity measurement can be found in table 3, where N corresponds to the number of Channels and the ${N-1 \choose 2}$ is a function to obtain all the possible combinations for N − 1 channels for 2 positions.

Table 3. Number of clock cycles required for the calculation of the different parameters.

	Clock cycles
Degree	$N{N-1 \choose 2}$+1
Degree (weighted)	$N{N-1 \choose 2}$+1
Clustering coefficient	$N{N-1 \choose 2}$+7
Transitivity	$N{N-1 \choose 2}$+6
Density	$N{N-1 \choose 2}$+2
Characteristic path length	N2+N+1
Eccentricity, radius, diameter	N2+N+1

In table 3, the measurement with the longer calculation time was the clustering coefficient. The reason for this is the extensive iterative processing done to calculate the geometric mean of all the possible triangles for a specific node. This triangle calculation could be performed in parallel if desired, instantiating as many modules as the number of nodes in the system. Nevertheless, without parallelization the design is still fast enough for real-time implementation, thus an increase in the speed is not necessary and will only increase the resources needed and the power consumption of the system. The total calculation time on MATLAB is 3.74 ms, which is about 28 times higher than ours (131 µs). It is worth mentioning is that the time required to perform the degree is almost as long as the Clustering coefficient because the degree is calculated during the same process. Indeed, the degree is a simple measurement that could be obtained faster than the actual time, but that implies creating another instance of the PLI matrix and extra logic to perform the calculation independently from the clustering coefficient. Since there is no benefit in calculating the degree faster, in this work it was opted to keep the degree calculation in conjunction with the clustering coefficient.

In this subsection, we compare the hardware savings achieved with our approximations compared to a one-to-one implementation of the algorithm. This is shown in table 4 in terms of 2-input NAND gates. In this table we converted all the adders, multipliers, dividers, and comparators complexities in terms of 2-input NAND gates, providing an easier way to compare these elements since each of them may require a different number of logic gates for its implementation. It is important to mention that on this table, the optimized degree is considered zero since it is implemented inside the clustering coefficient module. After all the optimizations the total saving comes to around 25% which translates in an equal reduction of area and power consumption of the final system compared to a one-to-one implementation.

After the functional test was performed and the system was performing as expected, the next step is to perform the synthesis of the system on an FPGA. The FPGA used for this system is the Stratix IV GX EP4SGX230K. The synthesis results for every module are shown in table 5.

The module with the highest resource usage was the PLI calculation, using 61% of the available resources on the board. In comparison to this, the proposed architecture for the calculation of the graph-theoretic parameters requires 9% of the total resources available at the board. Nonetheless, the calculation of the graph-theoretic parameters is slower than the PLI calculation, achieving a maximum frequency of 22.16 MHz. However, it is a well-known fact that although FPGA implementations provide design flexibility, they are inherently not optimized for area, speed and power compared to specialized ASIC implementation. To get a better understanding of hardware complexity, we performed an estimation of the logic gates required for the design by converting all the adders, multipliers, dividers, and comparators complexities in terms of 2-input NAND gates. We did not implement the FFT module here but implemented the full PLI and graph-theoretic parameter extraction circuits. The reason for not implementing FFT is that, there are several optimized implementations of it already proposed over the decades and any could be adopted for the current purpose depending upon the area and power budget. Therefore, for this estimation, we selected the implementation presented in [23] due to its low area and power consumption. The result is shown in table 6 in terms of 2-input NAND gate equivalent.

One important point for this design was to create a low power system. To demonstrate the low power capabilities of the devices a power estimation was obtained using the software provided by the manufacturer of the FPGA. This software guarantees a high level of accuracy when the software is provided with enough switching activity from the simulations. The switching activity was generated running the simulations previously mentioned in section 5.2. The EEG data is provided at the same sample rate at which it was acquired, and the frequency used in the system is 22.16 MHz which is the maximum frequency achieved. The dynamic power consumption of the system in this particular FPGA at a VCC voltage of 0.9 V is 51.84 mW. This total power consumption is important if this design is meant to be implemented as a portable continuous monitoring system for instance. Based on this power consumption and the current battery capacities in portable devices like mobile phones, with a capacity up to 4200 mAh [19], the proposed system could be energized continuously for days. It could be integrated with other systems for a further, more advanced processing of the signal as a standalone system. Despite this promising result, however, as mentioned earlier, it is a well-known fact that FPGAs are inherently not as efficient in terms of power consumption as an ASIC implementation where a design can be optimized from the point of view of power consumption. Therefore, we expect the system to show even better power performance when implemented in ASIC. To establish that, we performed an estimation of the power consumption in a 90 nm TSMC CMOS technology using Synopsis Design Compiler. In this estimation, we used the power reported in [23], which is an FFT module implemented in the same CMOS technology. Our results are presented in table 7. Based on the estimations, the power consumption of the complete system will be around 39.37 mW, which represents a power saving of around 24% from the 51.84 mW estimated for the FPGA implementation. Furthermore, it is important to note that FPGA uses a CMOS technology of 40 nm, less than half the size of the library used for the estimations. Thus, if a similar technology to the FPGA were to be used for the ASIC implementation, further power reduction can be expected.

In literature, only the authors from [11] have presented an architecture that performs a calculation of the graph connectivity measurements. The comparison between this and the proposed system is presented in table 8.

There are two major differences between these two systems. The first one is the inclusion of the module to calculate the PLI calculation in the proposed system. In [11], the authors did not perform the calculation of the instantaneous phase used for the calculation of the PLI because it was assumed that these phases were given as an input. In the proposed system, the input will be an EEG signal instead of phase values, and therefore includes the PLI calculation that was avoided in [11]. The second major difference is the way the graph-theoretic parameters are calculated. The whole architecture was redesigned with the main goal to reduce the number of arithmetic operations in the calculation. This not only allows to reduce the size of the final system but also is greatly reflected in the power consumption of the system. By reducing the number of arithmetic operations used, the total energy of the system is decreased as well. The reduction in resources used can be estimated by comparing the number of logic 2-input NAND gates that both architectures use. The number of gates reported in [11] is equal to 1116 K while in this work the estimate for the complete system came to around 937 K (table 6), about 16% less. Nevertheless, the architecture of [11] does not include the 19 channel PLI calculation and it is only doing the graph-theoretic parameter calculation. Thus, if we compare only the graph-theoretic parameter calculation of this work against [11], the 2-input NAND gates required for this architecture will be around 116 K, representing a reduction of about 89% in the total 2-input NAND gates used. Additionally, even though is not possible to compare the power between both architectures since the CMOS technology is different(130 nm) from ours(90 nm), it is expected that the reduction on logic elements previously discussed will have a similar impact in the total power consumption of the system. Even though there was a reduction in the number of arithmetic elements used for the operations, the proposed system is still able to perform the calculation in a similar time to that of the system presented in [11]. The difference in the calculation time between these two is seven clock cycles only. At the maximum frequency achieved, the total time required to calculate the graph-theoretic measurements in this design is equal to 7.89 µs compared to the 6.7 µs reported in [11]. Although the calculation time is longer in this design, this is a minor tradeoff that is compensated with fewer logic components and energy consumption.