Compressed sensing-enabled phase-sensitive swept-source optical coherence tomography

: 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 ﬁrst mathematically investigated the method, and veriﬁed it through computer simulations. We then conducted a vibrational frequency test and a ﬂow 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 ﬂow 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.


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][3][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][8][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][17][18][19][20][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

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 ∈ C M that has a sparse structure, from a linear observation y = Φx ∈ C P , where Φ ∈ C P×M is called the measurement operator, even when the number of observations is drastically lower than the size of the signal to recover (P M). The problem y = Φ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 Ψ ∈ C M×L , meaning that there exists a vector s ∈ C L such that x = Ψ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 estimatorx of the true signal x from the observation y by solving the following optimization problem: where Ψ * is the pseudo inverse of the operator Ψ. The quality of the estimationx of the ground truth 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]: where µ(Φ, Ψ), defined in [39], measures the incoherence between the matrices Φ and Ψ, and S is the sparsity level of x.

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 Φ ∈ {0, 1} P, M a random selection matrix, such that, for a calibration signal x, the resultant signal y = Φ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 Ψ = 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 F −1 then denotes the 2D-IDFT instead of the 1D-IDFT, and the selection matrix Φ belongs to {0, 1} N k P×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.

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 k is an integer between 2 and 8 with our set-up. Consequently, this method exploits the sparsity of the 2D-IDFT of the M × 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.

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: where 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 d is the optical path length difference between both mirrors , which is fixed for calibration channel. 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 µ(Φ, Ψ) 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: where S is an unknown constant, close to S.

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.

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). Depending on the size of the sub-sampling mask, the original B-scan will be segmented and randomly sub-sampled at the acquisition. Each segment is reconstructed by using NESTA algorithm and concatenated to form the reconstructed calibration B-scan.
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.

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

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,x between the original signal x and the reconstructed signalx, which is given by A r x,x close to 1 indicates a good reconstruction.

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.

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.

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. We then 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).

Ex vivo blood flow detection in swine heart
Epicardial ablation is an effective tool for the treatment of complex arrhythmias such as postmyocardial 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×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.

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 ofr x,x under different experimental conditions. Taking into account the randomness of the masks, five simulations were conducted for each sampling rate, and the averager x,x over these five simulations was computed.
The blue solid curve is the contour line that is given byr x,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,r x,x decreases exponentially with P/M. On the other hand, P/M could be reduced to ∼ 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 qualityr x,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 plottedr x,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 bringr x,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 log M/M, which fully agrees with the experimental results as shown in Fig. 3(d).

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.  The results of using pre-measured calibration (red dash lines) show excessive phase noise in both time domain and spectral domain, while the others present comparable performances to that of using the full calibration. It is worth noting that we could achieve a lower sub-sampling rate (P/M) by using a larger mask width (N k ) without compromising the performance. The measured noise floor for the proposed system is below 10 pm.

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

Doppler OCT
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. 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.

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.

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.

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

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. Fig. 8. The mis-synchronization due to the usage of two DAQ boards. (a) The raw phase angles obtained by CS recalibration over time and (b) their histogram. A fixed phase jump could be observed in (a), and the histogram in (b) shows two distinct peak with similar distribution. This phase jump is caused by the one pixel shift during the simultaneous triggering of the two DAQ baords. The numerically corrected phase angles (CS-based calibration) are plotted against that of (c) pre-measured calibration and (d) full calibration. The CS-based calibration shows a significantly superior phase stability to that of the pre-measured calibration and a comparable result to that of the full calibration.

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

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.

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.