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.