1. Introduction

Near-infrared spectroscopy (NIRS) is used to noninvasively characterize the absorption and scattering of light in centimeter thick tissues. This information can be used to assess tissue composition (hemoglobin, water, lipid) [1,2] and has been applied in clinical and preclinical research for noninvasive brain imaging [3,4] neonatal care [5], breast tumor imaging [1,6], muscle training [7], and critical care [8]. To date, clinically established NIRS methods are based on continuous-wave (CW) (i.e., time-independent) optical illumination, however time-resolved DOS in both the time- (TD-NIRS) and frequency-domains (FD-NIRS) can provide significant benefits. CW-NIRS only assesses changes in light intensity and therefore requires making assumptions about optical scattering that can introduce significant errors [9]. This limits the quantitative accuracy of CW-NIRS measurements and complicates direct comparisons of tissue measurements between individuals and longitudinally within individuals over long periods of time.

Alternatively, time-dependent methods such as FD-NIRS can separate light absorption from scattering, provide more uniform and deeper tissue sampling [10], measure fluorescence lifetime [11], and are widely recognized as providing substantially greater information content, accuracy, and 3D resolution [2,10,12]. Their quantitative features allow for absolute comparisons between individuals and longitudinal monitoring of subjects with exceptional accuracy and precision over days to years [13,14]. Scattering parameters can also be isolated and used to provide additional contrast [15]. However, despite these advantages, development of clinically-compatible FD- and TD-NIRS instrumentation has been impeded by their poor scalability for high density, high resolution imaging and is relatively complex and slow operation compared to CW-NIRS.

In this work we describe a new FD-NIRS platform that aims to overcome prior limitations by providing high-speed data acquisition and display of multi-frequency (1-400 MHz) and multi-wavelength (690, 785, 808, 850, 940, and 980 nm) FD-NIRS derived optical properties. Similar to recent work that has demonstrated high-speed (>10 Hz) FD-NIRS data acquisition [16–20], we employ an all-digital detection and signal generation scheme with a high-speed analog-to-digital converter (ADC) and a direct digital synthesizer (DDS) [21]. However, high-speed ADCs have the downside of generating extremely large datasets (in our system around 5.9 GB per second) which is difficult to stream or display in real-time. Because of this, to the best of our knowledge, no FD-NIRS systems have demonstrated real-time spatial mapping of chromophores at display frame rates >60 Hz (currently a standard refresh rate in displays [22] and mobile devices [23]). This is important for providing immediate lag-free feedback to users, particularly in a clinical setting. To remedy the shortcoming, we built-in a full digital streaming pipeline and control system inside the hardware which enables real-time chromophore display by using a field programmable gate array (FPGA) and co-processor (CPU), significantly reducing the data streaming output to around 48 bytes per view (all six wavelengths). We developed an ultrafast inverse model solution based on a novel $(O(log(n))$ lookup table approach that is about 1,000x faster than the typical Levenberg–Marquardt algorithm [24] and uses binary trees often deployed in k-nearest neighbor (knn) algorithms [25]. In addition, for maximum utility over a wide range of tissue optical properties, our system provides user-selectable single or multi-frequency FD-NIRS modulation.

The platform provides, for the first time, display of high-speed, high-density 2D quantitative chromophore concentration topography maps in real-time. Specifically, the system can display single frequency, single wavelength FD-NIRS amplitude and phase data at a maximum rate of 36,600 Hz, and provides tissue optical property and chromophore concentrations at a maximum of 17,891 and 10,211 Hz, respectively. Importantly, all processing is done on board without any external computation, greatly simplifying the display of chromophores and optical properties. Display rates for six wavelengths single frequency amplitude/phase capture, optical property inversion, and chromophore calculation are 2,560, 2,328 and 2,327 Hz, respectively. Optical property accuracy was compared to a traditional network analyzer-based system [1] and was found to be within 10% over a range of attenuating tissue-simulating phantoms.

These ultra high speeds can provide advantages in image resolution, accuracy and scalability. For example, a real-time high-speed system allows for higher spatial density imaging of tissue that lends itself to better assessment of tissue heterogeneity and overall resolution. In real-time, users can adjust imaging parameters such as modulation frequency selection and averaging to dynamically optimize image quality. Some tissue types and applications may only require one modulation frequency, while others benefit from a range of frequencies [26,27] to maximize accuracy. Furthermore, the system’s high-speed capability allows for scaling to hundreds of source-detector channels, while maintaining real-time display of optical properties and chromophores. We show data demonstrating use-case examples of the ultra-high-speed capabilities enabled by our approach including tumor simulating phantom images, as well as in vivo arterial occlusion and pulsatile measurements.

2. Methods and materials

The following section describes the developed system (Section 2.1–2.4), a brief summary of the optical property recovery methods used (Section 2.5), the characterization measurements undertaken to assess performance and instrument response (Section 2.6), and human in vivo measurement protocols are described in (Section 2.7).

2.1 System overview

Fig. 1 shows the overall system design and architecture. The system runs a single or multi-frequency sweep, electrically modulating each of the laser diode sources in a serial manner from an amplified direct digital synthesizer (DDS). Scattered light is collected by a 3mm diameter avalanche photodiode (APD), which is sampled by an analog to digital converter (ADC) at 250 mega-samples-per-second (MSPS). We currently employ sample sizes of 4096, though arbitrary sample sizes (limited by the 128 GB storage in this system) can be used at a cost of lower data collection speed, trading off accuracy with speed. An FPGA captures the ADC data codes and performs a Fourier transform (FT) to extract FD-NIRS phase and amplitude, calibrates the data based on a prior calibration phantom measurement to remove system response, and then calculates optical properties using a novel look-up table approach- all in hardware. Measurements are taken and displayed using a wireless device which communicates with the embedded hardware over peer-to-peer WiFi. The display code was written for Apple’s iOS in Swift, the firmware code was written in C++, while the FPGA elements were written in Verilog. The current version of the hardware, mostly based upon commercially-available evaluation boards (excluding a display), fits in a 12 x 12 x 3" volume.

 

Fig. 1. FD-NIRS system design and architecture. Terms: System on a chip (SoC), programming logic (PL), processing system (PS), nine degree of freedom (9-DOF), direct digital synthesizer (DDS), analog to digital converter (ADC), avalanche photodiode (APD).

Download Full Size | PDF

2.2 Signals and sources

RF from 1-400 MHz are generated by a DDS chip (Analog Devices AD9912) powered by a high-speed low phase noise 1 GHz clock (Crystek CCSO-914X-1000), important for precise signal generation. An attenuated portion of RF from the DDS is coupled (Mini-Circuits ZX30-9-4) to the ADC (Analog Devices 9613) as a reference signal, while the remaining is sent to a dynamic gain control chain where a low noise amplifier (LNA; Mini-Circuits TSY-172LNB+) is chained to a digitally controlled attenuator (Skyworks SKY12347-362) with 0 to 31.5dB of attenuation.

The system uses six light sources: five edge-emitting laser diodes (EEL) at 690, 785, 808, 850, and 980 nm (Thorlabs HL6738MG, L785P090, L808P200, L850P030, L980P030), and a vertical cavity surface emitting laser (VCSEL, Vixar) at 940 nm. The VCSEL was used due to its relatively low-cost compared to the EEL equivalent at that wavelength, and we have previously shown that VCSELs are suitable for use in FD-NIRS systems [28]. We designed an aluminum block heatsink to house the sources for temperature and power stability. Each laser is fiber coupled to a 400 $\mu$m multimode fiber, with all six fibers collected into a single ferrule. A single laser driver (Thorlabs MLD203CHB) is connected to two identical 1x3 DC switches (Analog Devices ADG802BRTZ) on the mainboard (Fig. 1) providing 6 DC outputs. Each output provides a DC bias and is connected to each source with a bias tee, while the RF for each source is switched with an SP8T RF switch (Analog Devices HMC253AQS24E).

2.3 Detection

Light is collected with a 3mm avalanche photodiode (Hamamatsu S2384) that directly contacts the sample surface. The detected signal is amplified +15 dB with an LNA (Mini-Circuits TSY-172LNB+) to allow for full use of the full-scale voltage (FSV) of the ADC, and can be bypassed for close source/detector separations so as not to exceed FSV of the ADC. Typical signal levels from the APD range from 0 to - 70 dBm, and thus fit within the ADC’s FSV range to maximize ADC resolution. The single ended signal from the APD is then converted into a differential signal and ac coupled to the ADC.

The ADC samples at 250 MSPS at 12-bits of resolution, dual channel, and with 1.75V FSV. The ADC outputs digital codes to the FPGA in the system on a chip (SoC). Modulation frequencies above 125 MHz are aliased (due to the 250 MSPS sampling rate), however since all modulation frequencies are known, aliased frequencies in overlapping Nyquist zones are removed [21].

2.4 Signal processing and control

12-bit data codes from the ADC encode the voltage levels for two channels: the measured APD channel and a reference DDS channel that serves as the amplitude and phase reference (shown in Fig. 1). The ADC 12-bit codes are processed in hardware using an FPGA/Co-processing chip (Xilinx Zynq 7020) facilitated by a commercial off the shelf (COTS) SoC platform (Krtkl Snicker-doodle Black). The system runs Linux in the processor system (PS), which then controls the data processing being done in the programming logic (PL). We developed a custom high speed TCP/IP data packet protocol that connects the mobile display device wirelessly to the system using AES encrypted peer-to-peer Wi-Fi.

The signal processing chain is shown in Fig. 2. In the pre-processing steps (done immediately after a user sets the measurement parameters in the GUI), the coefficients for the Goertzel algorithm (described in detail in the following paragraphs) are calculated on the device. Dual-channel raw ADC sample codes are captured by the FPGA and as each sample is received, they are pipelined to the Goertzel algorithm block which collects and converts the samples to real and imaginary values for each frequency used in a measurement as shown in Fig. 2. Of note is that with this pipeline method, the finished Fourier transform output of amplitude and phase values are ready within a single clock of the ADC sampling being finished (4 ns), allowing for near immediate use of the FT data. The measured signal is then calibrated with data from a calibration phantom to remove the internal system response.

 

Fig. 2. Signal processing block diagram.

Download Full Size | PDF

Raw digital data codes from the ADC are processed directly within the FPGA fabric. A typical approach [21] might use a fast Fourier transform (FFT) algorithm. However, because we need data from only one of the frequency bins at a time, we can significantly reduce computation time by using an optimized discrete Fourier transform (DFT). An algorithm invented by Goertzel [29] computes the $k$th DFT portion of a signal ${x[n]}$ of length $N$. Essentially, the Goertzel algorithm performs a computation of a single DFT coefficient. This can be advantageous over a full DFT for several reasons. When only a few values of the spectral components of a signal are required, rather than the whole spectrum, the algorithm can be significantly faster. Because one can consider the Goertzel as an infinite impulse response (IIR) filter [29], the computation can begin immediately at the arrival of a sample, a key to real-time processing. This is opposed to an FFT where the entire sample set is required first before computation. Additionally, there is no need to output data in bit-reverse order or build a re-order logic block in FPGA fabric. As such, an extremely low latency (single clock cycle) pipeline can be used to do the computation. Furthermore, FFT efficiency is dependent largely on the sample length, $N$, which is most effective as a power of two. For the Goertzel algorithm, any arbitrary N size is possible without adding computational complexity. while others have used the Goertzel algorithm for processing FD-NIRS signals [17,18], our approach differs in that the constants used in the Goertzel algorithm are pre-calculated and stored in registers in the FPGA for each frequency in a sweep, where others have only done single frequency calculations in an FPGA, or a multi-frequency Goertzel algorithm done less efficiently in software. Our approach allows for single-clock computation for each frequency used.

The Goertzel approach minimizes the number of calculations required for a real-valued input signal. For an input signal, $2N$ additions and $N$ multiplications are performed, resulting on the order of $3N$ +2 computations for each frequency. This approach is more advantageous as compared to an FFT when the number of frequency bins $K < 2log2N$. Memory requirements are also reduced when $K<4/7 N$ as compared to an FFT [29]. In our case $K$ is simply one (a single frequency bin). As an example, for a sample length of 65,536, the above approach achieves a 32x decrease in computation and a 37,000x decrease in required memory usage over an FFT.

During an FD-NIRS measurement the system will generate large quantities of raw time-domain data for a single measurement (16 kB at 4096 samples, $1 \lambda, 1f$) or around 8 gigabits-per-second (sampled at 250 MSPS). Streaming this much data in real time to an external device is not practical. Our approach reduces this data output by embedding the FT and recovery of optical properties into the hardware, reducing the data output to only 8 bytes per measurement ($1 \lambda, 1f$) which is easily transferred at a real-world 2 megabits-per-second wireless bandwidth.

2.5 Model inversion and chromophore estimation

Raw amplitude and phase data captured in the signal processing chain as described in Section 2.4 represent a convolution of attenuation and phase delay from the measured tissue and the instrument response. To remove instrument response, data is calibrated against a silicone phantom with known optical properties [12] and the data is stored in on-board memory. Pre-calculated inverse solutions to the p1 semi-infinite model [30,31] are then loaded into memory as lookup-tables for high-speed optical property recovery. There is a table for each frequency in the sweep and for each source detector separation (SDS) used. The values are stored as complex and each table spans $\mu _a$ and $\mu _s^{'}$ values. The tables are indexed into binary search trees to enable very fast searching on the order of $Olog(n)$ [25,32,33]. After the Goertzel algorithm calculates the measured real/imaginary data for each channel (in parallel to the device sampling from the ADC), the system completes a k-nearest neighbor (knn) search using the squared Euclidean norm (a common method well described in [32,34] ) for a $\mu _a$ , $\mu _s^{'}$ solution that matches the measured complex value, where k is the number of dimensions (in this case two, for the real and imaginary part). For multi-frequency sweeps, these $\mu _a$ , $\mu _s^{'}$ solutions are averaged to find final $\mu _a$ , $\mu _s^{'}$ values. The results of the knn lookup table method are covered in detail in the following sections. Chromophore concentrations are then calculated in firmware using known extinction coefficients [35–37].

2.6 Characterization methods

The measurements we performed to characterize the FD-NIRS platform performance are presented in this section. Unless specified otherwise, the following experimental conditions were used. All wavelengths of 690, 785, 808, 850, 940, and 980 nm were operated. We chose to modulate the sources at a single frequency of 70 MHz, though any frequency from 0.1-400 MHz could have been used. Measurements were collected in a reflection geometry at an SDS of 23 mm. 4096 samples were taken at a speed of 250 MHz with the ADC. Optical and chromophore properties were extracted using the method described in Section 2.5. The table size used for recovery was 2500 x 700 ($\mu _a, \mu _s^{'}$) as we found this table size brought the quantization error to below 1% for both $\mu _a$ and $\mu _s^{'}$. Optical power was set to 20 mW for all sources as measured from the fiber ferrule.

2.6.1 Noise characterization

The noise floor of the system was measured in two ways. First, a ’dark’ measurement was taken wherein the 3mm APD detector in the system was biased and optically blocked, while the rest of the system was energized. The electrical noise floor was also measured without the APD or laser source power supplies energized, which should result in a much lower noise floor as optical systems are typically detector limited. Samples were taken with the ADC at 4,096 and 65,536 samples for the dark measurement, and at 65,536 for the electrical measurements to accurately measure the extremely low noise floor of the system [38]. The gold standard vector network analyzer (VNA) based system used an intermediate frequency (IF) bandwidth of 1 kHz.

2.6.2 Accuracy and precision

A series of 8 solid silicon tissue simulating phantoms were prepared with varying $\mu _a$ from 0.002-0.021 mm$^{-1}$ (at 660 nm) and $\mu _s^{'}$ varied from 0.5-1.4 mm$^{-1}$ (at 660 nm) [39]. The solid phantoms contain four components: water-soluble nigrosin ink to control absorption properties, anatase titanium (IV) oxide scattering properties (both from Sigma-Aldrich, St. Louis, MO), p4 silicone rubber base and activator (Eager Polymers, Chicago, IL) in a process described in [39]. Optical properties for the system were recovered using the p1 semi-infinite model in the form of a lookup table as described in Section 2.5 and calibrated with a silicone phantom with known optical properties. The gold standard system was VNA based as described in [1] and was used to measure the phantoms with a self-calibrated multi-distance method (SDS ranging from 1.5 to 3.0 cm) [31,40]. To assess precision, 12 seconds of data (30,720 samples for all six sources; 2,560 Hz display rate, 70 MHz single frequency) was collected on a tissue simulating phantom (designed for $\mu _a$= 0.01 mm$^{-1}$ $\mu _s^{'}$= 1.0 mm$^{-1}$ at 660 nm). 12 seconds was chosen as an arbitrarily short period in order to observe whether optical power transients were present. For the longer measurement, we chose 1 hr because it would exceed the expected continuous tissue measurement time as well as be sufficient to observe any potential heating effects.

2.6.3 Inversion speed

Inversions per-second measurements were made on two different hardware schemes. The execution speeds were averaged over 10,000 look-ups, single threaded, and measured using on-board system clocks, accurate to around 1 $\mu$s. The ’Modern CPU’ system used in Section 3.4 was a modern workstation processor (AMD 3960X, 3.8 GHz), where as the ’SoC’ was the system on a chip (Zynq 7020, Xilinx).

2.6.4 1D Line scans and 2D imaging with tissue-simulating phantoms

To demonstrate high-density 1D line scans and 2D imaging, a series of solid silicone tissue simulating phantoms were prepared with a simulated tumor inclusion 1 mm below the surface with diameters of 1.0, 1.5, 2.0, and 3.0 cm (substrate: $\mu _a$= 0.005-0.008 $mm^{-1}$, $\mu _s^{'}$= 0.6-0.9 $mm^{-1}$; inclusion $\mu _a$= 0.013 $mm^{-1}$, $\mu _s^{'}$= 1.4 $mm^{-1}$ at 660 nm). A millimeter resolution measuring stick was placed on the phantom next to the probe and marks were placed at 0 cm and 6 cm on the phantom. The acquisition time was 21 seconds for each scan. The probe housing the optical ferrule and 3 mm APD was moved across the surface over the phantom and tumor inclusions at a constant rate of 0.3 cm/s. Acquisition was timed such that the start and stop of data coincided with the probe location at the start mark at 0 cm and the stop mark at 6 cm. 2D images of a 30 mm diameter inclusion phantom were created using high-density vertical and horizontal line scans with 1 cm spacing for both for comparison to typical FD-NIRS grid measurements [14]. The phantom was marked with a grid of 1 cm spaced marks in both the horizontal and vertical directions, with a metal measuring stick used as a straight edge guide to keep each measurement line scan straight along the grid lines, making a 7 x 7 cm grid. Similar to the 1D line scan, acquisition time was such that the start and stop of data acquisition coincided with the probe location at the start mark at 0 cm and the stop mark at 7 cm. All 6 sources were used. The line scans were combined into a single matrix, after which bi-cubic interpolation was employed to generate the image.

2.7 In vivo measurements

In vivo arterial occlusion measurements were obtained by placing a blood pressure cuff on the upper arm of a healthy 39 year old male volunteer. The human study was approved by the Institutional Review Board of the University of Notre Dame (#19-09-5578) and the participant provided written informed consent. Continuous chromophore measurements acquired with the system were displayed and recorded in real-time at 2.33 kHz with all 6 sources at a single modulation frequency of 70 MHz. Optical properties were recovered following the method outlined in 2.5. The volunteer began the experiment in a quiet room at rest for 20 minutes prior to measurement. The APD probe was allowed to reach an equilibrium temperature prior to the beginning of the experiment. Data was collected on the right forearm for 2 minutes at rest, with the volunteer breathing normal and relaxed (baseline). After 2 minutes, the cuff was manually inflated to 250 mmHg to restrict blood flow (occlusion). After 5 minutes of occlusion, the cuff was slowly released, after which data was recorded for another 5 minutes (recovery).

An extremely high-speed pulse measurement (5.1 kHz) was also taken on the right-wrist of the same volunteer. Tissue chromophore concentrations were displayed and recorded in real-time at a single modulation frequency of 70 MHz. Data was acquired for 35 seconds with two wavelengths (690 nm and 850 nm) to recover oxyhemoglobin ($HbO_2$) and deoxyhemoglobin ($HHb$) concentrations. Data was smoothed using a 100-pt moving average. For this measurement the SDS was set to 12 mm as the radial artery targeted was near to the surface of the skin. Prior to the measurement the volunteer was at rest for 20 minutes in a quite room. During the measurement the volunteer was breathing normally. The volunteer’s heart rate was simultaneously measured at 67 BPM via manual palpation.

3. Results

3.1 System characterization

The six light sources are capable of maximum optical powers ranging from 30-200 mW. The DC and RF switching integrated circuits can operate very fast (˜90 ns), but a typical laser settling time was found to be 1 $\mu$s for this system, limiting source switching speed to around 1 MHz. The DC laser driver exhibits excellent stability with 12 $\mu$A average root-mean-square current with $\mu$W power stability. The laser driver can deliver up to 200 mA with 3V compliance voltage and the DDS can deliver 0-13.7 mA, with a maximum 7dBm power into a 50 $\Omega$ load at 50 MHz. After the amplifier chain, the maximum RF output was just below 20 dBm.

The overall noise floor of the system was found to be -68 dBm at 4096 and 65,536 samples, suggesting it is the true noise floor. The noise floor data is shown compared to a VNA based gold standard system in Fig. 3 [1]. The electrical noise floor, measured without the detector was -93 dBm (not shown), illustrating that system sensitivity is limited while the source and detector circuitry is energized. The modulation efficiency during measurements was found to be between 90-95%, depending on the source.

 

Fig. 3. Noise floor of the FD-NIRS system at 4096 samples as compared to a reference VNA-based system.

Download Full Size | PDF

3.2 Accuracy

Figure 4 shows optical property recovery of a range of tissue simulating phantoms, comparing our system with a gold standard VNA system. We find that for nearly all absorption and reduced scattering values, the system was within 10% of the gold standard across all wavelengths, suggesting acceptable agreement between the systems. $\mu _s^{'}$ values had the most agreement with the reference system when values were near 0.75 mm$^{-1}$, while $\mu _a$ were in closer agreement with values lower than 0.013 mm$^{-1}$. The mean difference for both $\mu _a$ and $\mu _s^{'}$ were nearly identical to the VNA based gold standard.

 

Fig. 4. Bland-Altman plots of optical properties measured with the system as compared to the reference VNA based system. The center line shows the mean difference, while the top and bottom dashed lines denote 10% deviation from the mean difference. Nearly all absorption and reduced scattering differences fall within 10% of the VNA gold standard system.

Download Full Size | PDF

3.3 Broadband modulation

To demonstrate the system’s broadband modulation, a FD-NIRS measurement from 50-300 MHz in 1MHz steps was taken on a tissue-simulating phantom (Fig. 5). Data was filtered to remove frequencies surrounding the Nyquist zones at 125 MHz and 250 MHz [16]. The 850 nm amplitude and phase data are shown, though the other wavelengths had similar fits. The measured amplitude and phase were fit to the p1 semi-infinite diffusion approximation [41,42], and were found to be within 4.25% and 2.43% for $\mu _a$ and $\mu _s^{'}$, respectively, as compared to the gold standard system outlined in section 2.6.1.

 

Fig. 5. Broadband phase and amplitude data is shown for 850 nm with a 250 frequency sweep from 50-300 MHz on a tissue simulating phantom along with a p1 semi-infinite diffusion model fit.

Download Full Size | PDF

3.4 High-speed inversion

Table 1 shows the inversion rates of single and multi-frequency optical property recovery using the high-speed knn lookup table method outlined in Section 2.5. A low and high density lookup table for $\mu _a$ and $\mu _s^{'}$ was chosen to showcase the (ideally) $(O)log(n)$ [25] performance made possible with the knn method as values of n (table size) grows. As an example, a 7x increase in table size from 500x500 to 2500x700 resulted in only a 2% and 5% decrease in inversion/s speed for 1$f$ in the modern CPU and SoC, respectively. The knn lookup table algorithm running on the ’modern cpu’ (AMD 3960X, 3.8 GHz) is between 6-17x faster in inversion/s over the lower performance SoC (Zynq 7020, 866 MHz) in both single and multi-frequency inversion schemes, as was expected. Our results show that our knn inversion lookup table implementation on a modern CPU is between 2-3 orders of magnitude faster than prior approaches [19]. Furthermore, using the knn inversion method on the fully onboard SoC shows inversion rates to be 1-2 orders of magnitude faster than previously demonstrated [19]. We note that the SoC rates still results in significant improvement, even though the processor is less performant. We also note that for multi-frequency tables, the performance scales as $(O)(f * log(n))$ where $f$ is the number of frequency tables.

Table 1. High-speed knn lookup table inversions/s

View Table | View all tables in this article

3.5 High-speed acquisition

The device acquires high speed FD-NIRS measurements in a serial manner; each source is modulated sequentially and can be swept in modulation frequency from lowest to highest for multi-frequency measurements. The raw data capture speed of the 250 MSPS ADC dictates the maximum obtainable speed as we need to capture enough of the waveform to extract accurate amplitude and phase information. For example, at a 4096 sample rate, data acquisition takes around 16.4 $\mu$s, leading to an upper speed limit of 61,035 Hz (the fastest speed at that sample size). Smaller sample lengths can be used, though at the expense of phase resolution and thus accurate optical property recovery [21]. In practice, we found that 4096 is a good balance between data acquisition speed and accuracy of recovered optical properties.

Since our primary aim is to use the instrument to display data in real-time for use in a clinical system, we report our realized display rates for each stage of the data processing and display in Table 2. For a single-source/single-frequency measurement, the maximum rate for amplitude and phase measurement and display is 36,608 Hz, optical properties were found to be measured and displayed at 17,891 Hz, and chromophores at 10,211 Hz. Amplitude and phase from all six sources can be displayed at 2,560 Hz, optical properties at 1,615 Hz, and chromophores at 2,327 Hz. However, this can be further improved as there is a small delay introduced by the serial peripheral interface (SPI) communication driver associated with the switching hardware. A six-source, multi-frequency sweep from 50-300 MHz at 7 MHz steps (35 frequencies) results in amplitude/phase display at 175 Hz, with optical properties and chromophores displayed at 60.2 Hz. This speed could also be significantly improved if source multiplexing is used such as in [16,18].

Table 2. Display rates in Hz for each phase of the pipeline required to display four chromophore values ($HHb$, $HbO_2$, $H_2O$, and lipid).

View Table | View all tables in this article

If ultra high speed data display is not required, such as for single SDS configurations or applications requiring long-term continuous measurements, precision can be improved by averaging the data. Table 3 illustrates how averaging improves the phase precision at the expense of speed. The standard deviation decreases by more than 8.5x from 0.204 degrees (no averaging) to 0.024 degrees (30,000 samples averaged). The standard deviation and coefficient of variation (CV) for amplitude over the same conditions improved by 1.7x. This resulted in improved precision in recovered absorption and reduced scattering coefficients by 2 - 3x (no averaging vs 30,000 sample block average).

Table 3. Block averaging showing effect on phase precision. 70 MHz, 785 nm source (others wavelengths were similar).

View Table | View all tables in this article

3.6 Stability and precision

System precision and stability was assessed by analyzing data collected continuously over a short (12 s) and long (1 hr) periods. For the short-term precision assessment (Table 4), 12-seconds of high-speed data was collected on a tissue simulating phantom. The system displayed a coefficient of variance (CV) of less than 1% for all amplitudes and phase standard deviation typically between 0.21-0.33 degrees, though the 808 nm was slightly higher at 0.56 degrees due to its lower signal strength. A one-hour stability measurement was taken, wherein the same tissue simulating phantom was measured continuously at 2 second intervals. Optical properties were extracted using the method described in Section 2.5. The sources showed good optical property stability for $\mu _a$ and $\mu _s^{'}$ with an average CV of 3.66% and 1.92%, respectively, across all wavelengths as shown in Table 5. We note that the data sampled at 2 Hz was not averaged. If increased precision is necessary, the platform is capable of averaging optical properties for all six wavelengths over longer time periods, as was shown in Table 3.

Table 4. High-speed system stability during a 12-second measurement with a display rate of 2,560 Hz, for all six sources at 23 mm SDS and 70 MHz on a tissue simulating phantom ($\mu _a$= 0.01 $mm^{-1}$$\mu _s^{'}$= 1 $mm^{-1}$ at 660 nm).

View Table | View all tables in this article

Table 5. Standard deviation and coefficients of variation for optical properties over a 1-hr drift test across the six wavelengths used.

View Table | View all tables in this article

3.7 High-density, high-speed measurements of phantom inclusions

3.7.1 1D Line Scans

1D line-scan optical property measurements of tissue-simulating phantoms containing optical inclusions were taken as shown in Fig. 6 and the Visualization 1 video. Results for the 808 nm data is shown in Fig. 7, though the other five wavelengths had exhibited a similar shape. Data was acquired at full speed of 2.6 kHz, then 30 pt moving averaged, totaling 48,000 views and 288,000 total samples (for all six sources) for each line scan. Capture time for each measurement was 21 seconds. The $1/e^2$ values (a common beam diameter definition [43]) were calculated for each resulting curve for $\mu _a$ and $\mu _s^{'}$ at 1.0, 1.5, 2.0, and 3.0 cm diameters, resulting in errors of 81% (0.4%), 20.6% (0.6%), 8.5% (-1%) and 1.6% (3.3%), respectively. Measured inclusion widths of the subsurface inclusions for all wavelengths are shown in Table S1 in the supplemental material. Due to the extremely high density of data points, the curves show how smooth and information rich the data is when using high-speed data acquisition. The data suggests $1/e^2$ values of the high-spatial resolution reduced scattering information allows for an accurate estimate on the size of the underlying inclusion at these depths and SDS, but localization in the diffuse optical regime is generally a complicated function of SDS, inclusion depth, and optical properties. We also found that high resolution reduced scattering $1/e^2$ values better estimated the size of the underling inclusion when compared to low resolution (single cm) sampling. However, the choice of width metric is somewhat subjective [26]. Regardless, scattering information may give more accurate localization than absorption as others have suggested [10,44], but more inquiry is required.

 

Fig. 6. (a) A tissue simulating phantom (designed for $\mu _a$= 0.006 $mm^{-1}$, $\mu _s^{'}$= 0.9 $mm^{-1}$ at 660 nm) is impregnated with a 3.0 cm simulated tumor inclusion (designed for $\mu _a$= 0.01 $mm^{-1}$ $\mu _s^{'}$= 1.3 $mm^{-1}$ at 660 nm) 1 mm below the surface. (b) Six sources are fiber coupled to a handheld probe housing a 3 mm APD detector. (c) Absorption and reduced scattering values from the system are wirelessly transferred and displayed in real time at high speed (2,328 Hz at 70 MHz, 1$f$, 6$\lambda$) to a mobile device. The simulated tumor optical properties are visualized on the display as the probe is moved along the surface over the top of the tumor (see Visualization 1).

Download Full Size | PDF

 

Fig. 7. (a) 1D line-scans showing absorption and reduced scattering coefficient data taken of tissue simulating phantoms (substrate: $(\mu_a)= 0.005-0.008 mm^{-1}$, $\mu_s^{'}= 0.6-0.9 mm^{-1}$; inclusion $(\mu_a)= 0.013 mm^{-1}$, $\mu_s^{'}= 1.4 mm^{-1}$ at 660 nm), shown at 808 nm, single frequency (70 MHz). Simulated tumor inclusions were 1 mm below the surface in the center of the phantom on 10, 15, 20, and 30 mm diameter sizes. (b) Experimental setup diagram.

Download Full Size | PDF

3.7.2 2D spatial mapping

2D images of the 30 mm diameter inclusion phantom were created using high-density vertical and horizontal data (8 line scans in each direction) as shown Fig. 8(c). The high-density data set was captured at a 2.6 kHz data rate, at 22,000 samples per line scan for a total of 352,000 points per image, 2,112,000 total samples for all six sources. The image acquisition time was about 3 minutes. High-density images for all six wavelengths are shown in supplemental Fig. S1. For comparison Fig. 8(a) shows an image downsampled to 1 cm spacing (comprising of 49 points), which was the measurement spacing used in the ACRIN 6691 DOS imaging trial of breast cancer neoadjuvant chemotherapy response [14]. The 7 x 7 cm grid size was chosen to have 2 cm space on each side of the inclusion to show contrast between the background optical properties and the inclusion. Data from Fig. 8(a) was then interpolated to form Fig. 8(b). The lower spatial density acquisition image of Fig. 8(a) and Fig. 8(b) suggests valuable absorption and scattering information is missed, as compared to the high-density images in Fig. 8(c) where the 30 mm tumor topography is more accurately represented. This represents highest density 2D topographical FD-NIRS image created to date, with 43,100x times more spatial density than that used in ACRIN 6691 [14] and as reported more recently [19], about a 1,550x improvement in chromophore mapping speed. The line/stitching artifacts manifest in the images area a result of the bi-cubic interpolation of high resolution line scans in the x and y directions separated by 1 cm, creating a ’basket weave’ data pattern where the area in between each line scan contains interpolated data. Other interpolation methods will be explored in future work.

 

Fig. 8. (a) A 7 x 7 cm image of a tissue simulating phantom with raw data at 1 x 1 cm spacing (a typical spatial measurement for FD-NIRS systems) comprising 49 points. (b) Bi-cubic interpolation of (a). (c) High-density 2D spatial maps acquired with the system at 352,000 points per image, 2,112,000 total samples (for all six wavelengths), bi-cubic interpolation, 3 cm diameter inclusion at 850 nm and 70 MHz modulation frequency. The high-density image more accurately recovers the 3 cm inclusion shape.

Download Full Size | PDF

3.8 In vivo real-time display

To demonstrate the real-time processing and visualization abilities of the system, we performed an in vivo arterial occlusion measurement using a blood pressure cuff placed on the upper arm of a healthy volunteer. Amplitude and phase data taken with the system were recorded at 2.6 kHz using all 6 sources. Fig. 9 shows $HbO_2$ and $HHb$ concentrations recorded during the experiment. The data is displayed with a 1000-point (0.38s) moving average plotted for each chromophore. Cuff occlusion began at the 90 second mark, after which we observe a slight increase in $HbO_2$, likely due to slow manual cuff inflation. The data then shows a steady falling concentration of $HbO_2$, while $HHb$ concentration rises as expected due to the arterial occlusion. After 5 minutes of occlusion, the cuff was released and a sharp hyperemic response and rebound was observed. The SNR was found to be from 37 to 70 dB across all six wavelengths, with the lower SNR of 37 dB found at the 980 nm wavelength due to the strong absorption peak in water near that wavelength.

 

Fig. 9. An arterial occlusion experiment showing tissue $HHb$ (right) and $HbO_2$ (left) curves from the forearm of a healthy human subject. Lines indicate the start and end of the occlusion on the upper arm. Measurements were acquired at 2.6 kHz with a 1000-point moving average plotted for each chromophore.

Download Full Size | PDF

In a demonstration of the the system’s extremely high-speed sampling, measurements of chromophore data from two wavelengths (689 nm and 850 nm) on a healthy volunteer’s wrist are shown in Fig. 10. A pulsatile signal is clearly shown in the $HbO_2$ curve and was measured with the systems to be 67 beats-per-minute (BPM), matching manual palpation using a stopwatch. The majority of the high temporal resolution pulse waveforms show the dicrotic notch, and we also observe the respiratory envelope superimposed on the $HbO_2$ curve, estimated to be just under 10 seconds per breath. The inset in Fig. 10 displays $HbO_2$ changes displayed in a zoomed curve of a single pulse peak where high-density points are indicated by red markers.

 

Fig. 10. A pulsatile signal (~67 BPM) from a human wrist sampled at an extremely high-rate (5.4 kHz) is shown. The $HbO_2$ curve displays a respiratory envelope function due to normal subject breathing (~10 second cycles). Inset: $HbO_2$ changes are displayed in a zoomed curve of a single pulse peak where high-density points are indicated by red markers.

Download Full Size | PDF

4. Discussion and conclusion

Over 25 years of clinical research has shown that FD-NIRS based optical imaging can be a meaningful tool for understanding, diagnosing, and treating various human diseases. For example, the ACRIN 6691 multi-center study demonstrated that FD-NIRS based methods predict neoadjuvant chemotherapy response in breast cancer with performance similar or better than MRI [45]. For functional near-infrared spectroscopy (fNIRS), FD-NIRS can provide deeper sensitivity as compared to CW-fNIRS and therefore higher sensitivity to functional activation [10]. FD-NIRS has also shown promising clinical utility in breast cancer diagnosis, critical care, exercise physiology, among other applications. However, there is currently only one company producing FD-NIRS systems (ISS, Inc.), and none are FDA-cleared for clinical use. Given our experience, we suggest that there are three major limitations that have contributed to the slow clinical adoption of FD-NIRS. First, diffuse imaging methods typically exhibit low spatial resolution, thus providing data unlike standard medical imaging modalities. Second, the FD-NIRS research tools developed so far have been complex, slow, and prohibitively difficult to use. Finally, previous designs have lacked easy scalability necessary for high density and higher resolution imaging without greatly increasing the instrumentation footprint. The FD-NIRS platform we developed in this work was designed to overcome these outstanding issues through improved scalability, ease of use, and speed. Together, we believe that these improvements can lead to the higher spatial resolution, localization accuracy, sensitivity, and usability needed to accelerate clinical adoption.

Thus we demonstrated an all digital, scalable FD-NIRS platform capable of ultra-high speed chromophore processing and display utilizing an embedded system and external display. Digital ICs, RF and DC switches, filters and dynamic signal chain components were all specifically chosen for scalability, high-speed performance, and low noise. The DDS and RF circuitry allows for single or multi-frequency direct modulation of laser sources similar to VNA-based counterparts [1]. This design can enable, for example, a fully embedded FD-NIRS system in a single handheld device, which could increase FD-NIRS portability and accessibility. The system also showed comparable noise levels to VNA-based systems [1,46]. We characterized the accuracy of the system in comparison to a reference system and found it to be within 10% across all wavelengths, suggesting acceptable equivalence.

We showed that FD-NIRS amplitude/phase, optical properties, and chromophores can be acquired, processed entirely on an embedded platform, and displayed at speeds up to 36,608 Hz, 17,891 Hz, and 10,211 Hz respectively. To the best of our knowledge, this is the fastest FD-NIRS system demonstrated to date, especially in its capability to calculate and display optical properties. For comparison, Zhao et al. showed DNN based chromophore recovery and display at rates up to 27 Hz [18], Applegate et al. streamed FD-NIRS data to a display at 1.5 Hz [19], and the ISS OxiplexTS200 and ISS Imagent report 50-100 Hz rates. In comparison, the platform presented here can display optical property and chromophore data between 100x to 3,600x faster than these systems, as illustrated in Table 2. Our platform includes two key innovations to achieve this performance. The FPGA-based Goertzel implementation enabled amplitude and phase information to be extracted from the raw ADC time domain signal at just 1/2 of the actual sampling speed (first demonstrated in [47,48]). This strategy allows for a dramatic reduction in data storage (37,000x for 65,536 samples) and is a fundamental requirement to streaming data in real-time due to the large size of raw sampled ADC data. Second, we implemented an ultra-high-speed knn lookup table method into hardware, further reducing downstream computation and maximizing inversion/s speed, while also allowing chromophore computation to be done within the embedded SoC. Our knn lookup table method resulted in a 13x - 322x (in SoC) and 170x - 2,200x (modern CPU) improvement over previously reported literature [27]. The lookup table approach was chosen as it can be significantly faster (1,000x) than iterative linear optimization solutions typically used to find inverse model solutions without loss in accuracy [24]. Deep neural network (DNN) approaches can also provide fast inversion solutions, though at the expense of larger average errors than typical iterative solutions [18]. DNN solutions are also limited by their sensitivity to input parameters; a new DNN model may be required for each SDS, modulation frequency and may be more susceptible to noise [18].

A valid question to ask regarding high speed NIRS systems such as these is: why is it necessary to go so fast? We believe that these extremely high speeds provide significant advantages in scalability, accuracy, and imaging resolution. First, there is a desire among fNIRS researchers and users to cover the human head with greater area and density of optodes. It is estimated that 370 optodes are needed to cover the adult scalp for high-density fNIRS [49]. Our platform can easily scale to this number of optodes and thus provide quantitative, physiological information over the entire scalp in real-time. This could be of major clinical significance for patients that lack baseline measurements, such as those with brain trauma. The increased speed and sensitivity could also be used to further investigate the fast (<100ms) optical scattering-driven evoked-response optical signals that may be indicative of neural activity [50,51], but difficult to measure.

Second, high-speed FD-NIRS acquisition and processing gives the user significant control over modulation frequency choices and data averaging. These parameters can be optimized specifically for the desired application all while providing real-time display, similar to the various tissue modes of an ultrasound device. For example, some applications and tissue types can benefit from a range of modulation frequencies, while others require only one [26,27]. We also showed that the large amount of streaming data can be averaged to provide an order of magnitude improvement in phase precision at the expense of speed if desired.

Finally, this real-time high-speed FD-NIRS platform can be used to provide higher spatial density imaging of the tissue and lead to better assessment of tissue heterogeneity and overall resolution. A related optical imaging technique–spatial frequency domain imaging (SFDI)–can provide quantitative 2D spatial mapping in real-time, but has limited depth penetration (up to a few mm). In contrast the FD-NIRS system described here is a point-based imaging modality that can measure up to a few cm deep. To illustrate, we generated ultra-high-spatial-density 1D and 2D scans of phantoms with tumor-simulating inclusions at sizes from 1 - 3 cm. The results reveal a smooth, information-rich assessment of the phantom’s optical properties, as well as the ability to accurately differentiate the inclusion sizes. As an example of real-time high-speed optical property display, a 7 cm line scan in Fig. 6 took 21 seconds, recording 48,000 views and 288,000 samples (1$f,$6$\lambda$) demonstrating significantly faster imaging speed than other quantitative NIRS approaches [14,19]. The 2D quantitative optical property topography maps (Fig. 8) contain 352,000 points per image, representing 2,112,000 samples (for 6$\lambda$). This represents the highest density 2D topographical FD-NIRS image created to date, with 43,100x times more spatial density than the grids used in ACRIN 6691 [14]. Accurately assessing spatial heterogeneity can be significant in clinical applications; for example, it may be a marker in determining neoadjuvant chemotherapy response [52]. Importantly, the images can be constructed in real-time, allowing a user to re-take data if the acquisition was poor or to dynamically adjust the imaging region of interest.

In conclusion, we have demonstrated a new FD-NIRS platform with unprecedented scalability and speed that can lead to new compact and potentially wearable quantitative diffuse optical spectroscopy devices. These improvements are not only applicable to existing applications, such as fNIRS and breast cancer imaging, but can also open up new applications as the technology is further disseminated. Our future work entails miniaturization of the platform into handheld and wearable form factors, as well as application of the platform to acquire high resolution 2D and 3D diffuse optical imaging of tissue.