Biosignals - Basic signal processing

This post is by no means extensive, and emphasis is given to the conceptual understanding more than an in-depth treatment of the field. To begin we cannot escape from time series data or signals in research. Any time series data can be classified into continuous-time and discrete‐time signals. A discrete signal or discrete‐time signal is a signal that has been sampled at a certain interval from a continuous‐time. See below.

Fig-1: Example of a continuous signal (left) and discrete or digital signal (right)

What are Time and Frequency Domain?

In the time domain, any signal can be expressed as a combination of cosine and sine waveforms with differing periods and amplitudes. In the frequency domain, we reveal the frequency components of the same signal using some transformation.

A time-varying signal has different values across time, denoted as 𝓍(t). To reveal the frequency content of a signal, we use the Fourier Transform, where a signal can be expressed in the frequency domain, X(𝝎), with different frequency contents. The peak amplitude in the frequency domain will correspond to the harmonic components of the signal. This operation is invertible, that is, we can reconstruct the signal back to the time domain. The frequency is denoted either as 𝝎 (in radian unit) or  f (in hertz unit).

Briefly, it is also good to know that signals with different time‐domain characteristics have different frequency‐domain characteristics as follow:

       1) Continues‐time periodic signal → discrete non‐periodic spectrum

       2) Continues‐time non‐periodic signal → continuous non‐periodic spectrum

       3) Discrete non‐periodic periodic signal → continuous periodic spectrum

       4) Discrete periodic signal → discrete periodic spectrum.

The last transformation between time‐domain and frequency-domain is the most useful because discrete period signal is associated with both time‐domain and frequency, and that computers can only take finite discrete-time signals. 

What is Discrete Fourier Transform?

As the computer system cannot work with an almost-infinite continuous-time signal, we need to work on discrete non-period/period signals. A continuous signal 𝓍(t) has to be sampled at a certain rate F_s to obtain just a limited N number of samples, giving rise to a new discrete data 𝓍[n]. where n = 0, 1, 2, ..., N − 1; and 𝓍(n) is the signal amplitude at sample n, and N is the total number of samples. Suppose the discrete frequencies 𝝎_k = 2πk/N, then DFT operation is defined below to yield X[k], using an operation called the Discrete Fourier Transform (DFT). More technically, this DFT is a discrete or sampled version of the DTFT, Discrete-time Fourier Transform, a continuous function over a discrete time interval [−π,π].

One thing to note is that sampling a continuous signal must obey the Nyquist Theorem for the digital signal to be reliable.

Fig-2: An example of FFT in Matlab. The single-sided spectrum has double the amplitude of the FFT output.

The Fast Fourier Transform (FFT) is an efficient method for the evaluation of that operation using the Cooley-Tukey algorithm. It is very popular in various fields of signal processing being the default function in Matlab, Python, etc. The Matlab function returns a vector that appears to begin at the [0, F_s], where F_s is the sampling frequency. In reality, the discrete data can only contain up to F_s/2 based on Nyquist Theorem. We shall shift and adjust the mirroring to make the plot ranges from [−F_s/2, F_s/2]. The 'negative' frequencies come from the way the two-sided FFT is characteristically depicted because the Fourier Transform is mathematically defined on the interval [−∞,∞]. Refer to the figure below. Suppose there is a noisy sinewave of 2 sec, with a frequency of 3 Hz and sampled at 2000 Hz.

The Fourier Transform assumes that the signal is periodic and has infinite interval [−∞,∞]. When only a portion of data is analyzed, the data must be first truncated by a so-called window to preserve the frequency characteristics. An example of a window in digital signal processing (DSP) is a Hamming Window. Note: a window in the time domain is represented by a multiplication process, so it becomes a convolution in the frequency domain.

Fig-3: An example of stationary and non-stationary signals.

What is Power Spectral Density?

Power spectral density (PSD) function shows the strength of the energy as a function of unit frequency. In other words, it shows at which frequencies variations are strong or weak. The power spectrum is always applied to stationary signals, and has a real and non-negative value. To recap, if the statistical properties of the time series do not change over time, then that time series is said to be stationary. Look at an example below. Finally, the PSD of a fully uncorrelated random noise-like signal is pretty much flat over all frequencies, or 'white'. 

In practice, power spectral densities can be computed using a periodogram. This is done by applying the DFT of the signal and then computing the squared magnitude of the result. The spectrum amplitude and frequency are normalized by the sampling resolution employed to digitize the signal. Another popular method is Welch's method which gives a smoother PSD curve in exchange for reducing the frequency resolution. This method essentially uses 'frequency binning' of the original periodogram. First, the data segment is split up into m segments of a length L, but overlapping by d points. The overlapping segments are then windowed to taper the signal. This step is done because data segmentation yields discontinuity at both edges. Subsequently, the periodogram is calculated as in above and averages the spectra of these. If 𝒮_xx is the PSD of the signal from [−∞,∞], the one-sided version ranges 0 < f < ½ F_s, is usually denoted as G_xx.

Fig-4: An example of a periodogram of the same noisy signal as in the previous figure (using Matlab).

Fig-5: Useful block diagram showing how Welch's Method is achieved via FFT computation.

Note: In Matlab, pwelch function requires you to input nfft = number of FFT points. By convention, we set nfft to the power of 2 that is next above the N. For instance, if N = 1000 points, then nfft = 1024 (2^10). The number of overlap d is usually 50%, half of the window size L.

One more thing to share: if the signal is ergodic and non-windowed, this becomes the autocorrelation function (ACF) of that signal in the frequency domain (R_xx), but we won't discuss more mathematically. 

In layman's term, we say the autocorrelation is the degree of similarity of a signal with a delayed copy of itself over successive time interval. Informally, it is the similarity between a current observation and its past values. It can be used to understand patterns in the data, to assess periodicity, or randomness.

What is Cross Spectral Density?

By definition, the cross spectral density (CSD) measures the frequency components that are common between the signal 𝓍(t) and another signal 𝓎(t). Take note: if the PSD is related to the Fourier Transform of the autocorrelation function, the CSD are defined as the Fourier Transforms of the cross-correlation between two signals [𝓍(t) and 𝓎(t)]. The CSD is the distribution of power common to both signals per unit frequency. CSD uses the same computational methods as in PSD previously, e.g. Welch's method.

In layman's term, we say the cross-correlation is a measure of the similarity of the two signals as a function of the time lag relative to each other. In signal processing, both PSD and CSD are properties of the signal in the frequency domain.

What is Coherence?

Coherence (or spectral coherence) is a statistic to examine the relation between two signals or data sets in the frequency domain. Although in physics coherence can be defined in the spatial and temporal domain, in the case of time-series analysis and signal processing, it is always a function of frequency. Mathematically, coherence is a function of the power spectral densities, 𝒮_xx(f) and 𝒮_yy(f), and the cross power spectral density, 𝒮_xy(f), of signal 𝓍(t) and 𝓎(t). The magnitude of coherence varies in the interval [0, 1]. If coherence = 1, it means that both signals are perfectly correlated or linearly related each frequency. Conversely, if the value = 0 they are totally uncorrelated. To be valid, both signals have to be stationary within a period where the coherence is computed. 

Fig-6: An illustration of spectral coherence between two signals whose common frequencies are 100 and 200 Hz. In Matlab, mscohere function computes the coherence, where the magnitude can then be plotted.

Coherence has been used in neuroscience to assess the degree of similarity between two brain regions, which reflects the degree of functional coupling. The association between the brain and muscles, termed neuromuscular coupling, can be assessed by computing the coherence between the EEG (brain) and EMG (muscle) signals of interest.

Gaussian Mixture Models and EM algorithm

1. Gaussian Mixture Model (GMM)
In an earlier post, it was shown that we can partition our dataset into different subgroups or clusters, or how we hierarchically partition them step by step. One fundamental flaw of the k-means algorithm is the assumption that the dataset is more or less spherical. If it is not, then the algorithm will not do the job well (it will still partition, wrongly!) Another limitation of the k-means algorithm is that it is a hard-clustering algorithm, i.e. the clusters don't overlap. It does not provide any possibility or probability that the point can belong to any of the clusters.

A finite mixture model is a model comprised of an unspecified combination of multiple, but finite number of probability distribution functions.

A mixture model where the probability distribution functions are Gaussian or normal is called a Gaussian Mixture Model (GMM). It exactly attempts to provide a soft-clustering approach where clusters tend to overlap, i.e. there is some probability that a particular point belongs to more than one cluster. As an example, the probability density function of a 1D univariate dataset with 3 mixture components (K = 3) is below.

Fig-1: Probability density function (pdf) of a 1D dataset consisting of 3 mixture components. Note that the separation is not complete between adjacent clusters, and thus it reflects some probabilistic nature in it. [Ref 1]


So, a finite Gaussian mixture model means that we represent the dataset as being composed of K Gaussian functions, each has its own set of model parameters grouped as θ_k = { π, μ, Σ }:
     (1) A mean μ that defines its center.
     (2) A covariance matrix Σ that defines its width.
     (3) A mixing probability π or mixing weights.
In a univariate sense, a Gaussian or Normal distribution is modeled by the mean and variance (μ, σ²), satisfying the sufficient statistics. However, for a D-dimensional multivariate Gaussian distribution, we have a D × D covariance matrix Σ to represent the spread instead. In practice, the covariance matrix Σ can be further decomposed into several parameters that define certain geometric properties which at the same time provide certain constraints in the model for easy estimation (out of scope, but see Fig-3).

The mixing probability defines the probability that a certain point belongs to any clusters (Σπ = 1), which is essentially the cluster assignment problem in k-means. The larger this value is, the wider the shape of the Gaussian function will be, as it provides a higher chance that more points will be included under that function. A cluster assignment problem is estimated by the marginal probability distribution of a point in the mixture model, i.e. a weighted sum of the individual Gaussian function (defined by its three parameters above, Fig-1), where x is our dataset.

Fig-2: (a) Probability density function of 1D GMM; (b) Multivariate 2D gaussian probability density function; (c) An example of finite GMM consisting of 3 mixture components.

Fig-3: Comparison between k-Means and model-based GMM clustering algorithms. In the model-based method, we are not limited to spherical data and can estimate the shape of the distribution through an iterative process called EM algorithm. 


Now, how do we obtain the values for these parameters bearing in mind that we have K Gaussian functions? Of course, if we know the grouping sources beforehand (Fig-1), we can easily fit the Gaussian functions, compute the mean and standard deviation of each function (μ, σ²). In practice, we don't know the grouping sources, and the task can be even more challenging for the multidimensional Gaussian functions. Fortunately, we can estimate the relevant parameters with the help of the maximum likelihood estimation (MLE).

Note: Due to the fact that we need to determine model parameters, the soft-clustering approach using GMM is often called the model-based algorithms.

2. Expectation-Maximization Algorithm
Suppose we have a dataset x with n number of points in a GMM comprising K functions. The data is combined and the distributions are similar enough that it is not obvious to which distribution a given point may belong. To draw a sample from x, we select one of the components having a certain mixing probability π_k. Then with the help of a latent variable, z, for cluster assignment, we define a new probability (or in fact, likelihood) as a form of joint-probability distribution of the dataset x assuming that each data point is i.i.d.

Expectation-Maximization algorithm (EM) is an approach for maximum likelihood estimation with some latent variables. EM algorithm is a way to find maximum-likelihood estimates (MLE) for model parameters when the data is either incomplete, hidden, or multidimensional. The algorithm alternates between two steps:

• The first mode attempts to estimate the missing or latent variables called the estimation-step or E-step.
• The second mode attempts to optimize or maximize the parameters of the model to best explain the data called the maximization-step or M-step.

Nice references:
(1) Gaussian Mixture Model explained (more Maths involved!).
(2) Comparisons between k-means and EM algorithm
(3) Gaussian mixture modeling
# Gentle Introduction to Expectation Maximization
https://www.youtube.com/watch?v=REypj2sy_5U&ab_channel=VictorLavrenko

Diffusion Weighted Imaging

Although my PhD work involved functional MRI, I have not got a chance to work on diffusion-weighted imaging (DW). But more recently, more and more studies combine both MRI and DW together for more comprehensive measures.

Studies have shown that many developmental, aging and pathological processes influence the microstructural composition and architecture of the brain. As a primary consequence, water diffusion within the tissues will be altered by changes in the tissue microstructure and organization. Diffusion-weighted imaging (DW) becomes a popular tool to be used along with functional MRI. Another term closely associated with DW is diffusion tensor imaging (DTI). Whereas DWI is the raw data, containing diffusion maps in different directions, DTI is the process by which the raw data is transformed to calculate tensors (matrices that summarize the diffusion pattern in each voxel). In practice, DTI can be used for tractography as useful estimates of white matter fibers, albeit it does not capture or represent the actual myelinated fibers of the brain.

Water Diffusion and Tensor
Generally speaking, a diffusion model describes a transport phenomenon of a molecule from one location to another over time. In DW imaging, we are modeling the diffusion of water molecules, which are abundant inside a biological system. Water diffusion is very much influenced by the interactions with extracellular space, cellular membranes, and organelles. Cellular membranes hinder the diffusion of water, thus decreasing the mean squared displacement.

In fibrous tissues including white matter, water diffusion is relatively unrestricted along the fiber orientation. Conversely, it is highly restricted and hindered in the directions perpendicular to the fibers. In other words, the diffusion inside fibrous tissue is said to be directional-dependent or anisotropic. This phenomenon has been modelled using multivariate Gaussian distribution by Basser et al. (1994). In that model, a 3 × 3 variance-covariance matrix is used to capture diffusivity in a 3D space. The matrix is also known as a diffusion tensor D. The diagonal elements denote the diffusion variances along the X, Y and Z axes, and the off-diagonal elements are the covariance terms. For example, variability along X and Y directions are denoted by D_xy. Diagonalization of the diffusion tensor yields two outputs:

1. Eigenvalues λ, which describes the directions of the diffusion. Isotropic diffusion is represented with equal eigenvalues in the three directions. Conversely, the diffusion is anisotropic when the eigenvalues differ significantly.
2. Corresponding eigenvectors (e), apparent diffusivities along the axes of principal diffusion. The diffusion tensor can be visualized in Fig-1, with the eigenvectors defining the directions of the principal axes, with the radius defined by the eigenvalues. 

The solid shape shown in Fig-1 is also known as an ellipsoid, whose principal axes are in the direction of maximum diffusivity. In an isotropic condition (e.g. in CSF and gray matter), the ellipsoid will resemble a ball. However, in the neural tissues (e.g. white matter), the eigenvectors are more or less parallel to the tract orientation. Thus, diffusion tensors have been shown to be a sensitive tool to detect abnormality in the neural tissues.

Fig-1: (Top left) An illustration of fiber tracts with a certain orientation w.r.t scanner coordinate system. The fiber tissues exhibit directional dependence (anisotropy) on water diffusion. (Top right) The 3D diffusivity is modeled as a tensor, and visualized as an ellipsoid whose orientation is characterized by 3 eigenvectors (ϵ) and whose magnitude or length is characterized 3 eigenvalues (λ). The eigenvectors represent the major, medium, and minor axes of the ellipsoid and the eigenvalues represent the diffusivities in these three directions, respectively. (Bottom) This ellipsoid model is shown as a tensor, requiring a procedure known as matrix diagonalization. The major eigenvector (associated with the largest of the three eigenvalues), reflects the direction of maximum diffusivity, which, in turn, reflects the orientation of fiber tracts. Ref: Jellison et al, AJNR (2004).





On Image Acquisition
If we apply a certain MRI sequence that can manipulate this diffusion gradient, then we are able to see how the water diffuse in space. DW image acquisition is based on a specific but modified spin-echo sequence called the pulse gradient spin-echo (PGSE). In a way, it is a kind of a combination of spin-echo and gradient-echo sequence where there is a 90° and a 180° RF pulse, and a pair of diffusion gradient fields (a specific magnetic field with varying strength spatially) delivered on both sides of the 180° pulse. The purpose of this diffusion gradient is to allow the detection of phase differences in molecule spins.

Thus to produce DW imaging along the X-axis direction, for example, we apply strong magnetic field gradients along X. If molecules diffuse along X during the specified time interval, a signal attenuation will be observed compared to the signal without gradient.

For stationary water molecules, there will be no phase difference between before and after the gradient fields and no signal attenuation will occur. However, if there is a coherent flow in the direction of the applied gradient, a net phase difference will be introduced by the different amounts for each gradient field, resulting in signal attenuation. This phase difference is proportional to the gradient strength G, the duration of each gradient pulse 𝛿, and the time difference between the two gradient pulses Δ. It is well-known that PGSE is very sensitive to head motion. Because of this, a single shot echo pulse is used to acquire the signal readout in a very swift and efficient manner. The echo time (TE) is usually 100 msec in the DW sequence; with the repetition time (TR) between successive RF pulses to be rather long (6-7 sec).

Fig-2: Illustration of a pulse-gradient spin-echo MRI sequence used in the DW image acquisition (1 repetition).

The attenuated signal for moving molecules can be modeled as an exponential form containing the diffusivity and a b factor. This b factor depends only on the acquisition or MRI sequence parameters. Here the diffusivity is represented by the apparent diffusion coefficient (ADC) to indicate that the diffusion process is not free in the fibrous tissues, but restricted or modulated by many biological and anatomical factors. In the end, images are "weighted" by the diffusion process. This means that the signal is more attenuated the faster the diffusion and the larger the b factor is.

where S is the DW signal (or image), S₀ is the signal without any DW gradients, ADC is the apparent diffusion coefficient, and b factor is the diffusion-weighting that solely depends on the properties of the MRI sequence (gradient strength, time spacing). The optimum diffusion-weighting (b factor) for the brain is usually ~ 700 and 1300 s/mm² with 1000 s/mm² being the most common in use. From here, the 6 independent elements of the diffusion tensor D, i.e. 3 diagonal and 3 off-diagonal elements, may be estimated from the apparent diffusion coefficients using multiple linear least-squares methods or nonlinear modeling.

Fig-3: Illustration of different natures of anisotropy of the DW image inside the brain (from: FSL training slides).

Measurements tend to be quite sensitive to image noise, which can bias the anisotropy estimates. The accuracy of DTI measures (see below) may be improved by either: increasing the number of encoding directions or increasing the number of averages. Unfortunately, this will also affect the scan time during data collection. The image SNR can also obviously be improved by using larger voxels, although this will increase partial volume averaging of tissues, which can lead to errors in the diffusion tensor model. Moreover, image quality and thus spatial resolution may also depend on the application (clinical vs research). 

What is HARDI?

As stated above, the imaging process is usually obtained using b = 1000 s/mm2, N = 12 - 32 directions, parallel imaging with a very long repetition time (TR > 6 sec). DTI assumes water diffusion to have a Gaussian distribution. The estimation in 3D has six unknown coefficients to reconstruct, and thus a minimum of six diffusion-weighted images is required in addition to a baseline S0 image.

However, one great disadvantage of DTI imaging is dealing with crossing fibers. A technique called High Angular Resolution Diffusion Imaging or HARDI (Wedeen & Tuch, 2000) has the ability to discriminate multiple fiber populations crossing within the same voxel. To achieve that, HARDI requires the acquisition of > 50 gradient directions at a high b-value. In contrast, the traditional DTI only requires 6 directions at lower b-values. The higher angular resolution provides a more accurate representation of the 3D pattern of water diffusion within a voxel. All methods have in common the ability to provide the orientation of multiple white matter tracts within each voxel.

Two Main DTI Measures
In DTI, two most common measures are the mean diffusivity and anisotropy metric of the diffusion tensor. Mean diffusivity (MD, equivalent to ADC) is the average magnitude of diffusion in three main axes and quantified as the trace of the tensor divided by 3, tr(D)/3. Here, the trace is a sum of the diagonal elements of the tensor, which is equivalent to the average of the eigenvalues, λ. Many measures of anisotropy have been described, but the most widely used invariant measure of anisotropy is the fractional anisotropy (FA) (Basser and Pierpaoli, 1996). Both MD and FA are rotationally invariant.

In healthy people, WM regions have a wide range of FA values, ranging between 0.1 - 1.0, peaking at ~0.3, and much of this variation is due to crossing WM fibers (anisotropic!). Unfortunately, using FA should be handled with care. Although FA is a very sensitive and popular measure, it is not specific enough and does not capture the full shape of the diffusion tensor. Using this measure as a biomarker of WM integrity can be tricky. Some recent studies have suggested that the eigenvalues or their combinations demonstrate more specific relationships to white matter pathology. For example:
      (a) Radial diffusivity, D_R = (λ₂ + λ₃)/2, which appears to be modulated by WM myelin.
      (b) Axial diffusivity, D_A = λ₁, is more specific to axonal degeneration.
Some use cases of DTI in neuropathology is as follows. Example: demyelination could cause an increase in the radial diffusivity D_R, with little influence on the axial diffusivity measure. Increased tissue water in edema will increase the MD, whereas abnormal tissue growth may decrease the MD. In the case of ischemic stroke, FA and MD are found to increase and decrease during the acute phase. In the later chronic phase, FA and MD properties will reverse (decrease and increase).

Fiber orientation in the brain can be visualized using DW imaging. The strong assumption here is that the direction of maximum diffusivity in anisotropic voxels is an estimate of the major fiber orientation. Through some calculus operations, we can define the path taken by the white matter fibers. Applications of such tractography will be for future blog posts.

Motor Skill Learning
While there are ample studies examining the effect of motor learning on plasticity in the grey matter areas using functional neuroimaging, lesser studies have discussed changes in the white matter tracts as found using DWI. Results from those studies are not conclusive. Increase in FA values around the intraparietal sulcus has been observed after repeated juggling (Scholz, et al, 2009), and around the primary motor cortex related to motor adaptation (Landi, et al, 2011). Another study, however, reported significantly lower FA in both the left and the right CST in the musician group. Some recent studies have even made a claim to observe structural changes after a few sessions of motor skill learning (e.g. sequence tapping). 

Reference: Alexander AL, et al, Diffusion Tensor Imaging of the Brain (2007).