Recursive Filter With Exponential Kernel for Nonstationary Systems

Nonstationary stochastic systems in the Wiener-Kolmogorov sense have properties defined by their moments of probability, entropy, and distribution function. Filtering Theory, in general, describes indirectly, a stochastic system through the processes of parameter estimation and state identification. The objective of this article is to develop a Recursive Filter with an Exponential Kernel (RFEK) to reconstruct the response of a nonstationary stochastic system. To achieve this, first, a system viewed as a Black Box (BB) is analyzed. These systems are those whose internal dynamics are unknown, only their input-output is known from a set of responses measured with respect to a particular excitation. From these measurements and applying the proposed filter, a set of estimated parameters and identified states are obtained as a characterization of the system. Subsequently, a comparison is made between the filter output signal and the reference signal over time; that is, measuring their point-to-point convergence. The convergence of the stochastic reference makes it possible to indirectly observe its stability from a bounded estimation region. As a case study, bioelectric signals of the electroencephalographic (EEG) type are analyzed giving an improved approximation with respect to the Kalman filter results.


I. INTRODUCTION
In this paper, a Recursive Filter with an Exponential Kernel (RFEK) is proposed to create a tool to reconstruct nonstationary stochastic signals produced by Black Box (BB) systems. Within this type of signals are bioelectric signals such as Electrocardiograms (ECG), Electromyograms (EMG), Electrooculograms (EOG) and Electroencephalograms (EEG), these last one being the reference signals for the case study presented on this article.
Representing the signals without disturbances that alter the information is important because it will provide better input for different medical devices to have support for the prescriptions, and the analysis they make, reducing possible errors while introducing a preprocessing step. Therefore, the The associate editor coordinating the review of this manuscript and approving it for publication was Mauro Gaggero . RFEK looks to reconstruct signals achieving the greatest convergence between the reconstructed, and the reference EEG signal, having as future applications the insertion of this filter in the medical area.
The main contribution in this work is to perform a pointto-point parameter estimation, instead of performing it on average like traditional algorithms. This obtains a better description of the signal to be reconstructed. The innovative contribution of this filter is the exponential kernel since it corrects the signal during the estimation process.

A. ESTIMATION AND IDENTIFICATION
The filtering theory seeks to explain the internal condition of a system from two areas of development: the estimation of parameters, and the identification of internal and external states. Together, both areas allow a signal to be reconstructed, predicted, or filtered. Considering that given an excitation to the system a corresponding output is obtained, a transfer function that relates them can be defined, making it possible to propose equations to describe the response from its excitation [1].
The equations considered have associated gains that are related to the internal properties of the system from which the signals come. When describing these gains by means of some filtering technique, an estimation of parameters [2] is being carried out, and they are used to describe the internal dynamics of the system states [3].
One of the applications for these two areas is the Black Box (BB) system, which is a system, as represented in Figure 1, whose only measurable characteristics are the inputs and outputs, whilst the internal dynamics is unknown; thus, there is no direct strategy to model its exact behavior, but approximations through estimations and identification techniques. The complexity of these approximations increases when the system to analyze is stochastic.

B. STOCHASTIC SYSTEM
A stochastic system is one that has a non-deterministic behavior, whose current state is the combination of predictable actions with random elements. These elements have intrinsic calculable properties such as the distribution and probability moments, which allow observing the type and degree of stationarity of the system, as well as its rate of change measured by the correlation and variance between the random elements that integrate the measurable signals of the system [4].
The parameters and internal states of stochastic systems are commonly unknown as they are associated to responses with random dynamics, and since neither what happens within the system, nor the relationship between its variables are known, these systems are called BB. However, such dynamics can be approximated through its statistical properties using filtering techniques [5].
Estimation is a filtering technique that describes the interaction and dynamics of a system with stochastic characteristics [6]. In the Wiener-Kolmogorov sense, the estimation for stationary stochastic signals is carried out by means of techniques such as Least Mean Squares (LMS), Recursive Least Squares (RLS), the Kalman Filter (KF) [7]; or under the block-oriented Hammerstein-Wiener model [8], [9].
On the one hand, Volterra [10] developed a method that reconstructs the internal states of a non-linear system using several transformations [11]. For example, there is the Discrete Cosine Transform (DCT), which obtains coefficients whose base vectors depend on the order of the selected transformation, and not on the properties of the stochastic signal [12]; as well as Wavelets, which consider the fundamental wave and its level of decomposition to identify harmonics that contain noise, but do not estimate any characteristic parameter [13].

C. KERNEL FUNCTIONS
The filtering techniques use functions with different bases, which makes them approximations that relate the available information of the system based on transformations of the linear, polynomial, radial, or sigmoid type, for example. When these functions take the set of available system data as arguments seen as vectors, and operate them in a point-to-point way, then they are considered kernel functions [14]- [16].
Within the digital signal processing techniques, kernel estimation is used since it estimates parameters in a smooth way using the dot product, improving the convergence at the output and without increasing the computational complexity of the operations [17]- [19]. In this sense, the estimation can be used to describe the parameters of the transition matrix, that is, the relationship between the states and their evolution over time for non-linear systems [20].
According to [21] Multi-Output type adaptive kernel filtering algorithms have been proposed, which take advantage of the correlation between the input and the output of a system. A kernel method, as an extension of the Kalman filter for dynamic systems, is presented in [22], which is used to filter noise from non-linear signals and for prediction.
The kernel functions can be combined to give rise to hybrid configurations, such as the one presented in [23] where they present an algorithm for state estimation, merging different kernel functions for non-linear systems with unknown states, applied to neural networks and mentioning the use of global and local kernel functions.

II. FILTERING THEORY IN EEG SIGNALS
Among the bioelectrical signals that have different random variables, and that depend on time-varying parameters, are electroencephalographic (EEG) signals, which are measured using different means, allowing an example of a stochastic process output to be observed.
The importance of studying EEG signals lies in the fact that its description has allowed the analysis and classification of various diseases such as epilepsy, autism, depression, and Alzheimer [24]. For example, the electrodes make it possible to capture EEG signals, whose description is limited by varying distributions and their respective non-stationary probability moments.
The EEG signals, on the one hand, correspond to the composition of the n-neuron signals, which have associated weights that represent the intensity of interaction between each of them, having a high dynamic correlation, and whose values are unknown. On the other hand, these signals have multiple noises associated with the measurement instruments, which have an electric field between 40 and 100 µV, and an amplitude within a frequency band between 0 and 50 Hz. For VOLUME 10, 2022 the latter, the filtering theory proposes strategies to eliminate noise that can affect the interpretation of measurements.
Since both weights and their relationships are unobservable as they are within brain events, a filtering technique is required to describe their interaction and dynamics. Unfortunately, because EEG signals are stochastic with abrupt changes in amplitude and frequency, they exceed any of the aforementioned filtering theory schemes, requiring a different approach to the previous ones.
With this, a series of questions are associated: how to know that the signals received from an EEG are stable? Being a discrete signal, and given that the EEG signals are nonstationary, will it be possible to obtain their time-varying parameters? Is it possible to measure its stability with the rate of change between its estimated parameters? These questions must be resolved to improve the understanding and interpretation of the information flows provided by the brain [25].
Applying the filtering theory in the EEG signals implies seeking to reconstruct the reference signal through the coefficients obtained from an estimation process in such a way that the coefficients can be used in a model that represents said reference signal, even if the model is affected by unknown external noises [26], [27].
Estimation and identification with respect to the reference obtains the equivalent transition matrix to observe the stability of a signal. In general, the identification of non-stationary signals requires the estimation of variant parameters and depends on [28]: • The measurement of stochastic, or non-stationary inputs and outputs.
• The exponential transition matrix calculated by a filtering action.
• The variant profit, obtained from the correlation of identification errors and an innovation process.
In this sense, the estimation of variant parameters serves as a tool to describe what happens with the internal states of the system and with the level of instability presented, improving the training of signal identification systems [29]. In this way, the study of EEG signals, through the estimation and identification, improves the description of what happens internally in the stochastic processes of the brain [30].

III. INNOVATION PROPOSAL
In the time domain, it is possible to track the parameters of the model selected as identifier and describe the internal states that characterize a stochastic system, using the timevarying parameter estimator within the identifier, allowing convergence to a reference, which can be a non-stationary biological signal, through a process, in which different base and argument functions intervene, seen as kernel functions.
Due to the signals of interest for this case study are stochastic, with no defined or constant boundaries, and a quick change in their behavior or tendency is expected, an exponential Kernel is chosen for its properties, giving a faster convergence, and tracking respect to other functions.
The objective is to develop a technique of low computational complexity compared with the existing ones, which is based on the idea that an exponential kernel function can be modified and used as a complement to improve point by point the result of estimating existing strategies, such as the LS strategy; and since the EEG signals are of a nonstationary type [31], and due to the importance of their study, a set of EEG signals is considered as reference signals to characterize the performance of the proposed Recursive Filter with Exponential Kernel (RFEK). This is how the identification with an RFEK is proposed, seeking to converge in almost all points to a reference signal that requires the estimation of variant parameters to know its internal properties and evolution over time, using a kernel function configuration exponentially based, considering as arguments the results of the least squares estimation, different from those found in the state-of-the-art.
When applying the RFEK, its level of convergence and the similarity of its response with respect to the reference signal are evaluated. In this way, it is verified that both the estimation, as well as the stability and internal dynamics are valid, characterizing the signal with the least noise, using the innovation process to counteract it; this being possible since they have the same distribution function.

A. ARTICLE ORGANIZATION
This article is organized as follows. First, the concepts and antecedents that gives a better understanding of the work are introduced. Second, the mathematical model that represents the estimation and identification processes and its relationship with the LSM is described. Third, the results of applying the estimation and identification process to a set of EEG signals are shown. Next, a discussion of the results is presented, and finally, the advantages of the proposal are highlighted in the conclusions section.

IV. MATHEMATICAL MODEL
By considering that the LMS average least squares help to build the transition function represented by the relationship between the input and the output of the system through the estimation of parameters, and nonstationary conditions, in this work a recursive filter is proposed that performs the estimation of variant parameters in time through the identification of states by adding the kernel function (Ker), within the estimation [32].
In Figure 2, a block diagram that indicates the stages of the filter to characterize signals indirectly is presented, starting from the arrangement of an input signal and a system reference signal, both observables. According to the process, the input signal makes it possible to estimate the parameters A k that enter the identifier y k to generate the output signal that approaches the reference, and whose convergence difference gives rise to an error measured through the error function J k . The error function is recursive, using previous and the actual calculated errors, and it is updated in each iteration until the error between consecutive iterations is stationary. In each iteration of the filtering process, an adjustment of parameters is added, taking advantage of the information previously calculated, approximating the variant parameters in time, necessary to give a better follow-up of the reference, being in this adjustment the place where the innovation of the proposal is.
These two main sections of the proposed filter are presented in the following theorems 1 and 2. In the first, the signal identification process is described, while, in the second, the estimation of the time-varying parameters, considering the use of a kernel function that is shown.
Theorem 1 (Signal Identification Process): Let S be a BB system and let y k be an observable signal. If S is described by means of its internal states as (1) and its observable signal as (2), then there is an optimal recursive model in a sense of identification probability associated with S given by (3).
is the vector of coefficients and w k ∈ R [1×1] is the internal disturbance. Bounded in the domain N (µ, σ 2 < ∞).
where y k ∈ R [1×1] is the observable signal, C ∈ R [1×n] is the vector of coefficients, v k ∈ R [1×1] is the external disturbance and ε k ∈ R [1×1] is the measurement error.
whereÃ k ∈ R [n×n] is the unknown parameter matrix associated with the reference frame, andŷ k ,ŵ k ∈ R [n×1] are the identified output signals and noise, respectively. Proof: Let (4) be a model equivalent to the system S as in (1), delayed one instant of time.
where C + is the pseudo-inverse matrix.
By delaying (6) one instant of time, and replacing it in (5), gives the observable state as part of the recursion in (7).
The identification process (8) considers the signals of the system in a time k − 1 without the intervention of the current internal statex k .
Let ε k be the identification error of S, defined as the difference between the identified signal and the observable signal as in (9), so, the equivalence in (10) is given.
The identification error functional, J k = 0, to be minimal is defined by (11), and it is given in a recursive way by (12).
Evaluating the stochastic gradient with respect to the output y k that minimizes the functional of the error J k , there is the optimal recursive identification model as shown in (3).
Considering a kernel function Ker(v i ) as a linear operator that relates the variables v i , operating point to point, and maintaining the domain of the estimation and identification results, Ker : R [n×n] → R [n×n] , due to the rate of change required to follow the stochastic signals, the base selected for the kernel is of the exponential type, defined as in (13), where Ker(·), represents the kernel function, ε k is the identification error,Ã k the estimated parameters, sign (·) indicates the sign function, and e [·] indicates the exponential function. The use of (13) to correct the results from an average estimation, to approach the changes in the reference in an exponential rate, instead of a lineal one is as follows in Theorem 2. Theorem 2 (Estimation using the kernel function): Let S be a BB type system and letŷ k be the output signal identified as in (3). IfÃ k is the unknown parameter matrix associated with S obtained by minimizing the error functional ∇J k = 0 then VOLUME 10, 2022 there is an estimator with kernel Ker Ã k , ε k that corrects the estimate according to (14).
Proof: The identification error functional is defined by (11) and recursively expressed in (12). Developing the quadratic error there is (15).
Its stochastic gradient, with respect to the matrix of parameters ∇J k |˜Ã k is observed in (16).
Considering the kernel function in the estimator to update each of the parameter values, obtains (21), and considering Equation (14) within Equation (3) gives (22), With (22) being the filter that helps to reconstruct nonstationary signals by means of the identifier (3) and the parameter estimation (21). Thus, these expressions are applied to electroencephalograms (EEG), as an example of a non-stationary bioelectric signal, defining the performance of the proposal, as well as the performance with respect to existing filters in the reconstructions of said signals.

V. RESULTS: IDENTIFICATION OF EEG SIGNAL
To study the performance of the proposal described in the previous section, the estimation and identification method is applied to estimate the parameters of a set of electroencephalograms (EEG) signals obtained from a public database, which can be consulted in [33]. This database contains a set of EEG signals sampled from seven people in a baseline state while performing tasks related to human activities: multiplying, writing letters, counting, and moving.
Each subset of data considered for this work is made up of six EEG signals, sampled at 250 Hz for 10 seconds. The six responses correspond to channels C3, C4, P3, P4, O1, O2, which are illustrated in Figure 3. This nomenclature is in accordance with the 10/20 positioning system, which is an internationally recognized method to describe the location of the electrodes on the scalp [34]- [36].
For the validation of the proposal and obtaining results, only the six EEG signals of a single subject (Subject 1) associated with execution 1 (trial 1) of the ''counting'' activity are considered, which is a daily activity in the mailing office. In Figure 4 these signals are presented, showing the data corresponding to a period of 10 seconds, having a total of 25 00 points for each. The signals are shown separated and not superimposed to appreciate their variation over time and to distinguish the result in each channel.
From the estimation of parameters with the kernel (Ker) function and the identification of states developed in (22),  it is possible to know how to reconstruct the EEG signals, which have different variations in amplitude and frequency over time.
By applying the estimation and identification proposal to the six signals taken as reference and presented in Figure 4, the results shown in Figure 5 are obtained, which represents a comparison of the identification using an input signal with random properties, with respect to the original reference, showing the normalized signals to appreciate both the variation between signals, as well as the similarity in convergence.
Within the Figure 5, blue lines correspond to the reference signals, while the orange lines are the identified signals. All reference signals are presented normalized.
Each signal obtained has its distribution and stability according to the parameters that characterize it, and discretely these should not leave a region with an amplitude of [-1,1]. The parameters obtained from the identified signals of Figure 5 are observed in Figure 6.
Because the estimation approach is for a Black Box (BB), the parameters shown in Figure 6 do not represent a physical quantity of the reference system; therefore, their values are  interpreted as coefficients that are operated in scalar form at the input to the model over time. From the second moment of probability, as mentioned in Theorem 1, it is possible to have a measure of the rate of convergence between the reference and the identification with the estimated parameters.
Taking into account the reference EEG signals, applying the error functional described in (12), the level of convergence achieved is determined. The results of the functional error of the six signals are graphically presented in Figure 7. Table 1 shows, for each of the six EEG signals, seven values of those shown in Figure 7, allowing a numerical appreciation of the reduction of the error second moment of probability as the number of samples considered increases.
This is useful to corroborate that the error decreases indicating an improvement in convergence. When considering (12), it is observed that the values obtained from the six signals are of the order of thousandths, that is, ×10 −3 units.
To represent the second recursive probability moment of the six EEG signals in Figure 7, it is considered using a superscript i indicating which of the signals is being referred to, represented as J i k . Being the highest value J 5 k = 1.4 × 10 −3 . After obtaining the results with the RFEK proposal, then the same task is performed using the ARX, Hammerstein-Wiener and Kalman filter functions from the MATLAB R library modules, with the default initial setting; as well as the compare function to measure the percentage of approximation.
Having as reference the EEG signal number five (O1), the Figure 8 (above) shows the identification of the signal, and (below) the measure of the errors obtained with the four methods. The comparison results correspond to the case that had the highest response among the six considered signals.
As observed in Figure 8, the Kalman, ARK and RFEK filters have a good follow-up to the reference signal, being this last one with the highest percentage of approximation respect to the others.
It is important to remark that the initial setup for the other filters affects the results, being this a difficulty when the signals have variations or when the same method is applied to different references, making necessary a better preprocessing, and setting for these methods to achieve their best performance [34]- [37].

VI. DISCUSSION
By performing the parameter estimation described in (14) and with the identification according to (22) it is possible to reconstruct the non-stationary EEG bioelectric signal or any other type of signal with these characteristics. The kernel function application gives a punctual monitoring of the signal and not an average, thus achieving a greater convergence than with the ARX, Hammerstein-Wiener and Kalman filter methods as observed in Figure 8.
In Figure 7 it can be seen that the maximum value of the second recursive probability moment of the error corresponds to J k of signal 5, for sample k = 750 (1.4×10 −3 units), which will be described as ε. In the other samples, this value is no longer exceeded, indicating a reduction in the error, as shown in Table 1.
Some difficulties in the actual work are the signal acquisition to have data for validation, which is considered as a near future task. Thus, this proposal can help to improve the standard classification of signals from the different organs and senses of the human body, contributing as a tool that helps specialists in making decisions regarding a specific case.
The above, without only using tree relationships, hierarchies, or deep learning, which is developed for such purposes without having probabilistic models that analytically interact for various specialists to make decisions about a specific problem, and that leads to an answer with less divergence as is the case today.
Moreover, it is still necessary to go deeper into stability since it does not only have to do with the amplitude, but it is also necessary to consider the speed of change of the parameters. All this to be able to describe the different bioelectrical signals more precisely.
Finally, for further work it is also proposed examining and analyzing the stability and usability of this proposal in other engineering applications, as in determining fuzzy parameters or delayed hybrid systems as the one mentioned in [38].

VII. CONCLUSION
In this article, an RFEK for the identification of nonstationary signals was developed, thus having a tool that knows indirectly the internal dynamics of the Black Box (BB) reference system with non-stationary responses.
The estimation with the kernel function together with the identification, allowed filtering with a point-to-point convergence with respect to the six EEG signals.
In addition to seeing that the stability conditions were met in almost all the samples as they were limited in amplitude within the unit interval, it was observed that a relationship in frequency is established between them since they were of the periodic type in the face of a specific brain activity.
Future work will show if the speed of change concerns the stability of the EEG signal. In addition, the integration of fuzzy logic and other artificial intelligence techniques will be sought so that through them, a better description of the internal evolution of the brain will be achieved when faced with certain activities.
ROSAURA PALMA OROZCO received the Ph.D. degree in advanced technology from the CICATA, IPN. She is currently a Professor at the Escuela Superior de Cómputo, IPN. She is a member of the National System of Researchers. She has several publications in the area of stochastic linear filtering and recursive filtering, specifically in dynamic estimation. Her research interests include the computational and mathematical applications of filtering theory for analysis and study of complex systems represented in state space. She became a member of the Biophysical Society, in 2008.
KAREN ALICIA AGUILAR CRUZ received the B.S. degree in aeronautical engineering, the M.S. degree in computational engineering, and the Ph.D. degree in computing science from the Laboratory of Intelligent Systems for Automation at the Research in Computing Center (CIC-IPN), Instituto Politécnico Nacional, Mexico, in 2021. Her research interests include signal processing, identification, and estimation.
ROMEO URBIETA PARRAZALES was born in Tapanatepec, Oaxaca, Mexico, in 1947. He received the M.S. degree from the Instituto Politécnico Nacional (IPN), in 1990. From 1971 to 1987, he worked on diverse disciplines, such as electronics, automatic control design, and industrial processes fields for Mexican Government departments. Since 1977, he has been a Professor at the CIC, IPN, working on industrial intelligent control systems and stochastic filters using FPGA. From 1988 to 1996, he was a Professor at the CIDETEC, IPN, developing research on intelligent control and computing applied to industrial processes. He is the author of three books, and a dozen of articles in journals with strict arbitration.
JOSÉ LUIS FERNÁNDEZ MUÑOZ received the degree in physics from the Autonomous University of Zacatecas, the master's degree from the CINVESTAV, IPN, and the Ph.D. degree from the CICATA Legaria. He is currently working on research in applied physics and materials nanotechnology with the CICATA Legaria. He is a member of SIN level 1.