System and method for electrophysiological oscillatory rhythms extraction and characterization

Abstract

There is provided a signal processing method and a signal processing system for rhythm extraction. The method comprises acquiring, using an electrophysiological sensing device, an input signal comprising a rhythmic component and an arrhythmic component, performing a discrete wavelet decomposition of the input signal in a time-scale domain to generate a time-scale representation of the input signal, the time-scale representation comprising a plurality of wavelet coefficients each having a scale and a discrete time associated therewith, processing the plurality of wavelet coefficients to remove therefrom a plurality of arrhythmic features indicative of the arrhythmic component of the input signal, thereby obtaining a plurality of processed coefficients, reconstructing, using the plurality of processed coefficients, a rhythmic signal comprising a plurality of rhythmic features indicative of the rhythmic component of the input signal, and outputting the rhythmic signal.

Description

SYSTEM AND METHOD FOR ELECTROPHYSIOLOGICAL OSCILLATORY RHYTHMS EXTRACTION AND CHARACTERIZATIONCROSS-REFERENCE TO RELATED APPLICATIONS

[0001] The present application claims the benefit of United States Provisional Patent Application No. 63 / 684622 filed on August 19, 2024, the contents of which are hereby incorporated by reference.FIELD

[0002] The improvements generally relate to the field of signal processing, and more particularly to extraction of oscillatory rhythms from electrophysiological signals.BACKGROUND

[0003] Spectral parametrization methods have been recently introduced to describe the properties of rhythmic and arrhythmic processes in the spectral domain. However, these approaches fall short for disentangling rhythms from arrhythmicity regarding the time series. In other words, no existing method allows to correct for the impact of arrhythmicity on rhythmic estimates time series. This is however crucial as arrhythmicity can impact the detection and the morphology of individual oscillations, therefore affecting the estimation of their frequency, amplitude, waveform, and their cross-frequency coupling estimated transiently.

[0004] Therefore, there is a need for improvement.SUMMARY

[0005] In accordance with one aspect, there is provided a signal processing method for rhythm extraction. The method comprises acquiring, using an electrophysiological sensing device, an input signal comprising a rhythmic component and an arrhythmic component, performing a discrete wavelet decomposition of the input signal in a time-scale domain to generate a time-scale representation of the input signal, the time-scale representation comprising a plurality of wavelet coefficients each having a scale and a discrete time associated therewith, processing the plurality of wavelet coefficients to remove therefrom a plurality of arrhythmic features indicative of the arrhythmic component of the input signal, thereby obtaining a plurality of processed coefficients, reconstructing, using the pluralityof processed coefficients, a rhythmic signal comprising a plurality of rhythmic features indicative of the rhythmic component of the input signal, and outputting the rhythmic signal..

[0006] In at least one embodiment in accordance with any previous / other embodiment described herein, performing the discrete wavelet decomposition of the input signal comprises applying a plurality of first wavelet functions to the input signal to obtain the plurality of wavelet coefficients, applying a linear regression to the plurality of wavelet coefficients to estimate a scaling exponent p* indicative of the arrhythmic component of the input signal, and applying a plurality of second wavelet functions to the input signal to generate the time-scale representation of the input signal, the plurality of second wavelet functions parametrized based on the scaling exponent.

[0007] In at least one embodiment in accordance with any previous / other embodiment described herein, each first wavelet function and each second wavelet function is a fractional spline wavelet function.

[0008] In at least one embodiment in accordance with any previous / other embodiment described herein, each first wavelet function is a scaled and translated version of a reference wavelet function, further wherein each first wavelet function is parametrized by a regularity parameter a indicative of a smoothness of the reference wavelet function.

[0009] In at least one embodiment in accordance with any previous / other embodiment described herein, each second wavelet function is parametrized by a + p* / 2.

[0010] In at least one embodiment in accordance with any previous / other embodiment described herein, the regularity parameter has a value greater than 2.

[0011] In at least one embodiment in accordance with any previous / other embodiment described herein, the rhythmic signal is reconstructed based on a linear combination of the plurality of processed coefficients and the plurality of first wavelet functions parametrized by the regularity parameter.

[0012] In at least one embodiment in accordance with any previous / other embodiment described herein, acquiring the input signal comprises acquiring an Electroencephalography (EEG) signal using at least one EEG electrode.

[0013] In at least one embodiment in accordance with any previous / other embodiment described herein, acquiring the input signal comprises acquiring a Magnetoencephalography (MEG) signal using at least one magnetometer.

[0014] In at least one embodiment in accordance with any previous / other embodiment described herein, acquiring the input signal comprises acquiring a Near-infrared spectroscopy (NIRS) signal using at least one near-infrared detector.

[0015] In accordance with another aspect, there is provided a signal processing system for rhythm extraction. The system comprises a processing unit and a non-transitory computer-readable medium having stored thereon program instructions executable by the processing unit for receiving an input signal from an electrophysiological sensing device, the input signal comprising a rhythmic component and an arrhythmic component, performing a discrete wavelet decomposition of the input signal in a time-scale domain to generate a time-scale representation of the input signal, the time-scale representation comprising a plurality of wavelet coefficients each having a scale and a discrete time associated therewith, processing the plurality of wavelet coefficients to remove therefrom a plurality of arrhythmic features indicative of the arrhythmic component of the input signal, thereby obtaining a plurality of processed coefficients, reconstructing, using the plurality of processed coefficients, a rhythmic signal comprising a plurality of rhythmic features indicative of the rhythmic component of the input signal, and outputting the rhythmic signal.

[0016] In at least one embodiment in accordance with any previous / other embodiment described herein, the program instructions are executable by the processing unit for performing the discrete wavelet decomposition of the input signal comprising applying a plurality of first wavelet functions to the input signal to obtain the plurality of wavelet coefficients, applying a linear regression to the plurality of wavelet coefficients to estimate a scaling exponent p* indicative of the arrhythmic component of the input signal, and applying a plurality of second wavelet functions to the input signal to generate the timescale representation of the input signal, the plurality of second wavelet functions parametrized based on the scaling exponent.

[0017] In at least one embodiment in accordance with any previous / other embodiment described herein, each first wavelet function and each second wavelet function is a fractional spline wavelet function.

[0018] In at least one embodiment in accordance with any previous / other embodiment described herein, each first wavelet function is a scaled and translated version of a reference wavelet function, and each first wavelet function is parametrized by a regularity parameter a indicative of a smoothness of the reference wavelet function.

[0019] In at least one embodiment in accordance with any previous / other embodiment described herein, each second wavelet function is parametrized by a + p* / ? / 2.

[0020] In at least one embodiment in accordance with any previous / other embodiment described herein, the regularity parameter has a value greater than 2.

[0021] In at least one embodiment in accordance with any previous / other embodiment described herein, the program instructions are executable by the processing unit for reconstructing the rhythmic signal based on a linear combination of the plurality of processed coefficients and the plurality of first wavelet functions parametrized by the regularity parameter.

[0022] In at least one embodiment in accordance with any previous / other embodiment described herein, the program instructions are executable by the processing unit for acquiring the input signal comprising acquiring an Electroencephalography (EEG) signal using at least one EEG electrode.

[0023] In at least one embodiment in accordance with any previous / other embodiment described herein, the program instructions are executable by the processing unit for acquiring the input signal comprising acquiring a Magnetoencephalography (MEG) signal using at least one magnetometer.

[0024] In at least one embodiment in accordance with any previous / other embodiment described herein, the program instructions are executable by the processing unit for acquiring the input signal comprising acquiring a Near-infrared spectroscopy (NIRS) signal using at least one near-infrared detector.

[0025] Many further features and combinations thereof concerning embodiments described herein will appear to those skilled in the art following a reading of the instant disclosure.DESCRIPTION OF THE FIGURES

[0026] In the figures,

[0027] Figs. 1A to 1 H are plots illustrating the fractional spline wavelets used in the method described herein and the rhythmic spectroscopy achieved using the method described herein, in accordance with one embodiment;

[0028] Figs. 2A, 2B, 2C, and 2D are plots illustrating the rhythmic signal estimation in the wavelet domain performed using the method described herein, in accordance with one embodiment;

[0029] Fig. 3 shows plots illustrating the wavelet rhythmic signal estimation and multiple oscillations obtained using the method described herein, in accordance with one embodiment;

[0030] Figs. 4 and 5 illustrate a flowchart of a signal processing method for rhythm extraction, in accordance with one embodiment; and

[0031] Fig. 6 is a block diagram of an example computing device, in accordance with one embodiment.

[0032] It will be noted that throughout the appended drawings that like features are identified by like reference numerals.DETAILED DESCRIPTION

[0033] Described herein are signal processing methods and systems for accounting for arrhythmicity within the time domain. In particular, there is provided a wavelet-based method that accounts forthe arrhythmic background in an electrophysiological time series. As used herein, the term “wavelet analysis” refers to a signal processing technique based on a mathematical expansion using a wave-like oscillation of finite length (also referred to as a “wavelet” or “little wave”) having an amplitude that begins at zero, increases or decreases, then returns to zero (see Fig. 1 A described further below). In one embodiment, the systems and methods described herein were applied to a non-rapid eye movement (NREM) sleep’s electrophysiological time series. It should however be understood that, although reference is made herein to the systems and methods described herein being applied to NREM time series acquired based on Electroencephalography (EEG) recordingsindicative of electrical activity in a subject’s brain, it should be understood that other electrophysiological or optical signals may apply including, but not limited to, EEG signals, intracranial EEG signals, surface (or scalp) EEG signals, Magnetoencephalography (MEG) signals, and Near-infrared spectroscopy (NIRS) signals. It should also be understood that, while reference is made herein to using the systems and methods described herein for monitoring the function of a brain, the systems and methods described herein may be applicable for monitoring other organs or tissues.

[0034] Electrophysiological signals, such as EEG signals, typically comprise rhythmic and arrhythmic components. As used herein, the terms “rhythmic” and “rhythmicity” (when used with reference to a signal, such as an EEG signal indicative of electrical activity in a subject’s brain) refers to recurring patterns of the signal that exhibit a consistent frequency (e.g., cycles per second or Hertz) and shape. These regular and repeating (i.e. periodic) patterns may be described as oscillations which are sustained and predictable. Rhythmic components of EEG signals may be used to understand brain function since these oscillations can be associated with various states and activities (e.g., alertness, sleep, cognitive processes). In contrast, the terms “arrhythmic” (or “non-rhythmic”) and “arrhythmicity” as used herein refer to a signal which lacks stable recurring patterns. This can be seen as an aperiodic component of the signal which forms the background activity (or “noise”). Arrhythmic components of EEG signals may be associated with various neurological conditions or disorders and changes in their characteristics may allow understanding of brain function and pathology. For instance, alterations in the arrhythmic component of an EEG signal may be linked to cognitive decline.

[0035] As understood by those skilled in the art, brain rhythms from bioelectric activity reflect synchronization of neuronal assemblies and are often inventoried using the Fourier analysis of electrophysiological brain recordings. Alongside brain rhythms, new hypotheses suggest that arrhythmic brain activity plays complementary roles in NREM sleep functions. In particular, arrhythmic activity recruits a broad range of frequencies and is usually expressed as 1 / 13decay in the power spectral density (i.e. the power of the signal decreases as frequency increases, in contrast with the peaks associated with rhythmic activity). In other words, arrhythmicity is usually characterized as a spectral 1 / f slope (also referred to as a “scaling exponent”) and is expressed by dynamical scale-free fluctuations unrelated to specific oscillations. As used herein, the term “scale-free” refers to processes which do not have a typical or characteristic temporal scale and are betterdescribed by the exponent of a power-law function, in contrast to processes which have a characteristic temporal scale and which can be defined by their mean frequency and a narrow band spectral dispersion.

[0036] NREM collectively refers to sleep stages 1 to 3 during which there is generally little to no eye movement. Each stage exhibits distinct electroencephalographic characteristics shown and quantified by EEG recordings. NREM sleep comprises a first stage (or NREM1) which is the lightest and shortest stage, lasting for about 5% of the sleep cycle, with slow eye movement, a second stage (or NREM2) which is the longest stage (e.g., about 45% of the sleep cycle) where no eye movement occurs and sleep becomes deeper (the heart rate and body temperature decrease), and a third stage (or NREM3) which is the deepest stage (or slow-wave sleep, SWS). Pure arrhythmic EEG signal represents around 80% of the recording time in NREM2 and more than 50% in NREM3 sleep. In addition to this arrhythmic background activity, NREM sleep is characterized by cardinal rhythms critically involved in overnight information processing, such as sleep slow waves (SSW) (e.g., with a delta band of 0.5-4 Hz), theta bursts (e.g., with a delta band of 6-10 Hz), sleep spindles (e.g., with a delta band of 8-16 Hz), and sharp-wave ripples (e.g., with a delta band of 100-200 Hz). These rhythms are organized in complex wave sequences, with sporadic amplitude-phase coupling across frequencies, which occur locally and between remote brain regions. For instance, sleep spindles usually occur after SSWs’ down states, within the transition period towards their upstate. The loss of specificity of this phase coupling may contribute to impaired memory consolidation.

[0037] Arrhythmicity is expressed at the signal level as a ubiquitous desynchronized background activity that nonlinearly interferes with transient oscillatory events. This coexistence introduces noise in time-domain analyses, which decreases the robustness of phase and amplitude estimates and influences the detection of genuine rhythmic oscillations. Currently, no method allows to correct for the impact of arrhythmicity on rhythmic estimates time series. In other words, current spectral parametrization methods cannot address the interference between background and rhythms in the time domain.

[0038] To overcome the deficiencies of existing techniques, a wavelet-based method is described herein, which can be used to disentangle rhythmic time series from an arrhythmic background, allowing for time-domain analysis of rhythms free fromarrhythmicity. In one embodiment, the method described herein leverages a waveletbased scale-free estimation alongside signal denoising to whiten 1 / activity during NREM sleep time series. This allows for the extraction of a rhythmic time series that can be further analyzed. For instance, the synthetized rhythmic signal may be used to identify new biomarkers.

[0039] A wavelet expansion of the electrophysiological signals was considered to resolve the entanglement between rhythms and arrhythmic background at the signal level. This approach produces a rhythmic time series that preserves the oscillatory content of the signal. The method described herein denoises the signal from the arrhythmic background by performing the synthesis of the rhythmic signal sR(t) in the wavelet basis IJ from wavelet coefficientswJkobtained through the processing of the original signal s(t), as follows:

[0040] The waveletis referred to as a “transient little wave” (see Fig. 1A) which has a temporal scale = 2> and is timely located at t = k 2 where k and j are integers, as given by the following equation:

[0041] As used herein, the term “wavelet basis” refers to an orthogonal set of wavelet functions into which a signal, to which a discrete wavelet transform (DWT) is applied, is decomposed to represent the signal at different resolutions (or scales). DWT is preferably used in the embodiments described herein because it is more computationally efficient and suitable for signal processing than other wavelet transform techniques, such as continuous wavelet transform (CWT). As understood by those skilled in the art, wavelet functions are translated (i.e. shifted across time) and scaled (i.e. dilated or compressed) versions of a common reference function, referred to herein as the “mother wavelet”, whose type may depend on the features to be determined from the input signal. Thus, as used herein, the term “wavelet decomposition” refers to the fact that an input signal isdecomposed into different scales (also referred to as “resolutions”) by applying scaled and translated versions of the mother wavelet. The scales are related to the signal’s frequency while being inverse thereto, such that high scales correspond to low frequencies and low scales correspond to high frequencies. The wavelet functions therefore serve as so-called “wavelet filters” which are applied to the input signal in order to achieve the wavelet decomposition and obtain different representations of the input signal in different frequency bands, allowing for signal analysis at different scales and locations. As used herein, the term “wavelet coefficients” refers to the coefficients wjk(or w-k) that weight the contribution (i.e. the amplitude) of each wavelet function at each scale and discrete time location.

[0042] In equation (1), which defines the synthesis of the rhythmic signal, indices j and k respectively label the scale and the time location of the wavelet. As used herein, the term “wavelet synthesis” refers to the use of the wavelet transform to reconstruct a signal (e.g., the rhythmic signal). The parameter a characterizes the smoothness of the mother wavelet. This processing is only possible if the wavelet expansion of the original (e.g., EEG) signal (see equation (3) below) can account for the presence of an arrhythmic scaling exponent / ?*, as follows:

[0043] In one embodiment, the wavelet expansion of the original time series (see equation (3) above) and the synthesis of the rhythmic signal (see equation (1) above) rely on the use of two (2) distinctive fractional spline wavelet bases allowing to control for the differentiation of fractional order a + / ?* / 2 and a, as will be described further below.

[0044] The derivation of equations (1) and (3) above will now be described in further detail. While a uniform self-similarity cannot fully account for the scale-free dynamics of neural activity, it is assumed that self-similarity and second-order statistics of the signal’s fluctuations are suitable at the timescale of each epoch. The scale-free properties can however vary between epochs. This corresponds to a local linearization (additive) of the underlying arrhythmic and rhythmic processes. Thus, assuming a unique scaling exponentfor each four (4) second epoch, a robust estimator of the / ?* exponent is considered from a discrete wavelet decomposition of the signal, as follows:

[0045] As noted above, the smoothness of the mother wavelet ^(t) is controlled by the regularity parameter a. The little wave oscillates with a frequency fo that depends on the parameter a. Therefore, the main frequency of the wavelets at scale j in the expansion of equation (4) is = 2J / o. Fig. 1A (described further below) illustrates an example of this kind of wavelet and Fig. 1 C (described further below) illustrates the multiresolution timescale representation of a signal (e.g., an EEG signal) in such a wavelet basis. As used herein, the term “time-scale” (an expansion of a signal that exhibits properties that depend on different temporal scales at each instant) refers to the fact that the behavior of a given signal is examined across different time scales, where time and frequency information about the signal is provided simultaneously. As can be seen from Fig. 1 C and as will be described further below, the time-scale representation of the signal is a collection of wavelet coefficients at various discrete times and scales.

[0046] The log-log regression allowing to estimate the scaling exponents in the spectral domain log r( ) = - / ? log f (where r( ) is the power spectral density) can be reformulated in the wavelet domain as follows:

[0047] with the second-order statistics of the wavelet coefficients as follows:

[0048] In particular, the log-log regression involves plotting the wavelet coefficients (obtained by performing a first wavelet analysis in which first wavelet functions parametrized by the regularity parameter a are applied to the input signal) against the wavelet transform’s scale, and applying a linear regression to determine the slope of the log-log plot and thus allows to estimate the scaling exponent denoted by / ?*. Theregression (illustrated in Fig. 1 D described further below) is done across scales ranging from ji (high frequencies) to y2(low frequencies). The / ?* estimate (see equation (5)) in the time-scale domain is less likely to be biased by the presence of transient oscillations and specific peaks in the Fourier domain. This is because this regression relies on transient fluctuations, mostly dominated by the scale-free background in periods when oscillations are and are not present. Thus, the range from ji (high frequencies) to y2(low frequencies) should cover the spectral domain of interest with respect to the rhythms to be further extracted. It is also desirable for this estimate to be reliable with respect to the second order statistics involved in equation (6). There is a compromise between the number of scales to be used in equation (5) and the duration of the epochs in which the statistics are computed. The lowest scales (high frequencies) should not include dominant noisy activity, whereas the largest scales (low frequencies) are limited by epoch duration and edge effects. In summary, the approach described herein being intermediate between a fully time-resolved continuous wavelet transform (not suitable for signal processing) and a Fourier analysis (not suitable for time-resolved analysis), this tradeoff between accuracy of the scale-free model and the precision in time proves to be a benefit of the approach described herein based on the multi-scale representation of the signal.

[0049] In order to estimate the rhythmic signal, the wavelet expansion of the input signal is considered with a wavelet basis accounting for the scale-free property of the signal as characterized with / ?* as calculated in equation (5), as follows:

[0050] In particular, the wavelet coefficients wj,k (scale j, time location fc) are adjusted for the scale-free activity and are given by the scalar product (or correlation) between the w, = signal and the little wave, i.e.I ,k

[0051] In Equation (7), the summation over j ranges from 1 to / , where the parameter J refers to the highest scales (lower frequency) upon which the wavelet analysis operations described herein are performed to obtain the rhythmic signal, similarly to a high-pass filter. This largest scale is chosen so that the frequency range of this multiresolution decomposition covers the entire frequency range of the rhythmic activities of interest. In one embodiment, the scale runs from 1 to 9 and J is set to 8. It should however beunderstood that these parameters may be varied depending on the properties of the input data, such that other embodiments may apply.

[0052] The discrete wavelet coefficients displayed in Fig. 1C (described further below) illustrate the stationary gaussian process that governs the wavelet coefficients at each scale, of a signal mostly driven by a scale-free process.

[0053] The fractional spline wavelet bases have been chosen to handle the non-integer regularity of the original signal quantified with the parameter a and the following equality can be further shown:

[0054] where:

[0055] The relation (8) suggests that a whitening correction combined with a change of wavelet basis allows to synthetize a time-resolved fractional differentiation, able to reduce the arrhythmicity from the original signal, as follows:

[0056] To summarize, the rhythmic signal is defined by the following linear combination:

[0058] computed from the wavelet analysis of the original signal (see equation (7).

[0059] A denoising shrinkage has been introduced prior to the whitening correction to keep the significant largest wavelet coefficients and shrink (or set to zero) smaller waveletcoefficients which are indicative of noise. When performing the denoising shrinkage, the wavelet coefficients are compared to a threshold value and wavelet coefficients with absolute values below the threshold value (i.e. small wavelet coefficients) being modified and wavelet coefficients with absolute values above the threshold value (i.e. large wavelet coefficients) being retained. Each epoch will be analyzed with its appropriate wavelet basis parametrized with a specific / ?* whereas the synthesis (see equation (11)) is performed in a unique wavelet basis parametrized with a. In practice, the regularity of the rhythmic signal is accomplished with a > 2 (i.e. the regularity parameter a has a value greater than 2).

[0060] The wavelet analysis of an original signal and the synthesis of a rhythmic signal are illustrated in Figs. 1A, 1 B, 1 C, 1 D, 1 E, 1 F, 1 G and 2E. Fig. 1A shows a plot 102 of an example fractional spline wavelet that could be used for analysis of an EEG signal as illustrated in plot 104 of Fig. 1 B. Fig. 1 C shows a plot 106 of a time-scale decomposition of the EEG signal s(t) displayed in Fig. 1 B using the wavelet basis generated with the fractional spline wavelet of Fig. 1A. Each box 107 in the plot 106 of Fig. 1 C represents a wavelet coefficient (s, i / >(4)) at a particular scale ( / ) and discrete time ( / < 2 / ). The real time is indicated on the horizontal x-axis. To help the reading of the time-scale representations, visible boxes 107 account for 99% of the overall wavelet power of the signal. Grey intensity (globally normalized) indicates the amplitude of the largest wavelet coefficient. Fig. 1 D shows a plot 108 of the logarithmic wavelet power of the EEG signal of Fig. 1 B with respect to the dyadic scale parameter. The linear regression will provide the scaling exponent as given by equation (5). Fig. 1 E shows a plot 110 of a discrete wavelet representation of the rhythmic signal resulting from processing the EEG signal of Fig. 1 B. The lacunar representation exhibits wavelet coefficients at large scale (low frequencies) that could account for an oscillatory pattern localized in time (confirmed with an inspection of the EEG signal of Fig. 1 B).

[0061] Although the arrhythmic characteristic of the original discrete wavelet coefficients is attenuated in the discrete wavelet representations of the rhythmic time series (Fig. 1 E), the wavelet power of the rhythmic signal does not exhibit the details of the oscillatory content since its representation is driven by scales. However, the Fourier power spectrum of the synthetized rhythmic signals may be appropriate to extract the oscillatory components, as shown in Figs. 1 F, 1 G, and 1 H. Fig. 1 F shows a plot 112 of the rhythmic signal corresponding to the EEG signal illustrated in plot 104 of Fig. 1 B andsynthetized from the wavelet representation shown in Fig. 1 E. The value of the / ?*- exponent of the reduced or removed arrhythmic component is indicated in parenthesis on plot 112. From a collection of physiologically equivalent epochs (NREM3 for instance), the average ofthe spectral power densities from the rhythmic signal provides an inventory of rhythms referred to herein as “spectroscopy” (as illustrated in Fig. 1 H, described further below, for epochs collected in sleep). In addition to this spectral signature, the distribution of the / ^-exponent is obtained, which characterizes the arrhythmicity of the rhythm’s background can be seen from Fig. 1 G. Fig. 1 G indeed illustrates a plot 114 showing, from 300 epochs (NREM3, Anterior Cingulate), the probability distribution of the / ^-exponents. Fig. 1 H illustrates a plot 116 showing, from 300 epochs (NREM3, Anterior Cingulate), the rhythmic spectroscopy provided with the mean (curve 118).

[0062] Figs. 2A, 2B, 2C, and 2D further illustrate the rhythm’s spectroscopy obtained by Fourier analyzing the resulting rhythmic time series sR(t) of equation (1 1) above. In particular, to illustrate the method described herein, an alpha (10.5 Hz) oscillation in the background of a simulated arrhythmic time series was considered, as shown in plot 202 of Fig. 2A. In particular, plot 202 provides an example of a simulated input time series containing a mixture of a scale-free activity ( / ? = 2.2; arrhythmic background) and an alpha oscillation (four seconds 10.5 Hz rhythm). Plot 204 of Fig. 2B shows the Fourier spectra of the original time series s(t) displayed in plot 202, mostly dominated by the dominant arrhythmic background. The time series of plot 202 was processed using the method described herein to produce a rhythmic time series, as shown in plot 206 of Fig. 2C. In particular, the outputof the method described herein is shown in plot 206, revealing the rhythmic features indicative of the rhythmic component of the input time series synthesized from the wavelet’s coefficients filtered for the arrhythmic features indicative of the arrhythmic component of the input time series. Using the method described herein, the arrhythmic drift which dominated the original time series (see plot 202 Fig. 2A) is thus reduced or removed in the reconstructed signal (see plot 206 of Fig. 2C). The spectral density of the rhythmic component provided with the method described herein, as shown in plot 208 of Fig. 2D, exhibits a clear signature of the oscillation in comparison with the original Fourier spectrum (see plot 204 of Fig. 2B). In particular, the Fourier spectra of a singleepoch highlight the clear signature of the simulated oscillation that was added to the arrhythmic background, with a marginal arrhythmic contribution. This averaged across multiple epochs defining the rhythm’s spectroscopy for further analysis. In Fig. 2D,the black dotted line is the original spectra of the simulated oscillation, without the addition of an arrhythmic background.

[0063] It should be understood that, while reference is made herein to using a shrinkage approach to smooth spurious fluctuations, other embodiments may apply. For instance, in some embodiments, instead of using standard detection criteria into the rhythmic time series, matching pursuit methods may be applied to the rhythmic signal to better identify basic patterns corresponding to NREM rhythms. This may allow the development of detectors relaxing the need of prior detection criteria.

[0064] Another batch of simulations was conducted to assess the ability to address multiple rhythmic oscillations. With a similar strategy as for the neural mass simulations, a random set of 300 scale-free timeseries (mean exponent = 2.5 and variance = 0.01) was considered to which a combination of two oscillations (3 Hz and 13 Hz) were added in a 200 epochs subset (see plots 304 and 306 of Fig. 3). The simulation consisted of varying the amplitude of the 3 Hz oscillation only, keeping constant the amplitude of the 13 Hz oscillation. As previously, 60 samples were randomly selected from the pool of 300 epochs for analysis.

[0065] Plot 302 of Fig. 3 illustrates an example of a scale-free time series. Plot 304 of Fig. 3 illustrates an example of rhythmic oscillations to be added to plot 302. Plot 306 of Fig. 3 illustrates the spectral signature of the rhythmic component shown in plot 302. Plot 308 of Fig. 3 illustrates the response of the method described herein regarding the amplitude of the 3 Hz-oscillation with respect to varying amplitude of this oscillatory component. Plot 310 of Fig. 3 illustrates the response of the method described herein regarding the constant amplitude 13 Hz-oscillation with a varying amplitude of the 3 Hz- oscillation. Plot 312 of Fig. 3 illustrates the response of the method described herein regarding the scale-free exponent in a varying amplitude of the 3 Hz-oscillatory component.

[0066] As expected, the rhythmic amplitude estimation for the delta oscillation (amplitude out) scaled linearly with the amplitude of the simulated oscillation (amplitude in; see plot 308 of Fig. 3, b=0.87, p<0.0001 , R2=0.85). There was also no relationship between the amplitude of the simulated delta oscillation (amplitude in) and the estimated maximal sigma amplitude (see plot 310 of Fig. 3; b=-0.01 , p=0.27, R2< 0.01). The methoddescribed herein slightly overestimates the calculated exponent in presence of multiple oscillations, but this does not scale with the amplitude of the simulated delta oscillation (see plot 312 of Fig. 3, constant=2.7, p<0.001 , b=-0.03, p=0.31 , R2=0.001).

[0067] The simulations described above thus support the fact that rhythmic amplitude and scaling exponents outputs obtained using the method described herein are reliable to assess single or concurrent oscillations (sigma and delta), with minor caveats.

[0068] The method described herein is suitable for processing signals that exhibit well- defined scale-free properties, such as a consistent spectral slope ( / ?). The actual slope estimator used in the method described herein assumes a pure autosimilar signal. Intracranial EEG (iEEG) recordings are measured closer to the neural source, enabling a local measure of neural activity, in contrast to scalp EEG, which measures widespread cortical activity. As a result of this localized measure, scale-free dynamics are observed and quantified more clearly in iEEG data. Therefore, the current slope estimator used in the method described herein remains a reliable approach for iEEG recordings.

[0069] Nonetheless, iEEG recordings require invasive procedures, unlike scalp EEG, which is a non-invasive and widely used method in clinical settings. Thus, in some embodiments, it may be desirable to adapt the method described to scalp EEG recordings. Scalp EEG recordings typically exhibit a low signal-to-noise ratio (SNR) due to contamination from factors such as muscle activity (e.g., measured by electromyography, or EMG, signals) and eye movements. These elements can distort the estimation of the aperiodic component of the signal. Consequently, scalp recordings rarely adhere to a perfect power-law decay. To address this issue, the slope estimation method described herein above may be modified by relaxing the assumption of a perfectly auto-similar signal. Unlike the slope estimator used for iEEG, which applies scale-dependent weighting to reflect a power-law, the proposed modification does not rely on this behaviour. It is instead proposed to fit a linear model to the log-power across scales in the multiresolution decomposition of the EEG signals, making the slope estimation method more tolerant to deviations from perfect scale-free dynamics.

[0070] Fig. 4 and 5 illustrate flowcharts 400 and 500 of an exemplary signal processing method for rhythm extraction, in accordance with one embodiment. At step 402, an input signal (see plot 502 on Fig. 5) comprising a rhythmic component and an arrhythmiccomponent is acquired using an electrophysiological sensing device 403. The input signal is a time-series originally represented in the time domain. As described above, the input signal may comprise an electrophysiological or optical signal including, but not limited to, an EEG signal indicative of an electrical activity of a brain of a subject, a MEG signal, and a NIRS signal. The input signal may be acquired using an electrophysiological sensing device, such one or more EEG electrodes, one or more magnetometers, one or more near-infrared detectors, or other electrophysiological sensing devices configured to acquire electrophysiological signals related to a subject. For instance, the EEG electrode(s) may be placed on a surface of a subject’s brain (or within brain tissue) to acquire intracranial EEG signals). The EEG electrode(s) may alternatively be placed on the subject’s scalp to acquire surface (or scalp) EEG signal(s). The magnetometers) may be configured to acquire MEG signals) by measuring the magnetic fields produced by the electrical activity in the subject’s brain. The near-infrared detector(s) may be configured to acquire NIRS signals) by shining near-infrared light onto a tissue of the subject and measuring the amount of light absorbed or scattered by the tissue, thereby determining oxygenation and blood flow in the tissue.

[0071] Step 404 comprises performing a discrete wavelet decomposition of the input signal in a time-scale domain to obtain a discrete time-scale representation of the input signal. The time-scale representation comprises a plurality of wavelet coefficients. Step 404 is performed in the manner described herein above in relation to equation (3). In particular and as illustrated in Fig. 5, step 404 comprises sub-steps 504a and 504b. Step 504a comprises performing a first wavelet analysis (using a first wavelet basis) to estimate the scaling exponent / ?* (indicative of the arrhythmic component of the input signal) from first fractional wavelets filters (or functions) parametrized by a regularity parameter a indicative of a smoothness of the wavelet defined by the first wavelet analysis. In particular, the first wavelet functions are applied to the input signal to obtain the wavelet coefficients, and a linear regression is applied to the wavelet coefficients (in the manner described herein above) to estimate the scaling exponent. Step 504b follows step 504a and comprises performing a second wavelet analysis (using a second wavelet basis) to generate the discrete time-scale representation of the input signal using second fractional wavelet filters (or functions) parametrized by a + / ?* / 2. In particular, the second wavelet functions are applied to the input signal to generate the discrete time-scale representation. Plot 506 of Fig. 5 illustrates the resulting discrete time-scale representation.

[0072] Step 406 comprises using (i.e. processing) the wavelet coefficients obtained at step 404 to reduce arrhythmic features indicative of the arrhythmic component of the input signal. Step 406 is performed in the manner described herein above in relation to equation (1). Plot 508 of Fig. 5 illustrates one embodiment of the discrete time-scale representation of the input signal with the arrhythmic component removed.

[0073] Step 408 comprises synthesizing (i.e. reconstructing), using the coefficients processed at step 408, a rhythmic signal comprising rhythmic features indicative of the rhythmic component of the input signal. For this purpose, wavelet synthesis is performed with the processed coefficients, using the fractional filters parametrized by the regularity parameter a.

[0074] Step 410 comprises outputting the rhythmic signal, which is shown in plot 510 of Fig. 5. The rhythmic signal may be output in any suitable manner, such as rendered on a screen of a computing device. Plots 512a and 512b of Fig. 5 respectively illustrate the power spectral density of the input signal (plot 502) obtained at step 402, and the power spectral density of the rhythmic signal (plot 510) output at step 410.

[0075] In some embodiments, parts or all of the method described herein (e.g., method 400 of Fig. 4) are performed by a computing device 600, as illustrated in Fig. 6. The computing device 800 comprises a processing unit 602 and a memory 604 which has stored therein computer-executable instructions 606. The processing unit 602 may comprise any suitable devices configured to cause a series of steps to be performed such that instructions 606, when executed by the computing device 600 or other programmable apparatus, may cause functions / acts / steps described herein to be executed. The processing unit 602 may comprise, for example, any type of general-purpose microprocessor or microcontroller, a digital signal processing (DSP) processor, a CPU, an integrated circuit, a field programmable gate array (FPGA), a reconfigurable processor, other suitably programmed or programmable logic circuits, or any combination thereof.

[0076] The memory 604 may comprise any suitable known or other machine-readable storage medium. The memory 604 may comprise non-transitory computer readable storage medium, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. The memory 604 may include a suitable combination of anytype of computer memory that is located either internally or externally to device, for example random-access memory (RAM), read-only memory (ROM), electro-optical memory, magneto-optical memory, erasable programmable read-only memory (EPROM), and electrically-erasable programmable read-only memory (EEPROM), Ferroelectric RAM (FRAM) orthe like. Memory 604 may comprise any storage means (e.g., devices) suitable for retrievably storing machine-readable instructions 606 executable by processing unit 602.

[0077] The method described herein may be implemented, in part or entirely, in a high level procedural or object oriented programming or scripting language, or a combination thereof, to communicate with or assist in the operation of a computer system, for example the computing device 600. Alternatively, the method may be implemented in assembly or machine language. The language may be a compiled or interpreted language.

[0078] Embodiments of the method described herein may also be considered to be implemented by way of a non-transitory computer-readable storage medium having a computer program stored thereon. The computer program may comprise computer- readable instructions which cause a computer, or more specifically the processing unit 602 of the computing device 600, to operate in a specific and predefined manner to perform the functions described herein.

[0079] Computer-executable instructions may be in many forms, including program modules, executed by one or more computers or other devices. Generally, program modules include routines, programs, objects, components, data structures, etc., that perform particular tasks or implement particular abstract data types. Typically, the functionality of the program modules may be combined or distributed as desired in various embodiments.

[0080] In some embodiments, the systems and methods described herein may allow time-domain analyses of brain rhythms free from arrhythmicity, unlike conventional spectral methods which focus on properties of periodic peaks and aperiodic background. In comparison with spectral approaches, the method described herein estimates the spectral slope ( / ?), which is indicative of the arrhythmic component, in the time-scale domain and may therefore be less likely to be biased by the overlap between rhythmic and arrhythmic fluctuations at low frequencies. This time-scale approach also accountsfor the nonstationary nature of cerebral electrophysiological activity, which allows to perform reliable spectroscopy and time-domain analyses of brain rhythms.

[0081] In some embodiments, the systems and methods described herein may have the enhanced capability to characterize the hallmark features of NREM sleep rhythms. Furthermore, the systems and methods described herein may reveal rhythmic time series, highlighting key NREM sleep rhythms with improved sensitivity to their phase-amplitude coupling.

[0082] In some embodiments, the systems and methods described herein may provide improvements in terms of computational efficiency and computational speed compared to other computer-implemented signal processing techniques for rhythm extraction. In some embodiments, the method described herein is implemented on Graphics Processing Units (GPUs), allowing for a parallel computation across epochs, which provides significant performance improvements compared to existing techniques.

[0083] The above description is meant to be exemplary only, and one skilled in the art will recognize that changes may be made to the embodiments described without departing from the scope of the invention disclosed. Still other modifications which fall within the scope of the present invention will be apparent to those skilled in the art, in light of a review of this disclosure.

[0084] Various aspects of the systems and methods described herein may be used alone, in combination, or in a variety of arrangements not specifically discussed in the embodiments described in the foregoing and is therefore not limited in its application to the details and arrangement of components set forth in the foregoing description or illustrated in the drawings. For example, aspects described in one embodiment may be combined in any manner with aspects described in other embodiments. Although particular embodiments have been shown and described, it will be apparent to those skilled in the art that changes, and modifications may be made without departing from this invention in its broader aspects. The scope of the following claims should not be limited by the embodiments set forth in the examples but should be given the broadest reasonable interpretation consistent with the description as a whole.

Claims

WHAT IS CLAIMED IS:

1. A signal processing method comprising: acquiring, using an electrophysiological sensing device, an input signal comprising a rhythmic component and an arrhythmic component; performing a discrete wavelet decomposition of the input signal in a time-scale domain to generate a time-scale representation of the input signal, the timescale representation comprising a plurality of wavelet coefficients each having a scale and a discrete time associated therewith; processing the plurality of wavelet coefficients to remove therefrom a plurality of arrhythmic features indicative of the arrhythmic component of the input signal, thereby obtaining a plurality of processed coefficients; reconstructing, using the plurality of processed coefficients, a rhythmic signal comprising a plurality of rhythmic features indicative of the rhythmic component of the input signal; and outputting the rhythmic signal.

2. The method of claim 1 , wherein performing the discrete wavelet decomposition of the input signal comprises: applying a plurality of first wavelet functions to the input signal to obtain the plurality of wavelet coefficients; applying a linear regression to the plurality of wavelet coefficients to estimate a scaling exponent / ?* indicative of the arrhythmic component of the input signal; and applying a plurality of second wavelet functions to the input signal to generate the time-scale representation of the input signal, the plurality of second wavelet functions parametrized based on the scaling exponent.

3. The method of claim 2, wherein each first wavelet function and each second wavelet function is a fractional spline wavelet function.

4. The method of claim 2, wherein each first wavelet function is a scaled and translated version of a reference wavelet function, further wherein each first wavelet function is parametrized by a regularity parameter a indicative of a smoothness of the reference wavelet function.

5. The method of claim 4, wherein each second wavelet function is parametrized by a + p*!2.

6. The method of claim 4, wherein the regularity parameter has a value greater than 2.

7. The method of claim 4, wherein the rhythmic signal is reconstructed based on a linear combination of the plurality of processed coefficients and the plurality of first wavelet functions parametrized by the regularity parameter.

8. The method of claim 1 , wherein acquiring the input signal comprises acquiring an Electroencephalography (EEG) signal using at least one EEG electrode.

9. The method of claim 1 , wherein acquiring the input signal comprises acquiring a Magnetoencephalography (MEG) signal using at least one magnetometer.

10. The method of claim 1 , wherein acquiring the input signal comprises acquiring a Nearinfrared spectroscopy (NIRS) signal using at least one near-infrared detector.11 . A signal processing system comprising: a processing unit; and a non-transitory computer-readable medium having stored thereon program instructions executable by the processing unit for: receiving an input signal from an electrophysiological sensing device, the input signal comprising a rhythmic component and an arrhythmic component; performing a discrete wavelet decomposition of the input signal in a time-scale domain to generate a time-scale representation of the input signal, the timescale representation comprising a plurality of wavelet coefficients each having a scale and a discrete time associated therewith; processing the plurality of wavelet coefficients to remove therefrom a plurality of arrhythmic features indicative of the arrhythmic component of the input signal, thereby obtaining a plurality of processed coefficients; reconstructing, using the plurality of processed coefficients, a rhythmic signal comprising a plurality of rhythmic features indicative of the rhythmic component of the input signal; and outputting the rhythmic signal.

12. The system of claim 11 , wherein the program instructions are executable by the processing unit for performing the discrete wavelet decomposition of the input signal comprising: applying a plurality of first wavelet functions to the input signal to obtain the plurality of wavelet coefficients; applying a linear regression to the plurality of wavelet coefficients to estimate a scaling exponent / ^’indicative of the arrhythmic component of the input signal; and applying a plurality of second wavelet functions to the input signal to generate the time-scale representation of the input signal, the plurality of second wavelet functions parametrized based on the scaling exponent.

13. The system of claim 12, wherein each first wavelet function and each second wavelet function is a fractional spline wavelet function.

14. The system of claim 12, wherein each first wavelet function is a scaled and translated version of a reference wavelet function, further wherein each first wavelet function is parametrized by a regularity parameter a indicative of a smoothness of the reference wavelet function.

15. The system of claim 14, wherein each second wavelet function is parametrized by a + 72.

16. The system of claim 14, wherein the regularity parameter has a value greater than 2.

17. The system of claim 14, wherein the program instructions are executable by the processing unit for reconstructing the rhythmic signal based on a linear combination of the plurality of processed coefficients and the plurality of first wavelet functions parametrized by the regularity parameter.

18. The system of claim 11 , wherein the program instructions are executable by the processing unit for acquiring the input signal comprising acquiring an Electroencephalography (EEG) signal using at least one EEG electrode.

19. The system of claim 11 , wherein the program instructions are executable by the processing unit for acquiring the input signal comprising acquiring a Magnetoencephalography (MEG) signal using at least one magnetometer.

20. The system of claim 11 , wherein the program instructions are executable by the processing unit for acquiring the input signal comprising acquiring a Near-infrared spectroscopy (NIRS) signal using at least one near-infrared detector.

Images