2.1.

System Overview

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.

Fig. 1

System overview: (a) FD-NIRS source and (b) detector components integrated into our hybrid CW/FD-NIRS dynamic breast imager.

2.2.

Hardware

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.

Fig. 2

Schematic overview of the FD-NIRS instrument.

2.2.1.

Source component

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 $200\text{-}\mu \mathrm{m}$ 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 $600\mu \mathrm{s}$ 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.

Fig. 3

Complete source unit assembly: (a) source box contains clock generator, laser drivers, lasers, and galvo for optical multiplexing; (b) detector card.

2.3.

Demodulation Algorithm

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.¹

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 $N$ 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 $X[k]$ is computed as shown in Eq. (1), where $x[r]$ is the input signal, $W_N=e^{-j(2\pi /N)}$, and k denotes the frequency bin to be computed.

To simplify we turn to the Goertzel algorithm. As shown in Ref. 24, the $k$′th bin of the $N$-point DFT can be computed by feeding the signal into a system with impulse response $W_N^{-kn}u[n]$ which is initially at rest (i.e., all registers are initialized to zero). $u[n]$ is the unit step function. The desired result is then the $N$’th output value [see Eqs. (2) and (3)] of the system. See also flow graph in Fig. 4.

Fig. 4

Flow graph of the Goertzel algorithm.

Fig. 5

Flow graph of the Goertzel algorithm with an additional delay in the recursive loop.

The transfer function of the system of Fig. 5 is shown in Eq. (4) for the case when $l=2$. 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

Flow graph of simplified algorithm.

Fig. 7

Final implementation of the filter structure for the first modulation frequency.

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 $N=4·{10}^6$ 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.

3.1.

System Characterization

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 $<1.4\mathrm{pW}/\surd \mathrm{Hz}$, which is approaching the manufacturer specified noise floor of the APD module of $0.8\mathrm{pW}/\surd \mathrm{Hz}$.

Fig. 8

System characterization. (a) Typical plot of signal magnitude versus incident optical power. One second integration time per power level. Blue corresponds to the driven detector, while the signals from the three neighboring detectors are shown in red. Asymptotes are drawn in gray for the driven detector. The asymptotes intersect at 1.01 pW where the SNR is 1. (b) Typical plot of signal phase versus incident optical power. One second integration time per power level. The gray line marks the optical power where the amplitude SNR is 1.

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 $>100\mathrm{dB}$ (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 $>80\mathrm{dB}$. 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 $\sim 1.2\mu \mathrm{W}$ 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 $1.5\mu \mathrm{W}$ signal power.

We were not able to observe any interwavelength crosstalk nor amplitude to phase crosstalk up to the maximum specified input power of $1.5\mu \mathrm{W}$. The phase noise of the output signal is smaller than $6\mathrm{mrad}/\surd \mathrm{Hz}$ 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 $<1.5\%$, the phase $<3\mathrm{mrad}$ at an optical power of 5 nW. Our findings are summarized in Table 1.

Table 1

Summary of performance tests.

3.2.

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.²⁵ 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, $\sim 1\mathrm{mm}$ 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.

Fig. 9

Absorption stepping. (a) Recovered absorption coefficient versus ink concentration. The different lines represent different milk concentrations (the lowest line corresponds to the lowest milk concentration). (b) Recovered reduced scattering coefficient versus ink concentration. The different lines again represent different milk concentrations, with the lowest milk concentration corresponding to the lowest line. (c) The blue line shows the $y$-intercepts of linear fits to the curves in (a). Ideally, it would be constant at the level of water absorption ($0.0049{\mathrm{cm}}^{-1}$ at 685 nm). The red line shows the relative crosstalk from recovered absorption coefficient to recovered reduced scattering coefficient at different milk concentrations, derived from the slopes and mean magnitudes of the lines in (a) and (b).

Fig. 10

Scattering stepping. (a) Recovered reduced scattering coefficient versus milk concentration. The different lines represent different ink concentrations. (b) Recovered absorption coefficient versus milk concentration. The different lines again represent different ink concentrations, with the lowest ink concentration corresponding to the lowest line. (c) The blue line shows the $y$-intercepts of linear fits to the curves in (a). Ideally, it would be constant at zero. The red line shows the relative crosstalk from recovered reduced scattering coefficient to recovered absorption coefficient at different ink concentrations, derived from the slopes and mean magnitudes of the lines in (a) and (b).

The blue lines in Figs. 9(c) and 10(c) show the $y$-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 ($0.0049{\mathrm{cm}}^{-1}$), because at the $y$-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 $0.0049{\mathrm{cm}}^{-1}$. 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 $y$-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,²⁸ but we assume that the Imagent measurements have a precision of $\sim 5\%$ to 10%. We therefore optimized the assumed optical properties of the phantom within that range to minimize $y$-intercepts and crosstalk. The optimized calibration phantom optical properties were consistent between multiple titration datasets (data not shown).