## Abstract

Here we present a novel phase-sensitive swept-source optical coherence tomography (PhS-SS-OCT) system. The simultaneously recorded calibration signal, which is commonly used in SS-OCT to stabilize the phase, is randomly sub-sampled during the acquisition, and it is later reconstructed based on the Compressed Sensing (CS) theory. We first mathematically investigated the method, and verified it through computer simulations. We then conducted a vibrational frequency test and a flow velocity measurement in phantoms to demonstrate the system’s capability of handling phase-sensitive tasks. The proposed scheme shows excellent phase stability with greatly discounted data bandwidth compared with conventional procedures. We further showcased the usefulness of the system in biological samples by detecting the blood flow in *ex vivo* swine left marginal artery. The proposed system is compatible with most of the existing SS-OCT systems and could be a preferred solution for future high-speed phase-sensitive applications.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Phase-sensitive optical coherence tomography (PhS-OCT) is an emerging branch of optical coherence tomography (OCT), in which the depth-resolved phase angles are obtained in addition to the amplitudes from the reconstructed complex object functions. The retrieved phase angles could be used to sense picometer-level axial displacements via spectral domain phase microscopy (SDPM) [1], to provide better contrast in optical coherence tomography angiography (OCTA) [2–4], to detect the axial flow velocities in Doppler OCT (D-OCT) [5, 6], or to characterize the tissue stiffness through optical coherence elastography [7–9].

Towards the clinical translation of these functional OCT variants, the imaging speed of the system is highly valued in order to reject the motion artifacts [10, 11], facilitate lateral oversampling [12], or generate 4D time-lapse/time-resolved datasets [13, 14]. These make the swept-source (SS) implementation of OCT more favorable than its spectral domain (SD) counterpart due to the much higher A-line acquisition rate [15]. In fact, recent advancements in swept laser technologies have further widened the gap; the sweeping rate of the wavelength-swept laser source has been dramatically increased from 2 kHz to 28 MHz in the past decades [16–21].

However, the phase stability, a key performance specification, of SS-OCT is often lower than that of SD-OCT due to the following factors. Firstly, most wavelength-swept light sources used in SS-OCT utilize mechanical motions to rapidly tune the output wavelength, which leads to a poor repeatability of the wavenumber-versus-time scanning curve (*k*-*t* curve). Secondly, it is difficult to exactly synchronize the onset of the data acquisition with the beginning of the wavelength sweeping. The mis-synchronization will cause timing jitters. Both factors add multiplicative noises to the measured complex phase angle [22].

To address these issues, Yasuno *et al* first suggested that a calibration signal obtained in advance from a reference interferometer could be used to remap the OCT signal [23]. This method was later developed into a simultaneously recorded calibration channel by Gora *et al* [24], and further modified by Shangguan *et al* [25]. By using this simultaneously recorded calibration signal, both the scanning variability and the timing jitter could be corrected, and it is actively adopted by various research groups in recent years [14, 26–29].

Unfortunately, digitizing an extra channel will impose a hefty levy on the data bandwidth of the system. For example, a typical 200 kHz SS-OCT system could generate ∼800 MB of data every second, supposing that each A-line is digitized by 2,000 points at 16-bit precision. This figure will be quickly scaled up, if a faster laser is used (data rate larger than $>6$ GB/s needs to be handled for MHz light source [15]) or an oversampling scheme is adopted to improve the phase stability [22]. This inevitably forces SS-OCT designers to make a compromise among the imaging speed, data bandwidth, and the phase stability. In fact, the data transfer/record rate is now recognized as one of the major trade-offs in the SS-OCT system design [13].

To address this constraint, one straightforward approach is to compress the data channels utilizing Compressed Sensing (CS). The theory, which has been proposed in [30,31] and developed thereafter, states that a signal that has a sparse structure can be recovered with high precision from a vector of linear projections, even in the case where the size of the vector of projections is significantly smaller than the size of the original signal. A signal is said to have a sparse structure, if it can be written as the sum of very few vectors, in a given basis.

Previously, CS has been studied in the context of SD-OCT [32, 33], where the authors proved that it is possible to reconstruct the OCT images by using only a subset of the charge-coupled device (CCD) camera pixels. Their methods exploited the 1D-Fourier transform of the acquired OCT tissue, but did not take advantage of the 2D redundancy of the layers. Another intriguing approach has been proposed by Lebed *et al*, where OCT images are acquired by using an under-sampling grid [34]. However, their mathematical formulation is not related to CS, as requirements like matrix incoherence [31] are not met in their work. Other works have combined CS with SD-OCT, such as [35, 36], but their methods remain post-processing techniques instead of being real time. It is worth noting that all the previous CS researches are conducted on the basis of SD-OCT, which is not really useful in practice, since an addressable CCD camera is not commonly available and the data bandwidth of SD-OCT is not significant. The only exception is the work conducted by Mididoddi *et al*, in which they experimentally demonstrated a CS-based photonic time-stretch optical coherence tomography system [37]. However, their method is highly specialized and not compatible with most existing SS-OCT system.

In this manuscript, we proposed to apply CS on the simultaneously recorded calibration signal in PhS-SS-OCT: the calibration signal is randomly sub-sampled during the acquisition to reduce the data bandwidth, and later recovered based on the CS theory. The required sample numbers for exact recovery features a logarithmic growth when spectral oversampling is used. Moreover, the 2D redundancies of the calibration channel are also exploited to achieve the highest compression ratio. In this way, we could break the trade-off imposed by the digitization of an extra channel without losing the merits of possessing high phase stability. To our best knowledge, this is the first experimental implementation and demonstration of a CS-based acquisition in Fourier domain OCT systems. The proposed scheme is universal and could be applied to most existing SS-OCT system.

## 2. Compressed Sensing-enabled SS-OCT

#### 2.1. Compressed Sensing theory

The mathematical theory of CS has been introduced independently in 2006 by Candès, Romberg and Tao [30], and Donoho [31]. The goal of this theory is to provide tools and methods for the reconstruction of signals with sparse structures from a small collection of linear measurements. Mathematically, the problem of CS consists of recovering a signal $x\in {\u2102}^{M}$ that has a sparse structure, from a linear observation $y=\mathrm{\Phi}x\in {\u2102}^{P}$, where $\mathrm{\Phi}\in {\u2102}^{P\times M}$ is called the measurement operator, even when the number of observations is drastically lower than the size of the signal to recover ($P\ll M$). The problem $y=\mathrm{\Phi}x$ is strongly ill-posed in this context, while the prior information on the sparsity of the signal *x* is key. More generally, we consider that the signal *x* has a sparse representation in a known dictionary $\mathrm{\Psi}\in {\u2102}^{M\times L}$, meaning that there exists a vector $s\in {\u2102}^{L}$ such that $x=\mathrm{\Psi}s$ with very few non-zero coefficients. In this case, Ψ is called the sparsifying transform, and *s* is said *S*-sparse if only *S* of its coefficients are non-zero.

The CS theory shows that, under some constraints [30], it is possible to recover an estimator $\widehat{x}$ of the true signal *x* from the observation *y* by solving the following optimization problem:

*x*depends on several factors: both matrices Φ and Ψ, the level of sparsity of

*x*, and the intensity of noise in the observation.

In addition, note that, in the ideal case of noiseless acquisition, the minimum number *M* of measurements required for exact reconstruction is given by [38]:

*S*is the sparsity level of

*x*.

#### 2.2. Direct domain Compressed Sensing

The calibration signal in SS-OCT has two properties that motivated us to develop a CS-based technique, suited for faster SS-OCT acquisition:

First, the single calibration signal, corresponding to a single A-line, is a chirped sinusoidal function [24]. Such signal has a highly localized spectrum. In other words, these signals have sparse Fourier transform (see Fig. 1(b)), which is the best case scenario for direct domain CS acquisition, meaning that such signal can be reconstructed with a high precision from very few randomly acquired samples.

Second, any consecutive calibration signals (corresponding to consecutive A-lines, during the raster scanning of the whole sample) have a very similar appearance. Up to a few exceptions, where notable pixel shifts are observed, the pixel-wise difference of two consecutive calibration signals is equal to a constant. As a consequence, if we consider the whole stack of calibration signals as a calibration B-scan, the resulting image will present a strong redundancy along the B-scan direction (see Fig. 1(c)). This redundancy is exploited in this work, with CS reconstruction results improved when considering calibration signals as two dimensional images.

Most of state-of-the-art CS applications (such as MRI [40], Astrophysics [41], Radar [42], Holography [43]) exploit the natural sparsity of signals in the direct domain (Ψ being equal to the identity matrix, or to a wavelet transform) or of their gradient (Ψ being the total variation operator), while sampling the data in a very incoherent space, such as the Fourier or the Fresnel domain. Here, the direct domain refers to the space in which the data is acquired. We propose in this paper to design an alternative approach to CS, where the acquisition is done in the direct domain, and the sparsifying transform corresponds to the Fourier transform. Mathematically, we define $\mathrm{\Phi}\in {\{0,1\}}^{P,M}$ a random selection matrix, such that, for a calibration signal *x*, the resultant signal $y=\mathrm{\Phi}x$ is a random sub-sampling of the measures of *x*. In details, the matrix Φ contains at most 1 non-zero coefficient on each line, for a total amount of *P* non-zero coefficients, meaning that *y* is a selection of *P* samples among the *M* that define *x*. In addition, we define $\mathrm{\Psi}={\mathcal{F}}^{-1}$, which represents the inverse discrete Fourier transform (IDFT) operator. Then, the CS theory states that we can recover with high precision an estimator of *x* from the acquisition *y* by solving:

Note that, when *x* represents a group of *N _{k}* clock signals instead of only one clock signal, the approach is exactly the same. The operator ${\mathcal{F}}^{-1}$ then denotes the 2D-IDFT instead of the 1D-IDFT, and the selection matrix Φ belongs to ${\{0,1\}}^{{N}_{k}P\times {N}_{k}M}$.

Previous works in direct-domain CS include the spectral Compressive Sensing approach from [44] for general frequency-sparse signals, modeled as a combination of a small number of sinusoids. The spectral iterative hard thresholding (SIHT) algorithm using heuristic approximations was introduced in [44, 45]. However, such approach is not directly expandable to two dimensional objects, which is crucial in our approach, as we consider grouping consecutive calibration A-lines for improved results. For the sake of generality, we use the NESTA algorithm [46] to solve Eq. (3) in this work, as it can be easily adapted to both 1D and 2D cases.

#### 2.3. 1D vs 2D approach

Equation (3) gives an estimator of a calibration signal *x* from only a fraction of its samples, denoted as *y*. In this work, we propose two different approaches for the acquisition of *y*. While the mathematical formulation is the same, the algorithmic implementations are different.

The first option, which is standard in the domain, is to partially acquire one calibration A-line, and reconstruct an estimation of it through CS optimisation. This approach is one-dimensional, exploits the sparsity of the IDFT of each of the clock signals (see Fig. 1(b)), and can be solved using a multitude of state-of-the-art techniques, such as [44, 45] for instance. Then, the algorithm is run on the *N* calibration A-lines, for a complete reconstruction of the calibration B-scan.

The second option requires a technical improvement of the setup. As in this case, the acquired data *y* represents the sub-sampling of a group of *N _{k}* calibration A-lines, where

*N*is an integer between 2 and 8 with our set-up. Consequently, this method exploits the sparsity of the 2D-IDFT of the $M\times {N}_{k}$ portion of the calibration B-scan. Then, the algorithm is run on the $N/{N}_{k}$ groups of calibration A-lines, leading to the reconstruction of the complete calibration B-scan.

_{k}#### 2.4. Optimal sampling rate for CS-SS-OCT

Using the notations in the work of [22], we can model the clock signal as a chirped sine function that is sampled regularily over *M* data points:

*A*is the emission spectrum function, supposed to be constant,

*k*is the output wavenumber which is a function of time (sample number in discrete domain) , and

*z*is the optical path length difference between both mirrors , which is fixed for calibration channel.

_{d}Ideally, the *k*-*t* curve would be linear, leading to the calibration signal *I* being a perfect sinusoidal function. However, the actual function *k* of the laser is nonlinear (see Fig. 4(c)). The consequence of the nonlinear *k*-*t* curve is that the calibration signal consists in a sine wave with modulated amplitude, to which are added sine waves of lower frequencies. If we take the IDFT of the calibration signal, we observe two peaks corresponding to the highest frequency contained in the signal, of rather thin widths (less than 100 coefficients overall), as displayed on Fig. 1(b). Considering the other Fourier coefficients as noise terms, we deduce that the sparsity *S* of the clock signal is given by the number of coefficients contained at these two locations. As a consequence, since the mutual coherence $\mu \left(\mathrm{\Phi},\mathrm{\Psi}\right)$ is close to 1 in the case of direct-domain compressed sensing, the sub-sampling rate required for almost exact reconstruction of the clock signal is given by:

*S*

^{′}is an unknown constant, close to

*S*.

## 3. Methods

#### 3.1. OCT system

The SS-OCT system used in this study was developed on the basis of the one that was reported previously [47]. Two modifications are highlighted: (1) the original light source is replaced by a MEMS-based swept source (HSL-20-100-M-3-S, Santec, Japan), (2) the custom-made Michelson interferometer, used for calibration, is substituted by an integrated Mach Zehnder interferometer (INT-MZI-1300, Thorlabs, USA). The schematic of the system is shown in Fig. 2(a). The system have two output channels: (1) OCT channel and (2) calibration channel as illustrated in Fig. 2(b) and 2(c), respectively.

The center wavelength of the system is 1317.5 nm, the axial resolution is 16.8 *μ*m, and the lateral resolution is 8.77 *μ*m. The 6 dB coherence length is 6.5 mm, and the sensitivity of the system is 105.4 dB.

#### 3.2. Hardware implementation of sub-sampling

The schematic diagram of the proposed hardware-based signal sub-sampling is illustrated in Fig. 2(d). The OCT channel is digitized by one DAQ board (ATS 9373, AlazarTech, Canada) at full speed, while the simultaneously recorded calibration channel is digitized by another DAQ board of the same specifications but at a reduced rate via random sub-sampling. Specifically, a set of 65,536 registers, which is referred as “skip bitmap”, are used to specify which trigger event to be ignored during the acquisition by assigning binary values to the corresponding registers. Whenever a trigger event is skipped, no data will be generated. The skip bitmap could be made any arbitrary pattern and manipulated via the provided software API (ATS-SDK V7.1.4 (C++), AlazarTech, Canada).

For example, we can load a vector of 4,000 bits, which consists of alternating 1s and 0s, onto the skipping table. The DAQ board thus digitizes the signal when every other triggers are received during the first 4,000 clock periods. This process is repeated for the next 4,000 sampling points and so forth. We thus generate a sub-sampled calibration signal whose down-sampling factor is 2. The process is a segment-by-segment process as illustrated in Fig. 2(e).

The fully sampled OCT signal and the sub-sampled calibration signal are later transferred to the host computer via a Peripheral Component Interconnect Express (PCI-e) Gen. 3×8 bus.

It is worth noting that the two DAQ boards have to be carefully connected and synchronized to avoid extra jitter. Detailed discussions on this topic are provided later in Section 5.1.

#### 3.3. CS signal reconstruction

In order to reconstruct the calibration signal from the sub-sampled data, we solve the CS problem (Eq. (3)) using the NESTA algorithm [46], that has been shared by Becker *et al* on https://statweb.stanford.edu/ candes/nesta/. We implemented the algorithm using Matlab, on a PC workstation 2.93 GHz quad-core central processing unit (CPU), with 8 GB of RAM.

To be more precise, we used the ${\mathcal{l}}_{1}$ version of the proposed NESTA algorithm, adapted to the operators Φ and Ψ defined in Section 2.2. In addition, the regularization parameter *ε* has been set to 0 in our experiments, as we aim at exact reconstruction. Note that the denoising property of the regularization parameter is out of the scope of this manuscript.

Again, this reconstruction process is performed segment-by-segment as presented in Fig. 2(e). With these conditions, the reconstruction time of one calibration A-line is between 10 and 50 ms, depending on the sampling rate, and the sampling strategy.

#### 3.4. Experimental procedures

### 3.4.1. Reconstruction quality: numerical simulation

We first numerically studied the impact on reconstruction quality of different sub-sampling mask configurations including sub-sampling rate ($P/M$) and mask width (*N _{k}*).

Conventional SS-OCT calibration data were obtained at 200 Mega Samples/s (MS/s) and 1.6 GS/s for the numerical simulation. Random sub-sampling bitmaps were computer-generated and were used to digitally mask the calibration signal. The sub-sampled data was then reconstructed using the NESTA algorithm.

To evaluate the reconstruction quality, we use the correlation coefficient ${r}_{x,\widehat{x}}$ between the original signal *x* and the reconstructed signal $\widehat{x}$, which is given by

A ${r}_{x,\widehat{x}}$ close to 1 indicates a good reconstruction.

### 3.4.2. Phase stability

We then evaluated the scheme within its intended context by using the reconstructed calibration signal to remap and stabilize the corresponding OCT spectrum. A “standardized test” was conducted to measure the phase stability of the proposed system as described below.

We placed a microscope slide (1-mm thick, Microscope Slides, Fisherfinest, USA) under the sample arm, and blocked the reference arm. The interference pattern between the light reflected from the top surface and the bottom surface of the sample was obtained from the DAQ board $\#1$ at 800 MS/s. The sub-sampled calibration signal was recorded by DAQ board $\#2$ according to a predefined mask ($P/M=0.3$, *N _{k}* = 1), while the same calibration signal was fully digitized by DAQ board $\#1$ for comparison purposes. 1,000 M-scans were obtained at a fixed sample location.

We later reconstructed the calibration signal from its sub-sampled copy, and used this reconstructed signal to remap the OCT signal. The instantaneous phase angles over time at the peak location of the OCT A-lines were extracted and their standard deviation was calculated.

### 3.4.3. Proof-of-concept: vibrational frequency test

The system configuration is identical to that reported in Section 3.4.2. We used SDPM to measure the vibrational frequency of the sub-sampling-interval displacements of a sample to verify the system’s capability of conducting sensitive phase measurements. The test arrangement included a piezoelectric actuator (PZS001, Thorlabs, USA), which was driven by a a sinusoidal voltage signal from a function generator (AFG3022C, Tektronix, USA). The frequency of the driving sinusoidal is 10 kHz and the peak-to-peak voltage is 1V.

To showcase the performance of our proposed system, we used two other remapping schemes (Yasuno *et al* [23] and Shangguan *et al* [25]) on the same OCT dataset for comparison purpose. In the following sections, the terms pre-measured calibration and full calibration to will be used refer these two schemes, respectively.

### 3.4.4. Proof-of-concept: Doppler OCT

To further validate the proposed method, an experimental phantom was constructed to mimic blood flow. Intralipid emulsion (Sigma Aldrich, USA) was diluted to a concentration of $0.25\%$ in de-ionized water and stored within an intravenous fluid bag and tubing setup. The tubing was fitted into an irrigation pump (CoolFlow, Biosense Webster, USA) which allowed for precise manipulation of laminar flow rates. Imaging was performed over the short axis of the tubing for flow rates ranging between 1-3 mL/min. Wethen obtained Doppler OCT images in addition to optical micro-angiography (OMAG) [48]. Moreover, we compared resultant OMAG images by using two other remapping schemes (see Section 3.4.3).

### 3.4.5. *Ex vivo* blood flow detection in swine heart

Epicardial ablation is an effective tool for the treatment of complex arrhythmias such as post-myocardial infarction ventricular tachycardia (VT). In such procedures, care must be taken to avoid energy delivery proximal to coronary vasculature which can lead to vessel occlusion and thrombus formation [49, 50]. In this *ex vivo* experiment, we showcase the capability of the proposed system to detect the blood blow, which could potentially be used for real-time monitoring for the ablation.

Fresh swine hearts were obtained from Green Village Packing (Green Village, New Jersey, USA). Ventricular slabs containing epicardial vasculature were resected and cannulated using 900 *μ*m diameter flexible tubing coupled to an irrigation pump (CoolFlow, Biosense Webster, USA). Whole swine blood was then diluted at 1:10 ratio in phosphate buffered saline and perfused through the vessel at 5 mL/min. OCT volumes (2048×256×1024 pixels, 1.2×4.8×2.5 mm^{3}, long×width*es*depth) were acquired. The sub-sampling rate of the calibration channel was set to 30% and the mask width was 1. We then used the OMAG algorithm to visualize the blood flow.

## 4. Results

#### 4.1. Reconstruction quality evaluation

The sampling rate of the DAQ board was first set as 200 MS/s so that each A-line contained 1,000 sampling points ($M=1,000$). Result is in Fig. 3(a) show pseudo-color maps of ${\overline{r}}_{x,\widehat{x}}$ under different experimental conditions. Taking into account the randomness of the masks, five simulations were conducted for each sampling rate, and the average ${\overline{r}}_{x,\widehat{x}}$ over these five simulations was computed.

The blue solid curve is the contour line that is given by ${\overline{r}}_{x,\widehat{x}}=0.9$, which we take as a baseline in this study. It shows that for a mask that consists of only one A-line, a sub-sampling rate ($P/M$) $>40\%$ is necessary to achieve a result better than the baseline. Generally, ${\overline{r}}_{x,\widehat{x}}$ decreases exponentially with $P/M$. On the other hand, $P/M$ could be reduced to $\sim 15\%$ if mask width *N _{k}* is set to 8, thanks to the redundancy within B-scans.

We repeated the experiment for a different sampling rate (1.6 GS/s), and the result is plotted in Fig. 3(b). The same reconstruction quality ${\overline{r}}_{x,\widehat{x}}$ could be achieved at a lower sub-sampling rate/larger mask width: For example, by sub-sampling only 5% of the data (*N _{k}* = 8), the averaged correlation coefficient between the original calibration signal and the reconstructed one could be as high as $0.98$. To better convey the idea, we plotted ${\overline{r}}_{x,\widehat{x}}$ against the sub-sampling percentage for four combinations of sampling rates and mask widths in Fig. 3(c).

We further conducted experiments for three other sampling rates: 600 MS/s, 800 MS/s and 1 GS/s. We recorded the $P/M$ (*N _{k}* = 1) that is needed to bring ${\overline{r}}_{x,\widehat{x}}$ to the baseline (0.9) value at these sampling rates. Based on the prediction made in Eq. (5), sub-sampling rates should be linear with $\mathrm{log}\phantom{\rule{0.2em}{0ex}}M/M$, which fully agrees with the experimental results as shown in Fig. 3(d).

#### 4.2. Phase stability

An exemplary full calibration A-line is plotted in Fig. 4(a), while the sub-sampled acquisition and reconstruction are given in Fig. 4(b). The reconstructed calibration signal is almost identical to the fully sampled one in spite of a scaling factor, which is caused by the difference between the two DAQ boards.

The extracted *k*-*t* curve from the full and the reconstructed signals are illustrated side-by-side in Fig. 4(c). The error, which is four orders of magnitude smaller than the original signal, is also plotted in the same panel in red solid line.

In the last step, we computed the histogram of the phase angle distribution at the peak location. Results, using the full and the reconstructed calibration signal, are shown in Fig. 4(d) and 4(e), respectively. The standard deviation of the measured phase angles is equal to $4.53$ mrad (signal-to-noise ratio (SNR) = 52.65 dB) for full calibration and $4.49$ mrad (SNR = 50.86 dB) for reconstructed calibration. The phase stability of the proposed system is not compromised, while only 30% of the bandwidth for the extra calibration channel is used in this example.

It is worth noting that reconstructed results were numerically corrected before presentation. Details can be found in Section 5.3.

#### 4.3. Vibrational frequency test

Instantaneous phase angles extracted from the surface of the piezoelectric actuator over a small time frame are plotted in Fig. 5(a). Four different remapping strategies are used: pre-measured calibration, full calibration, CS-1 ($P/M=0.3,{N}_{k}=1$) and CS-2 ($P/M=0.15,{N}_{k}=8$). Both CS results are very similar to that of the full calibration case, while the pre-calibration result manifests apparent noises. We then conducted a discrete Fourier transform (DFT) on the acquired phase angles, and obtained the frequency response of the actuator which is illustrated in Fig. 5(b). The dominant frequency peak at 10 kHz matches that of the stimulation exactly, and is well preserved across all methods. However, if we take a closer look at the low frequency region of the spectra as shown in the inset, we observe the elevated background noise in the case where no simultaneous calibration is used (red dashed line). On the other hand, both reduced-data-rate schemes perform very well and are comparable to that of the conventional dual-channel scheme. The noise floor of the proposed scheme is below 10 pm, if we exclude the discrete spikes that are caused by the mechanical motion of the imaging platform.

#### 4.4. Doppler OCT

We set the flow velocity of the irrigation pump to be 2 mL/min. The Tygon tubing (Saint-Gobain, France) we used in the experiment has an inner diameter about 0.89 mm and an outer diameter of 1.56 mm. It should be noted that the tube was positioned at an angle of 86.05° to the vertical and the cross-sectional area of the tube is 87.11% smaller than that of the pump. The original OCT image of the phantom is shown in Fig. 6(a), and its DFT spectrum against the fast axis is shown in Fig. 6(e). The DFT spectrum shows a strong Doppler frequency shift in the inside region of the tube. After high-pass filtering the image, we obtain the OMAG image by following the steps described in Wang *et al* [48].

Similar to what we have done in previous subsection, three different remapping schemes (full calibration, pre-measured calibration, and proposed CS-based calibration) were used and the results are given in Fig. 6(b), 6(c), and 6(d), respectively. The CS-reconstructed results ($P/M=0.3$, *N _{k}* = 1) shows good image quality comparable to using the full calibration. On the contrary, pre-measured calibration presents multiple artifacts due to the phase noise.

We lastly computed the Doppler image by using the CS-reconstructed calibration, which is displayed in Fig. 6(f). To further showcase the velocity profile, the average of 4 adjacent Doppler A-lines as indicated by the bold red line in Fig. 6(f) was calculated and plotted in Fig. 6(g). The average measured flow speed (half of the maximum velocity) is 0.81 mm/s, which agrees with the prediction.

#### 4.5. Ex vivo blow flow detection in swine samples

The obtained OCT image and the calculated OMAG mapping are overlaid and presented in Fig. 7(a). Representative sectional views in both *y* (slow scanning axis)$-z$ and *x* (fast scanning axis)$-z$ planes are given in Fig. 7(b) and 7(c), respectively. The left marginal artery is segmented in OMAG image with high contrast. It should be noted that the artifacts presented in the lower boundary of the vessel is mainly due to the reduced SNR in those regions.

## 5. Discussion

Both simulations and experiments conclude that the proposed CS-based system is capable of conducting phase-sensitive SS-OCT measurements with reduced acquisition data rate. Hence, phase stability of the system is comparable to that of using full calibration. Here in this section, we provide some side notes to further help the readers to understand the results and their significances.

#### 5.1. Hardware-based framework suitable for future compression on OCT channel

The proposed hardware configuration could be extended to the OCT channel to randomly decimate the signal. A CS-based algorithm, which utilizes 1D sparsity of the A-lines, could be implemented either in real-time [51] or post-processing [33] to reconstruct the OCT images at a reduced data rate. Time-elapse imaging and Doppler OCT could be good candidate applications: considering that the scenes update slowly in either fast or slow axis, we could take advantage of the redundancy in those directions by combining scenes together and conducting a 2D random sub-sampling.

#### 5.2. Nyquist sampling versus Compressed Sensing

The calibration signal shown in Fig. 1 appears to be band-limited, and we theoretically should be able to reconstruct the original signal by just sampling at the theoretical Nyquist frequency. However, the actual calibration signal acquired at the Nyquist frequency suffers from aliasing caused by spectral leakage toward the high frequencies, which will eventually lead to compromised phase calibration of the system [22]. An over-sampling scheme is a plausible approach to alleviate this effect [22, 24, 29, 52]. Intuitively, a higher sampling rate could better localize the true starting time of the wavelength sweeping.

Furthermore, a much lower data rate than the Nyquist rate is reached with the proposed method if the 2D redundancy of the calibration signal is exploited, as presented in Fig. 3. The extra digitization channel implemented within this study was to experimentally show the equivalence of the compressed clock to the fully digitized clock. As shown within our study, by using a 1D mask only 30% of the clock signal will need to be recorded. For a 2D mask (8 A-lines), this is further reduced to 15%.

#### 5.3. Synchronization between the DAQ boards

One of the biggest challenges in this study was to synchronize the two DAQ boards in the experiment, so that both the sub-sampled signal and the fully digitized OCT signal are perfectly aligned. Misalignment between the two channels will cause an extra jitter that deteriorate the phase stability of the system. The reason why we are not using a single DAQ board to simultaneously perform the full digitization of OCT channel and the sub-sampling of the calibration channel is purely technical: the DAQ board we are using currently only supports one mode at a time.

In our setup, we have to first ensure that the two boards are of the exact same specifications and are connected to the same trigger signal through same length of electrical connections. But there is still a random jitter that is smaller than 1 clock period between the two boards due to the random starting phase of their crystal oscillators. In fact, if we remap the OCT signals by using reconstructed calibration signals as described in Section 4.2, we observe a fixed amount of jump in the phase as shown in in Fig. 8(a). The histogram of the uncorrected phase angles, which is given in Fig. 8(b), shows two peaks and two groups of distributions. One group corresponds to the reconstructed calibration that has no jitter, and one corresponds to the missed one clock period.

This jitter would not exist if we have a device that could support both modes at the same time. For example, a dual channel DAQ board with a customized field-programmable gate array (FPGA). Therefore, in this study, we performed the numerical correction by measuring the value between the two peaks and subtract it from the group that has the high value. The corrected results are shown in Fig. 8(c) along with the results by using pre-measured calibration. A side-by-side comparison with the full calibration is given in Fig. 8(d) for a short excerpt. To make the two comparable, we subtracted the random starting phase from the corrected phase.

#### 5.4. Optimal sampling rate

The discussion in Section 2.4 shows that the minimum sampling rate required for exact reconstruction of a clock signal by using direct domain CS is linear with $\mathrm{log}\phantom{\rule{0.2em}{0ex}}M/M$, with a constant that is close to the sparsity level of the signal according to Eq. (5). Prediction of the value of this constant *S*^{′} remains an open question, which will be the focus of our next work. Moreover, techniques other than direct domain CS such as dynamic compressed sensing algorithms [53] will be considered for comparison purpose as well.

#### 5.5. Towards real-time processing

As mentioned before in the Methods section, the CS reconstruction of the calibration signal is done in post-processing for the time being. Currently, the CPU-based processing requires tens of milliseconds to perform one A-line, which is orders-of-magnitude slower than the required processing rate for streaming. But thanks to the independent nature of the OCT dataset, it is possible to incorporate parallel computing to accelerate the processing [51]. Moreover, using a user-programmable Field-programmable gate array (FPGA) might further speed up the CS reconstruction [54]. For the next step, we will explore both pathways and strive to implement real-time OCT image reconstruction and visualization.

## 6. Conclusion

To summarize, we have proposed a hardware-based compressed sensing SS-OCT system, and we demonstrated within a biological sample *ex vivo* localization of flow within a coronary artery using phase based analysis. Three sets of proof-of-concept experiments were conducted and the proposed system is shown to be capable of conducting highly phase-sensitive tasks at a reduced data sampling rate, which is especially useful for ultrahigh-speed applications.

We have also theoretically studied the SS-OCT calibration signal under the framework of compressed sensing theory. The results could be used as a guideline for future system design and development.

## Funding

National Institutes of Health (NIH) (1DP2HL127776-01) (CPH), National Science Foundation (NSF) (CAREER Award 1454365) (CPH), Agence Nationale de la Recherche (ANR-10-INSB-04-06 France Bioimaging, ANR-10-LABEX-62-IBEID).

## Acknowledgments

The authors would like to thank Dr. Y. Gan and Dr. P. Spinicelli for very useful discussions. We acknowledge computing resources from Columbia University’s Shared Research Computing Facility project, which is supported by NIH Research Facility Improvement Grant 1G20RR030893-01, and associated funds from the New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR) Contract C090171, both awarded April 15, 2010.

## References

**1. **H. Y. Lee, P. D. Raphael, J. Park, A. K. Ellerbee, B. E. Applegate, and J. S. Oghalai, “Noninvasive in vivo imaging reveals differences between tectorial membrane and basilar membrane traveling waves in the mouse cochlea,” Proceedings of the National Academy of Sciences **112**, 3128–3133 (2015). [CrossRef]

**2. **R. K. Wang, S. L. Jacques, Z. Ma, S. Hurst, S. R. Hanson, and A. Gruber, “Three dimensional optical angiography,” Optics Express **15**, 4083–4097 (2007). [CrossRef] [PubMed]

**3. **D. M. Schwartz, J. Fingler, D. Y. Kim, R. J. Zawadzki, L. S. Morse, S. S. Park, S. E. Fraser, and J. S. Werner, “Phase-variance optical coherence tomography: A technique for noninvasive angiography,” Ophthalmology **121**, 180–187 (2014). [CrossRef]

**4. **A. Zhang, Q. Zhang, C.-L. Chen, and R. K. Wang, “Methods and algorithms for optical coherence tomography-based angiography: a review and comparison,” J. of Biomedical Optics **20**, 100901 (2015). [CrossRef]

**5. **Y. Zhao, Z. Chen, C. Saxer, S. Xiang, J. F. de Boer, and J. S. Nelson, “Phase-resolved optical coherence tomography and optical Doppler tomography for imaging blood flow in human skin with fast scanning speed and high velocity sensitivity,” Optics Letters **25**, 114–116 (2000). [CrossRef]

**6. **R. A. Leitgeb, R. M. Werkmeister, C. Blatter, and L. Schmetterer, “Doppler optical coherence tomography,” Progress in Retinal and Eye Research **41**, 26–43 (2014). [CrossRef] [PubMed]

**7. **S. Wang and K. V. Larin, “Shear wave imaging optical coherence tomography (SWI-OCT) for ocular tissue biomechanics,” Optics Letters **39**, 41–44 (2014). [CrossRef]

**8. **S. Wang and K. V. Larin, “Optical coherence elastography for tissue characterization: A review,” Journal of Biophotonics **8**, 279–302 (2015). [CrossRef]

**9. **X. Liang, A. L. Oldenburg, V. Crecea, E. J. Chaney, and S. A. Boppart, “Optical micro-scale mapping of dynamic biomechanical tissue properties, Optics Express **16**, 11052–11065 (2008). [CrossRef] [PubMed]

**10. **R. F. Spaide, J. G. Fujimoto, and N. K. Waheed, “Image artifacts in optical coherence tomography angiography,” Retina **35**, 2163–2180 (2015). [CrossRef] [PubMed]

**11. **L. An and R. K. Wang, “In vivo volumetric imaging of vascular perfusion within human retina and choroids with optical micro-angiography,” Optics Express **16**, 11438–11452 (2008). [CrossRef] [PubMed]

**12. **Y. K. Tao, A. M. Davis, and J. A. Izatt, “Single-pass volumetric bidirectional blood flow imaging spectral domain optical coherence tomography using a modified Hilbert transform,” Optics Express **16**, 12350–12361 (2008). [CrossRef] [PubMed]

**13. **W. Wieser, W. Draxinger, T. Klein, S. Karpf, T. Pfeiffer, and R. Huber, “High definition live 3D-OCT in vivo: design and evaluation of a 4D OCT engine with 1 GVoxel/s,” Biomedical Optics Express **5**, 2963–2977 (2014). [CrossRef] [PubMed]

**14. **S. Song, W. Wei, B.-Y. Hsieh, I. Pelivanov, T. T. Shen, M. O’Donnell, and R. K. Wang, “Strategies to improve phase-stability of ultrafast swept source optical coherence tomography for single shot imaging of transient mechanical waves at 16 kHz frame rate,” Applied Physics Letters **108**, 191104 (2016). [CrossRef]

**15. **T. Klein and R. Huber, “High-speed OCT light sources and systems [Invited],” Biomedical Optics Express **8**, 828–859 (2017). [CrossRef]

**16. **B. Golubovic, B. E. Bouma, G. J. Tearney, and J. G. Fujimoto, “Optical frequency-domain reflectometry using rapid wavelength tuning of a Cr${}^{4+}$:forsterite laser,” Optics Letters **22**, 1704–1706 (1997). [CrossRef]

**17. **R. A. Huber, M. Wojtkowski, and J. G. Fujimoto, “Fourier Domain Mode Locking (FDML): A new laser operating regime and applications for optical coherence tomography,” Optics Express **14**, 3225–3237 (2006). [CrossRef] [PubMed]

**18. **S. H. Yun, G. J. Tearney, J. F. de Boer, N. Iftimia, and B. E. Bouma, “High-speed optical frequency-domain imaging,” Optics Express **11**, 2953–2963 (2003). [CrossRef] [PubMed]

**19. **W. Wieser, B. R. Biedermann, T. Klein, C. M. Eigenwillig, and R. Huber, “OCT Multi-Megahertz: High quality 3D imaging at 20 million A-scans and 4.5 GVoxels per second,” Optics Express **18**, 14685–14704 (2010). [CrossRef]

**20. **I. Grulkowski, J. J. Liu, B. Potsaid, V. Jayaraman, C. D. Lu, J. Jiang, A. E. Cable, J. S. Duker, and J. G. Fujimoto, “Retinal, anterior segment and full eye imaging using ultrahigh speed swept source OCT with vertical-cavity surface emitting lasers,” Biomedical Optics Express **3**, 2733–2751 (2012). [CrossRef] [PubMed]

**21. **X. Wei, A. K. S. Lau, Y. Xu, K. K. Tsia, and K. K. Y. Wong, “28 MHz swept source at 1.0 μm for ultrafast quantitative phase imaging,” Biomedical Optics Express **6**, 3855–3864 (2015). [CrossRef]

**22. **Y. Ling, Y. Gan, X. Yao, and C. P. Hendon, “Phase-noise analysis of swept-source optical coherence tomography systems,” Optics Letters **42**, 1333–1336 (2017). [CrossRef] [PubMed]

**23. **Y. Yasuno, V. D. Madjarova, S. Makita, M. Akiba, A. Morosawa, C. Chong, T. Sakai, K.-P. Chan, M. Itoh, and T. Yatagai, “Three-dimensional and high-speed swept-source optical coherence tomography for in vivo investigation of human anterior eye segments,” Optics Express **13**, 10652–10664 (2005). [CrossRef] [PubMed]

**24. **M. Gora, K. Karnowski, M. Szkulmowski, B. J. Kaluzny, R. Huber, A. Kowalczyk, and M. Wojtkowski, “Ultra high-speed swept source OCT imaging of the anterior segment of human eye at 200 kHz with adjustable imaging range,” Optics Express **17**, 14880–14894 (2009). [CrossRef] [PubMed]

**25. **Z. Shangguan, Y. Shen, P. Li, and Z. Ding, “Wavenumber calibration and phase measurement in swept source optical coherence tomography (in Chinese),” Acta Physica Sinica **65**, 034201 (2016).

**26. **M. Singh, C. Wu, C.-H. Liu, J. Li, A. Schill, A. Nair, and K. V. Larin, “Phase-sensitive optical coherence elastography at 1.5 million A-Lines per second,” Optics Letters **40**, 2588–2591 (2015). [CrossRef] [PubMed]

**27. **S. Amarakoon, J. H. De Jong, B. Braaf, S. Yzer, T. Missotten, M. E. van Velthoven, and J. F. de Boer, “Phase-resolved doppler optical coherence tomographic features in retinal angiomatous proliferation,” American Journal of Ophthalmology **160**, 1044–1054 (2015). [CrossRef] [PubMed]

**28. **S. Kim, P. D. Raphael, J. S. Oghalai, and B. E. Applegate, “High-speed spectral calibration by complex FIR filter in phase-sensitive optical coherence tomography,” Biomedical Optics Express **7**, 1430–1444 (2016). [CrossRef] [PubMed]

**29. **W. Chen, C. Du, and Y. Pan, “Cerebral capillary flow imaging by wavelength-division-multiplexing swept-source optical Doppler tomography,” Journal of Biophotonics **11**, e201800004(2018). [CrossRef] [PubMed]

**30. **E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Information Theory **52**, 489–509 (2006). [CrossRef]

**31. **D. Donoho, “Compressed sensing,” IEEE Trans. Information Theory **52**, 1289–1306 (2006). [CrossRef]

**32. **N. Mohan, I. Stojanovic, W. C. Karl, B. E. A. Saleh, and M. C. Teich, “Compressed sensing in optical coherence tomography,” Proc. SPIE **7570**, 75700L (2010). [CrossRef]

**33. **X. Liu and J. U. Kang, “SD-OCT Compressive: the application of compressed sensing in spectral domain optical coherence tomography,” Optics Express **18**, 22010–22019 (2010). [CrossRef]

**34. **E. Lebed, P. J. Mackenzie, M. V. Sarunic, and M. F. Beg, “Rapid volumetric OCT image acquisition using Compressive Sampling,” Optics Express **18**, 21003–21012 (2010). [CrossRef] [PubMed]

**35. **W. Meiniel, Y. Gan, C. P. Hendon, J.-C. Olivo-Marin, A. Laine, and E. D. Angelini, “Sparsity-based simplification of spectral-domain optical coherence tomography images of cardiac samples,” in *Proceedings of 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI)*, (IEEE, 2016), pp. 373–376. [CrossRef]

**36. **W. Meiniel, Y. Gan, J.-C. Olivo-Marin, and E. Angelini, “A sparsity-based simplification method for segmentation of spectral-domain optical coherence tomography images,” Proc. SPIE **10394**, 1039406 (2017).

**37. **C. K. Mididoddi, F. Bai, G. Wang, J. Liu, S. Gibson, and C. Wang, “High-throughput photonic time-stretch optical coherence tomography with data compression,” IEEE Photonics Journal **9**, 1–15 (2017). [CrossRef]

**38. **E. Candès and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse Probl. **23**, 969 (2007). [CrossRef]

**39. **D. L. Donoho and X. Huo, “Uncertainty principles and ideal atomic decomposition,” IEEE Trans. Inf. Theory **47**, 2845–2862 (2001). [CrossRef]

**40. **M. Lustig, D. Donoho, and J. M. Pauly, “MRI Sparse: The application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine **58**, 1182–1195 (2007). [CrossRef]

**41. **J. Bobin, J.-L. Starck, and R. Ottensamer, “Compressed sensing in astronomy,” IEEE Journal of Selected Topics in Signal Processing **2**, 718–726 (2008). [CrossRef]

**42. **M. A. Herman and T. Strohmer, “High-resolution radar via compressed sensing,” IEEE Trans. Signal Process. **57**, 2275–2284 (2009). [CrossRef]

**43. **M. M. Marim, M. Atlan, E. Angelini, and J.-C. Olivo-Marin, “Compressed sensing with off-axis frequency-shifting holography,” Optics Letters **35**, 871–873 (2010). [CrossRef] [PubMed]

**44. **M. F. Duarte and R. G. Baraniuk, “Spectral compressive sensing,” Applied and Computational Harmonic Analysis **35**, 111–129 (2013). [CrossRef]

**45. **R. Baraniuk, V. Cevher, M. Duarte, and C. Hegde, “Model-based compressive sensing,” IEEE Trans. Inf. Theory **56**, 1982–2001 (2010). [CrossRef]

**46. **S. Becker, J. Bobin, and E. Candès, “NESTA: A Fast and Accurate First-Order Method for Sparse Recovery,” Journal on Imaging Sciences **4**, 1–39 (2011). [CrossRef]

**47. **Y. Ling, X. Yao, and C. P. Hendon, “Highly phase-stable 200 kHz swept-source optical coherence tomography based on KTN electro-optic deflector,” Biomedical Optics Express **8**, 3687–3699 (2017). [CrossRef] [PubMed]

**48. **R. K. Wang and L. An, “Doppler optical micro-angiography for volumetric imaging of vascular perfusion in vivo,” Opt. Express **17**, 8926–8940 (2009). [CrossRef] [PubMed]

**49. **J. Whitaker, H. Raju, C. Taylor, and C. A. Rinaldi, “Accelerated idioventricular rhythm after left atrial tachycardia ablation as a marker of acute coronary ischemia,” HeartRhythm Case Reports **1**, 99–102 (2015). [CrossRef] [PubMed]

**50. **E. M. Aliot, W. G. Stevenson, J. M. Almendral-Garrote, F. Bogun, C. H. Calkins, E. Delacretaz, P. D. Bella, G. Hindricks, P. Jaïs, M. E. Josephson, J. Kautzner, G. N. Kay, K.-H. Kuck, B. B. Lerman, F. Marchlinski, V. Reddy, M.-J. Schalij, R. Schilling, K. Soejima, and D. Wilber, “EHRA/HRS Expert Consensus on Catheter Ablation of Ventricular Arrhythmias: Developed in a partnership with the European Heart Rhythm Association (EHRA), a Registered Branch of the European Society of Cardiology (ESC), and the Heart Rhythm Society (HRS); in collaboration with the American College of Cardiology (ACC) and the American Heart Association (AHA),” EP Europace **11**, 771–817 (2009). [CrossRef]

**51. **D. Xu, Y. Huang, and J. U. Kang, “Real-time compressive sensing spectral domain optical coherence tomography,” Opt. Lett. **39**, 76–79 (2014). [CrossRef]

**52. **S. Moon and Z. Chen, “Phase-stability optimization of swept-source optical coherence tomography,” Biomed. Opt. Express **9**, 5280–5295 (2018). [CrossRef] [PubMed]

**53. **M. S. Asif and J. Romberg, “Sparse recovery of streaming signals using${\mathcal{l}}_{1}$-homotopy,” IEEE Trans. Signal Processing **62**, 4209–4223 (2014). [CrossRef]

**54. **A. Kulkarni and T. Mohsenin, “Accelerating compressive sensing reconstruction OMP algorithm with CPU, GPU, FPGA and domain specific many-core,” in Proceedings of 2015 IEEE International Symposium on Circuits and Systems (ISCAS), pp. 970–973.