Quantitative comparison of analysis methods for spectroscopic optical coherence tomography

Spectroscopic optical coherence tomography (sOCT) enables the mapping of chromophore concentrations and image contrast enhancement in tissue. Acquisition of depth resolved spectra by sOCT requires analysis methods with optimal spectral/spatial resolution and spectral recovery. In this article, we quantitatively compare the available methods, i.e. the short time Fourier transform (STFT), wavelet transforms, the Wigner-Ville distribution and the dual window method through simulations in tissue-like media. We conclude that all methods suffer from the trade-off in spectral/spatial resolution, and that the STFT is the optimal method for the specific application of the localized quantification of hemoglobin concentration and oxygen saturation. ©2013 Optical Society of America OCIS codes: (030.1640) Coherence; (070.4790) Spectrum analysis; (160.4760) Optical properties; (170.6510) Spectroscopy, tissue diagnostics. References and links 1. N. Bosschaart, T. G. van Leeuwen, M. C. G. Aalders, B. Hermann, W. Drexler, and D. J. Faber, “Spectroscopic low coherence interferometry”, Chapter 23 in Optical Coherence Tomography – Technology and Applications, W. Drexler and J.G. Fujimoto, eds. (Springer Berlin Heidelberg New York USA), 2nd edition (2013) 2. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Differential absorption imaging with optical coherence tomography,” J. Opt. Soc. Am. A 15(9), 2288–2296 (1998). 3. U. Morgner, W. Drexler, F. X. Kärtner, X. D. Li, C. Pitris, E. P. Ippen, and J. G. Fujimoto, “Spectroscopic optical coherence tomography,” Opt. Lett. 25(2), 111–113 (2000). 4. C. Xu, J. Ye, D. L. Marks, and S. A. Boppart, “Near-infrared dyes as contrast-enhancing agents for spectroscopic optical coherence tomography,” Opt. Lett. 29(14), 1647–1649 (2004). 5. F. E. Robles, C. Wilson, G. Grant, and A. Wax, “Molecular imaging true-colour spectroscopic optical coherence tomography,” Nat. Photonics 5(12), 744–747 (2011). 6. N. Bosschaart, M. C. G. Aalders, D. J. Faber, J. J. A. Weda, M. J. C. van Gemert, and T. G. van Leeuwen, “Quantitative measurements of absorption spectra in scattering media by low-coherence spectroscopy,” Opt. Lett. 34(23), 3746–3748 (2009). 7. N. Bosschaart, D. J. Faber, T. G. van Leeuwen, and M. C. G. Aalders, “In vivo low-coherence spectroscopic measurements of local hemoglobin absorption spectra in human skin,” J. Biomed. Opt. 16(10), 100504 (2011). 8. N. Bosschaart, D. J. Faber, T. G. van Leeuwen, and M. C. G. Aalders, “Measurements of wavelength dependent scattering and backscattering coefficients by low-coherence spectroscopy,” J. Biomed. Opt. 16(3), 030503 (2011). 9. J. Yi and V. Backman, “Imaging a full set of optical scattering properties of biological tissue by inverse spectroscopic optical coherence tomography,” Opt. Lett. 37(21), 4443–4445 (2012). 10. C. P. Fleming, J. Eckert, E. F. Halpern, J. A. Gardecki, and G. J. Tearney, “Depth resolved detection of lipid using spectroscopic optical coherence tomography,” Biomed. Opt. Express 4(8), 1269–1284 (2013). 11. R. N. Graf and A. Wax, “Temporal coherence and time-frequency distributions in spectroscopic optical coherence tomography,” J. Opt. Soc. Am. A 24(8), 2186–2195 (2007). 12. F. Robles, R. N. Graf, and A. Wax, “Dual window method for processing spectroscopic optical coherence tomography signals with simultaneously high spectral and temporal resolution,” Opt. Express 17(8), 6799–6812 (2009). 13. C. Xu, F. Kamalabadi, and S. A. Boppart, “Comparative performance analysis of time-frequency distributions for spectroscopic optical coherence tomography,” Appl. Opt. 44(10), 1813–1822 (2005). 14. R. Leitgeb, M. Wojtkowski, A. Kowalczyk, C. K. Hitzenberger, M. Sticker, and A. F. Fercher, “Spectral measurement of absorption by spectroscopic frequency-domain optical coherence tomography,” Opt. Lett. 25(11), 820–822 (2000). #195497 $15.00 USD Received 12 Aug 2013; revised 28 Sep 2013; accepted 15 Oct 2013; published 23 Oct 2013 (C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002570 | BIOMEDICAL OPTICS EXPRESS 2570 15. A. Wax, C. Yang, and J. A. Izatt, “Fourier-domain low-coherence interferometry for light-scattering spectroscopy,” Opt. Lett. 28(14), 1230–1232 (2003). 16. L. Cohen, “Time-frequency distributions – a review,” Proc. IEEE 77(7), 941–981 (1989). 17. Data tabulated from various sources compiled by S. Prahl, http://omlc.ogi.edu/spectra. 18. N. Bosschaart, M. C. G. Aalders, T. G. van Leeuwen, and D. J. Faber, “Spectral domain detection in lowcoherence spectroscopy,” Biomed. Opt. Express 3(9), 2263–2272 (2012). 19. J. Yi and X. Li, “Estimation of oxygen saturation from erythrocytes by high-resolution spectroscopic optical coherence tomography,” Opt. Lett. 35(12), 2094–2096 (2010).


Introduction
Spectroscopic optical coherence tomography (sOCT) and its closely related brother lowcoherence spectroscopy (LCS) are exciting extensions of low coherence interferometry (LCI) that enable the recovery of depth-resolved spectra from tissue [1].sOCT and LCS not only enhance the contrast of OCT images [1][2][3][4][5], but also allow for the quantification of depth-and wavelength resolved optical properties, such as the absorption [5][6][7], scattering and backscattering coefficient [8,9].These optical property spectra provide valuable information on the physiological status of tissue -for instance on the spatial distribution of hemoglobin concentration and oxygen saturation [7] and the identification of lipid in atherosclerotic plaques [10].
Methods that have been applied in practice for the extraction of spectral information from LCI signals are the short-time Fourier transform (STFT) [4,[6][7][8][9][10], wavelet transforms [3], the Wigner-Ville distribution (WV) [11] and the dual window method (DW) [5,12].In general, all methods aim to recover spectra from distinct depth intervals in the LCI data set.Important aspects in the analysis are 1) the spatial resolution at which spectra are recovered -i.e. the size of the depth interval; 2) the spectral resolution at which spectra are recovered, determining the resolvability of spectral features; and 3) the spectral recovery -i.e. the correct recovery of spectral amplitude and shape.The quality of the available analysis methods has been primarily evaluated in terms of the first two aspects [13], since optimization and minimization of the trade-off in spectral/spatial resolution is a well-known challenge in sOCT and LCS.However, the third aspect (spectral recovery) is equally, if not more important to evaluate when sOCT or LCS are used for the measurement of depthresolved optical properties and the quantification of chromophore concentrations from it.As we will show in this article, small changes in the shape and amplitude of the recovered absorption spectrum can lead to large errors in the derived chromophore concentrations.The same problem may occur when spatially resolved spectra of the scattering and backscattering coefficient are used to derive the micro-and nano-scale organization of the imaged tissue [8,9].
In this article, we therefore quantitatively compare the performance of the available analysis methods in terms of all three aspects, focusing specifically on the measurement of depth-resolved absorption spectra of hemoglobin.Hereto, we perform simulations on tissue simulating media with multiple reflectors and physiologically relevant concentrations of oxyand deoxyhemoglobin.In order to quantify the spectral recovery, we analyze the mean squared difference (χ 2 ) between the input and recovered absorption coefficient spectra, and we evaluate the quantification of the total hemoglobin concentration and oxygen saturation.

Theory
In a standard LCI system -with sample field E S originating from the sample arm and reference field E R originating from the reference arm -the general expression of the LCI interferogram is written as: with wavenumber k = 2π/λ, wavelength λ and optical path length difference 2z, such that z/n is the assigned depth location in the tissue with refractive index n.For simplicity, we assume n = 1 in the remaining part of this paper.We start with OCT data in spatial (z-) domain, either directly acquired ("time domain", TD-OCT) or the result of signal processing after Fourierdomain (FD-OCT) detection.We assume real-valued representation of the OCT interferograms (i.e.including the modulation term).Together, the wavenumber k and optical path length difference 2z form the fundamental Fourier pair in LCI data analysis: where ℑ denotes the Fourier transform.Since the wave number k is directly related to wavelength λ, wavelength dependent spectra S(λ) can be obtained from the backscattered LCI signal.
The goal of sOCT and LCS is to obtain information in the λ-and z-domain simultaneously, with high resolution in both λ and z.Due to wavelength-dependent scattering and absorption by the different structures in tissue, the spectral content of the LCI signal changes with depth.If we would directly apply the Fourier transform on the LCI signal (Eq.( 2), either the depth, or wavelength varying information will be lost -depending on time, or spectral domain signal acquisition, respectively.Instead, depth-resolved spectroscopic information can be obtained using time-frequency (TF) analysis methods (or, for the present context: spectral/spatial analysis).The four accepted methods for spectral-spatial analysis in sOCT and LCS will be described next.We point out that these methods can be performed equivalently on OCT data that is acquired in the spatial domain (described here) or in the spectral domain [14,15].

Short time Fourier transform (STFT)
The short-time Fourier transform (STFT) is the most widely applied spectral/spatial analysis method: where w(z,Δz) is an analysis window confined in space around z with spatial width Δz, for example a Gaussian function.The multiplication with a relatively short window effectively suppresses the signal outside the analysis point z ± ½Δz.Physically, the STFT can be considered as the result of passing a signal through an array of band-pass filters with linearly increasing center frequency and constant bandwidth which is inversely proportional to Δz.Thus, there is an inherent trade-off between spectral and spatial resolution.A window with short spatial width Δz will localize the signal well in space but will have reduced k-resolution; conversely a signal with long spatial width will be less well localized in space with the benefit of increased spectral resolution.For a Gaussian window, a spatial domain width Δz will yield a spectral (wavenumber) resolution of Δk = 1/(2Δz), or equivalent wavelength resolution Δλ: In practice, starting with a real-valued time-domain A-scan, a spatial window w z with center at z and width Δz is chosen.This window is multiplied with the A-scan and subsequently Fourier transformed with z→k (Eq.( 3).Repeating for every z results in a complex valued STFT(k,z; w z ) spectrogram with z-axis equal to that of the time-domain A-scan, and k-axis from k = 0 to k = π/2δz where δz is the sampling increment of the A-scan.The resolution in the spatial domain is given by de window width Δz; the resolution in the spectral domain is Δk = 1/2Δz assuming a Gaussian window (Eq.(4 the complex spectrogram STFT(k,z; w k ).The spectral range equals that of the interference spectrum; the spatial range runs from z = 0 to z = π/2δk where δk is the spectral domain sampling interval.For further processing the amplitude of the complex spectrogram is taken.

Wavelet transforms
The wavelet transform was introduced to partially overcome the trade-off in spectral/spatial resolution of the STFT, by adjusting the window size to the frequency being considered.The basic difference between the wavelet transform and the STFT is that the spatial width and the spectral bandwidth of the wavelet are both changed, while its shape remains the same.Physically, the wavelet transform can also be regarded as an array of band-pass filters with constant relative bandwidth with respect to the center frequency.As a consequence, the wavelet transform uses short spatial windows at high frequencies, and long spatial windows at low frequencies.Whereas the trade-off in spectral/spatial resolution remains, it now depends on frequency: the spectral (resp.spatial) resolution becomes poorer (resp.better) with increasing analysis frequency.
In the practice of sOCT, wavelet transforms have only been applied using Gaussian windows with fixed bandwidth [3], which is essentially the same as the STFT.

Wigner-Ville distribution (WV)
Bi-linear spectral/spatial distributions do not suffer from the resolution trade-off between both domains.The most important member of this class is the Wigner-Ville distribution [14]: The Wigner-Ville distribution is the Fourier transform of an autocorrelation measure of the signal i D (z).Whereas in conventional autocorrelation computations the result is a function of lag (z-z') only, here the dependence on z is maintained.In concordance with the Wiener-Khinchine theorem, the straightforward physical interpretation of the WV is as a localized power spectral density of the detector signal.The result of Eq. ( 5) is complex; for further processing the absolute value is taken.The drawback of the WV lies in its quadratic nature: whereas the STFT of two signals X and Y yields STFT x (k,z) + STFT y (k,z), the WV will contain interference terms i.e.WV x+y (k,z) = WV x (k,z) + WV y (k,z) + 2Re[WV x,y (k,z)].Even though the interference term contains information on the separation of X and Y (i.e. in space), when it overlaps with the signal terms, interpretation of the WV becomes challenging.In practice, the signal is analyzed using a window w, similar to the STFT.This Pseudo-WV effectively smoothes the spectral/spatial distribution, suppressing the interference terms.However, a short spatial smoothing window will be wide in frequency, leading to a good spatial resolution but bad spectral resolution and vice versa.It is possible to add a degree of freedom by considering a separable smoothing function Δ(k,z) = w k (k,Δk)w z (z,Δz) that allows independent control in both depth and frequency of the smoothing applied to the WV; Δk and Δz denote the width of the window functions in the spectral and spatial domain, respectively.The STFT compromise between spatial and spectral resolution is now replaced by a compromise between the joint spectral/spatial resolution and the level of suppression of the interference terms.

Dual window method (DW)
Robles et al. showed the pseudo WV can also be obtained by STFT analysis using two kdomain window sizes Δk 1 >>Δk 2 [12] (or equivalently: two z-domain window sizes Δz 1 <<Δz 2 ).The result STFT 1 (k,z) × STFT 2 (k,z) is mathematically equivalent to a pseudo WV with window widths of w k (k,Δk = Δk 2 /2) in the spectral domain and window w z (z,Δz = 1/(2Δk 1 )) in the spatial domain.As for the WV, any remaining interference terms may carry information, for example on average spatial scatterer separation.Robles et al state that the DW method is superior to the STFT in terms of the trade-off in spectral/spatial resolutioni.e. the spatial resolution equals the width of the narrow spatial window, and the spectral resolution equals the width of the narrow spectral window.When applied in practice, interference terms between nearby reflectors (both located within one long spatial window) appear as frequency modulations on the recovered spectra.These modulations can be removed by low-pass frequency domain filtering [12].

Methods
We quantitatively compare the methods described in Section 2 in terms of spatial resolution, spectral resolution and spectral recovery through computational simulations in Labview 2011.The STFT and the practical implementation of wavelet transforms are essentially the same.The same holds for the (pseudo) WV distribution and the DW method.Therefore, we limit our comparison to the STFT and the DW method.Our simulations are performed in the context of time domain data acquisition.Results are evaluated in both the time and frequency domain; therefore the outcome will be analogous for spectral domain data acquisition.

Simulations
We start by generating a spatial domain OCT A-scan with 4096 data points and a maximum depth z max of 375 µm.In the A-scan, we introduce two point reflectors with a spacing of D = 25 µm: R1 at 220 µm and R2 at 245 µm (Fig. 1(a)), representing the walls of a blood vessel.The spectral reflection S R1 (k) of R1 (Fig. 1(b): actual) is defined by a Gaussian 'source' spectrum, centered around k 0 = 11.5 µm −1 (corresponding to wavelength λ 0 = 2π/k 0 = 546nm) on a wavenumber range k = 8-15 µm −1 (λ = 419-785nm) with a full width half max of k FWHM = 2.35 µm −1 (λ FWHM = 113nm).The spacing between R1 and R2 is simulated as an absorbing blood volume or vessel, with a total hemoglobin concentration of [tHb] = 150 g/L, an oxygen saturation of SO 2 = 85% and zero scattering.We can therefore define the spectral reflection S R2 (k) of R2 (Fig. 1(c): actual) as the source spectrum S R1 (k) from R1, multiplied with the influence of absorption from the blood volume: , where the absorption coefficient spectrum: is obtained from the absorption spectra of deoxygenized (Hb) and oxygenized (HbO 2 ) hemoglobin in g/L [17].The factor 0.01 accounts for the definition of SO 2 as a percentage.Subsequently, we convolute the point reflectors with an OCT interferogram that has a coherence length corresponding to the spectral content at R1 and R2 (l c ≈1.3 µm) and calculate the actual A-scan envelope (Fig. 1(a): actual).
To compare the performance of the STFT to the DW method, we define a medium STFT window w M with Δz M = 22 µm in the spatial domain -as we used in our previous work [6][7][8] to quantify hemoglobin absorption in human skin -and a short (w S ) and a long (w L ) window for the DW method -comparable to those used by Robles et al [5] to quantify hemoglobin absorption in an animal model.Table 1 summarizes the window sizes in all domains and the resulting wavelength resolutions (Eq.( 4).Since the wavenumber resolution is linear in k, the resulting wavelength resolution is higher at the shorter wavelengths compared to the longer wavelengths.According to the definition in Section 2, the resolutions for the STFT are 22 µm in depth and 6.9 nm (at λ 0 ) spectrally, whereas the resolutions for the DW method are 5 µm in depth and 2 nm (at λ 0 ) spectrally.For both methods, the spectral/spatial (i.e.time-frequency) distributions are calculated and A-scans are obtained by taking the cross section through these distributions in the depth direction (Fig. 1a).The spectra at R1 and R2 are obtained by taking the cross section through the spectral/spatial distributions in the wavenumber direction at the corresponding depth positions (Fig. 1(b) and Fig. 1(c)).Figure 1(a) also shows the resulting A-scans and recovered spectra for the STFT's with the individual windows w L and w S .
Since the long window w L of the DW method probes both R1 and R2 simultaneously, a frequency modulation is visible in the DW spectra due to the interference between both reflectors.When taking the Fourier transform of the DW and STFT spectra in Fig. 1(c), the modulation frequency appears exactly at the reflector spacing of 25 µm for the DW (Fig. 1(d)).Hence, this modulation can be removed efficiently by applying an (ideal) low-pass filter with a cut-off frequency of 7 µm on the Fourier transformed spectra.Note that this cutoff frequency is higher than the previously reported cut-off frequency of 3 µm for the DW method [12], since frequencies lower than 7 µm contain contributions from the features of the hemoglobin absorption spectrum that need to be preserved.The filtered spectra are shown in Fig. 1(e) and Fig. 1(f) for R1 and R2, respectively.
The absorption coefficient µ a of the simulated blood volume between R1 and R2 is recovered using Beer's law: with S R1 (k) and S R2 (k) the simulated spectra at R1 and R2 from both methods.

Simulation cases
We will describe four simulation cases for which we evaluate the recovery of µ a : Case 1: Simulation of the two reflectors R1 and R2, as described in Section 3.1.The aim of this case is to recover µ a in an 'ideal' situation, i.e. without the influence of other reflectors.
Case 2: Simulation of three reflectors, with a third reflector R3 introduced at depths larger than R2.R3 is simulated as a strong reflector without any absorption between R2 and R3 (magnitude S R3 = 3.5 × magnitude S R2 ).The distance between R3 and R2 is varied from 0 -100 µm.In this case, we also investigate the influence of noise on the recovery of µ a for 20 µm and 100 µm separation distance.
Case 4: Identical to case 2, but with moderate power for R3 (S R3 = S R2 ) and varying size of the long DW window w L (Δz L = 5 -75 µm).

Quantification of performance
For our simulation cases, we define three measures to quantify the performance of both methods in terms of how well they recover the µ a of the simulated blood volume: 1) The normalized chi-squared χ 2 value of µ a recovery.The χ 2 value of µ a is obtained by: where µ a,act is the actual (input) µ a and µ a,sim is the recovered µ a by either the DW, or the STFT method.All χ 2 values are obtained between 475 and 650 nm, hence λ 1 = 475 nm and λ N = 650 nm.Although the wavelength resolutions differ between the DW and STFT method (Table 1), N is equal for both methods in the 475-650 nm range, due to zeropadding in the calculation of the spectral/spatial distributions.Equation (6) shows that the lower χ 2 , the better the recovery of µ a .
2) The effective spatial resolution.It is expected that χ 2 will increase when a third reflector R3 (at a distance >R2) approaches R2 towards a distance that is similar to the spatial resolution of both methods, since S R2 (k) cannot be recovered independently of S R3 (k) in that case.We define the effective spatial resolution as the distance between R2 and R3 where χ 2 deviates more than 15% from the χ 2 without the presence of R3.(The allowed deviation can be chosen dependent on application).
3) The correct recovery of the actual total hemoglobin concentration [tHb] and oxygen saturation (SO 2 ) from µ a .The [tHb] and SO 2 are recovered from the simulated µ a spectra through a non-linear least squares chromophore fit of the model of Eq. ( 7) to the simulated µ a with fit parameters [tHb] and SO 2 .
Fig. 2. Recovered absorption spectra µ a from the blood volume between R1 and R2 for the DW and STFT method (Case 1).The recovered µ a spectrum from the blood volume is shown in Fig. 2 for both methods, using the filtered spectra of R1 and R2.Compared to the STFT, the spatial recovery is much better for the DW method (Fig. 1(a)).However, the recovered spectra by the DW method deviate stronger from the actual spectrum than the spectra recovered by the STFT (Figs. 1(e), 1(f) and Fig. 2), which is reflected in the values for χ 2 .Table 2 summarizes the performance of the STFT and DW method for the spectral recovery of R1 and R2 and the µ a in between.Lowpass filtering of the spectra enhances the χ 2 values for the DW method, but it remains higher with 1-2 orders of magnitude than the χ 2 values for the STFT method.Since the STFT method does not suffer from the influence of the interference between both reflectors for this window size and reflector separation, no difference is found between the χ 2 values for the unfiltered and filtered spectra of R1 and R2 within the presented significance (Table 2).The slight increase in the χ 2 value of µ a after filtering is attributed to the minor loss of frequencies containing hemoglobin-specific absorption features above 7 µm (Fig. 1(d)).Whereas the recovered [tHb] is only slightly lower than the actual [tHb] of 150 g/L for both methods (STFT error: 3.2 g/L, DW error: 7.0 g/L), the actual SO 2 of 85% is severely underestimated by both the STFT (with 18%) and the DW method (with 40%).This underestimation results from the flattening of the HbO 2 absorption peaks due to spectral smoothing by the finite window sizes.Figure 3 shows the dependency of the fitted SO 2 (Fig. 3(a)), fitted [tHb] (Fig. 3(b)), and χ 2 of µ a (Fig. 3(c)) on the actual (input) value of the SO 2 .Since the HbO 2 peaks become sharper for increasing SO 2 values, the amount of peak flattening increases with SO 2 .As a consequence, the underestimation of the fitted SO 2 and [tHb] increases, as well as the χ 2 of µ a .Due to the consistent linear behavior of the fitted SO 2 and [tHb] on the actual SO 2 , the fitted values in practice can be recalculated to the actual values if this relation is known -e.g. from a simulation.For the STFT method, this correction can also be applied directly by smoothing of the µ a,HbO2 and µ a,Hb spectra in the chromophore fit with the spectral resolution of the STFT [18].

Case 2: three reflectors, varying separation
This case describes the influence of a third reflector R3 (S R3 = 3.5 × S R2 ) on the determination of the µ a of the blood volume between R1 and R2.It is expected that R3 will influence this determination if R3 approaches R2 within the spatial resolution of the analysis methods (within 22 µm for the STFT, and 5 µm for the DW). Figure 4 shows a situation where R3 is located at a distance of 20 µm from R2 (Fig. 4(a)).Although this is well outside the spatial resolution of the DW method, the recovered µ a between R1 and R2 by the DW method is highly influenced by the presence of R3 (Fig. 4(b)).On the other hand, the recovered µ a by the STFT remains practically unaffected and agrees well with the actual (input) spectrum of µ a .Also for this case, the A-scan recovery of the DW method is superior to the STFT.Nevertheless, the spectral performance of the DW method is expectedly compromised by poorer spectral recovery when the large reflector R3 is probed simultaneously with R2 by the long spatial DW window w L .
We added a movie file to Figs. 4(a) and 4(b) (Media 1) that shows the effect of the position of R3 (300 to 245 µm) on the recovered µ a between R1 and R2. Figure 5 shows the influence of R3 on the determination of the µ a of the blood volume between R1 and R2 for several more separations of R2 and R3 (345 to 245 µm).5(c) on the fitted SO 2 .For all parameters, a deviation from the baseline value is observed at separations smaller than ~16 µm for the STFT, and smaller than ~30 µm for the DW method.Although the χ 2 of µ a for the DW method first decreases for separations <~30 µm (indicating better spectral recovery), the fitted values of the [tHb] and SO2 are also affected at these separations -indicating worse spectral recovery.This is caused by the introduction of frequency modulations on the recovered µ a (see Fig. 4, Media 1) for separations <~30 µm, which 'deceptively' results in a better χ 2 of µ a .For this case, the effective spatial resolution (defined as a deviation in χ 2 of >15% from the χ 2 without the presence of R3) of both methods is therefore 16 µm for the STFT, and 30 µm for the DW method -compared to the theoretical resolutions of 22 µm and 5 µm, respectively.
To investigate the influence of noise on the χ 2 of μ a recovery, we added uniform white noise to the simulated signals as illustrated in Fig. 4(c) and Fig. 4(d).We varied noise powers, represented as a signal-to-noise ratio, here defined as: amplitude R1 20 log , stand.dev. of noise Figure 6(a) shows the result for a separation between R2 and R3 of 100 μm, e.g. more than the large DW window size w L .In this case, the influence of R3 is negligible.As can be expected, χ 2 increases with decreasing SNR for both methods.Performance of the DW method based on the metric χ 2 is better than STFT under low SNR conditions.When R3 is positioned within the large DW window (separation 20 μm; Δz S < 20 μm < Δz L ), the performance of the STFT is expectedly better for low noise conditions (Fig. 6(b)).For low SNR, the STFT and the DW method appear to perform equally in this situation.In the practice of sOCT and LCS, the influence of noise on the measurement of optical property spectra is reduced by temporal and/or spatial averaging of the acquired spectra [7].

Case 3: three reflectors, varying separation and power
Since the spectral power of R3 is likely to influence the effective spatial resolution, we repeated Case 2 for a range of spectral powers of R3 (S R3 = α × S R2 , with α = 0.25-3.5). Figure 7 shows the resulting effective spatial resolutions as a function of the power ratio α between R3 and R2.As expected, the effective spatial resolution becomes smaller (i.e.better) for lower power reflected from R3.When the power of R3 is very small compared to the power of R2 (α = 0.25), the effective spatial resolution of the DW method (8 µm) is close to the theoretical spatial resolution of 5 µm.Although the effective spatial resolution of the DW is better than that of the STFT for this power ratio α, the spectral recovery for the STFT method remains better with a χ 2 of 0.15, compared to a χ 2 of 3.7 for the DW.   7. Effective spatial resolution of the DW and STFT method as a function of the power ratio between R3 and R2 (Case 3).The effective spatial resolution is defined as the separation between R2 and R3 (see Fig. 4a), where the χ 2 of µ a between R1 and R2 deviates more than 15% from the χ 2 of µ a without the presence of R3.

Case 4: three reflectors, varying separation and window sizes
Since the effective spatial resolution of the DW method is larger than its theoretical spatial resolution (defined by the short spatial window Δz S ), we expect that the long spatial window (Δz L ) influences the effective spatial resolution.To test this hypothesis, we repeated Case 2 with moderate power for R3 (S R3 = S R2 ) and varying sizes of Δz L . Figure 8(a) shows that the effective spatial resolution of the DW method decreases as a function of the size of Δz L and only comes close to the theoretical spatial resolution of 5 µm when Δz L ≈Δz S .For window sizes Δz L smaller than approximately 45 µm (~2*Δz M ), the effective resolution of the DW method is better than that of the STFT.However, the χ 2 of µ a at this effective resolution is 3.5 for the DW method (i.e. at a separation between R2 and R3 of ~45 µm), whereas the χ 2 of µ a for the STFT is only 0.12 for the same separation between R2 and R3 -indicating better spectral recovery for the STFT.Comparable to Fig. 5(a), Fig. 8(b) displays the χ 2 of µ a as a function of the separation between R2 and R3 for the STFT method and for a selection of long spatial window sizes (Δz L = 5, 7.5, 15 and 75 µm) of the DW method.The χ 2 of the DW method is better (i.e.smaller) than that of the STFT method when the separation between R2 and R3 is smaller than ~7 µm, for all window sizes Δz L .However, the spectral recovery in this case is poor for both methods with a χ 2 of µ a around 10 and higher, resulting in deviations from the actual χ 2 of µ a comparable to the DW result in Fig. 4(b).
Note that the results presented in this case will depend on the spectral power of R3: for higher powers of R3, the effective spatial resolution will be larger (i.e.worse) at equal window size Δz L (see Case 3).

Discussion
In this article, we evaluated the methods for spectral/spatial (i.e.time-frequency) analysis that are used in the practice of sOCT and LCS.We showed that the STFT and wavelet transforms are essentially the same if the latter is applied using Gaussian windows with fixed bandwidths.Also the (pseudo) Wigner Ville distribution and the DW method are identical, as was shown in Ref [12].We therefore limited our evaluation to the comparison of the STFT to the DW method in terms of spectral/spatial resolution and spectral recovery.As was emphasized before by Xu et al [13], the optimal spectral/spatial analysis for sOCT and LCS depends on the application and imaging scheme that is applied.In this article, we focused on the quantification of optical property spectra (µ a ) and the application of localized measurements of hemoglobin concentrations [tHb] and oxygen saturation (SO 2 ) in tissue.Besides the spectral/spatial resolution, the spectral recovery of the analysis method is an important aspect for this application, since small changes in spectral shape and amplitude can lead to large errors in the parameters of interest.
From the simulation of the two reflectors enclosing a blood volume in Case 1, we learned that the STFT performs better than the DW method in terms of spectral recovery, based on the χ 2 values of the recovered spectra in Fig. 1(e), 1(f) and Fig. 2 (Table 2).As a consequence, the [tHb] and SO 2 of the blood volume that were derived from the DW spectra show more deviation from the actual values than the same parameters derived from the STFT spectra.Both the STFT and the DW method underestimate the actual SO 2 of the simulated blood volume (Fig. 3(a)), due to recovery of the spectra with a lower spectral resolution than the reference spectrum that was used for quantification of the SO 2 .Fortunately, this underestimation can be compensated for by spectral smoothing of the reference spectrum prior to the fitting procedure.
Interference terms occurred in both the STFT and DW spectra when more than one reflector was probed within one spatial window.Since the long spatial window Δz L of the DW method is larger than the STFT spatial window Δz M , the DW method suffers more from this artifact than the STFT.Low-pass filtering is an effective solution to this artifact, but is limited to a cut-off frequency of 7 µm for applications involving hemoglobin absorption, due the features of its absorption spectrum.Hence, when reflectors are spatially separated by less than 7 µm, a reliable recovery of hemoglobin absorption is not possible by any method if window sizes exceed 7 µm.We expect that this problem is of minor influence in real tissue, since multiple reflectors at varying separations will average out this effect.
Although the spatial recovery of the DW method -with a theoretical resolution that equals the size of the short spatial window Δz S -seems to be better than that of the STFT (Fig. 1(a)), the recovered spectra by the DW method are influenced by reflectors that are located outside the theoretical spatial resolution (Fig. 4).Hence, the effective spatial resolution of the DW method is worse than its theoretical spatial resolution.We explain this difference by the influence of the large spatial window Δz L on the effective spatial resolution (Fig. 8(a)) -i.e.Δz L and Δz S are not entirely decoupled in the spectral/spatial recovery.On the other hand, the effective spatial resolution of the STFT was better than its theoretical spatial resolution Δz M .This can be explained by the Gaussian shape of the STFT window, which favors those reflections that are probed in its center and attenuates those reflections that are probed in its tails.Note that the effective spatial resolution is a qualitative measure of resolution, since it depends on the relative power of the reflectors (Fig. 7), noise power and choice of maximum allowed χ 2 of µ a deviation of baseline.It is intended to facilitate comparisons between the STFT and DW methods, and is indicative of the spatial range over which spectral information is obtained.
We performed our simulations with window sizes that have been used in the practice of localized measurements of [tHb] and SO 2 in tissue with LCS (using STFTs: Δz M = 22 µm, as in Ref [7]) and with sOCT (using the DW method: Δz S = 5 µm and Δz L = 75 µm, as in Ref  [5]).In this article, we showed that the effective spatial resolution of the DW method becomes better for smaller window sizes of Δz L .For a power ratio of 1 between the second and third reflector (which is a reasonable scenario in tissue given the small spatial variations in refractive index), the effective spatial resolution of the DW method is better than that of the STFT when Δz L ≈2Δz M (Fig. 7(a)).However, in this case the spectral recovery (χ 2 of µ a ) for the DW method only becomes better than that of the STFT when Δz L ≈7 µm (i.e.Δz L ≈1/3Δz M, or Δz L ≈1.5Δz S , Fig. 8(b)) -although the spectral recovery for neither of the two methods is good enough for the quantitative interpretation of µ a in this situation.Hence, the DW method can perform better than the STFT, but only if the spatial distribution of tissue reflectivity is relatively homogeneous, and not for the condition that Δz L >> Δz S .
Adding noise to our simulations decreases the quality of spectral recovery for both the STFT and DW method (Fig. 6), but it does not alter our conclusions on the comparison between the STFT and DW method.In the practice of sOCT and LCS, temporal and/or spatial averaging of the acquired spectra reduces the influence of noise [7].One parameter that we did not include in our simulations is the influence of tissue scattering.Tissue scattering will lead to higher signal attenuation, which influences the quality of spectral recovery -similar to the influence of noise.Therefore, tissue scattering may influence the choice of optimal window size(s) [19].Also here, temporal and/or spatial averaging of the acquired spectra will enhance signal quality [7].In the presence of scattering, both sOCT and LCS will measure the total attenuation coefficient (µ t = µ a + µ s , with µ s the scattering coefficient) rather than the absorption coefficient only.Several models exist to accurately separate the scattering and absorption contributions to µ t and from that, the tissue chromophore concentrations [7,10].
In summary, the results of our simulations show that the STFT is a reliable and relatively straight forward method to use for the quantitative spectral/spatial analysis of optical properties in sOCT and LCS.Its results are predictable as a function of window size w M (Δz M ,Δk M ).Although we investigated its performance only for one size of w M , we showed that hemoglobin absorption spectra can be recovered accurately with this window within the expected spectral and spatial resolutions.Other sizes for w M may be more effective in other situations, e.g. for other wavelength ranges and other chromophores.This can be checked with simulations similar to those presented in this paper.The DW method requires more careful optimization in terms of window sizes, as we showed for the long spatial window size Δz L .Reducing the size of the small spatial window Δz S will further increase the quality of spatial (A-scan) recovery until Δz S equals the coherence length of the source spectrum.However, this will lead to a further decrease in spectral recovery.Since spectral shape is preserved quite well by the DW method (Fig. 4(b)), it is excellently suited for contrast enhancement and color coding of high resolution OCT images.However, based on our simulations, it is less well suited for the quantitative recovery of optical property spectra, because the spectral amplitude is less well preserved (Fig. 4(b)).
All simulation software can be obtained by contacting the authors, and from our group's website: http://www.biomedicalphysics.org.

Conclusion
In conclusion, we compared the performance of the available spectral/spatial analysis methods for sOCT and LCS.All methods suffer from a trade-off in spectral/spatial resolution.The effective spatial resolution of the STFT was found to be better than its theoretical spatial resolution, whereas the effective spatial resolution of the DW method was found to be worse than its theoretical spatial resolution.In general, for the application of localized measurements of hemoglobin concentrations and oxygen saturation, the STFT is the most optimal method.

Fig. 1 .
Fig. 1.Simulation of two reflectors (R1 and R2) for the DW and STFT method (Case 1): (a) recovered A-scans.(b) recovered spectra from R1, no filtering.(c) recovered spectra from R2, no filtering.(d) Fourier transform on the recovered spectra from R2, revealing an interference term at a distance of 25 µm (separation of R1 and R2) for the DW method.(d) recovered spectra from R1, low-pass filtered at 7 µm.(e) recovered spectra from R2, low-pass filtered at 7 µm.The actual, and STFT-recovered spectra overlap in figures b, c, e and f.

Figure 1
Figure1shows the results of the simulation with two reflectors R1 and R2, representing two vessel walls (25 µm separation) enclosing an absorbing blood volume ([tHb] = 150 g/L, SO 2 = 85%).The black line in Fig.1(a) corresponds to the actual A-scan envelope; the black lines in Figs.1(b), 1(c), 1(e) and1(f) to the actual reflection spectrum (i.e. the simulation inputs).The recovered µ a spectrum from the blood volume is shown in Fig.2for both methods, using the filtered spectra of R1 and R2.Compared to the STFT, the spatial recovery is much better for the DW method (Fig.1(a)).However, the recovered spectra by the DW method deviate stronger from the actual spectrum than the spectra recovered by the STFT (Figs.1(e), 1(f) and Fig.2), which is reflected in the values for χ 2 .Table2summarizes the performance of the STFT and DW method for the spectral recovery of R1 and R2 and the µ a in between.Lowpass filtering of the spectra enhances the χ 2 values for the DW method, but it remains higher with 1-2 orders of magnitude than the χ 2 values for the STFT method.Since the STFT method does not suffer from the influence of the interference between both reflectors for this window size and reflector separation, no difference is found between the χ 2 values for the unfiltered and filtered spectra of R1 and R2 within the presented significance (Table2).The slight increase in the χ 2 value of µ a after filtering is attributed to the minor loss of frequencies containing hemoglobin-specific absorption features above 7 µm (Fig.1(d)).

Fig. 3 .
Fig. 3.For both the DW and the STFT method: (a) fitted oxygen saturation SO 2 .(b) fitted total hemoglobin concentration [tHb].(c) χ 2 of µ a as a function of actual SO 2 for the case of two reflectors (Case 1).Error bars represent the 95% confidence intervals of the fit parameters.

Fig. 4 .
Fig. 4. Simulation of three reflectors (R1, R2 and R3) for the DW and STFT method (Case 2) with R2 and R3 separated at 20 µm: (a) recovered A-scans.(b) recovered absorption spectra µ a from the blood volume between R1 and R2.(c) recovered A-scans in the presence of noise.(d) recovered absorption spectra µ a from the blood volume between R1 and R2 in the presence of noise.Media 1 shows for panel a and b (without noise) how R3 approaches R2 and its effect on µ a for both methods.

Fig. 5 .
Fig. 5.For the case of three reflectors (Case 2): (a) χ 2 of µ a .(b) fitted total hemoglobin concentration [tHb].(c) fitted oxygen saturation SO 2 as a function of the separation between R2 and R3.Error bars represent the 95% confidence intervals of the fit parameters.

Figure 5 (
Figure 5(a) depicts the influence of R3 on the χ 2 of µ a , Fig. 5(b) on the fitted [tHb] and Fig. 5(c) on the fitted SO 2 .For all parameters, a deviation from the baseline value is observed at separations smaller than ~16 µm for the STFT, and smaller than ~30 µm for the DW method.Although the χ 2 of µ a for the DW method first decreases for separations <~30 µm (indicating better spectral recovery), the fitted values of the [tHb] and SO2 are also affected at these separations -indicating worse spectral recovery.This is caused by the introduction of frequency modulations on the recovered µ a (see Fig.4, Media 1) for separations <~30 µm, which 'deceptively' results in a better χ 2 of µ a .For this case, the effective spatial resolution (defined as a deviation in χ 2 of >15% from the χ 2 without the presence of R3) of both methods

Fig.
Fig.7.Effective spatial resolution of the DW and STFT method as a function of the power ratio between R3 and R2 (Case 3).The effective spatial resolution is defined as the separation between R2 and R3 (see Fig.4a), where the χ 2 of µ a between R1 and R2 deviates more than 15% from the χ 2 of µ a without the presence of R3.

Fig. 8 .
Fig. 8. (a) Effective spatial resolution of the DW and STFT method as a function of Δz L for the large spatial window w L (Δz L ,Δk L ) of the DW method (Case 4).(b) χ 2 of µ a for the STFT method and for the DW method, with a selection of sizes for Δz L .

Table 1 . Definition of simulated window sizes and theoretical resolutions
Δz: spatial window size (i.e.depth resolution), Δk: spectral window size, N w : number of data points in Δz or Δk, N total total number of data points in one spectrum or A-scan, Δλ: wavelength resolution.Note: these values do not account for the refractive index of the medium.If Δz is defined in tissue with group refractive index n, Δz will be a factor n larger in air, hence Δλ will be a factor n smaller than presented in this table.