
\documentclass[12pt, a4paper, twoside, openany]{report}

% ==========================================
% PACKAGES
% ==========================================
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[english]{babel}
\usepackage{graphicx}
\usepackage{amsmath, amssymb, amsfonts}
\usepackage{geometry}
\geometry{a4paper, top=3cm, bottom=3cm, left=3.5cm, right=3cm}
\usepackage{fancyhdr}
\usepackage{titlesec}
\usepackage{cite}
\usepackage{setspace}
\usepackage{algorithm}
\usepackage{algorithmic}

\usepackage[rgb]{xcolor} 

\usepackage[a-2b,mathxmp]{pdfx} 

% ==========================================
% HYPERREF & ANNOTATION FLATTENING SETUP
% ==========================================
\hypersetup{
    pdfpagemode={UseOutlines},
    pdfstartview={FitH},
    bookmarksopen=true,
    bookmarksnumbered=true,
    % 3. FIX: ELIMINATE INVISIBLE BOX ANNOTATIONS COMPLETELY
    colorlinks=true,      
    linkcolor=black,      
    filecolor=black,
    urlcolor=black,       % Temporary change to black to guarantee no uncalibrated color spaces leak
    citecolor=black,
    pdfa=true
}

\setstretch{1.15} % Line spacing

% ==========================================
% DOCUMENT START
% ==========================================
\begin{document}

\begin{titlepage}
    \centering
    % Ensure you have a file named 'polito_logo.png' in the same folder
    \includegraphics[width=0.5\textwidth]{polito_logo.png}\par\vspace{1cm}
    
    {\LARGE \textbf{POLITECNICO DI TORINO}\par}
    \vspace{0.5cm}
    {\Large ICT for smart societies\par}
    \vspace{1cm}
    {\Large Master's Degree Thesis\par}
    \vspace{1.5cm}
    
    {\huge \textbf{Kalman Filters to Denoise EEG Signals: A Time-Varying Autoregressive Approach for Event-Related Potential Enhancement}\par}
    \vspace{1.5cm}
    
    \begin{flushleft}
        {\large \textbf{Supervisors:}}\par
        \large Prof. Roberto GARELLO\par
        \large Prof. Vito DE FEO\par
    \end{flushleft}
    
    \vspace{1cm}
    \begin{flushleft}
        {\large \textbf{Candidate:}}\par
        \large Hadi GHAVIPEYKAR SHIRINSOO\par
    \end{flushleft}
    
    \vfill
    {\large \textbf{July 2026}\par}
\end{titlepage}

\pagenumbering{Roman}

% ==========================================
% ABSTRACT
% ==========================================
\chapter*{Abstract}
\addcontentsline{toc}{chapter}{Abstract}

Brain-Computer Interfaces (BCIs) use Event-Related Potentials (ERPs)---which are measurable brain responses to different sensory or motor events---to translate neural activity into actionable features for communication and control. Inside these ERPs, the Readiness Potential (RP) is a very important biomarker to identify motor intentionality and movement preparation. However, extracting these fast neural markers from raw electroencephalographic (EEG) data requires a strong artifact suppression and advanced state-space tracking. 

This thesis investigates and tries to optimize the BCI pre-processing architecture of BVibes, a specialized signal processing toolbox developed at the University of Essex to which this project contributes. We aim to improve this toolbox by adding a Time-Varying Autoregressive (TVAR) Kalman smoother. Rather than just applying existing algorithms one after the other, this work analyzes the actual mathematical interactions between frequency filtering, blind source separation, spatial filtering, and state-space estimation across different clinical datasets.

Through our tests, we noticed that using a 10 Hz zero-phase FIR filter together with standard Infomax ICA performs much better than Infinite Impulse Response (IIR) filters and hybrid Wavelet/SOBI methods (FORCe) for offline RP extraction. More importantly, this research shows that applying the Kalman filter on continuous data before the artifact correction actually disrupts ICA, because it inadvertently "Gaussianizes" the physiological spikes. We also simulated different non-linear filters (Extended and Unscented Kalman Filters) and found that they do not offer any real mathematical advantage over the standard Linear Kalman Filter for TVAR tracking, meaning we can save a lot of computational time.

A major contribution of this work is finding an antagonistic relationship between traditional Jitter Compensation techniques and state-space tracking. While shifting single-trial epochs in time can improve a simple average, it breaks the temporal phase continuity of the background noise, causing the Kalman tracking to fail. Because of this, our optimized architecture keeps Jitter Compensation disabled. Furthermore, the good noise suppression achieved by the FIR filter and the post-epoching KF-TVAR loop allows to isolate clean topographies of the Lateralized Readiness Potential (LRP), making the complex Spherical Spline Interpolation Spatial Filtering (SISF) mostly redundant. In the end, this optimized and lighter pipeline provides classifiers with better features, helping to improve the overall reliability and accuracy of BCI applications.




% ==========================================
% TABLE OF CONTENTS
% ==========================================
\tableofcontents
\listoffigures
\listoftables


% ==========================================
% ACKNOWLEDGMENTS
% ==========================================
\chapter*{Acknowledgments}
\addcontentsline{toc}{chapter}{Acknowledgments}

First and foremost, I would like to express my sincere gratitude to my supervisors, Prof. Roberto Garello and Prof. Vito De Feo. I deeply thank them for giving me the invaluable opportunity to work on such an engaging and challenging Brain-Computer Interface project, and for their continuous guidance, expertise, and support throughout the completion of this Master's degree at Politecnico di Torino.

A special thank you goes to Cristina Del Prete, with whom I collaborated closely for almost the entirety of this project. Her patience, brilliant insights, and constant availability were a tremendous help to me, and I am incredibly grateful for her daily support and friendship during this journey.

Finally, I would like to thank my family and friends for their endless encouragement, sacrifices, and emotional support during my university studies. This milestone is dedicated to all of you.

\cleardoublepage
\pagenumbering{arabic} % Resets page counter to 1 and uses 1, 2, 3...

% ==========================================
% CHAPTER 1
% ==========================================
\chapter{Introduction and Background}

\section{Introduction}
Patients suffering from severe brain injuries often experience the sad condition of being conscious, but they are completely unable to interact or communicate with the outside world \cite{BonomiThesis}. Usually, these severe injuries lead to a period of coma. Individuals who wake up from a coma frequently survive in a non-responsive or vegetative state, and this makes clinical diagnosis very difficult. In these situations, the lack of communication prevents doctors from having a complete understanding of the patient's actual state of consciousness. Unfortunately, this often leads to misdiagnoses that can compromise the choice of the best medical treatments.

To obtain more reliable assessments, Brain-Computer Interfaces (BCI) represent a very promising strategy \cite{Kostoglou2024}. BCIs offer a non-invasive method that continuously analyzes electroencephalographic (EEG) and electromyographic (EMG) signals. The main goal of these interfaces is to extract useful information about the patient’s intention to move, which helps clinicians to distinguish between voluntary, and involuntary motor tasks \cite{Luck2014}. 

However, raw EEG data is basically a noisy mixture of neural, physiological, and environmental signals. The brain's responses to these cognitive events—known as Event-Related Potentials (ERPs)—are deeply buried under this noise \cite{Mohseni2007}. Specifically, the electrical signals of motor intention, like the Readiness Potential (RP), are often masked by ocular artifacts, muscle twitches, and the normal background activity of the brain. Because of this, my thesis focuses on advanced pre-processing algorithms and filtering techniques to clean the raw EEG signals. The final objective is to reliably extract these ERPs using advanced mathematical tracking, specifically using a Kalman filter combined with a Time-Varying Autoregressive (TVAR) model, to generate better features for BCI classification \cite{Khan2006, Tarvainen2001}.

In this project, rather than just doing a simple software implementation, we systematically investigate and optimize the whole BCI pre-processing chain. Through extensive simulations on datasets with different motor tasks, we compare multiple pipeline configurations. We look at optimizing frequency filters using zero-phase FIR filters \cite{Widmann2012} and compare traditional Infomax Independent Component Analysis (ICA) \cite{Bell1995} against newer hybrid Wavelet/SOBI algorithms (like FORCe) \cite{Daly2015} for offline artifact suppression. 


\section{Background and Literature Review}

\subsection{Neurological Definitions of Consciousness}
To better understand this clinical challenge, it is necessary to define consciousness neurologically. According to foundational neuroscience, it involves two primary components: \textit{arousal} (phenomenal consciousness or wakefulness, clinically marked by spontaneous eye opening) and \textit{awareness} (focal consciousness, describing the subjective experience and the ability to follow commands) \cite{Koch2018}. Disorders of consciousness (DOC) are broadly classified into Coma, Vegetative State (VS), and Minimally Conscious State (MCS). Furthermore, conditions like Locked-In Syndrome (LIS) leave patients fully aware and vigilant but suffering from quadriplegia, effectively trapping their active minds in unresponsive bodies \cite{Koch2018, DelPreteThesis}. 

\subsection{The Mystery and Definitions of Consciousness}

\subsubsection{The Philosophical and Neurological Divide}
The concept of consciousness has historically been a topic of profound philosophical debate, broadly divided between Dualism (the separation of the immaterial mind and the material body, championed by Descartes) and Materialism (the belief that mental states are entirely identical to scientifically provable brain states) \cite{DelPreteThesis}. However, explaining the exact physical processes that give rise to subjective feelings and sensations—often referred to as \textit{qualia}—remains an elusive challenge. Philosopher David Chalmers famously defined this as the ``hard problem of consciousness,'' separating it from the ``easy problems'' such as the ability to discriminate stimuli or the mechanical difference between wakefulness and sleep \cite{Chalmers2007, BonomiThesis}.

\subsubsection{Neuronal Correlates of Consciousness (NCC)}
To bypass the purely philosophical ``hard problem,'' modern neuroscientists focus on the Neuronal Correlates of Consciousness (NCC), defined as the minimal neuronal mechanisms jointly sufficient for any one specific conscious experience \cite{DelPreteThesis}. The NCC is not a single unified entity but is composed of several distinct fragments \cite{BonomiThesis, DelPreteThesis}:
\begin{itemize}
    \item \textbf{Level of Vigilance (Arousal):} The state of alertness and wakefulness, clinically marked by spontaneous eye opening. It is primarily modulated by the reticular formation in the brain stem.
    \item \textbf{Contents of Consciousness (Awareness):} The subjective experience, thoughts, and feelings occupying the mind at a given moment, managed by the specific and associative areas of the cerebral cortex.
    \item \textbf{Selective Attention:} The mental concentration focused on a given stimulus, controlled by the posterior parietal and prefrontal cortices.
    \item \textbf{Self-Memory:} The continuous comparison between current sensorial inputs and past experiences to form a cohesive self-identity.
\end{itemize}

\subsubsection{Leading Neuroscientific Theories: GNW and IIT}
Currently, two primary theories attempt to model how these neurological fragments combine to form consciousness \cite{DelPreteThesis}:
\begin{enumerate}
    \item \textbf{Global Neuronal Workspace (GNW):} Proposed by Bernard J. Baars, this theory states that consciousness arises only when sensory information is globally broadcasted to multiple cognitive systems simultaneously. If a signal remains localized in a specific sensory or motor system, it remains unconscious.
    \item \textbf{Integrated Information Theory (IIT):} Developed by Giulio Tononi, IIT posits that consciousness is the intrinsic causal power of a complex system. It relies on two properties: \textit{Integration} (every conscious experience is a unified, indivisible scene) and \textit{Differentiation} (a conscious state is a highly specific selection from billions of possible states). IIT introduces a mathematical measure of consciousness, $\Phi$ (Phi); a system with zero $\Phi$ is entirely unconscious, whereas the human brain's massive interconnectivity yields an extremely high $\Phi$ \cite{Tononi1998, DelPreteThesis}.
\end{enumerate}

\subsubsection{Clinical Disorders of Consciousness}
Understanding these neurological definitions is critical when treating severe traumatic brain injuries. When the Neuronal Correlates of Consciousness are disrupted, patients fall into various Disorders of Consciousness (DOC). These are broadly classified into Coma (neither aroused nor aware), Vegetative State (VS, where patients regain arousal and sleep-wake cycles but show no focal awareness), and Minimally Conscious State (MCS, where minimal but reproducible awareness is present) \cite{Koch2018}. Furthermore, conditions like Locked-In Syndrome (LIS) leave patients fully aware and vigilant (high $\Phi$ and intact GNW) but suffering from severe quadriplegia, effectively trapping their active minds in unresponsive bodies \cite{Koch2018, DelPreteThesis}. 


\subsection{State-of-the-Art in BCI Signal Processing}
Because purely behavioral assessments fail to distinguish between these states, BCI tracking is utilized. Grounded in the principle that consciousness is closely tied to motor cognition, monitoring the brain's preparatory electrical activity before an imagined movement can verify if intentionality remains preserved \cite{Luck2014}. 

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{libet_timeline.png}
    \caption{Timeline of motor intention depicting the Readiness Potential preceding conscious awareness and the actual muscle contraction .} % \cite{DelPreteThesis}
    \label{fig:libet_timeline}
\end{figure}


Extracting these signals requires advanced multivariate time-series modeling. Table \ref{tab:lit_review} summarizes the primary state-of-the-art methodologies currently employed in BCI applications to denoise and track cerebral physiologic signals.

\begin{table}[htbp]
    \centering
    \renewcommand{\arraystretch}{1.3}
    \caption{Literature Review: Multivariate Time-Series Models for BCI and Cerebral Physiologic Signals .} % \cite{Vakitbilir2025, Kostoglou2024, Khan2006, Tarvainen2001}
    \vspace{0.2cm}
    \resizebox{\textwidth}{!}{
    \begin{tabular}{|p{3cm}|p{4.5cm}|p{6.5cm}|}
        \hline
        \textbf{Author \& Year} & \textbf{Model / Algorithm} & \textbf{Primary Contribution / Findings} \\
        \hline
        \hline
        \textbf{Tarvainen et al. (2001)} & Time-Varying ARMA (TVAR) & Demonstrated that using a Kalman smoother eliminates the temporal tracking lag inherent in adaptive algorithms like RLS. \\
        \hline
        \textbf{Khan \& Dutt (2006)} & EM-Kalman Smoother (EMKS) & Introduced the Expectation-Maximization algorithm to automatically optimize hyperparameters ($Q, R$) without manual tuning. \\
        \hline
        \textbf{Ming (2012)} & Sequential Monte Carlo (SMC) & Applied non-linear, non-Gaussian particle filters to successfully isolate single-trial ERP dynamics from severe background artifacts. \\
        \hline
        \textbf{Kostoglou et al. (2024)} & Genetic Algorithm (GA) KF-TVAR & Proved that extracting the time-varying trace of the Kalman covariance matrix significantly improves BCI classification accuracy. \\
        \hline
        \textbf{Vakitbilir et al. (2025)} & VAR \& Fractional Integration & Reviewed state-space applications for cerebral modeling, highlighting the Kalman filter's superior capability in capturing transient ERD/ERS events. \\
        \hline
    \end{tabular}}
    \label{tab:lit_review}
\end{table}

Historically, standard pipelines employ Independent Component Analysis (ICA) specifically the \textit{runica} algorithm based on information-maximization \cite{Bell1995}—to separate artifacts. Slow drifts are handled via robust detrending frameworks \cite{DeCheveigne2018}. However, tracking the highly transient cognitive shifts associated with motor execution requires the transition from non-parametric static filters to dynamic state-space models like the Kalman filter \cite{Khan2006}.


\section{Clinical and Real-World Applications of Advanced EEG Processing}

\subsection{The Core Application of the Proposed BCI Framework}
The primary application of the methodology developed in this thesis is the creation of a non-invasive, low-cost Brain-Computer Interface (BCI) designed as an objective diagnostic tool. Following severe traumatic brain injuries, patients may emerge from a coma but survive in a Vegetative State (VS), a Minimally Conscious State (MCS), or Locked-In Syndrome (LIS). In these conditions, patients may be fully awake and aware, but profound paralysis prevents them from interacting with the outside world \cite{BonomiThesis, Koch2018}. Because standard clinical assessments rely entirely on observed physical behaviors, they frequently lead to devastating clinical misdiagnoses.

By applying an advanced preprocessing pipeline, this project rigorously cleans heavily corrupted EEG signals.


\subsection{The Clinical Significance of Event-Related Potentials (ERPs) and the Readiness Potential (RP)}

Event-Related Potentials (ERPs) represent the macroscopic neural responses to specific sensory, cognitive, or motor events. They are typically extracted from continuous EEG through rigorous time-locked averaging \cite{Cohen2014, Mohseni2007}. In clinical neuropsychology, ERPs are utilized as broad, objective biomarkers to assess various stages of brain processing. Because different ERP components map to distinct cognitive functions—such as visual processing or auditory attention—they are invaluable for assessing the overall functional integrity of the brain in non-responsive patients, entirely bypassing the need for physical communication \cite{Mohseni2007, Luck2014}.

Within this broader spectrum of ERPs, Motor-Related Cortical Potentials (MRCPs) are specifically tied to movement. The most clinically significant of these for Brain-Computer Interfaces is the Readiness Potential (RP) \cite{DelPreteThesis}. RPs are slow, negative-shifting electrical potentials generated in the motor cortex that arise up to 2 seconds before an actual voluntary muscle contraction occurs \cite{Luck2014}. While general sensory ERPs prove that a patient can passively process external stimuli, RPs act as direct electrophysiological biomarkers of active motor planning and conscious intent.

If a paralyzed patient is asked to imagine or attempt to move their hand, the appearance of an ERP—specifically the RP component—in their EEG trace objectively proves they heard the command, understood it, and consciously formulated a motor plan \cite{DelPreteThesis}. Furthermore, this project specifically isolates the Lateralized Readiness Potential (LRP). By applying a double-subtraction method across the pre-motor cortex channels (e.g., Fc3 and Fc4), the LRP verifies that conscious motor programming remains fully preserved in the undamaged brain hemisphere of hemiplegic patients \cite{BonomiThesis}.

\subsection{Broader Real-World Applications}
Beyond the specific diagnosis of consciousness in paralyzed patients, the advanced mathematical tools utilized in this project represent the cutting edge of biomedical signal processing. Currently, these exact state-space tracking and artifact rejection processes are deployed across various medical and engineering domains \cite{Vakitbilir2025}:
\begin{itemize}
    \item \textbf{BCIs for Communication and Prosthetics:} Algorithms like the Informed Kalman Filter (IKF) and FORCe are used to decode motor intentions and Event-Related Desynchronization (ERD) in real-time. This allows severely paralyzed individuals to steer robotic arms, navigate motorized wheelchairs, or type using P300 spellers \cite{Daly2015, Enshaeifar2016}.
    \item \textbf{Epilepsy Diagnosis:} Kalman filters and TVAR models are heavily utilized to predict background brain rhythms and isolate non-stationary anomalies. This drastically enhances the Signal-to-Noise Ratio (SNR) of hidden epileptic spikes, aiding in accurate seizure-onset mapping and reducing false positives in ambulatory monitoring \cite{Vakitbilir2025}.
    \item \textbf{Neurocritical Care:} In Intensive Care Units (ICUs), multivariate time-series models track continuous streams of cerebral blood flow, intracranial pressure, and tissue oxygenation in Traumatic Brain Injury (TBI) patients. These models reveal dangerous breakdowns in the brain's autoregulatory mechanisms, guiding life-saving clinical interventions \cite{Vakitbilir2025}.
    \item \textbf{Psychiatric Disease Diagnosis:} Advanced time-frequency EEG analyses are increasingly used to discover biomarkers for mental illnesses (such as schizophrenia, autism, and depression). Researchers use these machine-learning-driven EEG features to build massive databases aimed at predicting which pharmacological treatments will be most effective for individual patients.
\end{itemize}




\section{Fundamentals of Electroencephalography (EEG)}

\subsection{Neural Origins and the Equivalent Current Dipole}
Electroencephalography (EEG) is a prominent, non-invasive technique used to acquire brain potentials reflecting human cognitive and motor behavior. Discovered by Hans Berger in 1924, it was proven that the synchronous electrical activity of the brain could be registered by applying electrodes on the scalp and amplifying the resulting signal \cite{DelPreteThesis}. 

The physiological origin of the EEG signal primarily lies in the postsynaptic potentials (PSPs) generated by pyramidal cells contained within the folded sheet of the cerebral cortex \cite{BonomiThesis}. When an excitatory neurotransmitter binds to receptors on the apical dendrite of a pyramidal cell, it causes a positive ionic current to flow into the cell, leaving a net negative charge outside the membrane \cite{BonomiThesis, DelPreteThesis}. This creates a microscopic dipole. Because pyramidal cells are oriented parallel to each other and perpendicular to the cortical surface, the electromagnetic fields originated from thousands of adjacent neurons undergo constructive interference \cite{BonomiThesis}. 

The superimposition of these individual microscopic dipoles from thousands of simultaneously activated neurons forms a single equivalent current dipole \cite{BonomiThesis}. The voltage recorded on the surface of the scalp will be positive on one side of this equivalent dipole and negative on the other, separated by a neutral zero-voltage line. However, brain networks operate across multiple spatial scales. At the microscopic scale, neural activity does not produce an electric field powerful enough to be recorded from the scalp. It is only at the mesoscopic and macroscopic scales—where neural populations span several cubic millimeters to centimeters—that EEG can successfully detect these dynamics \cite{BonomiThesis}.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.7\textwidth]{dipole_model.png}
    \caption{Schematic of pyramidal cells and the formation of a single equivalent current dipole from the cerebral cortex .} % \cite{BonomiThesis}
    \label{fig:dipole}
\end{figure}

\subsection{EEG Recordings and Frequency Bands}
The standard clinical procedure for EEG recording utilizes the International 10-20 system, which places electrodes at proportional distances (10\% or 20\% of the total nasion-inion or left-right ear distance) to standardize measurements \cite{DelPreteThesis}. EEG waveforms measured at the scalp typically exhibit small amplitudes ranging from 10 $\mu$V to 500 $\mu$V and are heavily attenuated by the volume conduction of the skull and brain tissue \cite{BonomiThesis, DelPreteThesis}. 

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.5\textwidth]{10_20_system.png}
    \caption{The International 10-20 system for placing EEG electrodes on the scalp .}% \cite{BonomiThesis}
    \label{fig:10_20_system}
\end{figure}

Beyond amplitude, EEG signals are heavily classified by their frequency spectrum (typically 0.1 Hz to 100 Hz), which directly correlates to different states of consciousness and cognitive behaviors \cite{DelPreteThesis}:
\begin{itemize}
    \item \textbf{Delta ($\delta$) waves (0.5--4 Hz):} High-amplitude waves associated with deep, dreamless sleep.
    \item \textbf{Theta ($\theta$) waves (4--8 Hz):} Associated with light sleep, deep relaxation, and meditation.
    \item \textbf{Alpha ($\alpha$) waves (8--13 Hz):} Present during physical and mental relaxation with eyes closed.
    \item \textbf{Mu ($\mu$) waves (7.5--12.5 Hz):} Originating in the motor cortex, these waves are physically suppressed when a voluntary movement is executed or intended, making them a crucial marker for Brain-Computer Interfaces.
    \item \textbf{Beta ($\beta$) waves (13--30 Hz):} Associated with conscious actions, active thinking, problem-solving, and concentration.
    \item \textbf{Gamma ($\gamma$) waves (30--100 Hz):} Related to perception, hyper-vigilance, and full awareness \cite{DelPreteThesis}.
\end{itemize}

\subsection{Event-Related Potentials and the Readiness Potential}
Because raw EEG is highly noisy and consists of a mixture of neural sources, sensory reactions, and physiological artifacts, specific neurocognitive responses cannot be easily observed in a single continuous trace. To extract meaningful information, the Event-Related Potential (ERP) technique relies heavily on averaging multiple temporal epochs strictly aligned to a specific trigger or event (the "onset") \cite{BonomiThesis, DelPreteThesis}.

In the context of motor cognition, the ERPs occurring in relation to movements are known as Motor Related Cortical Potentials (MRCPs). The most important MRCP for this thesis is the \textbf{Readiness Potential (RP)}, or \textit{Bereitschaftspotential} (BP), which acts as a reliable biomarker of the intentionality behind motion. The RP features a slow negative cortical trend that precedes the actual voluntary muscle contraction by up to 2 seconds \cite{Luck2014, DelPreteThesis}.
 It can be sub-divided into two main components \cite{BonomiThesis}:
\begin{enumerate}
    \item \textbf{Early BP:} A bilaterally symmetrical negative slope originating in the pre-supplementary motor area (pre-SMA). It reflects the subconscious preparatory state and general attention required before a movement.
    \item \textbf{Late BP:} A steeper negative slope originating from the primary motor cortex that reflects the conscious intention to act. Its amplitude increases with the complexity and precision required by the specific movement \cite{BonomiThesis}.
\end{enumerate}


\subsection{The Volume Conduction Problem}
Even if all artifacts are perfectly removed, BCI systems must confront \textbf{volume conduction}. When a postsynaptic potential (PSP) occurs in the cerebral cortex, it generates a microscopic dipole, which in turn creates an instantaneous electric field and a flow of charged particles (ions and electrons) \cite{BonomiThesis}. 

However, electricity does not simply travel in a straight line between the two poles of the neural dipole. The human head consists of multiple layers of varying electrical resistance. While the brain tissue and cerebrospinal fluid (CSF) are highly conductive, the skull represents a highly resistive barrier. As the electrical current encounters the high resistance of the skull, the signal is forced to disperse laterally to find less resistive paths to the surface \cite{BonomiThesis}. 

Consequently, the electrical signal originating from a localized neuronal population spreads across the entire scalp. The scalp essentially becomes divided into three broad macro-regions: a positively charged area, a negatively charged area, and an infinitesimally narrow neutral band separating them \cite{BonomiThesis}. This creates several severe analytical challenges:
\begin{enumerate}
    \item \textbf{Loss of Spatial Resolution:} Each superficial electrode does not merely record the neural activity directly beneath it. Instead, it measures a complex, weighted, and attenuated sum of components originating from neuronal generators spread all over the brain. The distribution of an ERP component becomes extremely broad, making it difficult to pinpoint the active motor cortex regions \cite{BonomiThesis}.
    \item \textbf{The Inverse Problem:} The phenomenon of volume conductance renders the mathematical mapping of scalp topographies back to their original neural sources "ill-posed" or "undetermined." Because the signal disperses so widely, multiple internal dipole configurations can produce the exact same voltage topography on the scalp. For example, a single highly-active deep dipole can yield a topography identical to two moderately-active superficial dipoles \cite{BonomiThesis}. 
    \item \textbf{Destructive Interference:} If multiple underlying dipoles possess opposing orientations, their electrical fields will destructively interfere as they propagate through the skull. This can result in significant motor intention activity completely canceling out before it ever reaches the EEG cap \cite{BonomiThesis}.
\end{enumerate}

Ultimately, overcoming volume conductance is the primary motivation for introducing spatial filters (such as the Laplacian and Spherical Spline Interpolation filters) into the BCI pre-processing chain. By estimating the second spatial derivative of the voltage distribution, these filters counteract the spatial blurring of the skull, isolating the localized brain components required to accurately classify the Lateralized Readiness Potential (LRP).

\section{Standard Preprocessing Pipelines}

Because raw electroencephalographic signals are inherently corrupted by physiological artifacts and severely dispersed by the physical limitations of the skull, a rigorous, multi-stage pre-processing pipeline is absolutely required before any meaningful Event-Related Potentials (ERPs) can be extracted. Standard pipelines, such as those implemented in the EEGLAB and BVibes environments, systematically employ Blind Source Separation (BSS), detrending algorithms, and Spatial Filters to clean the signal sequentially \cite{DelPreteThesis, BonomiThesis}.

\subsection{Independent Component Analysis (ICA)}
To isolate and remove physiological artifacts—such as eye blinks and muscular twitches—without discarding entire epochs of valuable data, Independent Component Analysis (ICA) is extensively utilized. ICA is a specific class of Blind Source Separation (BSS) that exploits the temporal statistical independence of the recorded mixed signals to separate them into their original, underlying sources \cite{DelPreteThesis}. This approach is fundamentally analogous to the ``cocktail party problem,'' where a system must isolate a single voice in a noisy room using multiple spatially distributed microphones.

The mathematical foundation of ICA relies on finding an unmixing matrix $W$. Assuming the signal recorded at the scalp is $X$, the ICA algorithm decomposes this signal into a matrix of Independent Component (IC) activations $U$:
\begin{equation}
    U = WX
\end{equation}
Conversely, the inverse of the unmixing matrix, $W^{-1}$, acts as the mixing matrix, allowing the reconstruction of the scalp data through multiplication with the IC activations ($X = W^{-1}U$) \cite{DelPreteThesis}. 

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{ica_flowchart.png}
    \caption{Schematic flowchart of the ICA decomposition and back-projection process .} % \cite{DelPreteThesis}
    \label{fig:ica_flowchart}
\end{figure}

For EEG signals, ICA successfully operates under several physiological assumptions: (1) the dynamics can be modeled as a collection of statistically independent brain processes; (2) the volume conduction is instantaneous with negligible propagation delay; and (3) the probability distribution of the underlying analog sources is non-Gaussian, which is accurate since active neuronal populations generate sparsely active, non-normal distributions \cite{DelPreteThesis}. Pre-processing pipelines often utilize the `runica` algorithm (developed by Bell and Sejnowski), which relies on a neural network employing a natural gradient to maximize the joint entropy andminimize the mutual information among output units \cite{Bell1995}.

\subsection{Automated Artifact Detection Criteria}
Once the signal is decomposed, automated algorithms evaluate specific features of each Independent Component in both the temporal and frequency domains to identify whether it is of cerebral or artifactual origin \cite{DelPreteThesis}:
\begin{itemize}
    \item \textbf{Correlation with EOG:} If Electrooculography (EOG) reference channels are available, Pearson's correlation coefficient $\rho$ is calculated. Components showing a maximum linear dependence with the EOG are immediately flagged as ocular artifacts.
    \item \textbf{Kurtosis and Skewness:} Kurtosis measures the sharpness of a probability distribution's peak. A distribution with a high kurtosis (leptokurtic) indicates the presence of sudden, sharp spikes. Components heavily corrupted by eye blinking, cardiac rhythms, and intense muscular twitches typically exhibit a leptokurtic distribution. Skewness measures the symmetry of the distribution. 
    \item \textbf{Entropy:} Evaluated as $H = - \sum p(x) \log_2(p(x))$, entropy acts as a statistical measure of randomness. Unpredictable or erratic artifacts present large entropies, whereas genuine brain activity typically displays more structured probability values.
    \item \textbf{Average Log Band Power and Fit-Error:} Muscular components often present high spectral energy in the 25 Hz to 40 Hz band. A fit-error evaluates the deviation of a component's spectrum from a prototypical curve; low fit-errors in high-frequency bands identify EMG noise.
    \item \textbf{Focal Components:} Evaluated via the maximum z-score of the mixing matrix $W^{-1}$, this highlights components that are exclusively confined to a single faulty electrode (e.g., electrode pop artifacts). 
\end{itemize}
Once these artifactual components are flagged, they are set to zero. The EEG signal is then back-projected ($X = W^{-1}U$), effectively removing the artifact's influence without sacrificing the underlying neural data.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{robust_detrending.png}
    \caption{Visual comparison of standard detrending versus robust trial-masked detrending in preserving the genuine EEG baseline .} % \cite{DelPreteThesis}
    \label{fig:detrending}
\end{figure}

\subsection{Spatial Filtering for Volume Conductance}
After temporal noise is removed via ICA and detrending, the pipeline must address the spatial blurring caused by volume conductance. Spatial filters, often referred to as surface Laplacian filters, estimate the second spatial derivative of the instantaneous voltage distribution. By doing so, they heavily emphasize localized radial neuronal sources while attenuating broad, dispersed signals \cite{BonomiThesis}. 

There are several families of spatial filters employed in BCI applications:
\begin{enumerate}
    \item \textbf{Small Laplacian Spatial Filter (SLSF):} This filter acts as a high-pass spatial filter by subtracting a weighted average of the immediate neighboring electrodes (typically 4 or 8 neighbors) from a reference electrode. The equation is defined as $v_i = c_i - \sum_{j=1}^{N} \alpha \cdot c_j$, where $\alpha = 1/N$. This effectively removes DC spatial noise, isolating localized high-spatial-frequency activity \cite{BonomiThesis}.
    
    \item \textbf{Large Laplacian Spatial Filter (LLSF):} Utilizing the next-closest set of neighboring electrodes, the LLSF operates identically to the SLSF but is tuned to preserve middle-high spatial frequencies, making its focal point slightly wider \cite{BonomiThesis}.
    
    \item \textbf{Common Average Reference Spatial Filter (CARSF):} In this variation, the voltage of all other electrodes in the entire montage ($N=33$ for a 34-electrode cap) is averaged and subtracted from the reference. While simple, it retains a mix of low and middle spatial frequencies, often resulting in blurred topographies that fail to precisely isolate the pre-motor cortex \cite{BonomiThesis}.
    
    \item \textbf{Spherical Spline Interpolation Spatial Filter (SISF):} The most mathematically rigorous approach, the SISF projects the 3D scalp surface onto a computational unit-sphere \cite{Perrin1989, BonomiThesis}. To correct 2D planar distortion, the system estimates a cosine distance matrix among all pairs of electrodes. Filter matrices $G$ and $H$ are computed using Legendre polynomials ($P_n$) based on these cosine distances:
    \begin{equation}
        G_{ij} = \frac{1}{4\pi} \sum_{n=1}^{\text{leg.order}} \frac{(2n+1) \cdot P_n(\cos(d_{ij}))}{(n(n+1))^m}
    \end{equation}
    The filter's behavior is explicitly tuned by the integer $m$ (which shifts behavior from a low-pass to a high-pass spatial filter) and the Legendre order $n$ (which defines the spatial spherical harmonics). When optimally tuned (e.g., $m=2, n=10$), the SISF generates highly detailed topographies, vastly improving the Signal-to-Noise Ratio (SNR) and allowing for strict observation of the Lateralized Readiness Potential (LRP) in paralyzed patients \cite{BonomiThesis}.
\end{enumerate}

\begin{table}[htbp]
    \centering
    \renewcommand{\arraystretch}{1.3}
    \caption[Comparison of Spatial Filtering Techniques]{Comparison of Spatial Filtering Techniques utilized in BCI preprocessing to counteract volume conduction \cite{Perrin1989, BonomiThesis}.}

    \vspace{0.2cm}
    \resizebox{\textwidth}{!}{
    \begin{tabular}{|p{3.5cm}|p{4.5cm}|p{5.5cm}|}
        \hline
        \textbf{Spatial Filter} & \textbf{Mathematical Mechanism} & \textbf{Performance on Readiness Potentials} \\
        \hline
        \hline
        \textbf{CARSF} (Common Average Reference) & Subtracts the average voltage of all 33 other montage electrodes. & Blurs the topography; retains mixed low/middle frequencies; poor LRP isolation. \\
        \hline
        \textbf{SLSF} (Small Laplacian) & Subtracts weighted average of the 4 to 8 immediate nearest neighbors. & Highly local; removes DC spatial noise; excellent for strict focal pre-motor isolation. \\
        \hline
        \textbf{LLSF} (Large Laplacian) & Subtracts weighted average of the next-outer ring of neighbors. & Preserves middle-high spatial frequencies; slightly wider focal region than SLSF. \\
        \hline
        \textbf{SISF} (Spherical Spline Interpolation) & Estimates 2nd spatial derivative via Legendre polynomials on a unit sphere. & Most rigorous; dramatically improves SNR; provides the cleanest topographies if temporal noise is removed first. \\
        \hline
    \end{tabular}}
    \label{tab:spatial_filters}
\end{table}

Ultimately, these standard preprocessing steps act as the foundational layer upon which advanced Time-Varying Autoregressive (TVAR) modeling and Kalman smoothing algorithms can be successfully implemented.


\subsection{The Challenge of Pipeline Unification}
While the individual algorithms described above—FIR filtering, Infomax ICA, robust detrending, Spatial Filtering (SISF), and Woody's Jitter Compensation—are widely accepted as standard BCI practices, simply stacking these mathematical processes sequentially does not guarantee optimal signal extraction. In fact, advanced state-space estimators like the Kalman filter rely heavily on specific statistical properties of the data, such as phase continuity and non-Gaussian variance.

A major challenge addressed in this thesis is that applying these standard pre-processing techniques blindly can actively destroy the mathematical foundations required by the Kalman filter, leading to signal divergence. Therefore, the following chapters are dedicated to deeply analyzing the structural interactions between these algorithms to empirically define a unified, conflict-free Time-Varying Autoregressive Kalman architecture.

% ==========================================
% CHAPTER 2
% ==========================================
\chapter{The Kalman Filter: Theory and State-Space Modeling}

\section{State-Space Representations}

Within the significant toolbox of mathematical tools utilized for stochastic estimation from noisy sensor measurements, the Kalman filter is one of the most prominent and widely adopted \cite{Laaraiedh2009}. Published by Rudolf E. Kálmán in 1960, the algorithm estimates unknown variables from a series of measurements observed over time, producing estimates that are significantly more accurate than those based on a single measurement \cite{WikiKF}. 

In the context of biomedical time series, signals such as the electroencephalogram (EEG) are inherently non-stationary stochastic processes driven by complex, underlying neural dynamics. A state-space model provides an elegant mathematical framework to capture these dynamics by emphasizing the existence of ``hidden'' or unobservable states that directly influence the measured variables \cite{Ming2012, Vakitbilir2025}. Fundamentally, state-space models operate through a two-fold process: the state equation, which describes how the system evolves over time, and the observation equation, which details how the unobservable states contribute to the observed signals recorded at the scalp \cite{Vakitbilir2025}. The Kalman filter can be regarded as analogous to a hidden Markov model (HMM), with the key difference that these hidden state variables exist in a continuous space rather than a discrete state space \cite{WikiKF}.

\subsection{The Linear Stochastic Difference Equations}
The Kalman filter is built upon linear dynamic systems discretized in the time domain. It assumes that the true state of the system at time step $k$ evolves from the previous state at $k-1$ according to the following linear stochastic difference equation, known as the \textbf{State Equation} or transition model \cite{WikiKF, Laaraiedh2009}:
\begin{equation}
    x_k = F_k x_{k-1} + B_k u_k + w_k
\end{equation}
where:
\begin{itemize}
    \item $x_k \in \mathbb{R}^n$ is the true, hidden state vector of the system at time $k$.
    \item $F_k$ is the $n \times n$ state transition matrix that maps the previous state $x_{k-1}$ to the current state, assuming no driving function or process noise.
    \item $u_k \in \mathbb{R}^l$ is the optional control-input vector.
    \item $B_k$ is the $n \times l$ control-input matrix that maps the control vector $u_k$ to the state space.
    \item $w_k$ represents the process noise, which accounts for unmodeled dynamics or inherent physiological system uncertainty.
\end{itemize}

Simultaneously, the filter relies on the \textbf{Observation Equation} (or measurement model), which relates the hidden state to the actual noisy measurements collected by the sensors (i.e., the EEG electrodes) \cite{Laaraiedh2009, WikiKF}:
\begin{equation}
    z_k = H_k x_k + v_k
\end{equation}
where:
\begin{itemize}
    \item $z_k \in \mathbb{R}^m$ is the observed measurement vector at time $k$.
    \item $H_k$ is the $m \times n$ observation matrix that maps the true state space into the observed measurement space.
    \item $v_k$ represents the observation noise, which includes sensor inaccuracies, volume conduction blurring, and environmental artifacts.
\end{itemize}

\subsection{Statistical Assumptions and Covariance Matrices}
A cornerstone of the standard Kalman filter's optimality lies in its assumptions regarding the noise parameters. Both the process noise $w_k$ and the measurement noise $v_k$ are assumed to be mutually independent, white, zero-mean multivariate Gaussian processes \cite{Laaraiedh2009, WikiKF}. Mathematically, their probability distributions are defined as:
\begin{equation}
    w_k \sim \mathcal{N}(0, Q_k)
\end{equation}
\begin{equation}
    v_k \sim \mathcal{N}(0, R_k)
\end{equation}
where $Q_k$ is the process noise covariance matrix and $R_k$ is the measurement noise covariance matrix \cite{WikiKF, Laaraiedh2009}. 

In practical implementations, $Q_k$ and $R_k$ might change with each time step; however, they are often assumed to be constant matrices ($Q$ and $R$) to ease the computational load \cite{Laaraiedh2009}. Optimizing these covariance hyperparameters is notoriously difficult but highly necessary for biomedical applications. If $Q_k$ is poorly estimated, the filter might diverge or fail to track rapid transient changes in brain states; conversely, a poorly estimated $R_k$ might cause the filter to overly trust noisy sensor data, failing to adequately smooth the artifactual physiological spikes typical of EEG signals \cite{Khan2006, WikiKF}.

When applied specifically to EEG signal analysis, this generic state-space representation is frequently adapted. Rather than representing physical kinematics (such as position or velocity), the hidden state vector $x_k$ often encapsulates the time-varying autoregressive (TVAR) coefficients \cite{Khan2006}. In this configuration, the observation matrix $H_k$ is replaced by a regressor vector containing past lags of the EEG time-series, allowing the Kalman filter to optimally track and estimate the dynamically changing spectral properties of brain networks \cite{Khan2006, Vakitbilir2025}.

\subsection{The Requirement of Temporal Phase Continuity}
A fundamental, yet often overlooked, requirement of the discrete-time state-space representation is the strict preservation of temporal phase continuity. The state transition equation ($x_k = A x_{k-1} + w_{k-1}$) and the continuous updating of the error covariance matrix ($P_k$) assume that the chronological step from $k-1$ to $k$ represents a genuine physiological progression. 

In biomedical signal processing, if the temporal alignment of the data is artificially manipulated—for instance, by shifting single-trial epochs in time to align specific peaks—the underlying phase continuity of the background stochastic noise ($w_{k-1}$ and $v_k$) is violently broken. When the Kalman equations process phase-shifted data, the recursive covariance tracking diverges because the mathematical history of the noise floor no longer matches its actual physiological progression. Therefore, maintaining strict chronological epoch boundaries is theoretically mandatory for stable state-space estimation.


\section{Standard Kalman Filter Derivation}

The Kalman filter is fundamentally a recursive estimator. This means that to compute the estimate for the current state, the algorithm only requires the estimated state from the immediately previous time step and the current measurement, completely eliminating the need to store the entire historical sequence of observations \cite{WikiKF}. 

Functionally, the Kalman filter is implemented as a predictor-corrector algorithm that operates in two distinct, alternating phases: the \textit{Predict} phase (time update) and the \textit{Update} phase (measurement update) \cite{Laaraiedh2009}. 

\subsection{The Predict Phase (Time Update)}
During the predict phase, the filter utilizes the state estimate from the previous time step to produce an \textit{a priori} estimate of the state at the current time step. It is termed \textit{a priori} because it is generated before the new sensor observation $z_k$ has been incorporated \cite{WikiKF}.

The filter predicts the current state vector $\hat{x}_{k|k-1}$ by applying the state transition matrix $F_k$ and any known control inputs $B_k u_k$:
\begin{equation}
    \hat{x}_{k|k-1} = F_k \hat{x}_{k-1|k-1} + B_k u_k
\end{equation}

Simultaneously, the algorithm must project the estimated error covariance matrix $P$ forward in time. This matrix acts as a quantitative measure of the estimated accuracy (or uncertainty) of the state prediction:
\begin{equation}
    P_{k|k-1} = F_k P_{k-1|k-1} F_k^T + Q_k
\end{equation}
Here, $Q_k$ is the process noise covariance matrix. It is added during the prediction step to account for the injection of new uncertainty due to unmodeled dynamics or inherent physiological noise driving the system between time steps \cite{WikiKF, Laaraiedh2009}.

\subsection{The Update Phase (Measurement Update)}
Once the new noisy measurement $z_k$ is observed from the sensors (e.g., the EEG electrodes), the update phase incorporates this information to refine the \textit{a priori} prediction into a highly accurate \textit{a posteriori} state estimate \cite{WikiKF}.

In order to guarantee numerical stability and preserve the positive-definiteness of the state error covariance matrix $P_{k|k}$ during multi-trial offline computations, the standard update equation is substituted with the robust \textit{Joseph Stabilized Form}. This mathematical reformulation prevents accumulation of round-off errors and is expressed as:
\begin{equation}
    P_{k|k} = (I - K_k H_k) P_{k|k-1} (I - K_k H_k)^T + K_k R_k K_k^T
\end{equation}
By incorporating this symmetric structure, the filter remains bound to physically realistic bounds, eliminating the divergence caused by precision limitations in localized continuous brain volatility.

The first step is to calculate the measurement \textbf{innovation} or pre-fit residual, $\tilde{y}_k$. This vector represents the discrepancy between the actual sensor measurement and the measurement the model expected to see:
\begin{equation}
    \tilde{y}_k = z_k - H_k \hat{x}_{k|k-1}
\end{equation}
The covariance of this innovation, which defines the uncertainty of the prediction mapped into the measurement space, is computed as:
\begin{equation}
    S_k = H_k P_{k|k-1} H_k^T + R_k
\end{equation}
where $R_k$ is the measurement noise covariance matrix \cite{WikiKF}.

Next, the algorithm calculates the \textbf{Optimal Kalman Gain}, $K_k$. The Kalman filter is mathematically constructed as a Minimum Mean-Square Error (MMSE) estimator. Consequently, the Kalman gain is rigorously derived by taking the matrix derivative of the trace of the \textit{a posteriori} error covariance and setting it to zero. This results in the weight matrix that strictly minimizes the a posteriori error:
\begin{equation}
    K_k = P_{k|k-1} H_k^T S_k^{-1}
\end{equation}
Intuitively, the Kalman gain acts as a dynamic weighting factor. If the sensor noise covariance $R_k$ is exceptionally high, the gain approaches zero, and the filter heavily trusts the internal model's \textit{a priori} prediction. Conversely, if the predicted error covariance $P_{k|k-1}$ is high, the filter heavily trusts the new measurement $z_k$ \cite{WikiKF}.

The \textit{a posteriori} state estimate is then finalized by adjusting the \textit{a priori} state with the scaled innovation:
\begin{equation}
    \hat{x}_{k|k} = \hat{x}_{k|k-1} + K_k \tilde{y}_k
\end{equation}
Finally, the \textit{a posteriori} estimate covariance matrix must be updated to reflect the increased confidence in the state gained from the new measurement. In its standard (and computationally cheapest) format, this is expressed as:
\begin{equation}
    P_{k|k} = \left( I - K_k H_k \right) P_{k|k-1}
\end{equation}
While this standard form is nearly always utilized in practice, it is only mathematically robust if the optimal gain is used. In applications suffering from numerical instability or severe round-off errors, an alternative formula known as the Joseph form is often preferred to guarantee that the covariance matrix remains symmetric and positive-definite \cite{WikiKF, Laaraiedh2009}.



\section{Variants: EKF, UKF, and RTS Smoothing}

The standard Kalman filter's optimality is strictly bounded by the assumption that the system dynamics are entirely linear and driven by Gaussian noise. However, realistic biomedical signals, particularly electroencephalograms (EEG), are generated by highly complex neural mechanisms that frequently exhibit non-linear and non-stationary behaviors \cite{Ming2012, Vakitbilir2025}. General non-linear and non-Gaussian models admit no precise closed-form filtering solutions. Consequently, to apply the state-space framework to complex biological systems, several advanced variants and approximations of the Kalman filter have been developed \cite{Ming2012, WikiKF}.

\subsection{Extended Kalman Filter (EKF)}
In the Extended Kalman Filter (EKF), the state transition and observation models are no longer constrained to be linear functions of the state. Instead, they can be represented by non-linear, but strictly differentiable, functions \cite{WikiKF}:
\begin{equation}
    x_k = f(x_{k-1}, u_k) + w_k
\end{equation}
\begin{equation}
    z_k = h(x_k) + v_k
\end{equation}
Here, the non-linear function $f$ is used to compute the predicted state from the previous estimate, while $h$ is used to compute the predicted measurement from the predicted state. 

Because these non-linear functions cannot be directly applied to the covariance matrices without distorting the Gaussian noise assumptions, the EKF performs a step-wise linearization. At each time step, the algorithm computes a matrix of partial derivatives—known as the \textbf{Jacobian}—evaluated around the current predicted state \cite{WikiKF}. The EKF essentially linearizes the non-linear system locally using a first-order Taylor series expansion, allowing the standard Kalman equations to be applied to the Jacobians. While historically popular, the EKF suffers from severe limitations: if the system is highly non-linear, the local linearization becomes highly inaccurate, causing the filter to diverge quickly \cite{WikiKF}.

\subsection{Unscented Kalman Filter (UKF)}
To address the severe linear approximation errors of the EKF, the Unscented Kalman Filter (UKF) was developed. When the transition and observation models ($f$ and $h$) are highly non-linear, the UKF accurately captures the posterior mean and covariance to the 3rd order (Taylor series expansion) for any non-linearity, vastly outperforming the EKF \cite{WikiKF}. 

Instead of explicitly calculating complex Jacobians, the UKF utilizes a deterministic sampling technique known as the \textbf{Unscented Transform (UT)}. The UT selects a minimal, carefully chosen set of sample points around the current mean, called \textit{sigma points}. These sigma points perfectly capture the true mean and covariance of the state distribution \cite{WikiKF}. 

These deterministic sigma points are then individually propagated directly through the true non-linear functions $f$ and $h$. The new predicted mean and covariance are then empirically reconstructed from the weighted average of these transformed points. This approach completely removes the requirement to analytically calculate Jacobians, making it highly robust for complex physiological modeling where exact derivatives are impossible to find \cite{WikiKF}.

\subsection{Rauch-Tung-Striebel (RTS) Smoother}
The standard Kalman filter, EKF, and UKF are all strictly causal, recursive filters; they estimate the current state utilizing only past and present measurements. While this is necessary for real-time, online Brain-Computer Interfaces (BCI), many clinical EEG analyses are performed offline on fully recorded datasets (such as extracting Event-Related Potentials). If future measurements are available, they can be used to dramatically improve the accuracy of past state estimates \cite{Khan2006, WikiKF}. 

This non-causal estimation is achieved through fixed-interval smoothing, with the \textbf{Rauch-Tung-Striebel (RTS) Smoother} being the most prominent optimal algorithm \cite{WikiKF}. The RTS smoother operates in a two-pass mechanism:
\begin{enumerate}
    \item \textbf{Forward Pass:} A standard Kalman filter is run forward through the entire time series data. The \textit{a priori} and \textit{a posteriori} state estimates ($\hat{x}_{k|k-1}$, $\hat{x}_{k|k}$) and their respective covariances ($P_{k|k-1}$, $P_{k|k}$) are saved in memory for every time step \cite{WikiKF}.
    \item \textbf{Backward Pass:} Starting at the final time step $N$, the algorithm proceeds backward in time, calculating the optimally smoothed state $\hat{x}_{k|N}$ using the following recursive retrodiction equations \cite{WikiKF}:
    \begin{equation}
        \hat{x}_{k|N} = \hat{x}_{k|k} + C_k (\hat{x}_{k+1|N} - \hat{x}_{k+1|k})
    \end{equation}
    \begin{equation}
        P_{k|N} = P_{k|k} + C_k (P_{k+1|N} - P_{k+1|k}) C_k^T
    \end{equation}
    where the smoothing gain $C_k$ is computed as:
    \begin{equation}
        C_k = P_{k|k} F_{k+1}^T P_{k+1|k}^{-1}
    \end{equation}
\end{enumerate}

In the context of EEG processing, the Kalman smoother is heavily utilized to estimate time-varying autoregressive (TVAR) coefficients. Because the backward pass incorporates all available data in the epoch, it effectively removes the inherent tracking lag associated with standard adaptive algorithms (like Recursive Least Squares or standard Kalman filtering). This makes the RTS smoother an exceptionally powerful tool for precisely isolating and tracking transient neural dynamics, such as Event-Related Desynchronization (ERD) and Synchronization (ERS) within the motor cortex \cite{Khan2006, Tarvainen2001}.

\subsection{Theoretical Justification for the Linear Standard in TVAR Modeling}
While non-linear variants like the Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF) are exceptionally powerful for highly dynamic systems (such as ballistic trajectory tracking or non-linear robotics), their application to Brain-Computer Interfaces must be carefully justified. 

In the context of this thesis, the Kalman filter is utilized specifically to estimate the parameters of a Time-Varying Autoregressive (TVAR) model. In this specific mathematical formulation, the state transition matrix ($A$), which defines how the AR coefficients evolve from time step $k-1$ to $k$, is defined as an Identity matrix ($I$). 

Because an Identity matrix dictates a strictly linear, one-to-one state transition, the system entirely lacks the severe non-linear properties that the Jacobian matrices of the EKF or the deterministic sigma-points of the UKF were designed to solve. Consequently, from a purely theoretical standpoint, applying an EKF or UKF to an Identity-based TVAR model yields no mathematical advantage over the standard Linear Kalman Filter, while introducing massive, unnecessary computational overhead. This theoretical hypothesis will be empirically tested and validated during the pipeline optimization in Chapter 4.


\section{State-of-the-Art in Biomedical Signals}

In the field of electroencephalography (EEG) and advanced Brain-Computer Interfaces (BCIs), the Kalman filter has evolved from a simple noise-reduction tool into a foundational mechanism for tracking complex, unobserved neural dynamics \cite{Vakitbilir2025}. Because the brain is a highly dynamic system, static filtering approaches fail to capture the transient nature of cognitive processes. Consequently, state-space representations are now heavily utilised to predict future physiological states, integrate multi-modal data (such as fusing EEG with electromyography or intracranial pressure signals), and isolate event-related brain activity from intense physiological noise \cite{Vakitbilir2025}.


\begin{table}[htbp]
    \centering
    \renewcommand{\arraystretch}{1.3}
    \caption[Comparison of Multivariate Time-Series Models]{Comparison of Multivariate Time-Series Models for Cerebral Physiologic Signals \cite{Vakitbilir2025}.}

    \vspace{0.2cm}
    \resizebox{\textwidth}{!}{
    \begin{tabular}{|p{2.5cm}|p{6cm}|p{6cm}|}
        \hline
        \textbf{Model} & \textbf{Primary Advantages} & \textbf{Limitations / Challenges} \\
        \hline
        \hline
        \textbf{VAR} & Captures dynamic linear interactions and temporal relationships across multiple EEG channels. & Requires strictly stationary data; highly sensitive to the chosen lag length ($p$). \\
        \hline
        \textbf{TVAR} & Captures transient events (like ERD/ERS) by allowing parameters to vary over time. & Risk of overfitting; computationally intensive; requires advanced parameter estimation (e.g., Kalman). \\
        \hline
        \textbf{VARFI / VARFIMA} & Handles fractional differencing to accurately model long-range memory and dependencies. & High computational cost to estimate fractional parameters; interpretability issues. \\
        \hline
        \textbf{Kalman Filter} & Optimal recursive real-time state estimation; highly adaptable to non-stationary environments. & Assumes linearity and Gaussian noise; heavily dependent on initial covariance hyperparameter tuning. \\
        \hline
    \end{tabular}}
    \label{tab:multivariate_models}
\end{table}

\subsection{Artifact Suppression and Real-Time Denoising}
A primary application of the Kalman filter in biomedical signal processing is the continuous suppression of physiological artifacts. Standard algorithms often struggle to separate high-amplitude Ocular Artifacts (OA) from genuine frontal cortex activity. The Kalman filter addresses this by recursively estimating the true state of the brain signal while treating the artifact as a component of the measurement noise $v_k$ \cite{Ferdous2017, Vakitbilir2025}. While some comparative offline studies indicate that Wiener filters can achieve a lower Mean Squared Error (MSE) when the global statistical properties of the entire signal are known a priori \cite{Ferdous2017}, the Kalman filter remains vastly superior for online BCI systems. Its recursive predictor-corrector nature allows it to seamlessly adapt to sudden variations without requiring the entire data history in memory \cite{Ferdous2017, Vakitbilir2025}.
Furthermore, recent comparative studies have evaluated the Kalman filter against other deterministic approaches such as the Discrete Wavelet Transform (DWT) and Empirical Mode Decomposition (EMD) for EOG artifact rejection \cite{Salis2013}. While Wavelet transforms (e.g., using a db6 mother wavelet) offer excellent time-frequency localization, they rely on fixed basis functions that can inadvertently remove genuine cerebral activity if the neural signal heavily overlaps with the artifact's frequency band. The Kalman filter, conversely, dynamically adapts to the non-stationary covariance of the system, achieving comparable Signal-to-Noise Ratio (SNR) improvements to Wavelets without suffering from the severe mode-mixing issues inherent in EMD \cite{Salis2013, Enshaeifar2016}.

\subsection{Tracking Spectral Dynamics: ERD/ERS and TVAR Models}
Beyond simple denoising, the state-of-the-art heavily relies on the Kalman filter to track the time-varying spectral properties of the brain. When a subject performs or imagines a motor task, the brain exhibits Event-Related Desynchronization (ERD) in the Mu and central-Beta rhythms, followed by Event-Related Synchronization (ERS) \cite{Khan2006}. To quantify this, the EEG signal is modelled as a Time-Varying Autoregressive (TVAR) process, where the AR coefficients act as the hidden state variables estimated by the Kalman filter \cite{Tarvainen2001, Kostoglou2024}. 

Historically, TVAR coefficients were tracked using adaptive algorithms such as Least Mean Squares (LMS) or Recursive Least Squares (RLS) \cite{Tarvainen2001}. However, these algorithms suffer from a severe "tracking lag"—a temporal delay between the actual cognitive shift and the algorithm's detection of it \cite{Tarvainen2001, Khan2006}. If the forgetting factor $\lambda$ in an RLS algorithm is set too high, the algorithm adapts too slowly; if set too low, the estimates become violently unstable \cite{Tarvainen2001}. By employing a Kalman smoother—which utilises both past and future measurements in a two-pass mechanism—this tracking lag is completely eliminated, yielding perfectly aligned and highly stable estimates of the Readiness Potential and ERD transitions \cite{Tarvainen2001, Khan2006}.

\subsection{Hyperparameter Optimisation Strategies}
The most significant hurdle in applying the Kalman filter to biomedical signals is accurately setting the model hyperparameters: the process noise covariance ($Q$), the measurement noise covariance ($R$), and the initial error covariance ($P_0$) \cite{Khan2006, Kostoglou2024}. If $Q$ is improperly tuned, the filter will either diverge or fail to track rapid shifts in motor intention.

To resolve this, modern BCI frameworks decouple the Kalman filtering from manual tuning by employing fully automatic, data-driven optimisation techniques:
\begin{itemize}
    \item \textbf{Expectation-Maximization (EM):} The EM algorithm is an iterative statistical technique used to find the maximum likelihood estimates of the state-space parameters. By processing a training dataset of EEG trials, the EM algorithm automatically finds the optimal covariance matrices without requiring tedious manual experimentation, vastly improving the smoothness of the ERD estimates \cite{Khan2006}.
    \item \textbf{Genetic Algorithms (GA):} Recent BCI optimisations employ global optimisation techniques like Genetic Algorithms to minimize the Normalized Mean Squared Error (NMSE) across multiple EEG channels simultaneously. Furthermore, by tracking the temporal evolution of the trace of the Kalman covariance matrix $P(n)$, researchers can extract an additional, highly discriminative feature for machine learning classifiers, as sharp changes in the trace directly correlate to transitions in mental states \cite{Kostoglou2024}.
\end{itemize}

\subsection{Non-Linear and Non-Gaussian Advancements}
Despite the mathematical elegance of the Kalman smoother, the foundational algorithm assumes that both the state transitions and the driving noise are strictly linear and Gaussian. However, biological signals, particularly single-trial Event-Related Potentials (ERPs), often exhibit impulsive noise (non-Gaussian) and complex, non-linear neural dynamics \cite{Ming2012}. 

Because general non-linear, non-Gaussian models admit no precise closed-form filtering solutions, recent advancements have shifted towards Sequential Monte Carlo (SMC) methods \cite{Ming2012}. State-of-the-art frameworks now apply techniques such as the Rao-Blackwellized Particle Filter (RBPF) to partially observed diffusion models \cite{Ming2012}. By modelling the EEG dynamics with non-Gaussian stochastic volatility—such as the Cox-Ingersoll-Ross process—these advanced particle filters can capture abrupt, transient changes that Gaussian models inappropriately blur, demonstrating remarkable Signal-to-Noise Ratio (SNR) improvements (e.g., up to 12.46 dB in simulated evoked potentials) and successfully isolating single-trial ERP dynamics from severe background artifacts \cite{Ming2012}.


% ==========================================
% CHAPTER 3
% ==========================================
\chapter{Methodological Framework: Integrating Kalman Filters with EEG}

\section{Introduction to EEG Signal Processing}
Monitoring the spectral characteristics of brain signals provides critical insights into the underlying processes responsible for their generation. In Brain-Computer Interface (BCI) applications, extracting highly accurate temporal features is fundamental in decoding neural activity, as it provides a reliable means to differentiate between various motor tasks or mental states \cite{Kostoglou2024}. 

However, raw electroencephalographic (EEG) data is inherently problematic to process. It is characterized by massive non-stationarity, meaning that the statistical properties of the signal—such as its mean, variance, and frequency composition—change rapidly over time \cite{Cohen2014}. Traditional frequency-domain analyses, such as the Fast Fourier Transform (FFT), assume that the signal is stationary over the evaluated time window. Applying standard Fourier analyses to highly non-stationary Event-Related Potentials (ERPs), such as the Readiness Potential (RP), inevitably results in a severe loss of temporal precision, blurring the exact onset of cognitive motor intentions \cite{Cohen2014, Mohseni2007}. To overcome these limitations, advanced parametric modeling combined with dynamic state-space estimation must be utilized.

\section{Time-Varying Autoregressive (TVAR) Modeling}
To track the dynamic, non-stationary nature of the brain, parametric spectral estimation methods are preferred over non-parametric Fourier methods. Specifically, Autoregressive (AR) modeling provides a highly efficient mathematical framework for describing EEG time series. 

In a standard Autoregressive model, the current value of the EEG signal, $y(t)$, is predicted as a linear combination of its previous values, plus a white noise error term. The standard AR model of order $p$ is defined as:
\begin{equation}
y(t) = \sum_{i=1}^{p} a_i y(t-i) + e(t)
\end{equation}
where $a_i$ represents the autoregressive coefficients, $p$ represents the model order (the number of past data points used to predict the current point), and $e(t)$ is the stochastic white noise driving the system.

While the standard AR model provides excellent frequency resolution, it assumes that the coefficients $a_i$ remain constant. Because the brain's baseline activity changes during the execution of a motor task (such as the Event-Related Desynchronization observed during movement), a standard AR model will fail to capture these transient cognitive shifts \cite{Khan2006}. 

To account for this, the model must be upgraded to a Time-Varying Autoregressive (TVAR) framework. In a TVAR model, the coefficients are not static; they are allowed to evolve continuously over time to perfectly track the shifting neurophysiological dynamics:
\begin{equation}
y(t) = \sum_{i=1}^{p} a_i(t) y(t-i) + e(t)
\end{equation}
The critical challenge in BCI signal processing then becomes: how do we accurately and recursively estimate these time-varying coefficients $a_i(t)$ in the presence of massive physiological noise? The mathematically optimal solution to this dynamic tracking problem is the Kalman filter.

\section{The Kalman Filter as a Dynamic State-Space Estimator}
The Kalman filter is an optimal recursive algorithm designed to estimate the internal state of a linear dynamic system from a series of noisy measurements. It is particularly powerful because it minimizes the mean square error of the estimated parameters, making it the ideal engine for tracking the TVAR coefficients of an EEG signal \cite{Tarvainen2001, Salis2013}.

To apply the Kalman filter to our TVAR model, the system must first be formalized into a discrete-time state-space representation. This representation consists of two primary equations: the State Transition Equation and the Measurement Equation.

\subsection{State-Space Formulation}
Let the state vector $\mathbf{x}_k$ represent the true, hidden AR coefficients at time step $k$. Thus, $\mathbf{x}_k = [a_1(k), a_2(k), \dots, a_p(k)]^T$. The evolution of these cognitive states over time is modeled as a random walk, defined by the \textbf{State Transition Equation}:
\begin{equation}
\mathbf{x}_k = \mathbf{A} \mathbf{x}_{k-1} + \mathbf{w}_{k-1}
\end{equation}
where $\mathbf{A}$ is the state transition matrix that maps the previous state to the current state, and $\mathbf{w}_{k-1}$ represents the process noise, assumed to be zero-mean multivariate Gaussian noise with covariance matrix $\mathbf{Q}$. In the specific case of modeling EEG TVAR coefficients, the state transition matrix $\mathbf{A}$ is defined as an Identity matrix ($\mathbf{I}$), assuming the coefficients evolve smoothly via the stochastic driver $\mathbf{w}$.

The relationship between these hidden cognitive states and the actual electrical voltages measured by the EEG electrodes on the scalp is defined by the \textbf{Measurement Equation}:
\begin{equation}
z_k = \mathbf{H}_k \mathbf{x}_k + v_k
\end{equation}
where $z_k$ is the observed EEG voltage at time $k$, $\mathbf{H}_k$ is the observation matrix (which contains the past $p$ values of the EEG signal: $[y(k-1), y(k-2), \dots, y(k-p)]$), and $v_k$ is the measurement noise—representing muscular artifacts, 50 Hz powerline interference, and volume conduction—with a variance of $R$.

\subsection{The Recursive Tracking Algorithm}
The Kalman filter operates in a continuous, two-step recursive loop: the Prediction (Time Update) step and the Correction (Measurement Update) step.

\subsubsection{Step 1: Prediction (A Priori Estimation)}
Before a new EEG sample is recorded, the filter predicts the current cognitive state based solely on the mathematical history of the system. 
\begin{equation}
\hat{\mathbf{x}}_{k}^{-} = \mathbf{A} \hat{\mathbf{x}}_{k-1}^{+}
\end{equation}
\begin{equation}
\mathbf{P}_{k}^{-} = \mathbf{A} \mathbf{P}_{k-1}^{+} \mathbf{A}^T + \mathbf{Q}
\end{equation}
Here, $\hat{\mathbf{x}}_{k}^{-}$ is the \textit{a priori} state estimate, and $\mathbf{P}_{k}^{-}$ is the \textit{a priori} error covariance matrix, which represents the mathematical uncertainty of the prediction.

\subsubsection{Step 2: Correction (A Posteriori Update)}
Once the actual EEG voltage $z_k$ is measured by the amplifier, the filter calculates the difference between what it predicted and what it actually observed. This difference is called the Innovation (or residual). The filter then computes the \textbf{Kalman Gain} ($\mathbf{K}_k$), which dictates how much the filter should trust the new measurement versus its own prediction:
\begin{equation}
\mathbf{K}_k = \mathbf{P}_{k}^{-} \mathbf{H}_k^T (\mathbf{H}_k \mathbf{P}_{k}^{-} \mathbf{H}_k^T + R)^{-1}
\end{equation}
If the measurement noise $R$ is extremely high (e.g., during a jaw clench or eye blink), the Kalman Gain drops near zero, and the filter ignores the noisy data, relying entirely on its smooth internal prediction. Conversely, if $R$ is low, the filter heavily incorporates the new data.

Finally, the state estimate and the error covariance are updated (\textit{a posteriori}) to reflect the assimilated information:
\begin{equation}
\hat{\mathbf{x}}_{k}^{+} = \hat{\mathbf{x}}_{k}^{-} + \mathbf{K}_k (z_k - \mathbf{H}_k \hat{\mathbf{x}}_{k}^{-})
\end{equation}
\begin{equation}
\mathbf{P}_{k}^{+} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}_k) \mathbf{P}_{k}^{-}
\end{equation}
By continuously looping through these equations at 128 Hz or 512 Hz, the KF-TVAR model acts as a highly aggressive, intelligent temporal smoother that dynamically suppresses physiological noise while perfectly preserving the low-frequency morphological curves of the Readiness Potential.


\subsection{Data-Driven Localized Variance Estimation for Hyperparameters}
In the state-space formulation of the TVAR model, the performance of the Kalman filter heavily depends on how we initialize its hyperparameters: the process noise covariance ($Q$), the measurement noise covariance ($R$), and the initial state error covariance ($P_0$). Historically, these values were found through tedious manual trial-and-error, or by utilizing computationally heavy machine-learning solvers like the Expectation-Maximization (EM) algorithm.

To optimize the pipeline for Brain-Computer Interfaces—where computational speed is critical—we bypassed these heavy iterative algorithms. Instead, we implemented a data-driven localized variance estimation technique directly tied to the temporal gradients of the EEG signal.

First, we set the Time-Varying Autoregressive model order to $p=2$. This provides enough spectral resolution for the low-frequency Readiness Potential while keeping the computational load very low. The state vector is thus configured as $x_k = [a_1(k), a_2(k)]^T$.

Rather than guessing the initial uncertainty, we calculate the baseline variance of the localized EEG trial dynamically ($\sigma^2_{local} = \text{var}(y)$). This localized variance drives the initialization of all Kalman matrices:
\begin{equation}
P_0 = \mathbf{I}_{p \times p} \cdot \sigma^2_{local}
\end{equation}

The process noise covariance matrix ($Q$), which tells the filter how rapidly the AR coefficients are allowed to change, is tuned based on these signal temporal gradients. Through empirical testing, we applied a scaling factor of $10^{-2}$ to the localized variance:
\begin{equation}
Q = \mathbf{I}_{p \times p} \cdot (10^{-2} \cdot \sigma^2_{local})
\end{equation}

Finally, the measurement noise covariance is fixed to the standard deviation of the expected artifact noise ($R = \sigma_{noise}^2$). By utilizing this localized, data-driven method, the Kalman filter optimally scales its tracking responsiveness to the exact noise floor of the individual patient's trial. This completely removes the need for slow offline training datasets and guarantees stable, high-speed convergence, making it highly suitable for future real-time applications.


\section{Theoretical Justification for the Linear Standard}
While non-linear variants like the Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF) are exceptionally powerful for highly dynamic systems (such as ballistic trajectory tracking or non-linear robotics), their application to Brain-Computer Interfaces must be carefully justified. 

In the context of this thesis, the Kalman filter is utilized specifically to estimate the parameters of the TVAR model. As defined in Equation 3.3, the state transition matrix ($\mathbf{A}$), which dictates how the AR coefficients evolve from time step $k-1$ to $k$, is strictly defined as an Identity matrix ($\mathbf{I}$). 

Because an Identity matrix dictates a strictly linear, one-to-one state transition, the system entirely lacks the severe non-linear properties that the Jacobian matrices of the EKF or the deterministic sigma-points of the UKF were designed to solve. Consequently, from a purely theoretical and mathematical standpoint, applying an EKF or UKF to an Identity-based TVAR model yields absolutely no mathematical advantage over the standard Linear Kalman Filter. Utilizing non-linear variants in this specific BCI architecture only introduces massive, unnecessary computational overhead without improving the Signal-to-Noise Ratio.


\subsection{Phase 1: Frequency Filtering and Phase Distortion}
The first stage of signal conditioning involves removing frequency bands that are physiologically irrelevant to motor execution. Historically, standard pipelines utilize Infinite Impulse Response (IIR) Butterworth filters. However, IIR filters introduce severe non-linear phase distortion \cite{Cohen2014}. Phase distortion means that different frequency components within the EEG signal are delayed by different amounts of time. For a precisely time-locked signal like the Readiness Potential, non-linear phase delays will artificially smear the peak of the waveform, making it impossible to determine the exact millisecond of cognitive motor onset.

To rectify this, the optimized pipeline replaces IIR filters with zero-phase Hamming-windowed sinc FIR (Finite Impulse Response) filters (specifically utilizing EEGLAB's \texttt{pop\_eegfiltnew} architecture). By applying the filter forward and then backward across the data, the FIR filter perfectly neutralizes any temporal delays, maintaining absolute phase continuity \cite{Cohen2014}. Furthermore, to heavily suppress the 20 Hz to 150 Hz spectral energy associated with myogenic (EMG) muscle twitches, the low-pass boundary is aggressively restricted to 10 Hz, effectively neutralizing powerline interference and high-frequency stochastic noise before the signal reaches the state-space estimator.

\subsection{Phase 2: Independent Component Analysis (ICA) for Artifact Rejection}
Even after aggressive 10 Hz filtering, large-amplitude physiological artifacts—such as Electrooculographic (EOG) eye blinks—remain deeply embedded in the low-frequency spectrum ($0.1 - 5$ Hz) and overlap directly with the Readiness Potential. Because these signals share the same frequency band, temporal filters cannot separate them. 

Instead, the pipeline utilizes Infomax Independent Component Analysis (ICA). ICA is a blind source separation algorithm that linearly unmixes the multidimensional scalp EEG signals into statistically independent spatial components. The mathematical foundation of Infomax ICA relies on maximizing joint entropy and isolating non-Gaussian distributions \cite{Cohen2014}. Ocular blinks and muscular spikes are highly leptokurtic (they possess a sharp, non-Gaussian probability distribution). By evaluating the kurtosis and skewness of the unmixed components, the algorithm identifies which components represent eye blinks, zeros out their activation matrices, and back-projects the cleaned cerebral data to the scalp.

\subsection{Phase 3: The Critical Placement of the KF-TVAR Layer}
A fundamental methodological breakthrough introduced in this thesis is the strict positioning of the Kalman filter \textit{after} single-trial epoching, rather than applying it to the continuous EEG stream. 

If the Kalman state-space smoother is applied to continuous, un-epoched data prior to Phase 2, its mathematical variance-suppression mechanism actively dampens the sharp, leptokurtic spikes of eye blinks. By smoothing these spikes, the Kalman filter inadvertently "Gaussianizes" the artifacts. Because ICA strictly requires non-Gaussian distributions to separate independent sources, applying the Kalman filter first actively destroys the unmixing matrix of the ICA algorithm, resulting in a catastrophic failure of artifact rejection.

By epoching the data first—extracting precise temporal windows (e.g., $-4500$ ms to $+3000$ ms) strictly aligned to the EMG motor onset—the Kalman filter can be applied safely. Operating on independent epochs allows the state-space equations to track the intra-epoch variance without disrupting the global non-Gaussianity required by earlier spatial separation layers.

\subsection{Phase 4: The Redundancy of Spherical Spline Spatial Filtering}
Historically, EEG signals suffering from severe volume conductance required Spatial Filters—such as the Spherical Spline Interpolation Spatial Filter (SISF) or the Surface Laplacian—to calculate the second spatial derivative of the voltage distribution and isolate local neuronal sources \cite{Cohen2014}. Volume conduction occurs because the electrical fields generated by the cortex spread tangentially as they travel through the resistive layers of the skull and scalp, causing significant spatial blurring across neighboring electrodes.

However, a critical methodological discovery of this framework is the redundancy of this spatial layer when paired with advanced temporal tracking. Because the combination of the aggressive 10 Hz zero-phase FIR filter and the post-epoching KF-TVAR layer efficiently eradicates the vast majority of background and volume-conducted stochastic noise in the temporal domain, the SISF has minimal residual spatial noise to suppress. Consequently, complex spatial filtering can be bypassed in the optimized pipeline to preserve computational efficiency without sacrificing topographical clarity.

\subsection{Phase 5: The Antagonistic Role of Jitter Compensation}
Standard BCI pipelines frequently utilize Jitter Compensation techniques (such as Woody's method or RIDE) to temporally align single-trial epochs based on peak latencies prior to averaging. While shifting single-trial epochs forward or backward in time is highly beneficial for simple arithmetic averaging, this thesis establishes that jitter compensation is fundamentally antagonistic to state-space modeling. 

As established in Equation 3.3 and 3.6, the Kalman filter relies on the continuous temporal evolution of the covariance matrix ($\mathbf{P}_k$) and the phase continuity of the background stochastic noise ($\mathbf{w}_{k-1}$). Shifting single-trial epochs in time violently breaks the chronological temporal phase of the micro-oscillations. When the Kalman filter processes these artificially shifted, phase-broken epochs, the mathematical history of the noise floor no longer matches its actual physiological progression, causing the AR tracking to diverge and heavily degrade the signal. Therefore, maintaining strict chronological epoch boundaries—and strictly disabling Jitter Compensation—is theoretically mandatory for stable state-space estimation.

\section{The Finalized Pipeline Architecture}
Based on rigorous mathematical testing and theoretical validation across the benchmark datasets, the standard BVibes sequence was significantly restructured. The finalized sequence of operations, yielding the highest Signal-to-Noise Ratio and the cleanest Readiness Potential morphology, strictly defines the integration of the KF-TVAR model while purposefully omitting conflicting spatial and temporal shifting algorithms. The complete logical flow is formalized in Algorithm \ref{alg:pipeline2}.

\begin{algorithm}[htbp]
\caption{Optimized KF-TVAR Preprocessing Pipeline (Pipeline 2)}
\label{alg:pipeline2}
\begin{algorithmic}
\REQUIRE Raw EEG Signal sampled at $f_s = 512$ Hz
\STATE \textbf{Resample} EEG signal to $128$ Hz to optimize model order limitations.
\STATE \textbf{Temporal Filtering:} Apply \texttt{pop\_eegfiltnew} (Zero-phase FIR filter, $0.05$ Hz - $10$ Hz) to eliminate EMG noise and non-linear phase distortion.
\STATE \textbf{Bad Channel Rejection:} Remove channels with repetitive abnormal variance (preserving essential motor channels C3, C4, Cz).
\STATE \textbf{Artifact Detection \& Correction:} Execute Infomax ICA (\texttt{runica}) utilizing non-Gaussian component isolation.
\STATE \textbf{Channel Interpolation:} Spherical spline interpolation of rejected peripheral channels.
\STATE \textbf{Epoching:} Extract epochs strictly aligned to EMG onset (Jitter Compensation explicitly disabled to preserve temporal phase continuity).
\STATE \textbf{Kalman Filtering (KF-TVAR):} Apply linear state-space smoothing on the independent epoched trials to track dynamic cognitive state transitions.
\STATE \textbf{Bad Epoch Rejection:} Remove trials containing residual high-variance artifacts.
\STATE \textbf{ERP Estimation:} Average trials chronologically to extract the high-SNR Readiness Potential.
\RETURN Cleaned ERP and Topographical features for BCI classification.
\end{algorithmic}
\end{algorithm}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\textwidth]{pipeline1.png}
\par\vspace{0.15cm}
\makebox[0.45\textwidth][c]{a)} \hfill \makebox[0.45\textwidth][c]{b)}
\caption{Initial Pre-processing Pipelines. a) The standard unoptimized baseline pipeline utilizing basic artifact correction and spatial filtering without any state-space Kalman tracking. b) The preliminary architecture incorporating Kalman filtering, but still utilizing phase-breaking jitter compensation.}
\label{fig:pipeline1}
\end{figure}


\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.85\textwidth]{pipeline2.png}
    \caption{Optimized Pipeline 2. The Linear Kalman filtering layer is applied post-epoching, while redundant spatial filters and phase-breaking jitter compensation algorithms are strictly disabled to preserve state-space tracking mechanics.}
    \label{fig:pipeline2}
\end{figure}%

% ==========================================
% CHAPTER 4
% ==========================================
\chapter{Implementation, Experimental Framework, and Results}

\section{Introduction to the Experimental Framework}
Moving from raw, multi-channel EEG recordings to clean Readiness Potentials (RPs) requires a very optimized signal processing architecture. As we discussed in the previous chapters, the fact that brain signals are non-stationary, combined with a lot of physiological and environmental noise, makes traditional averaging techniques quite limited \cite{Mohseni2007}. The main goal of this thesis is to engineer a pipeline that can isolate these small cognitive intentions, even in severely paralyzed patients.

This chapter explains the step-by-step optimization of the pre-processing pipeline. Instead of just assuming that standard techniques work well together, we tested every stage of the signal processing chain—like temporal filtering, blind source separation, state-space Kalman estimation, spatial filtering, and jitter compensation. After defining our mathematical evaluation criteria, we validate our optimized Time-Varying Autoregressive Kalman Filter (KF-TVAR) architecture on multiple real clinical datasets, covering voluntary, and totally involuntary motor tasks.


\section{Mathematical Evaluation Criteria}
To objectively quantify the success of the proposed signal processing techniques across different datasets, qualitative visual inspection must be strictly supported by robust mathematical parameters. The evaluation of Slow Cortical Potentials (SCPs), such as the Readiness Potential, requires specific metrics that capture both the suppression of high-frequency noise and the isolation of the low-frequency motor intention without penalizing the natural morphology of the brainwave.

\subsection{Localized Onset Signal-to-Noise Ratio ($SNR_{onset}$)}
A successfully extracted Readiness Potential (RP) must show a clearly positive Signal-to-Noise Ratio. The most important part of the SNR evaluation is strictly around zero ($t=0$). Why? Because the actual Readiness Potential represents the exact moment of motor execution, which peaks precisely at the onset. Looking at a general, basic SNR calculation over the entire epoch can be highly misleading because it dilutes the real performance of the pipeline by including non-motor background periods.

Instead of using a basic overall power ratio, we evaluate the mathematical mean of the continuous instantaneous SNR array strictly within a 500 ms analytical window around the motor onset. To do this, we simply define a sample window of $\pm 500$ ms based on the sampling frequency of our data. We then locate the exact index of the motor onset ($t=0$) and extract all the continuous SNR values from 500 ms before to 500 ms after the movement. Finally, we calculate the average of just these specific values.

Mathematically, this translates to extracting the Localized Onset SNR ($SNR_{onset}$) as:
\begin{equation}
SNR_{onset}(ch)=\frac{1}{N_{win}}\sum_{k=t_{start}}^{t_{end}}SNR_{continuous}(ch,k)
\end{equation}
where $t_{start}$ and $t_{end}$ represent the sample boundaries of our $\pm1000$ ms window, and $N_{win}$ is the total number of samples in that range. 


We use this specific targeted approach to accurately compare the baseline non-Kalman pipeline with the optimized Kalman version. By strictly isolating this tight window, we guarantee that our metric captures how well the pre-processing pipeline decoupled the true neural transition exactly when the patient moved, without being tricked by unrelated SNR fluctuations occurring seconds before or after the event.




\subsection{Residual Oscillation Index ($OI_{res}$)}
Standard evaluation metrics utilized in biomedical literature, such as the Normalized Root Mean Square Error (NRMSE) or the Percentage Root Mean Square Difference (PRD) \cite{Salis2013}, are generally designed for high-frequency, highly repetitive signals like Electrocardiograms (ECGs). Applying standard derivative boundaries over raw EEG averages systematically penalizes the natural physiological curvature and turning-null points of Slow Cortical Potentials. 

To selectively isolate and evaluate background non-cerebral oscillations without distorting the underlying RP morphology, a trend-decoupling framework is implemented to calculate the Residual Oscillation Index ($OI_{res}$). For each electrode channel $ch$, the macroscopic intentional motor trend is captured via a low-order polynomial estimation matrix:
\begin{equation}
y_{trend}(ch, t) = \sum_{m=0}^{N} \beta_m t^m
\end{equation}
The high-frequency residual noise component, $y_{res}(ch, t)$, is isolated by subtracting this slow behavioral baseline from the continuous signal configuration:
\begin{equation}
y_{res}(ch, t) = RP_{data}(ch, t) - y_{trend}(ch, t)
\end{equation}
The localized residual oscillation index is subsequently computed by taking the absolute second-order discrete temporal derivative strictly across the isolated noise matrix:
\begin{equation}
OI_{res}(ch) = \sum_{t} \left| \frac{\partial^2 y_{res}(ch, t)}{\partial t^2} \right|
\end{equation}
Because this metric is a discrete sum of second-order temporal derivatives used strictly for relative baseline comparison, it is quantified in arbitrary units (a.u.). This mathematically guarantees that a reduction in the $OI_{res}$ reflects the authentic suppression of microvoltal myogenic (EMG) artifacts and Gaussian noise floors.


\section{Structural Optimization of the Pre-Processing Pipeline}
Before applying the architecture to the broader clinical database, the individual components of the signal processing chain were independently optimized to prevent the propagation of mathematical errors through the pipeline.

\subsection{Frequency Filtering and Artifact Rejection}
Initial temporal filtering utilized a standard Infinite Impulse Response (IIR) Butterworth approach with a wide passband of $0.05$ Hz to $50$ Hz. Muscular artifacts (EMG) caused by jaw clenching or neck tension encapsulate massive spectral energy in the 20 Hz to 150 Hz band. The 50 Hz low-pass boundary inadvertently allowed a massive amount of this high-frequency myogenic noise to leak into the final average, causing severe high-frequency oscillations in the Readiness Potential. 

To rectify this, the pipeline was updated to utilize a zero-phase Hamming-windowed sinc FIR filter with a strict cutoff at $10$ Hz. Unlike IIR filters, the zero-phase FIR filter introduces absolutely no non-linear phase distortion, significantly improving the temporal alignment of the final RP peak while successfully eradicating the high-frequency muscular interference.

\begin{figure}[htbp]
    \centering
    \begin{minipage}{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RP_SNR_with_old_filter.png}
        \caption{RP and SNR extraction using the outdated 50Hz IIR filter, showing heavy residual EMG oscillations ($OI_{res}$).}
        \label{fig:old_filter}
    \end{minipage}\hfill
    \begin{minipage}{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RP_SNR_with_new_filter.png}
        \caption{Significant morphological improvement utilizing the optimized zero-phase 10Hz FIR filter.}
        \label{fig:new_filter}
    \end{minipage}
\end{figure}

To permanently remove electrooculographic (EOG) artifacts, the standard pipeline employs Infomax Independent Component Analysis (ICA). In pursuit of a more automated approach, the FORCe algorithm (Fully Online and Automated Artifact Removal for Brain-Computer Interfacing)—which utilizes a hybrid Wavelet-SOBI decomposition \cite{Daly2015}—was tested. 

Despite its theoretical advantages in real-time environments, offline comparative analyses revealed that FORCe did not improve the $OI_{res}$ or the average SNR in highly corrupted clinical data. In fact, standard Infomax ICA proved vastly superior at permanently isolating and eliminating ocular artifacts without blurring the underlying low-frequency motor intent.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.7\textwidth]{FORCe_RP_SNR.png}
    \caption{RP and SNR metrics demonstrating that the FORCe hybrid SOBI algorithm failed to outperform standard Infomax ICA in offline Readiness Potential extraction.}
    \label{fig:force_test}
\end{figure}

\section{Optimizing State-Space Estimation: The Kalman Architecture}
The core contribution of this thesis is the integration of the Time-Varying Autoregressive (TVAR) Kalman filter. Because the Kalman filter acts as a recursive state-space estimator, its mathematical interaction with other processing blocks must be rigorously defined.

\subsection{Optimal Structural Placement of the Filter}
The structural placement of the filter heavily dictates pipeline success. Applying the Kalman filter on continuous data \textit{prior} to artifact correction severely degraded performance. The Kalman smoothing mechanism actively suppresses variance; when applied to continuous raw data, it inadvertently "Gaussianizes" leptokurtic physiological spikes (like eye blinks). Because Independent Component Analysis (ICA) strictly relies on the non-Gaussianity of sources to achieve blind separation, pre-smoothing the data actively destroys ICA's separation matrix.

Conversely, embedding the Kalman filter \textit{after epoching} perfectly isolated the single-trial dynamics. The state-space estimator tracks the intra-epoch variance independently, yielding the highest SNR and the cleanest frequency spectra without disrupting the artifact rejection layers.

\begin{figure}[htbp]
    \centering
    \begin{minipage}{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RP_SNR_Kalman_before_artifact.png}
        \caption{Degraded SNR when Kalman smoothing is erroneously applied to continuous data before ICA.}
        \label{fig:kf_before}
    \end{minipage}\hfill
    \begin{minipage}{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RP_SNR_Kalman_after_artifact.png}
        \caption{Significant SNR improvement when Kalman smoothing is placed securely after single-trial epoching.}
        \label{fig:kf_after}
    \end{minipage}
\end{figure}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{Spectrum_kalman_after_epoching.png}
    \caption{Frequency spectrum stabilization confirming the successful suppression of broad-band noise when the Kalman filter is applied strictly after epoching.}
    \label{fig:spectrum_epoching}
\end{figure}

\subsection{Evaluating Kalman Variants: KF vs. EKF vs. UKF}
Simulations were executed to determine if non-linear filter variants could improve the tracking of the TVAR coefficients. The standard linear Kalman Filter (KF) was compared against the Extended Kalman Filter (EKF)—which utilizes Jacobian matrices to linearize non-linear models—and the Unscented Kalman Filter (UKF)—which uses deterministic sigma-point sampling for highly non-linear distributions. 

The empirical results definitively proved that there are virtually no morphological or statistical differences between the three variants. Because the defined state transition matrix ($A$) mapping the AR coefficients is an Identity matrix, the system lacks the severe non-linear temporal transitions required to benefit from EKF or UKF architectures. Thus, the standard linear KF is mathematically optimal and computationally superior, avoiding unnecessary algorithmic complexity.

\begin{figure}[htbp]
    \centering
    \begin{minipage}{0.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{GF99_with_kalman.png}
        \caption{Standard Linear KF.}
        \label{fig:gf99_kf}
    \end{minipage}\hfill
    \begin{minipage}{0.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{GF99_with_EKF.png}
        \caption{Extended KF (EKF).}
        \label{fig:gf99_ekf}
    \end{minipage}\hfill
    \begin{minipage}{0.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{GF99_with_UKF.png}
        \caption{Unscented KF (UKF).}
        \label{fig:gf99_ukf}
    \end{minipage}
\end{figure}

\subsection{The Antagonistic Impact of Jitter Compensation}
A crucial empirical discovery in this research was the antagonistic relationship between Woody's Jitter Compensation and state-space Kalman tracking. Traditionally, jitter compensation aligns epochs based on peak latency to improve the amplitude of the standard mathematical average. 

However, shifting single-trial epochs in time misaligns the underlying phase consistency of the background stochastic micro-oscillations. The Kalman filter relies on the continuous temporal evolution of the covariance matrix ($P_k$). When the filter processes these artificially shifted, phase-broken epochs, the state-transition equations are disrupted, causing the AR tracking to diverge or excessively smooth the target onset. As shown in Figure \ref{fig:jitter_kalman_fail}, applying Jitter Compensation \textit{with} the Kalman filter severely degrades the signal morphology. Therefore, to preserve true state-space dynamics, Jitter Compensation must be strictly disabled when utilizing the KF-TVAR approach.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.7\textwidth]{GF99_with_jitter_Compensation_with_kalman.png}
    \caption{Degradation of the Readiness Potential when Jitter Compensation is mistakenly combined with Kalman state-space tracking, breaking the temporal phase continuity.}
    \label{fig:jitter_kalman_fail}
\end{figure}

\subsection{The Redundancy of Spatial Filtering (SISF)}
Similarly, comparative results indicated that applying the Spherical Spline Interpolation Spatial Filter (SISF) yielded nearly identical, if not slightly dampened, SNR values compared to the spatially unfiltered data. Because the aggressive 10 Hz FIR filter and the post-epoching Kalman smoother already eradicated the vast majority of the background and volume-conducted temporal noise, the spatial filter had minimal residual noise to suppress. While SISF remains theoretically valid for high-density topographical mapping, it is redundant for localized RP tracking and can be bypassed to reduce computational load.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.7\textwidth]{GF99_with_spetial_filtering_jitter_Compensation_kalman.png}
    \caption{RP extraction including Spatial Filtering (SISF). The spatial filter provides no significant enhancement when the Kalman filter is already aggressively active.}
    \label{fig:spatial_redundancy}
\end{figure}

\section{In-Depth Analysis on the GF89 Benchmark Dataset}
To rigorously quantify the effects of the optimized pipeline, a deep analysis was performed on the GF89 benchmark dataset. Table \ref{tab:gf89_benchmark} details the regional telemetry across the scalp, proving the specific enhancements generated by the Kalman filter on the primary motor pathways compared to the global average.

\begin{table}[htbp]
\centering
\small
\renewcommand{\arraystretch}{1.15}
\caption{Quantitative performance telemetry on the GF89 dataset: Regional analysis tracking the Motor Pathways (FC1, FC2, C3, Cz, C4) versus the Global 34-channel average before and after state-space Kalman integration.}
\vspace{0.2cm}
\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Region} & \textbf{Kalman Layer} & \textbf{Average SNR (dB)} & \textbf{Residual $OI_{res}$} \\ \hline \hline
\textbf{Global (34 Channels)} & Disabled  & -7.71  & 3.78 \\
                              & Enabled   & -5.83  & 0.51 \\ \hline
\textbf{Motor Pathways}       & Disabled  & -2.62  & 4.29 \\
                              & Enabled   & +4.60  & 0.51 \\ \hline
\end{tabular}
\label{tab:gf89_benchmark}
\end{table}

As observed in the regional analysis, the unfiltered baseline for the motor cortex was buried under a negative Signal-to-Noise Ratio (-2.62 dB). By strictly tracking the true TVAR state, the Kalman filter suppresses the stochastic spikes and successfully pulls the focal motor intentionality into a positive threshold (+4.60 dB). Simultaneously, the Residual Oscillation Index ($OI_{res}$) collapses from 4.29 down to 0.51, mathematically verifying the aggressive eradication of high-frequency myogenic artifacts.

\begin{figure}[htbp]
\centering
\begin{minipage}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{GF99_without_kalman.png}
\caption{GF89: Unfiltered Raw RP showing significant negative-dB noise floors.}
\end{minipage}\hfill
\begin{minipage}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{GF99_with_kalman.png}
\caption{GF89: Stabilized with KF-TVAR. Motor SNR elevates into a positive threshold.}
\end{minipage}
\end{figure}

\section{Validation Across Multi-Task Clinical Datasets}
To ensure the pipeline is robust and generalizable, the finalized architecture (128 Hz resampling, 10 Hz FIR filtering, Infomax ICA, post-epoching Linear Kalman smoothing, NO jitter compensation, NO spatial filtering) was executed across additional clinical datasets spanning varying levels of patient consciousness and signal corruption.

\subsection{Standard Motor Tasks: GF89, AL85, MI99, and SS89}
Across all standard motor task datasets, the application of the post-epoching Kalman loop yielded exceptionally clean, non-oscillating Readiness Potentials.


\begin{figure}[htbp]
\centering
\begin{minipage}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{AL85_without_kalman.png}
\caption{AL85: Unfiltered Raw RP containing overlapping, messy motor topographies.}
\end{minipage}\hfill
\begin{minipage}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{AL85_with_kalman.png}
\caption{AL85: Stabilized with KF-TVAR, revealing pristine signal divergence.}
\end{minipage}
\end{figure}

\begin{figure}[htbp]
\centering
\begin{minipage}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{MI99_without_kalman.png}
\caption{MI99: Unfiltered Raw RP.}
\end{minipage}\hfill
\begin{minipage}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{MI99_with_kalman.png}
\caption{MI99: Stabilized with KF-TVAR.}
\end{minipage}
\end{figure}

\begin{figure}[htbp]
\centering
\begin{minipage}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{SS89_without_kalman.png}
\caption{SS89: Unfiltered Raw RP. High-frequency oscillations mask the onset.}
\end{minipage}\hfill
\begin{minipage}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{SS89_with_kalman.png}
\caption{SS89: Stabilized with KF-TVAR. $OI_{res}$ successfully collapses.}
\end{minipage}
\end{figure}

\subsection{Complex Motor Tracking: TC99}
The TC99 dataset captures complex motor reactions. These movements naturally have a much lower physiological amplitude and a higher baseline variance compared to standard executed movements.


\begin{figure}[htbp]
    \centering
    \begin{minipage}{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{TC99_without_kalman.png}
        \caption{TC99: Unfiltered Raw RP, showing the target motor intention buried in a negative SNR floor.}
    \end{minipage}\hfill
    \begin{minipage}{0.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{TC99_with_kalman.png}
        \caption{TC99: Stabilized with KF-TVAR, successfully elevating the signal into a positive threshold.}
    \end{minipage}
\end{figure}


\begin{table}[htbp]
\centering
\small
\renewcommand{\arraystretch}{1.2}
\caption{Comprehensive performance telemetry across multiple clinical datasets comparing the Unfiltered Baseline against the Optimized Linear KF-TVAR Pipeline (No Jitter Compensation, No Spatial Filter).}
\vspace{0.2cm}
\resizebox{0.95\textwidth}{!}{
\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Dataset} & \textbf{Measurement} & \textbf{Baseline (No Kalman)} & \textbf{Optimized Linear KF} \\
\hline
\hline
\textbf{GF89} & Motor Cortex SNR (C3/C4) & -2.6 dB & \textbf{+4.6 dB} \\
              & Pre-Motor $OI_{res}$ (FC1) (a.u.) & 4.29     & \textbf{0.51} \\
\hline
\textbf{AL85} & Motor Cortex SNR (C3/C4) & +8.5 dB  & \textbf{+21.4 dB} \\
              & Pre-Motor $OI_{res}$ (FC2) (a.u.) & 3.10     & \textbf{0.41} \\
\hline
\textbf{MI99} & Motor Cortex SNR (C3/C4) & +12.1 dB & \textbf{+28.5 dB} \\
              & Pre-Motor $OI_{res}$ (FC1) (a.u.) & 4.05     & \textbf{0.55} \\
\hline
\textbf{SS89} & Motor Cortex SNR (C3/C4) & +9.4 dB  & \textbf{+19.8 dB} \\
              & Pre-Motor $OI_{res}$ (FC2) (a.u.) & 3.88     & \textbf{0.49} \\
\hline
\textbf{TC99} & Motor Cortex SNR (C3/C4) & -2.1 dB  & \textbf{+11.3 dB} \\
\hline
\textbf{DS\_uni} & Motor Cortex SNR (C3/C4) & -3.8 dB & \textbf{+1.2 dB} \\
                 & Pre-Motor $OI_{res}$ (FC1) (a.u.) & 5.14 & \textbf{0.82} \\
\hline
\end{tabular}
}
\label{tab:comprehensive_results}
\end{table}

\section{Multi-Dataset Statistical Validation}
To summarize the overall mathematical stabilization achieved by the pipeline, a Grand Average statistical validation was performed across all tested clinical datasets. The telemetry was isolated to the primary motor network (FC1, FC2, C3, CZ, C4) and evaluated using the Wilcoxon Rank-Sum test to determine the statistical significance of the Kalman filter's enhancements. 

The comprehensive channel-by-channel grand averages, alongside their standard deviations (SD) and resulting $p$-values, are detailed in Table \ref{tab:grand_average_statistics}.

\begin{table}[htbp]
\centering
\small
\renewcommand{\arraystretch}{1.3}
\caption{Grand Average Statistical Validation Benchmark across the Motor Network. The table reports the Mean $\pm$ Standard Deviation for SNR and the Mean Residual Oscillation Index. Significance is calculated via Wilcoxon Rank-Sum test ($p < 0.05$).}
\vspace{0.2cm}
\resizebox{\textwidth}{!}{
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Channel} & \textbf{Baseline SNR (dB)} & \textbf{Kalman SNR (dB)} & \textbf{Baseline $OI_{res}$ (a.u.)} & \textbf{Kalman $OI_{res}$ (a.u.)} \\
\hline
\hline
\textbf{FC1} & 0.26 $\pm$ 10.29 & \textbf{3.96} $\pm$ 4.51 & 5.09 & \textbf{1.04} \\
\textbf{FC2} & 1.68 $\pm$ 9.52 & \textbf{2.03} $\pm$ 4.62 & 6.99 & \textbf{1.54} \\
\textbf{C3} & 3.20 $\pm$ 11.66 & \textbf{5.34} $\pm$ 7.00 & 4.88 & \textbf{0.99} \\
\textbf{CZ} & 1.51 $\pm$ 10.95 & \textbf{4.61} $\pm$ 5.53 & 5.57 & \textbf{1.06} \\
\textbf{C4} & 6.92 $\pm$ 12.43 & \textbf{13.84} $\pm$ 7.46 & 6.16 & \textbf{1.23} \\
\hline
\textbf{TOTAL (MEAN)} & 2.71 $\pm$ 10.97 & \textbf{5.95} $\pm$ 5.82 & 5.74 & \textbf{1.17} \\
\hline
\hline
\multicolumn{5}{|c|}{\textbf{Wilcoxon Rank-Sum Significance Test}} \\
\hline
\textbf{$p$-value} & \multicolumn{2}{c|}{$p = 0.15079$ } & \multicolumn{2}{c|}{\textbf{$p = 0.00794$ }} \\
\hline
\end{tabular}
}
\label{tab:grand_average_statistics}
\end{table}

\begin{table}[htbp]
\centering
\small
\renewcommand{\arraystretch}{1.3}
\caption{Quantified Improvements of the Kalman Filter: Enhancements in SNR and reductions in Standard Deviation (variance) and Residual Oscillation Index (myogenic artifacts) across the Motor Network.}
\vspace{0.2cm}
\resizebox{0.95\textwidth}{!}{
\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Channel} & \textbf{SNR Improvement (+dB)} & \textbf{STD Reduction (-dB)} & \textbf{$OI_{res}$ Reduction (-a.u.)} \\
\hline
\hline
\textbf{FC1} & +3.70 dB & -5.78 dB & -4.05 \\
\textbf{FC2} & +0.35 dB & -4.90 dB & -5.45 \\
\textbf{C3}  & +2.14 dB & -4.66 dB & -3.89 \\ % FIXED: Changed from -2.14 dB to +2.14 dB to align with master results
\textbf{CZ} & +3.10 dB & -5.42 dB & -4.51 \\
\textbf{C4} & +6.92 dB & -4.97 dB & -4.93 \\
\hline
\textbf{TOTAL (MEAN)} & \textbf{+3.24 dB} & \textbf{-5.15 dB} & \textbf{-4.57} \\
\hline
\end{tabular}
}
\label{tab:grand_average_improvements}
\end{table}

As visualized in Figure \ref{fig:diff_datasets_stats}, the application of the optimized Kalman filter produces a massive, statistically significant suppression of the Residual Oscillation Index ($p = 0.00794$). The global $OI_{res}$ collapses from a highly erratic baseline of 5.74 down to a tightly bounded 1.17. This mathematically confirms that the state-space tracking loop successfully and consistently eradicates high-frequency myogenic artifacts across all patient profiles. 

Regarding the Signal-to-Noise Ratio, the Grand Average demonstrates a positive mean amplification from +2.71 dB up to +5.95 dB. However, due to the massive inter-subject variance inherent in evaluating healthy individuals alongside severely hemiplegic patients, the raw SNR amplification does not reach strict statistical significance ($p = 0.15079$). Despite this, the Kalman filter provides a crucial stabilizing effect: the standard deviation of the global SNR is cut perfectly in half, dropping from $\pm 10.97$ dB in the baseline down to just $\pm 5.82$ dB. This proves that the Kalman architecture effectively standardizes highly volatile clinical EEG data, creating a much more reliable feature baseline for Brain-Computer Interface classification.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.95\textwidth]{Different_data_set_results.png}
    \caption{Grand Average statistical validation mapping. (Top) SNR comparison across the motor network showing a general positive trend and a massive reduction in inter-subject variance (tighter error bars). (Bottom) Residual Oscillation Index showing the statistically significant ($p < 0.01$) collapse of high-frequency artifacts across all channels.}
    \label{fig:diff_datasets_stats}
\end{figure}

% ==========================================
% CHAPTER 5: CLEAN & HUMANIZED VERSION
% ==========================================
\chapter{Conclusions, Clinical Implications, and Future Work}

\section{Research Objectives and Clinical Context}
Objective consciousness assessment in non-responsive patients with severe traumatic brain injuries remains a critical frontier in modern neurology. Patients surviving in vegetative or minimally conscious states, as well as those presenting Locked-In Syndrome (LIS), frequently retain distinct profiles of cognitive awareness that cannot be identified via standard behavioral metrics. This systematic clinical decoupling routinely induces diagnostic misclassifications, compromising therapeutic efficacy and long-term neurorehabilitation strategies \cite{BonomiThesis}.

Brain-Computer Interfaces (BCIs) address this limitation by bypassing physical motor pathways to monitor cortical electrical activity directly. The core objective of this dissertation was to design and evaluate an optimized preprocessing framework capable of isolating non-stationary Event-Related Potentials (ERPs) from highly corrupted scalp EEG records to systematically reveal preserved motor intentionality.

\section{Methodological Contributions}
This work systematically demonstrates that embedding a Linear Kalman Filter within a Time-Varying Autoregressive (TVAR) structure provides superior tracking capabilities for transient biological time-series. Rather than treating pipeline components in isolation, this research mapped the underlying mathematical and statistical conflicts within the legacy BVibes architecture, leading to several key pipeline updates:

\subsection{Phase Preservation via Zero-Phase FIR Architectures}
Infinite Impulse Response (IIR) configurations introduce non-linear phase distortion that systematically shifts and smears precisely time-locked waveforms like the RP. Replacing these models with a 10~Hz zero-phase Hamming-windowed sinc FIR filter neutralized these tracking delays, isolating low-frequency trends while effectively eliminating peripheral myogenic noise boundaries.

\subsection{Statistical Artifact Separation}
Comparative analysis confirmed that standard Infomax Independent Component Analysis (ICA) outperforms hybrid SOBI-wavelet techniques (such as FORCe) for offline neural feature preservation. Infomax ICA effectively exploits non-Gaussian source distributions, enabling clean isolation of ocular and high-amplitude physiological spikes without dampening low-frequency motor waveforms.

\subsection{Optimal Structural Placement of State-Space Operators}
A primary architectural finding of this thesis involves the relative sequence of spatial decomposition and state-space filtering blocks. Applying Kalman smoothing directly to continuous data streams prior to blind source separation attenuates local signal variances and inadvertently "Gaussianizes" leptokurtic artifact profiles. Because ICA relies explicitly on non-Gaussianity to resolve independent components, this sequence induces a collapse of the spatial unmixing matrix. Restructuring the pipeline to execute the Linear Kalman loop strictly after single-trial epoching preserves the statistical properties required by earlier artifact rejection layers while maximizing localized variance smoothing.

\section{Synthesis of Empirical Findings}

\subsection{Linearity vs. Non-Linearity in Parametric Parameter Estimation}
Cross-comparison of Linear Kalman Filters against non-linear variants (EKF and UKF) demonstrated that advanced linearization frameworks offer no statistical or morphological advantages for TVAR tracking. Because the transition matrix mapping the evolution of AR coefficients is defined as an Identity matrix ($A = I$), the parameter space remains strictly linear. Utilizing complex non-linear approximations adds significant computational latency without improving the structural Signal-to-Noise Ratio.

\subsection{Phase Disruption via Temporal Jitter Shifting}
This research establishes a fundamental conflict between state-space temporal estimation and traditional epoch alignment techniques like Woody's Jitter Compensation. While temporal shifting improves classical cross-trial arithmetic averages by aligning peak latencies, it destroys the phase continuity of the background stochastic noise. Since the recursive updates of the Kalman error covariance matrix ($P_k$) assume chronological noise progression, this phase disruption causes model divergence, distorting the true motor onset geometry. Consequently, jitter compensation must be systematically omitted when using state-space tracking.

\subsection{Topographical Spatial Spatial Optimization}
Empirical evaluations showed that combining post-epoching Linear Kalman filters with zero-phase FIR conditioning achieves exceptional noise suppression, isolating clear topographies of the Lateralized Readiness Potential (LRP). This spatial performance renders heavy spatial filters, such as the Spherical Spline Interpolation Spatial Filter (SISF), structurally redundant for localized motor pathway tracking, allowing for a lighter architecture with lower processing requirements.

\section{Clinical Applications and Diagnostic Performance}
The finalized pipeline was validated across several clinical databases spanning different task constraints and levels of physical injury. For voluntary datasets (GF89, AL85, MI99, SS89), the architecture systematically collapsed high-frequency noise bands, as synthesized in Table~\ref{tab:comprehensive_results}. 

In the TC99 database, where target intentions are routinely masked by elevated baseline variances, the framework successfully recovered the underlying RP morphology, shifting the localized onset signal-to-noise ratio from a negative baseline to a positive threshold of $+11.3$~dB.

The ultimate clinical diagnostic utility of this framework was validated using the involuntary,  \texttt{DS\_uni} dataset. Under completely involuntary conditions, the lack of physical muscle execution removes the typical structural reference triggers ($t = 0$), causing standard threshold pipelines to fail. By adapting baseline correction parameters to accommodate the prolonged cognitive transitions typical of severe brain injuries, the optimized post-epoching Kalman tracking successfully exposed the latent LRP profile. This empirical confirmation proves that the proposed state-space architecture can objectively identify and map trapped cognitive intent within an entirely paralyzed physiological framework.

\section{The Final Proposed Architecture}
The finalized BCI architecture is streamlined and highly effective: it resamples data to 128 Hz, applies a 10 Hz zero-phase FIR filter, utilizes Infomax ICA for artifact removal, and executes Linear Kalman Smoothing strictly on independent, un-shifted single-trial epochs.

\section{Limitations of the Current Framework}
\subsection{High-Density Topographical Dependency of ICA}
While Infomax ICA is highly effective, it heavily depends on having a sufficient number of electrodes. If applied to low-density clinical caps (e.g., 8 or 16 channels), the algorithm may struggle to mathematically separate the independent sources.

\subsection{Constraints of Offline Processing}
Currently, the pipeline requires the entire dataset to be recorded before processing begins. Features like zero-phase filtering (which looks forward and backward in time) prevent this current exact setup from being used instantaneously. 

\section{Future Work and Extensions}
\subsection{Translation to Real-Time Online BCIs}
The natural next step for this research is to move from offline analysis to a real-time, online Brain-Computer Interface. This will require replacing the non-causal FIR filters with low-latency causal filters and adapting the ICA into an adaptive, online spatial filter.

\subsection{Expansion into Advanced Parametric Domains: VARFI and VARMA}
Future iterations could expand the state-space estimation into the Vector Autoregressive Fractionally Integrated (VARFI) or VARMA domains \cite{Vakitbilir2025}. Incorporating fractional integration would allow the system to account for long-range memory dependencies in the brain, yielding deeper insights into the functional connectivity of damaged neural networks.

\subsection{Application to Broader Neurological Disorders}
Beyond hemiplegia and comatose states, this highly sensitive temporal tracking could be applied to track tremors in Parkinson's disease or detect early-onset neural decoupling in other severe neurodegenerative disorders.

\section{Concluding Remarks}
Combining advanced state-space mathematics with clinical neuroscience gives us a powerful tool to understand the damaged human brain. This thesis proves that by strictly following the mathematical rules of state-space modeling, preserving temporal phase continuity, and removing conflicting algorithms, we can mathematically bypass physical paralysis. By successfully extracting the hidden, conscious intentions of severely paralyzed patients, this optimized Kalman architecture moves the medical community one step closer to ending the tragedy of diagnostic misclassification in Disorders of Consciousness.


% ==========================================
% BIBLIOGRAPHY
% ==========================================
\cleardoublepage
\addcontentsline{toc}{chapter}{Bibliography}
\bibliographystyle{ieeetr}
\bibliography{references} 

\end{document}

