Frequency domain near-infrared spectroscopy (FD-NIRS) is a well-established technique for measuring tissue optical properties.1,2 The FD-NIRS instrumentation delivers radio-frequency (RF) intensity-modulated light and measures the amplitude and phase of the diffusely re-emitted light either in reflection or transmission geometry. These measurements enable the calculation of absolute optical absorption and scattering properties of tissue. Such data acquired at multiple near-infrared wavelengths allow quantification of tissue chromophore concentrations, such as oxy- and deoxy-hemoglobin, water, and lipids.3,4 Numerous applications of FD-NIRS techniques have been reported, including breast,5,6 brain,7,8 and muscle 9,10 tissues to name a few.
With spatially distributed source and detector probes, FD-NIR technology has also been used in diffuse optical tomography (DOT) systems, primarily in the area of breast imaging for both diagnosis and therapy-monitoring applications.11–14 In this context, FD-NIR imagers provide valuable information regarding the spatial variation of absolute chromophore concentrations either by RF measurements per se12 or enhanced with additional continuous wave (CW) source and detector sets.14,15 However, unlike CW optical systems which have achieved high imaging rates through a combination of frequency and time multiplexing,16–19 FD imagers still require fairly long acquisition times (up to 60 s) to obtain a full set of tomographic data. The lengthy measurement is a result of sequential cycling through both laser wavelengths and output fibers, where sufficient integration time is needed for each data point to overcome the inherently lower signal-to-noise ratio (SNR) of RF measurements. The need to improve these long acquisition times has become more important recently with the increased use of dynamic breast optical tomography to characterize the breast tissue response to external stimuli20,21 which requires obtaining a full dataset at least every few seconds.
Here, we present a multiwavelength FD-NIR imager design that combines a fast galvo-based fiber multiplexer on the source side with direct analog-to-digital conversion (ADC) on the detection side. The direct sampling method, enabled by advances in high speed ADCs, offers two significant advantages compared to the commonly used homo- and heterodyne I/Q demodulation schemes.1 On the one hand, at least 4 to 8 modulation frequencies can be used simultaneously in conjunction with discrete Fourier transform (DFT) demodulation before detector saturation or decreased noise performance become a concern. This allows operation of multiple laser sources without the need for multiple, frequency specific I/Q circuits, as required in homodyne demodulation; on the other hand, by placing the ADC close to the photodetector module the RF analog circuit design is simplified and there is much less chance for RF leakage, interference, and channel-to-channel crosstalk.
While the direct sampling approach has been increasingly used in the telecommunications field, it has only recently been implemented for diffuse optical spectrometry by our group22 as well as researchers at U.C. Irvine.23 This article contains, to our knowledge, the first report of a direct digital full-sampling FD-NIR breast optical tomography device. We describe the hardware design and show system characterization data and titration measurements. Our results demonstrate performance that matches or exceeds other current FD-NIR instrumentation at significantly higher data acquisition rates. The robust modular design of the system can be easily expanded to larger configurations with respect to both source wavelengths, as well as source and detector arrays.
The FD-NIR instrument described in this paper has been designed to integrate into a multimodal optical imaging–digital breast tomosynthesis breast imager capable of static and dynamic DOT with mammographic guidance. The complete system is shown in Fig. 1 and contains both CW and frequency domain subunits. The two components of the 2 wavelength, 24 source, and 20 detector FD-NIR subunit are on the upper side of the cart. The hardware units are linked by control and signal routing cables. The lower housing marked (a) in the figure contains the source component, while the upper housing marked (b) contains the detection component, both described below.
A schematic overview of the FD-NIR system can be seen in Fig. 2. The source unit provides dual wavelength (685 and 830 nm) simultaneous illumination by modulating two lasers at two separate frequencies (frequency division multiplexing), while the detector unit provides multichannel simultaneous multifrequency demodulation enabled by full digital sampling.
The 50-MHz reference clock for the instrument is generated by a temperature compensated crystal oscillator (Connor-Winfield D75J). A 2.7-GHz clock signal is derived from this reference by a phase-locked loop (PLL; Analog Devices AD9517). The laser modulation frequencies, as well as the sampling clock for the ADC, are then derived from the PLL output by integer division. Final division by two for the laser modulation frequencies is only done on the laser driver boards to keep the signals as localized as possible to simplify shielding of electromagnetic interference. Distribution of the clock signal to all 20 detector cards is accomplished by two low-voltage differential signaling (LVDS) clock fanout buffers (Analog Devices ADCLK854). The circuit design has been carefully optimized to minimize clock jitter and ensure detector card synchronization.
The light sources consist of one 50-mW 685-nm laser diode (Opnext HL6750MG) modulated at 67.5 and a 50-mW 830-nm laser diode (Opnext HL8338MG) modulated at 75 MHz. The beams of both diodes are combined by a dichroic mirror (Semrock FF757-Di01) and then launched via a two-dimensional (2-D) scanning galvanometer (ThorLabs GVS002) into multimode silica fibers. Twenty-five fibers are arranged in an array so that the illuminated fiber can be selected by steering the galvanometer to the appropriate angle. Additionally, the beam can be steered into a beam dump, which allows the system to acquire the dark signal. Using an off-the-shelf galvanometer instead of a traditional optical switch as an optical multiplexer offers a cost-effective solution to allow dynamic adjustment of scanning patterns and to accommodate more source locations without needing additional hardware. The galvanometer is controlled by its own dedicated microcontroller (Atmel ATmega2560) to ensure fast switching between fibers. It also contains the necessary position calibration procedures. The complete source unit assembly is shown in Fig. 3(a). When switching between neighboring source fibers, as it is typically the case during a measurement, we achieve a switch time of from trigger to illumination of the fiber. Stabilization of the signal is achieved within 2 ms. In the worst case transition we achieve a switch time of 2 ms, and the signal is completely settled within 3 ms.
After light passes through tissue, it is collected by fused silica fiber bundles and routed to the photodetectors. Each of the 20 detection channels is built on its own printed circuit board [see Fig. 3(b)] to allow for easy replacement and future expandability. On each of the detector boards an avalanche photo diode (APD) module (Hamamatsu C5331-04) detects the optical signal. The module also contains the first gain stage in the form of a low noise trans-impedance amplifier, as well as the necessary high voltage supply to bias the APD. The bias voltage is temperature compensated to keep the gain of the APD constant.
After the APD module, the signal is further amplified by a high-speed, low noise current feedback op-amp (Analog Devices AD8000). Then the signal is filtered with a bandpass filter with a 63- to 77-MHz passband (maximum attenuation 1.5 dB) and rejection below 50 MHz and above 95 MHz (Mini-Circuits SXBP-70+) in order to block other interference (e.g., in our breast tomography system, light coming from the CW-NIR instrument, which could otherwise saturate the ADCs). In the next step, the single ended signal is converted into a differential signal by a transformer. Then it is fed into a differential amplifier (Linear Technology LTC6409) which serves as the last gain stage, as well as ADC buffer and anti-aliasing low-pass filter. Each section of the analog signal chain is individually shielded to prevent interchannel crosstalk between neighboring detector cards.
Finally, the signal is sampled 180 million times per second (180 MSPS) with 16-bit resolution by a high-speed ADC (Linear Technology LT2209). The ADC is directly connected to its own dedicated low-cost field-programmable gate array (FPGA, Xilinx Spartan-6 LX9) which demodulates the signal.
All 20 detector cards as well as the control card are connected by a common backplane. A microcontroller (Atmel ATmega2560) collects the data from all detector cards via a serial peripheral interface bus with LVDS. The data are then sent on to the PC via USB for further signal processing and data recording.
As mentioned in the previous section, the analog signal is sampled 180 million times per second at 16-bit resolution. From this raw data stream, the signals from the 685 and 830 nm lasers have to be extracted. In a traditional analog instrument, this would be accomplished using homo- or heterodyne detection with a mixer for down conversion and then a slow ADC for sampling the in-phase and quadrature signals.1
In our instrument, however, this has to be done digitally. One possibility would be to use a fast Fourier transform (FFT) algorithm. However, a fairly large amount of frequency bins is required to keep the width of each individual bin small enough, and thus the memory and computational requirements would be very high, and a cost-effective solution for a 20 channel instrument is not possible.
The FFT algorithm computes all bins simultaneously while for our application we only need to obtain the values from two bins. Therefore, another possibility is to directly implement a DFT, where we compute the two bins we are interested in individually, without the need of having to compute all the other bins. In our design, we achieve a frequency bin width comparable to a 4 million point FFT at a much lower computational cost. We have also further optimized the DFT algorithm, as described below, to avoid multiplication in the high frequency section.
The standard DFT is computed as shown in Eq. (1), where is the input signal, , and k denotes the frequency bin to be computed.
To simplify we turn to the Goertzel algorithm. As shown in Ref. 24, the ′th bin of the -point DFT can be computed by feeding the signal into a system with impulse response which is initially at rest (i.e., all registers are initialized to zero). is the unit step function. The desired result is then the ’th output value [see Eqs. (2) and (3)] of the system. See also flow graph in Fig. 4.Fig. 5. Note that some nonrecursive elements have been added to compensate for the added delay and keep the overall transfer function unchanged.
The transfer function of the system of Fig. 5 is shown in Eq. (4) for the case when . The structure can be further simplified by multiplying the numerator and denominator with a common factor as shown in Eqs. (5) and (6) (similar to Ref. 24):Fig. 6. The complex multiplication in the recursive loop is replaced with a real multiplication at the expense of an additional term in the nonrecursive part. The modulation frequencies of the instrument can be chosen freely within a certain range. Thus, if we choose the modulation frequency of one laser as corresponding to , the cosine becomes zero, and thus the structure can be further simplified as shown in Fig. 7. The recursive loop is now reduced to a simple adder, for a real input sequence the numbers stay real in this part of the computation. A similar structure can be found for the case of and the modulation frequency of 75 MHz.
Overall this novel algorithm results in minimal resource requirements in the FPGA. Specifically, for a 4 million point DFT computation only one 36 bit adder is required, and the adder can have a latency of up to four clock cycles, which allows us to employ pipelining in its implementation. The complex multiplications in the nonrecursive structure can all be done in real time, e.g., in MATLAB®, since only the last four results of the recursive structure are needed, and thus the amount of data to be transferred and the amount of computation are fairly low.
We implemented the above algorithm on the FPGA of each detector card. For each modulation frequency, we calculate two DFT’s with 50% overlap and length of in real time. This results in an output data rate of 90 Hz. All 180 million samples per second the ADC acquires are actually used for calculation of the output data.
For this implementation only four 36 bit adders and some control logic are required. Everything fits even in the smallest FPGA model of Xilinx’s low cost Spartan-6 line, with room to add additional frequencies if required later on. On the PC side, MATLAB® has to perform up to five complex multiplications, and up to five complex additions per sample per channel, or up to 18,000 complex multiplications and additions per second overall to decode and record the amplitude and phase raw data in real-time, a task any recent computer can easily handle.
Figure 8 shows a typical plot of signal magnitude versus incident optical power of our system, obtained by sending the modulated light through neutral density filters of different attenuation. The optical powers given are corrected for modulation depth, which is typically 85%. The intersection of the gray asymptotes drawn in the figure indicates the optical power corresponding to an SNR of one. We typically measure a noise equivalent power of , which is approaching the manufacturer specified noise floor of the APD module of .
In the same figure, the responses of the neighboring, nondriven detectors are plotted in red. As can be seen, even at high optical input power into the main, driven detector, the neighboring, nondriven detectors show no signal. The channel separation on the detector side is thus (five orders of magnitude). However, on the laser source side there is some measurable crosstalk between sources because of the proximity of the fibers in the optical multiplexer. We measure a channel separation of . The overall system is thus limited by the optical multiplexer and shows a channel separation over 80 dB.
The ADC converter saturates when a signal of at 830 nm is fed to the APD. This is an order of magnitude higher than the performance of our commercial CW-NIR system (TechEn CW6), and it is also much higher than the signals encountered in the transmission type measurements used in our breast scanner. Together with the noise floor of 1.4 pW this results in an instantaneous dynamic range of more than 115 dB. At 685 nm, the dynamic range is even slightly larger because saturation is only reached at signal power.
We were not able to observe any interwavelength crosstalk nor amplitude to phase crosstalk up to the maximum specified input power of . The phase noise of the output signal is smaller than at 100 pW input power. Stability over 10 h was measured after leaving the instrument on for 1 h in a climate controlled room. The measured amplitude changed by , the phase at an optical power of 5 nW. Our findings are summarized in Table 1.
Summary of performance tests.
|Phase noise||at 100 pW|
|Maximum input power|
|Amplitude stability||at 10 h, 5 nW|
|Phase stability||3 mrad at 10 h, 5 nW|
Tissue-Like Phantom Optical Property Recovery
To test the ability of the instrument to recover absolute optical properties, we performed several titration experiments using a mixture of water, whole milk bought at the local grocery store, and black India ink (Higgins 44204). In particular, we measured the optical properties of 200 mixtures with milk concentrations ranging from 7% to 37%, and ink concentrations from 0 to 45.6 ppm. The upper concentration limits correspond to maximal optical properties typically observed in patients.25 The measurements were performed in multidistance reflectance mode,26,27 with source detector separations of 15, 20, 25, and 30 mm, and analyzed using the semi-infinite form of the diffusion approximation.
We used 10 equally spaced milk concentrations between 7% and 37%. For each milk concentration, we took measurements over a series of 20 absorption steps corresponding to 0 to 45.6 ppm of ink. Each phantom started as a 1500-mL mixture of milk (7% and 37%) and water (to 100%; no ink was present in the phantom at baseline). The probe was held in place by a holder on the surface of the liquid, submerged to ensure good optical contact. Ink was then added to the milk and water mixture by a computer-controlled syringe pump. Measurements were taken in between each ink adding step. To ensure proper mixing of all three components, we used a magnetic stirrer, and waited 30 s after each ink addition. This experiment was repeated 10 times, for each level of milk concentration. The required calibration measurements were taken before and after the whole measurement sequence.
The results can be seen in Figs. 9 and 10. We only show the results for 685 nm, as the salient features are similar at 830 nm. The leftmost graphs (a) show that the recovered absorption and scattering coefficients change linearly with ink and milk concentrations, respectively. The closeness of the lines also shows that the recovered absorption coefficients are largely independent of milk concentration, and that recovered scattering coefficients are largely independent of ink concentrations. The middle graphs (b) show again the independence of recovered reduced scattering coefficients with changing ink concentration, and recovered absorption coefficients with changing milk concentrations. We observe that our results are close to the ideal.
The blue lines in Figs. 9(c) and 10(c) show the -intercept of a linear fit to each of the curves in Figs. 9(a) and 10(a), respectively. Ideally the blue curve in Fig. 9(c) would be constant at the level of water absorption (), because at the -intercept the ink concentration is zero, and thus the only remaining absorber is water. We can see that the reported values are nearly independent of milk concentration, but slightly above the ideal value of . One possible explanation is contamination of the mixture with another absorber.
The blue line in Fig. 10(c) extrapolates the recovered scattering coefficient at zero milk concentration for different ink concentrations. We would expect this value to be zero, because without milk the mixtures do not contain any scatterer. The results are fairly close to the ideal, and show only a small positive correlation with ink concentration. The slight increase in scattering coefficient with increased ink could be explained by scattering from the ink particles.
The red line in Fig. 9(c) shows the slope of the linear fit to the lines in Fig. 9(b), normalized by the slopes of the linear fits to the lines shown in Fig. 9(a). Ideally the slopes of the lines in Fig. 9(b) would be zero, representing zero crosstalk from absorption to scattering coefficients. The red line in Fig. 10(c) shows the reciprocal measurement. We can see that the crosstalk generally remains below 10%.
Of note, we observed that the values of -intercepts and crosstalk depend strongly on the assumed optical properties of the calibration phantom. We measured these properties with our commercial ISS Imagent FD-NIRS system,28 but we assume that the Imagent measurements have a precision of to 10%. We therefore optimized the assumed optical properties of the phantom within that range to minimize -intercepts and crosstalk. The optimized calibration phantom optical properties were consistent between multiple titration datasets (data not shown).
In the previous sections, we showed that we successfully built a direct sampling multiwavelength FD-NIR system that performs close to the noise limit of the built in APD modules, while being significantly faster than previous FD-NIR DOT devices due to frequency division multiplexing, and improvements in the speed of the optical multiplexing system. The modular nature of the instrument would easily allow for additional detectors or sources to be fitted if required. The employed direct digital sampling and digital demodulation is independent of the particular detector used. Therefore, if an even lower noise floor is required, photo multiplier tubes could be used instead of APDs, e.g., at the disadvantage of higher cost, potentially smaller dynamic range and the introduction of amplitude–phase crosstalk.29
We also developed an RF signal decoding algorithm which employs only adders in the high speed section of the filter structure, thus allowing the use of low cost FPGAs. If required, it would easily be possible to increase the number of decoded frequencies to at least four, thus allowing a four wavelength instrument, or alternatively in a two wavelength system increase the speed by having two active RF sources at any time.
If even more frequencies need to be decoded, or if frequency sweeping is required, a digital implementation of a more traditional demodulation scheme, e.g., homodyne or heterodyne detection would have to be implemented. For this purpose, a larger and more expensive FPGA might be necessary.
Our system characterization data demonstrate an instantaneous dynamic range of with 60 dB SNR at a 90-Hz bandwidth, matching or exceeding other reported FD-NIRS diffuse optical spectrometers and imagers.4,23,30–32 The combination of these two parameters enables fast acquisition of several channels of data in several seconds of transmission-geometry breast optical tomography that in turn allows monitoring of tissue hemodynamics driven by intrinsic and extrinsic factors. Furthermore, the fast switching time of the galvo multiplexing system requires us to discard only two samples per transition between illuminated source positions. Thus if one aims for a 50% overall duty cycle (individual sources will have obviously a lower duty cycle), all the data from all possible combinations of sources and detectors can be acquired in 1.16 s. In our particular application, we are aiming for a duty cycle above 90%, and hence chose settings that result in a 10-s cycle time and 94% duty cycle.
We built a high-performance multiwavelength FD-NIR imager that achieves high temporal resolution while matching the performance of other slower instruments in the field. The receiving side of the instrumentation uses a modular design that can be easily expanded in the future to increase the number of detection channels. Adding wavelengths and source positions is also fairly easy, albeit in the case of source positions at the penalty of decreased duty cycle. The devices have been integrated into our next generation dynamic breast optical tomography system and, together with a fast CW-NIR imager offers high speed imaging of breast tissue hemodynamics.
This work has been funded by National Institutes of Health grants R01-CA142575, R01-CA097305, R01-CA187595, and R00-EB011889. The authors thank Joe Stadtmiller, Robert Dewsnap, Ron Altman, William Johnson, and Arthur DiMartino at TechEn Incorporated for collaborative efforts in developing the instrument.
Bernhard B. Zimmermann is currently a PhD candidate in the Department of Electrical Engineering and Computer Science at Massachusetts Institute of Technology, working on diffuse optical imaging techniques and instruments for breast cancer detection and therapy monitoring, as well as brain perfusion monitoring. He received his BS and MS degrees in Electrical Engineering at the Swiss Federal Institute of Technology in Zürich, Switzerland.
Qianqian Fang is currently an assistant professor in the Department of Bioengineering at Northeastern University. He received his PhD degree in biomedical engineering in 2005 at Dartmouth College, and postdoctoral training at Massachusetts General Hospital (MGH) between 2005 and 2009. He became an instructor in 2009 and assistant professor of radiology in 2012 at MGH. His current research interests include multimodality imaging, translational optical tomography and spectroscopy, point-of-care devices, and high-performance computational imaging methods.
David A. Boas is the director of the Martinos Optics Division in the Department of Radiology at Massachusetts General Hospital, and a professor of radiology at Harvard Medical School. He received his BS degree in physics at RPI and PhD degree in physics at the University of Pennsylvania. He is the founding president of the Society for Functional Near Infrared Spectroscopy, founding Editor-in-Chief of Neurophotonics, and the recipient of the 2016 Britton Chance Award in Biomedical Optics.
Stefan A. Carp is currently an assistant professor of radiology at Harvard Medical School and a Member of Massachusetts General Hospital Martinos Center Optics Division. He received BS degree in chemistry and chemical engineering from MIT and a PhD degree in chemical engineering from the University of California, Irvine. His research interests include the development of optical imaging tools for breast cancer management, as well as novel techniques for tissue perfusion and oxygen metabolism monitoring.