Open Access
Add to BibSonomy
Abstract
Here we show a method of the lock-in amplifying near-infrared signals originating within a human brain. It implies using two 90-degree rotated source-detector pairs fixed on a head surface. Both pairs have a joint sensitivity region located towards the brain. A direct application of the lock-in technique on both signals results in amplifying common frequency components, e.g. related to brain cortex stimulation and attenuating the rest, including all components not related to the stimulation: e.g. pulse, instrumental and biological noise or movement artefacts. This is a self-driven method as no prior assumptions are needed and the noise model is provided by the interfering signals themselves. We show the theory (classical modified Beer-Lambert law and diffuse optical tomography approaches), the algorithm implementation and tests on a finite element mathematical model and in-vivo on healthy volunteers during visual cortex stimulation. The proposed hardware and algorithm complexity suit the entire spectrum of (continuous wave, frequency domain, time-resolved) near-infrared spectroscopy systems featuring real-time, direct, robust and low-noise brain activity registration tool. As such, this can be of special interest in optical brain computer interfaces and high reliability/stability monitors of tissue oxygenation.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Light penetrating tissues is highly diffused (random scattering direction after a few millimetres of travel) and as such, methods of focusing the diffused light itself or shaping the origin of carried information are in high demand [1]. In applications of an adult brain imaging, based on the near-infrared spectroscopy (NIRS), light carries mixed information on haemodynamic of grey matter and superficial head tissues. This effect limits significantly potential of clinical and/or neurophysiological applications of the functional NIRS technique. Thus, research focuses on removing the extracerebral information from measured signals [2–4] as this can lead to false positives and false negatives in functional NIRS (fNIRS) [5].
An approach frequently used to reduce the extracerebral contamination component utilizes an observation that the detected light penetration depth increases with the source-detector distance on a head surface. The extracerebral information can be removed with a linear regression using signals at short (e.g. ≤15 mm) source-detector distances in post processing [6] or/and by the high-density diffuse optical tomography (HD-DOT) approach [7–10] where depth/spatial discrimination is based on analysis of signals registered at thousands (over 10 thousand) of overlapping source-detector pairs at a varying spatial separation. The time-resolved NIRS approach offers a depth discrimination as well [11]. A measured distribution of time of flight of photons (DTOF) provides the so-called late and early photons [12–14] that travel through paths of different lengths. Hence, with varying probability of penetrating at different depths. Further, statistical moments of the measured DTOFs [15] reveal varying sensitivity in depth [16], depending on the moment order. Thus, allowing for the penetration depth discrimination [17]. A selective measurement focused towards the deeper head structures was recently proposed in [18,19]. This approach uses a spatial differencing of in-line source-detector pairs with slightly shifted detectors to shape the total resulting sensitivity as focused towards deeper (brain) tissues. One may use this method to design a desired and narrowed sensitivity region of an optical measurement [20,21].
We put forward a following research hypothesis focused on potential reduction of the influence of extracerebral tissue layers on the optical signals measured on the surface of a head in the reflectance geometry: a brain originating component of an optical signal can be amplified using interference methods where spatial sensitivity distributions of interfering signals have a common volume within the brain region of interest. Thus, the contrast to noise ratio of the haemodynamic response appearing in the brain cortex can be increased.
A method of ultrasonically encoded optical focusing was presented in [1,22] where a small tissue volume is stimulated with a focused ultrasound wave. The wave frequency becomes visible in an optical signal passing through the stimulated tissue volume. The optical signal amplitude oscillating at the ultrasound wave frequency can originate only within the small stimulated tissue volume. However, this method is difficult to apply in brain imaging as the skull creates a barrier hard to penetrate for the ultrasound.
Here we propose to utilize the natural, biological variability present in optical signals to limit the volume from which the information origins. The general concept can be, at some degree, compared to the mentioned above idea of encoding diffusely scattered light with an ultrasound wave focused inside a tissue [22]. The tissue displacement induced by an ultrasound modulation is here replaced with a brain functional stimulation occurring within a region of the grey matter. We propose to interfere (lock-in) optical signals having spatial sensitivity profiles overlapping within the brain cortex region of interest. The absorption sensitivity profile quantitates the volume light penetrates for a given source-detector pair [8] and reflects sensitivity of the measurement to a change in the tissue absorption. Considering signals measured at two source-detector pairs with overlapping sensitivity volumes, the signals’ interference strength will follow haemodynamic changes originating at the same time within the overlapping sensitivity volume.
Here we show principles of the method we call lock-in functional near-infrared spectroscopy. The lock-in operation is applied to a geometry of two overlapping source-detector pairs where the measurement sensitivity is at the geometry centre and focused within the tissue, towards the brain. We present how to apply the lock-in functional NIRS into the modified Beer-Lambert law and the diffuse tomography reconstruction. The concept is tested against finite element numerical model and using in-vivo brain signals registered during visual cortex stimulation. The research suggests the lock-in method might be of interest in brain computer interface applications. By design, it provides brain activity in real-time, is focused within the grey matter, with high contrast to noise ratio, suppressed biological noise and movement artefacts.
2. Methods
2.1 Motivation
We consider a measurement geometry at two 90-degree rotated source-detector pairs as shown in Fig. 1. We will call the pairs arms or reference arms, interchangeably. A change in light intensity attenuation ΔA measured on the surface depends on the spatial location r of the local absorption perturbation (change) Δμa(r). The relation ΔA/Δμa(r) = J(r) = J, called the Jacobian (sensitivity) is expressed in length units and represents the mean photons path length while traveling through the absorption perturbation volume. The sensitivity spatial distribution has the so-called banana shape as in Fig. 1(a)-b.
Fig. 1. Sensitivity volumes of light attenuation measurements to changes in absorption for the proposed geometry of two 90-degree rotated source-detector pairs (30 mm source-detector distance). Results of FEM analysis (details of the model presented in method section) for a two-layer model where the upper layer is 13 mm thick. Sensitivities (mean partial path lengths) for pair #1 (a) and pair #2 (b), the sum of both (c) and the square root of the element-wise product (d) are within the same colour scale.
Here, we have two sensitivity volumes J1 and J2 as shown in Fig. 1(a)-(b) for both pairs. Let us consider two scenarios where a new sensitivity profile is shaped: (i) by superimposing J1 and J2 or (ii) deriving the square root of the element-wise product of J1 and J2 (the square root of the Hadamard product $\sqrt {{{\textbf J}_1} \circ {{\textbf J}_2}\; } $). The superposition of both sensitivities is presented in Fig. 1(c). This operation does not lead to a solution of the problem as the superimposed sensitivity is less selective (more non-unique regions) and covers large volume of the medium. The goal is the opposite. As shown in Fig. 1(d), the square root of the Hadamard product approaches zero outside of the volume in which sensitivities J1 and J2 intersect. Therefore, we focus further analysis on the second scenario where the element-wise product of sensitivities is used.
2.2 Theory
A classical lock-in amplification mechanism (equivalent of a homodyne detection followed by a low-pass filtering) is used to decode from a signal amplitude and/or phase at a given carrier frequency and is widely used in high-bandwidth communication, e.g. it forms the basis of Wi-Fi networks (the quadrature modulation/demodulation). A great benefit of the lock-in approach is that the amplitude/phase at the carrier frequency can be detected even if it is under the noise floor (orders of magnitude). Here we propose to use the lock-in amplification where the signal to decode from and the signal template (carrier) are given by the two source-detector pairs as shown in Fig. 1. Hence, components present in both signals at the same time will amplify. This is true for the brain components originating at the same time within the common sensitivity region as in Fig. 1(d). Furthermore, it is expected other ‘non-periodic’ signals such as biological noise and/or movement artefacts will be attenuated. Where the biological noise is considered as all that happens outside of the brain (i.e. not the cortex haemodynamic response). The digital lock-in amplification will be carried out within a moving window of length T. As such, the filtering operation will have a low-pass nature with a cut-off frequency around 1/T. Therefore, any fast changes such as movement artefacts or instrumental noise will be removed preserving amplified common components (within the stop-band, e.g. the heartbeat). The moving window length expressed in samples Nw should give best results when Nw is a power of 2 as the fast Fourier transform is applied at one point of the algorithm. The moving window width in seconds will be T = Nw/Fm, where Fm is the measurement sampling frequency. E.g. in this research we use T = 4.57s = 128/28Hz, Fm = 28 Hz. The moving window width is set to locate the heartbeat frequency in the stop band (about 4–5 heartbeats per window).
Here we use the procedure of digital lock-in filtering adapted to our needs. The amplitude I of locked-in light intensities I1 and I2 as measured for the two pairs (Fig. 1) is calculated as:
where the non-shifted component isHow to apply the Modified Beer-Lambert law
The locked-in signal I calculated using Eq. (1) cannot directly replace the light intensity when the modified Beer-Lambert law is used. The signal I can be applied directly into the modified Beer-Lambert law equation if the mean path length of light is doubled in the equation. The signal I is built as lock-in amplifications of two signals I1 and I2. Therefore, it relates to photons that travel twice through the medium. Let us introduce the mathematical formalism.
First, we narrow the problem to a single frequency component, where signals as measured with both reference arms are pure sinusoids having amplitudes I1 and I2 and the phase shift between them φ. In the classical homodyne detection approach, where the reference signal is a cosine at a frequency f with unit amplitude (I2 = 1) and zero phase shift (φ = 0), Eq. (2) gives I0 = I1 and Eq. (3) results in I90 = 0. In our case, the amplitude I2 is not the unity and φ might not be 0. Therefore, Eq. (2) and (3) yield I0 = I1I2cos(φ) and I90 = I1I2cos(φ + π/2). Thus, Eq. (1) gives $I = \sqrt {I_0^2 + I_{90}^2} = {I_1}{I_2}$. A change in the attenuation ΔA of the locked-in signal induced by e.g. a brain stimulation episode (measured as I) referred to a rest/reference period (measured as Iref) can be expressed as follows:
Equation (4) shows that the locked-in signal attenuation represents sum of attenuations at both reference arms. As a rule of thumb, amplitudes of the light intensities I1 and I2 measured at both arms will attenuate due to the brain activation almost equally. As such, ln(I1,ref/I1) ≈ ln(I2,ref/I2). Therefore, we can further approximate the attenuation change ΔA of the locked-in signal as:
The modified Beer-Lambert law bonds the attenuation change of a measured light intensity I1 with tissue absorption change Δμa as follows ln(I1,ref/I1) = Δμa$\langle$l$\rangle$, where $\langle$l$\rangle$ is the mean path of the light diffusively traveling through a tissue. Therefore, the modified Beer-Lambert law applied to the lock-in of signals as measured at both reference arms is:
2.3 Algorithm implementation
We provide an example of MATLAB implementation in the Algorithm 1 table. The lock-in filtering operation can be applied in a ‘sliding window’ manner, where the window of the Nw samples in width ‘slides’ through signals by one sample. For further convenience we share the full MATLAB implementation in Code 1 [23] and Code 2 [24], which uses an example of in-vivo data as measured in this research and available in Dataset 1 [25].
To serve the reader convenience the modified Beer-Lambert law data processing pipeline is:
2.4 Testing methods
FEM analysis
We use a simple two-layer head model as this geometry helps to investigate, separate and understand signal originating within extracerebral and brain tissues. The testing model is a slab of 120×120×50 mm (length, width and height) divided into 1,535,567 uniformly distributed tetrahedral elements (261,433 nodes) with an average volume of 0.45 ± 0.12 mm3 and 1.67 ± 0.29 mm internode distance. The extracerebral tissue thickness (grey matter surface location) is set to 13 mm following statistics carried out on 24 structural images registered with the magnetic resonance imaging technique as shown in [26]. The locked-in source-detector pairs are placed at the model centre on the top side where the inter-optode distance can be 15 mm or 30 mm as shown in Fig. 2. The exact fixing/positioning of the source-detector pairs on the surface is not relevant here as the medium properties change only in depth (layers). However, we keep the 90-degree rotated source-detector pairs organization for a better visualization.
Fig. 2. FEM model used for testing the lock-in signal processing principle (a). Combinations of signals to lock-in (b). ‘long’ refers to 30 mm source-detector separation and ‘short’ is 15 mm. Superficial layer thickness is 13 mm.
The brain tissue optical parameters were set to match haemoglobin changes that can be observed at the visual stimulation [27] oscillating around the rest state absolute haemoglobin concentrations of CHbO2 = 37 μM (oxy-haemoglobin HbO2) and CHb = 27 μM (deoxy-haemoglobin Hb). The extracerebral tissue haemoglobin content fluctuation is considered here as a biological noise fluctuating around the extracerebral layer rest state absolute values of CHbO2 = 33 μM and CHb = 24 μM. All rest state values of haemoglobin concentrations were set following [28,29]. The extracerebral tissue biological noise was measured in-vivo at a short (15 mm) source-detector separation. A representative signal was measured on a healthy volunteer during visual stimulation and it was observed, signals measured at the 15 mm separation do not show visual stimulation components. Final model parameters are shown in Fig. 3.
Fig. 3. Finite elements 2-layer head model haemodynamic parameters. Biological noise in the extracerebral tissue as measured in-vivo in this study at 15 mm source-detector separation.
We use common NIRS wavelengths of 750 nm and 850 nm to match our in-vivo studies. Following [30] the water content is set to 78% within both layers and the reduced scattering is approximated according to the Mie scattering theory μ′s = aλ-b, where light wavelength λ is expressed in μm. The scattering amplitude a (related to scattering objects concentration) and power b (related to their size) is set to a = 0.654 /mm−1 and b = 0.926 /-. That gives the reduced scattering coefficient at 750 nm and 850 nm as μ′s = 0.85 mm−1 and μ′s = 0.76 mm−1 respectively.
Time courses of intensities at both wavelengths are calculated using the NIRFASTer [31] software which is a GPU accelerated version of the NIRFAST [32] package designed to solve the diffusion equation using the finite element method.
In-vivo visual cortex stimulation
The in-vivo data were collected on six healthy volunteers following protocols adopted from our previous study [27]. Please refer to this research for details on subject position, head-screen distance, etc. A flickering (at 8 Hz) black and white, round checkerboard was displayed on a computer screen to stimulate the visual cortex. The protocol was set as: ≥30 s of rest followed by 10 epochs of 40 s, built as 20 s of stimulation and 20 s of rest. The visual stimulation is synchronized with the fNIRS measurement through an external PC registering sync signals under the same sampling clock with a DAQ and a LabVIEW software.
Sources and detectors were arranged within a mesh fixed to a head surface as shown in Fig. 4. Optodes were fixed over the left side of the visual cortex of a healthy volunteer. There is a 4–5 mm free space between the optodes holder surface and the fibre tips. This space will be filled with hair. The fibre tips create a ‘hair brush’ that can brush its way through the hair to the skin. The reasoning for the fixing system size and location was to maximize number of source-detector pairs registering visual stimulation signals. Optical bundles of 2.5 mm active diameter were used. Source points combine two wavelengths of 750 nm and 850 nm using a bifurcated bundle of optical fibres.
Fig. 4. Sources and detectors arranged within a mesh (a) fixed on a head surface (b) over the left side of visual cortex of a healthy volunteer. Panel (a) shows all possible combinations of pairs.
The high-density diffuse optical tomography measurement system is in-house developed and capable of high dynamic range measurements (up to 7 orders of magnitude). This is enabled by in-house developed avalanche photodiode detection and LED light source modules. The system is briefly shown in [33]. It uses a frequency and temporal encoding strategies to simultaneously detect data at up to 96 sources × 96 detectors × 2 wavelengths. The encoding and decoding are carried out internally using an in-house-developed FPGA (field-programmable gate array) module. For the source-detector arrangement as in Fig. 4 only a portion of the system capabilities is used to register just 64 signals (4 sources × 8 detectors × 2 wavelengths). The measurement frequency for the setup as in Fig. 4 (all pairs) was set to 28 Hz. Data are transferred in real-time through the Ethernet interface (TCP/UDP protocol) to an in-house-developed dynamic library linked to a convenient MATLAB routine (the *.mex file) providing a MATLAB plugin for easy, real-time visualization, processing and data storage.
In-vivo motion artefacts testing
One of the healthy subjects participated in motion artefacts test. Fibres were fixed on the left side of the forehead. The subject was responsible to hold/support all fibres bundled using his right hand. The protocol was resting – walking – resting – jumping – resting with the following timing: resting – sitting for 60 s, walking for 5 s and sitting for the next 25 s, jumping for 5 s followed by resting – sitting for 25 s.
Lock-in amplification strategies
The proposed lock-in fNIRS method makes use of two source-detector pairs at a large distance with the common sensitivity region overlapping within the brain cortex. However, for testing purposes, we use additional three lock-in filtering patterns: both pairs are at the large distance (30 mm) but there is no sensitivity region overlapping, two collinear pairs – one at the large distance and the second one at the short one (15 mm) and finally both pairs are at the short distance. All four scenarios are summed up in Fig. 5. The colinear long-short pair at Fig. 5(c) corresponds to the commonly used concept of NIRS signals collection at long and short distances.
Fig. 5. Data lock-in amplification strategies: (a) the lock-in method – both pairs at large distance (30 mm) with overlapping sensitivity region, (b) both pairs at large distance but no overlapping sensitivity region, (c) collinear pairs – one pair is at large distance and the second one is at the short one (15 mm), (d) both pairs are at the short distance.
The procedure as in the Algorithm 1 table is applied to all lock-in amplification strategies as shown in Fig. 2 and Fig. 5. The lock-in filtering operation is applied in a sliding-window manner. The 128 sample-long window, used at the 28 Hz system measurement acquisition frequency, gives one lock-in result from a 4.6 s long signal. Haemoglobin concentration changes were calculated using the modified Beer-Lambert law at 2 wavelengths (750 nm and 850 nm). Assuming extinction coefficients of haemoglobin as compiled by Scott Prahl and available at [34]. The mean optical path lengths are taken from simulations carried out on the physical phantom as shown in Fig. 2.
3. Results
3.1 FEM analysis
Figure 6 shows time courses of data simulated and recovered for combinations of overlapping and non-overlapping measurements at 30 mm and 15 mm source-detector distances as proposed in Fig. 2. Please refer to the ‘Testing methods’ section and Fig. 3 for details on how the data was generated. To quantify the recovery of the haemodynamic response, we have calculated the fitting error for the amplitude and the offset as the fitted function is periodic. The amplitude fitting error is expressed in percent. The offset fit error is expressed in μM and represents difference between mean values of the recovery and ground truth.
Fig. 6. Time courses of ground truth changes in haemoglobin concentrations (thick lines) and results of recovery for overlapping and non-overlapping geometries. The top row of the time series shows results of the typical modified Beer-Lambert law analysis at long (30 mm) and short (15 mm) source-detector distances. The light attenuation panel shows the entry data for the analysis. The amp. error represents the amplitude fit error in percent and the off. error shows error in offset between both periodic signals. Data for the short distance does not provide the brain haemodynamic response and the recovered haemoglobin changes approach zero. The bottom row of panels shows combinations of lock-in filtered data as proposed in Fig. 2(b).
As shown in Fig. 6, the lock-in approach fits best to the ground truth. The offset fit drops significantly for the lock-in operation. For sake of comparison, the modified Beer-Lambert law analysis for single pairs at long (30 mm) and short (15 mm) source-detector distances were filtered with a moving average filter of 128-samples window (4.6 s) to match the ‘cut-off’ frequency of the lock-in filter. The lock-in operation on signals at long and short source-detector separations (with and without the brain activation components) allows to detect the brain activation and remove the DC drift. Furthermore, this lock-in combination is still beneficial as the brain activation is recovered. Finally, data for the short distance does not provide the brain haemodynamic response and the recovered haemoglobin changes approach zero.
The frequency analyse of the temporal data is shown in Fig. 7. Haemoglobin signals decomposed into frequency components reveal the nature of the lock-in operation.
Fig. 7. The frequency analysis of simulated temporal data from Fig. 6. Comparison of the lock-in and the classical data processing approach. Please see the removal of the DC component and frequencies around the visual stimulation using the lock-in filter.
The frequency analysis shows that the lock-in filter attenuated the DC component (drifts) and frequencies around the visual stimulation (see Fig. 6). The lock-in filter preserves (amplifies) the biological components located within the stop-band. The moving average filtering attenuates all. The lock-in operation preserved the respiratory rate and heart rate, despite the 4.6 s long filtering window. This is not possible with a common moving average filtering.
The proposed approach works within a moving window that differs in size by about two orders of magnitude as compared to the classical homodyne detection. The classical lock-in needs 10 or more periods of a frequency component to detect its amplitude. It is a single number taken from the 10 epochs. The proposed method works within the window size of a fraction of the component period. E.g. let us consider the sampling frequency of Fm = 28 Hz, brain cortex activation epoch of 40 s (1/40 = 0.025 Hz) and the proposed window width of Nw = 128 samples. The Nw samples give the window length of about 11% of the brain activation epoch at which we want to increase the contrast to noise ratio. The classical homodyne detection (lock-in) would need around 400 seconds of the signal (10 epochs of the brain activation). The proposed method is not designed to directly detect amplitudes of desired components. The signal I as in Eq. (1) oscillates with frequency components common for both light intensities I1 and I2. These components are amplified. Finally, the method attenuates (but not rejects) all frequencies higher than the window reciprocal (1/T ≈ 0.219 Hz).
3.2 In-vivo visual cortex stimulation
Representative time courses of data registered in-vivo in two subjects are shown in Fig. 8,9. Figure 8 represents an example of a high-quality data in which the changes related to visual stimulation are clearly visible. Figure 9 shows an example of a less pronounced visual stimulation response. Data as registered on other subjects represent the high-quality data example. The signals of light intensity changes registered at all possible source-detector combinations were of good quality with the coefficient of variation at the largest source-detector distance (34 mm) calculated as <2% at both wavelengths. The variation is mainly related to the pulse wave. We show results as registered for the four signal processing scenarios as shown in Fig. 5. As a reference, source-detector pairs #1 and #2 used for the lock-in procedure were processed separately in the classical manner. Furthermore, a zero-phase moving average filter of the same window size as the moving lock-in window (128 samples equal 4.57s) was applied. This approach allowed to show the difference between the proposed method and one of the simplest/robust filtering technique.
Fig. 8. Good quality data example. Time courses of changes in concentrations of oxy- (red lines) and deoxy-haemoglobin (blue lines) registered in-vivo in a healthy volunteer under visual stimulation for varying lock-in signal strategies (see Fig. 4 and 5). Results for pair #1 and #2 show the average filtering (thick lines) overlaying direct photodetector signals (thin black/grey lines). For high-resolution see the supplemental material.
Fig. 9. Example of a less pronounce visual response. Time courses of changes in concentrations of oxy- (red lines) and deoxy-haemoglobin (blue lines) registered in-vivo in a healthy volunteer under visual stimulation for varying lock-in signal strategies (see Fig. 4 and 5). Results for pair #1 and #2 show the average filtering (thick lines) overlaying direct photodetector signals (thin black/grey lines). For high-resolution see the supplemental material.
Results as in Fig. 8,9 show that the lock-in filtering approach in the experiment conditions is always beneficial, regardless of the spatial location of interfered signals combination used. Under the condition that at least one interfering arm registers the brain haemodynamic response. Furthermore, the real benefit of the interference approach reveals in the noisy input data scenario (Fig. 9). Additionally, the lock-in filtering preserves the respiratory rate and pulse components present in both interfering arms. Which is not achievable for the averaging method. This is shown with the frequency analysis in Fig. 10. Interfering signals at short distances (with no brain haemodynamic response) is neutral. It filters data preserving the biological components, e.g. Mayer waves, respiratory rate or pulse.
Fig. 10. An example of the lock-in method frequency characteristics as compared to a common moving average filtering. The moving average filtering is a zero-phase filter (MATLAB ‘filtfilt’ function). Time data shown in the first column in Fig. 8 (quality data) were used. The oxygenated haemoglobin spectrum is shown.
Mathematically, the proposed lock-in fNIRS method is a homodyne (wide-band homodyne) detection followed by a low-pass filter with the cut-off frequency equal to the reciprocal of the lock-in window width. An example of the filter characteristics is shown in Fig. 10.
The lock-in procedure preserves biological data even after the low-pass filtering (the sliding window size in Fig. 10). The simple averaging low-pass filtering cuts out all in the stop-band and does not change the components in the pass-band. The proposed lock-in method has a zero-phase shift character as shown in the phase characteristics (angle of the complex transfer function of the filter). This preserves temporal data shape and time relations as the frequency components phase shift in the pass-band is 0 or 2π radians (the filter phase in Fig. 10). The phase shift within the stop-band becomes random. However, for the biological components preserved within the stop-band by the lock-in operation, the phase shift of the respiratory rate and pulse approaches π.
The next testing step it to increase the lock-in filter moving window size to span on over 2 visual stimuli epochs. Thus, the visual stimulation frequency will be located within the lock-in filter stopband. The fast Fourier transform is used to calculate orthogonal signals as in Eq. (3). Therefore, the window width in samples is not required to be the power of 2 but makes the algorithm more efficient in terms of execution speed. The lock-in procedure selectively amplifies all common frequency components (drifts, brain response, Mayer waves, pulse, etc.) in the entire spectrum. The amplification is followed by the attenuation in the stopband (frequencies above the reciprocal of the filter window width in seconds). However, frequencies of e.g. respiration and heartbeat are attenuated but preserved and can be still detectable. Such approach might be of interest in high noise data processing where stimulation was repeated numerous times. Conventional filtering removes all components in the stop-band as shown in Fig. 11.
Fig. 11. The lock-in method with large lock-in sliding window size. Data equal to the first column in Fig. 8 (quality data).
We define a local parameter, contrast to noise ratio (CNR). The brain haemodynamic response (peak at the visual stimulation frequency fstim = 1/40 = 0.025 Hz) related to the noise floor around the brain response. Frequency spectra of haemoglobin were calculated from 400 s of signals registered during the 10 epochs of the visual stimulation. The CNR gives a measure of the brain stimulation ‘detectability’. Here, the CNR noise floor is calculated within the span of 0.5fstim = 0.0125 Hz and 2fstim = 0.05 Hz. In opposite to the signal to noise ratio, where the noise floor is global usually. The biological noise (Mayer waves – 0.1 Hz, pulse – 1 Hz, respiratory rate – 0.25 Hz) is outside of the CNR noise floor calculation range. The CNR for six subjects and four 90-degree shifted lock-in pairs per subject (a total of 24 measurements) is shown in Fig. 12. Results are sorted by the CNR in spectrum of the oxygenated haemoglobin for the lock-in method.
Fig. 12. The CNR in oxygenated (a) and deoxygenated (b) haemoglobin for six subjects and four 90-degree shifted lock-in pairs per subject (a total of 24 measurements). Results are sorted by the CNR of the oxygenated haemoglobin for the lock-in method. Measurements where the lock-in method provides the CNR gain are highlighted with circle markers (big green and smaller orange). Three quality zones in CNR are marked with colours. The bottom one (CNR < 2) marks measurements with non-detectable brain haemodynamic response to the visual stimuli.
The CNR in Fig. 12 is divided into three quality groups, where the CNR thresholds are 2 and 4. The brain haemodynamic response to the visual stimuli is non-detectable for the CNR < 2. Further, measurements where the lock-in procedure reveals CNR improvement are highlighted with round markers: big green when the lock-in has the best CNR and smaller orange when the lock-in approach has better CNR than the mean of the CNR for both reference arms processed with the average filter. As shown, the lock-in procedure performs best for the ‘low quality’ data where the CNR is between 2 and 4. An interesting pattern for the deoxygenated haemoglobin can be observed in Fig. 12(b). The CNR in deoxy-haemoglobin drops significantly if the CNR in oxygenated haemoglobin approaches the ‘low quality’ region (CNR < 4). Further, if the visual stimulation is visible in deoxygenated haemoglobin signal, the CNR is usually high (between 5 and 8).
The motion artefact testing results are shown in Fig. 13. All locked-in pairs show the same results. Therefore, a single pair (with the strongest artefacts registered) is shown for a clear visual presentation.
Fig. 13. Results of the in-vivo motion artefacts testing. First 40 s of flat rest signal not shown to zoom in into the motion artefact region.
As shown in Fig. 13, the proposed lock-in approach successfully filters the ‘non-periodic’ motion artefacts. The biological noise (such as the pulsatile component) is attenuated but preserved. The average filtering rejects pulsation and ‘over-smoothes’ the time course losing the details. The same window size of 1.14 s (32 samples @ 28 Hz sampling rate) was used for the lock-in and average filtering.
4. Discussion
The conventional lock-in filtering is based on the homodyne signal mixing i.e. detection at a single frequency is employed. Here we mix two biological signals, similar to each other. The reference (oscillator) signal has many frequency components. Therefore, the homodyne detection converts to a broad-band homodyne i.e. homodyne detection at multi-frequencies present in the reference (oscillator) signal. The homodyne mixing step is followed by a low-pass filtering with the cut-off frequency linked to the lock-in window size sliding through the analysed signal. The lock-in procedure preserves components present in both interfering signals at the same time. Further, the components close in frequency to the visual stimulation are attenuated, increasing the contrast to noise ratio. Components within the stop-band like biological signals but present in one reference arm only (instrumental noise or movement artefacts) are attenuated. However, biological signals within the stop-band present in both reference arms (e.g. respiratory and pulse) are preserved. If the brain-origin activation component (in the pass band) is present in one of the two signals only, the lock-in fNIRS is still advantageous by applying the low-pass filtering. However, the resulting amplitude of the visual response will be slightly attenuated.
The proposed lock-in fNIRS as introduced here is not designed to remove the extracerebral components but to amplify brain components and increase the contrast to noise ratio. This technique leads to increase in probability of detecting a brain activity. Furthermore, the lock-in fNIRS signal can be further processed with techniques of extracerebral information removal, if desired. A comprehensive guide on such methods is compiled in [35] where authors show quantitative comparison of a number of methods e.g. spatial and temporal filtering, regression, component analysis or short-separation regression. Model of the extracerebral component can be built of a signal at single short source-detector distance (e.g. 15 mm), an average of signals at such separations or an average biased by a presence of Mayer wave oscillations [36]. In the proposed lock-in fNIRS technique, the extracerebral component model can be set to the lock-in signal for two short separations (e.g. 15 mm and 15 mm).
Additionally, the proposed lock-in fNIRS is recognized as extendable to post-processing analysis in two ways. First, signals (measured intensities or even calculated haemoglobins) can be locked-in with the stimulation model (sinus of the main frequency of a stimulation, e.g. 1/40 Hz in our in-vivo experiment or even the paradigm pattern). This will result in the amplitudes of changes in haemoglobin concentrations related to the brain stimulation expressed in μM and the phase shift. The shift will give temporal lag between the visual stimulation and brain response that can be expressed in absolute time units. This might be of interest if comparison of fNIRS lock-in signals at different spatial locations is performed.
The lock-in fNIRS is not limited to continuous wave measurements of light intensity. It can be directly transferred to other NIRS modalities like frequency or time domains. The lock-in procedure can be applied to the light wave phase shift or parameters of the measured distribution of time of flight of photons, e.g. mean time of flight, variance or intensity in time windows of the distribution curve as measured in the proposed here x-shaped geometry.
Using the proposed x-shaped measurement geometry, the two source-detector pairs employed in the lock-in fNIRS are at the same distance and consequently at similar electronics’ gain and noise levels. Therefore, there is no considerable issue on the detectors dynamic range as it can be an issue for a long and short source-detector pair combination.
The broad-band homodyne detection as presented in Algorithm 1 table can be supplement by calculation of the phase shift between both locked-in signals as: φ = arctan(I90/I0). This phase should be irrelevant as signals are similar. In our in-vivo data φ is always stable, below 0.01 deg and oscillates around 0 or a bias. Sometimes, there is a constant offset in the phase between used wavelengths of 750 nm and 850 nm. This shift does not impact the fNIRS results. However, this effect might need further investigation and it is expected that this is hardware/algorithm implementation specific. No visual stimulation traces were found in temporal fluctuations of all observed phase shifts. However, the Mayer waves, respiration and pulse are still visible in the power density spectrum of φ. We predict that NIRS source-detector pairs fixed around a human body and locked-in as proposed here (Algorithm 1) can trace spatial/temporal propagation of Mayer or pulse waves if amplitudes of the desired frequency components of φ are traced. Therefore, the lock-in NIRS technique can measure velocity of such waves that can be further correlated with the blood pressure as it was carried out in arteries [37] or serve as a tool to assess cardiovascular abnormalities [38].
We have investigated several approaches to the problem of building a dense mesh of sources and detectors paired in the proposed lock-in fNIRS regime that enables a three-dimensional tomography reconstruction as the HD-DOT does. A future work in this field is needed to explore opportunities for incorporating this algorithm to the high-density grid geometries as the combination may provide enhanced spatial/depth specificity. As for now, a mesh of 90-degree rotated pairs of sources and detectors working in the lock-in measurement regime fits well to the topography brain imaging approach where measured values can be spatially assigned to the brain surface. Please follow the Appendix on theory how to apply the proposed method to the diffuse optical tomography reconstruction.
5. Conclusions
We show that a brain originating component present in an optical signal can be amplified (and contrast to noise increased) with the lock-in (broad-band homodyne detection) method where spatial sensitivity distributions of interfering signals are common in the brain region of interest.
The lock-in fNIRS method can be implemented with common photodetectors and LEDs as there is no demanding dynamic range requirements and by default, signals undergo homodyne (broad-band homodyne) detection followed by low-pass filtering. As the lock-in detection relies on two biological signals, the lock-in filter characteristic will vary constantly adapting to current brain activation strength, biological noise, movement artefacts, etc.. The algorithm implementation and complexity are suitable for low resource implementations. Combined hardware and software features form a real-time, direct, robust and low-noise brain activity registration tool. This fits well into brain computer interface applications (including mental communication [39] or entertainment, computer gaming, etc.) or all sorts of fNIRS monitoring devices. We show a potential path to make a (f)NIRS device simple and reliable with a potential for robust implementation in a home use environment or entertainment like a bio-feedback in computer games.
Appendix
How to apply the diffuse tomography reconstruction
In the diffuse tomography reconstruction approach, the equation of light diffuse transport is linearized and inverted using the perturbation theory. This operation generates a set of linear equations U = -Jδμa that links changes in M measurements UM×1 carried out on e.g. a tissue surface induced by small changes in absorption $\mathrm{\delta }\boldsymbol{\mathrm{\mu}}_\textrm{a}^{N \times 1}$ within the tissue divided into N volumes (layers, voxels, finite element mesh nodes, etc.). When the Rytov approach is used in the linear perturbation, a change in the m-th measurement is set to Um = ln(Um/Uref,m) where ‘ref’ indicates value before the absorption change. Please refer to extensive sources of knowledge on the optical diffuse tomography in general [40] and a comprehensive report on how to derive the mathematical formalism of the optical tomography using analytical forms or the finite element model [41,42].
If U is considered in terms of the attenuation A of the light radiant flux Φ (the intensity/power having the W unit), A = ln(Φref/Φ) and the tomography problem converts to A = Jδμa. Please notice the cancelled out minus sign as the reference signal in the attenuation is in the nominator position and the Rytov formula uses reference in the denominator. The Jacobian matrix JM×N represents mean optical path lengths (expressed in the unit of length) given for M measurements and within N tissue volumes. The J matrix can be calculated numerically (or even measured [43]) element-wise as follows: Jm,n = ∂Am/∂μa,n (change in attenuation at m-th measurement induced by a small change in absorption within the n-th tissue element) or imaged with a time-gated CCD (charge-coupled device) camera using the adjoint fields method [44].
Now, let us consider the Hadamard product JH = J1∘J2 as in Fig. 1(d). An element of the JH matrix indexed by m and n equals to:
For the symmetrical geometry of pair #1 and pair #2 as in Fig. 1 the following condition is met:
Entries of the Jacobian matrices have the logarithmic nature and decay with the distance. The inverse problem of diffuse optical tomography uses the Tikhonov regularization to calculate a pseudoinverse of a Jacobian matrix [45,46]. As a rule of thumb, the regularization parameters discard the Jacobian entries below 10% of its maximum value. Values below the 10% introduce reconstruction noise. Figure 1 colour scale is set between the max of Jacobians and the 6% of the max. Therefore, considering Eq. (8) the attenuation changes at both arms induced within the nonzero JH elements are close to each other ∂A1,m ≈ ∂A2,m.
Equation (7) yields:
And finally in the matrix form, in terms of measurements at both arms (pair #1 and pair #2):
The left-hand side of Eq. (10) and Eq. (9) reduces to the attenuation as A1 ≈ A2 and the square root of the JH goes towards 0 if J1 and J2 do not intersect as in Fig. 1(d). This gives a substantial benefit as the sensitivity in the superficial layer (close to sources and detectors) cancels out.
Funding
National Science Centre Poland (2012/05/B/ST7/01162, 2016/21/N/ST7/03117, 2020/39/D/ST7/03425); Nalecz Institute of Biocybernetics and Biomedical Engineering (ST213/2021).
Disclosures
The authors declare no conflicts of interest related to this article.
Data availability
A MATLAB script with the implementation of the proposed method and a sample of in-vivo collected data to run it are available in Code 1 [23] and Code 2 [24] and Dataset 1 [25].
Supplemental document
See Supplement 1 for supporting content.
References
1. X. Xu, H. Liu, and L. V. Wang, “Time-reversed ultrasonically encoded optical focusing into scattering media,” Nat. Photonics 5(3), 154–157 (2011). [CrossRef]
2. M. Caldwell, F. Scholkmann, U. Wolf, M. Wolf, C. Elwell, and I. Tachtsidis, “Modelling confounding effects from extracerebral contamination and systemic factors on functional near-infrared spectroscopy,” NeuroImage 143, 91–105 (2016). [CrossRef]
3. M. D. Pfeifer, F. Scholkmann, and R. Labruyère, “Signal Processing in Functional Near-Infrared Spectroscopy (fNIRS): Methodological Differences Lead to Different Statistical Results,” Front. Hum. Neurosci. 11, 641 (2018). [CrossRef]
4. D. Milej, A. Abdalmalak, A. Rajaram, and K. St. Lawrence, “Direct assessment of extracerebral signal contamination on optical measurements of cerebral blood flow, oxygenation, and metabolism,” Neurophotonics 7(04), 045002 (2020). [CrossRef]
5. I. Tachtsidis and F. Scholkmann, “False positives and false negatives in functional near-infrared spectroscopy: issues, challenges, and the way forward,” Neurophotonics 3(3), 031405 (2016). [CrossRef]
6. R. B. Saager and A. J. Berger, “Direct characterization and removal of interfering absorption trends in two-layer turbid media,” J. Opt. Soc. Am. A 22(9), 1874–1882 (2005). [CrossRef]
7. K. M. Bergonzi, T. M. Burns-Yocum, J. R. Bumstead, E. M. Buckley, P. C. Mannion, C. H. Tracy, E. Mennerick, S. L. Ferradal, H. Dehghani, A. T. Eggebrecht, and J. P. Culver, “Lightweight sCMOS-based high-density diffuse optical tomography,” Neurophotonics 5(03), 1 (2018). [CrossRef]
8. H. Dehghani, B. R. White, B. W. Zeff, A. Tizzard, and J. P. Culver, “Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography,” Appl. Opt. 48(10), D137–143 (2009). [CrossRef]
9. A. T. Eggebrecht, S. L. Ferradal, A. Robichaux-Viehoever, M. S. Hassanpour, H. Dehghani, A. Z. Snyder, T. Hershey, and J. P. Culver, “Mapping distributed brain function and networks with diffuse optical tomography,” Nat. Photonics 8(6), 448–454 (2014). [CrossRef]
10. E. Vidal-Rosas, H. Zhao, R. Nixon-Hill, G. Smith, L. Dunne, S. Powell, R. Cooper, and N. Everdell, “Evaluating a new generation of wearable high-density diffuse optical tomography technology via retinotopic mapping of the adult visual cortex,” Neurophotonics 8(02), 025002 (2021). [CrossRef]
11. A. Torricelli, D. Contini, A. Pifferi, M. Caffini, R. Re, L. Zucchelli, and L. Spinelli, “Time domain functional NIRS imaging for human brain mapping,” NeuroImage 85, 28–50 (2014). [CrossRef]
12. J. Selb, J. Stott, M. Franceschini, A. G. Sorensen, and D. Boas, “Improved sensitivity to cerebral hemodynamics during brain activation with a time-gated optical system: analytical model and experimental validation,” J. Biomed. Opt. 10(1), 011013 (2005). [CrossRef]
13. D. Orive-Miguel, L. Hervé, L. Condat, and J. Mars, “Improving Localization of Deep Inclusions in Time-Resolved Diffuse Optical Tomography,” Appl. Sci. 9(24), 5468 (2019). [CrossRef]
14. D. Milej, A. Abdalmalak, A. Rajaram, A. Jhajj, A. Owen, and K. St. Lawrence, “Incorporating early and late-arriving photons to improve the reconstruction of cerebral hemodynamic responses acquired by time-resolved near-infrared spectroscopy,” J. Biomed. Opt. 26(05), 056003 (2021). [CrossRef]
15. A. Liebert, H. Wabnitz, J. Steinbrink, H. Obrig, M. Moller, R. Macdonald, A. Villringer, and H. Rinneberg, “Time-resolved multidistance NIRS of the adult head: intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons,” Appl. Opt. 43(15), 3037–3047 (2004). [CrossRef]
16. H. Wabnitz, D. Contini, L. Spinelli, A. Torricelli, and A. Liebert, “Depth-selective data analysis for time-domain fNIRS: moments vs. time windows,” Biomed. Opt. Express 11(8), 4224–4243 (2020). [CrossRef]
17. A. Sudakou, L. Yang, H. Wabnitz, S. Wojtkiewicz, and A. Liebert, “Performance of measurands in time-domain optical brain imaging: depth selectivity versus contrast-to-noise ratio,” Biomed. Opt. Express 11(8), 4348–4365 (2020). [CrossRef]
18. P. Sawosz and A. Liebert, “Method to improve the depth sensitivity of diffuse reflectance measurements to absorption changes in optically turbid medium,” Biomed. Opt. Express 10(10), 5031–5041 (2019). [CrossRef]
19. S. Fantini, G. Blaney, and A. Sassaroli, “Transformational change in the field of diffuse optics: From going bananas to going nuts,” J. Innovative Opt. Health Sci. 13(01), 1930013 (2020). [CrossRef]
20. G. Blaney, A. Sassaroli, and S. Fantini, “Design of a source–detector array for dual-slope diffuse optical imaging,” Rev. Sci. Instrum. 91(9), 093702 (2020). [CrossRef]
21. G. Blaney, A. Sassaroli, and S. Fantini, “Dual-slope imaging in highly scattering media with frequency-domain near-infrared spectroscopy,” Opt. Lett. 45(16), 4464–4467 (2020). [CrossRef]
22. Y. Liu, P. Lai, C. Ma, X. Xu, A. A. Grabar, and L. V. Wang, “Optical focusing deep inside dynamic scattering media with near-infrared time-reversed ultrasonically encoded (TRUE) light,” Nat. Commun. 6(1), 5904 (2015). [CrossRef]
23. S. Wojtkiewicz, “MATLAB implementation of the lock-in functional NIRS - main script,” figshare (2022), https://doi.org/10.6084/m9.figshare.19166456
24. S. Wojtkiewicz, “MATLAB implementation of the lock-in functional NIRS - support function,” figshare (2022), https://doi.org/10.6084/m9.figshare.19166462
25. S. Wojtkiewicz and K. Bejm, “Dataset for the MATLAB implementation of the lock-in functional NIRS - in-vivo data measured under visual stimulation condition,” figshare (2022), https://doi.org/10.6084/m9.figshare.19166468
26. M. Doulgerakis, A. T. Eggebrecht, and H. Dehghani, “High-density functional diffuse optical tomography based on frequency-domain measurements improves image quality and spatial resolution,” Neurophotonics 6(03), 1–14 (2019). [CrossRef]
27. K. Bejm, S. Wojtkiewicz, P. Sawosz, M. Perdziak, Z. Pastuszak, A. Sudakou, P. Guchek, and A. Liebert, “Influence of contrast-reversing frequency on the amplitude and spatial distribution of visual cortex hemodynamic responses,” Biomed. Opt. Express 10(12), 6296–6312 (2019). [CrossRef]
28. H. Auger, L. Bherer, E. Boucher, R. Hoge, F. Lesage, and M. Dehaes, “Quantification of extra-cerebral and cerebral hemoglobin concentrations during physical exercise using time-domain near infrared spectroscopy,” Biomed. Opt. Express 7(10), 3826–3842 (2016). [CrossRef]
29. G. Giacalone, M. Zanoletti, D. Contini, R. Re, L. Spinelli, L. Roveri, and A. Torricelli, “Cerebral time domain-NIRS: reproducibility analysis, optical properties, hemoglobin species and tissue oxygen saturation in a cohort of adult subjects,” Biomed. Opt. Express 8(11), 4987–5000 (2017). [CrossRef]
30. A. Torricelli, A. Pifferi, P. Taroni, E. Giambattistelli, and R. Cubeddu, “In vivooptical characterization of human tissues from 610 to 1010 nm by time-resolved reflectance spectroscopy,” Phys. Med. Biol 46(8), 2227–2237 (2001). [CrossRef]
31. .“GitHub - nirfaster/NIRFASTer: Open source software for multi-modal optical molecular imaging” (2020), retrieved https://github.com/nirfaster/NIRFASTer.
32. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Meth. Engng. 25(6), 711–732 (2009). [CrossRef]
33. S. Wojtkiewicz, P. Sawosz, M. Kacprzak, A. Gerega, K. Bejm, R. Maniewski, and A. Liebert, “Towards optical tomography of an adult human head,” in Biomedical Optics 2016, OSA Technical Digest (online) (Optical Society of America, 2016), OM4C.2.
34. S. Prahl, “Tabulated Molar Extinction Coefficient for Hemoglobin in Water” (2020), retrieved https://omlc.org/spectra/hemoglobin/summary.html.
35. H. Santosa, X. Zhai, F. Fishburn, P. Sparto, and T. Huppert, “Quantitative comparison of correction techniques for removing systemic physiological signal in functional near-infrared spectroscopy studies,” Neurophotonics 7(03), 035009 (2020). [CrossRef]
36. D. Wyser, M. Mattille, M. Wolf, O. Lambercy, F. Scholkmann, and R. Gassert, “Short-channel regression in functional near-infrared spectroscopy is more effective when considering heterogeneous scalp hemodynamics,” Neurophotonics 7(03), 035011 (2020). [CrossRef]
37. Y. Ma, J. Choi, A. Hourlier-Fargette, Y. Xue, H. U. Chung, J. Y. Lee, X. Wang, Z. Xie, D. Kang, H. Wang, S. Han, S.-K. Kang, Y. Kang, X. Yu, M. J. Slepian, M. S. Raj, J. B. Model, X. Feng, R. Ghaffari, J. A. Rogers, and Y. Huang, “Relation between blood pressure and pulse wave velocity for human arteries,” Proc. Natl. Acad. Sci. 115(44), 11144–11149 (2018). [CrossRef]
38. J. Poleszczuk, M. Debowska, W. Dabrowski, A. Wojcik-Zaluska, W. Zaluska, and J. Waniewski, “Subject-specific pulse wave propagation modeling: Towards enhancement of cardiovascular assessment methods,” PLoS ONE 13(1), e0190972 (2018). [CrossRef]
39. A. Abdalmalak, D. Milej, L. C. M. Yip, A. R. Khan, M. Diop, A. M. Owen, and K. St. Lawrence, “Assessing Time-Resolved fNIRS for Brain-Computer Interface Applications of Mental Communication,” Front. Neurosci. 14, 105 (2020). [CrossRef]
40. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010). [CrossRef]
41. S. R. Arridge, “Photon-measurement density functions. Part I : Analytical forms,” Appl. Opt. 34(31), 7395–7409 (1995). [CrossRef]
42. S. R. Arridge and M. Schweiger, “Photon-measurement density functions. Part 2: Finite-element-method calculations,” Appl. Opt. 34(34), 8026–8037 (1995). [CrossRef]
43. M. Kacprzak, A. Liebert, P. Sawosz, N. Zolek, and R. Maniewski, “Time-resolved optical imager for assessment of cerebral oxygenation,” J. Biomed. Opt. 12(3), 034019 (2007). [CrossRef]
44. P. Sawosz, M. Kacprzak, W. Weigl, A. Borowska-Solonynko, P. Krajewski, N. Zolek, B. Ciszek, R. Maniewski, and A. Liebert, “Experimental estimation of the photons visiting probability profiles in time-resolved diffuse reflectance measurement,” Phys. Med. Biol. 57(23), 7973–7981 (2012). [CrossRef]
45. W. Fan, H. Dehghani, and A. Eggebrecht, “Investigation of effect of modulation frequency on high-density diffuse optical tomography image quality,” Neurophotonics 8(04), 045002 (2021). [CrossRef]
46. G. Perkins, A. Eggebrecht, and H. Dehghani, “Quantitative evaluation of frequency domain measurements in high density diffuse optical tomography,” J. Biomed. Opt. 26(05), 056001 (2021). [CrossRef]
