1. Introduction

There is growing interest in using optical coherence tomography (OCT) for highly sensitive vibrational measurements in the auditory system. This is motivated by the ability to spatially resolve the vibratory response with subnanometer sensitivity. The detailed functional images generated have been used as a tool to probe fundamental auditory biomechanics and investigated as potential clinical diagnostics. Most commonly, a pure-tone sine wave stimulus is applied and the vibration of the tissue of interest at the stimulus frequency is measured. However, multi-tones and/or broad-band stimuli can also be used to probe non-linear biological phenomenon. These applications and more promise to drive the continued development of this nascent field.

Similar approaches have been published under a variety of names such as phase-sensitive optical coherence tomography for vibrometry [1], optical coherence tomographic vibrography [2], and volumetric optical coherence tomography vibrometry (VOCTV) [3]. They all utilize the interferometric phase of a low-coherence interferometer to extract periodic small-scale motion in the sample of interest. It is the low-coherence counterpart to the well-known technique, Laser Doppler Vibrometry (LDV). LDV has been used for many years in a variety of applications including materials testing [4,5] and biomechanical measurements [6–8]. The primary advantage of OCT based vibrometry is the spatial resolution. An LDV only reports the vibratory response of the dominant reflector; as a consequence in many applications a strongly-reflective object is placed on the object of interest. OCT based vibrometry can spatially resolve the position of the vibrating structure deep within the sample with a resolution of 10s of microns and sensitivity to vibrations with amplitudes under 10 pm [9].

A clear understanding of the fundamental noise statistics is key to interpreting results and pushing the sensitivity limit of this technology for scientific and clinical applications. While various aspects have been investigated and reported in the literature there has not been a detailed and complete derivation with consistent notation and assumptions. The sensitivity for vibrometry is fundamentally linked to the noise statistics of the OCT signal, hence we start from the spectral interferogram and derive the noise statistics through the entire processing chain.

Before beginning we offer a brief review of the literatures surrounding this topic. Most of the work is not directly linked to OCT based vibrometry, but rather other related techniques that make use of the interferometric phase in OCT. The noise statistics for optical coherence tomography have been derived previously [10–12] to demonstrate the advantages of spectral interferometry over time-domain interferometry for OCT. These were derived for shot-noise, defining the fundamental limit to OCT sensitivity. Shortly thereafter, researchers began to develop techniques which exploited the phase-sensitivity inherent to spectral interferometry. In the context of spectral domain phase microscopy, Choma et al. [13] noted that the displacement sensitivity is proportional to (2SNR)^−1/2, where SNR is the signal to noise ratio and we have corrected the factor of 2 multiplier. Joo et al. [14] showed a similar result along with experimental verification of the distribution. Both results are for the time-domain phase. The noise statistics in the frequency domain were in part derived by Szkulmowski et al. [15] who noted the Rician distribution that arises for the magnitude of the phase for flow velocity estimation. Many of the relevant statistical arguments we will make have also been addressed in the context of Magnetic Resonance Imaging [16].

In this paper the effect of additive noise on vibrometry was analyzed by considering the noise statistics through the entire processing chain with results for both the time-domain and frequency-domain phase. These are shown to be consistent with both numerical simulation and experiments with a swept-source OCT system.

2. Derivation of noise statistics

We will start by assuming all of the noise is additive since the main noise sources such as shot-noise, thermal noise, and relative intensity noise (RIN) are additive [17–19]. This is slightly more general, but consistent with prior derivations which assumed additive shot-noise. We further assume that the additive noise term, n_k(k, t), is independent, identically distributed (i.i.d) white Gaussian noise with zero mean and a standard deviation of σ_k. In order to determine the effect of additive noise on the extraction of the vibratory response, it is necessary to understand how the additive noise in the k-domain propagates through each mathematical operation, starting with the spectral interferogram H(k, t) and ending with the magnitude, |ϕ(z, f)|, and phase, ϕ(z, f), of the vibratory response. A flow chart showing the processing steps along with the statistical results, e.g. the mean and standard deviation, is shown in Fig. 1. We define signal to noise ratio (SNR) for the OCT system in the conventional way, signal power divided by the variance of the noise, i.e. SNR = (s/n)² or on a dB scale SNR = 20log(s/n), where s is the signal amplitude and n is the standard deviation of the noise.

 

Fig. 1. A flow chart showing the processing steps along with the statistical results such as distribution, mean (µ) and standard deviation (σ). Symbols: k, wavenumber; t, time; z, optical pathlength difference, related to the tissue depth (Δz) by z = 2nΔz, where n is refractive index; f, frequency; H(k, t), a time series of spectral interferogram; i${\mathfrak{F}}$_k(·), inverse Fourier transform along the k dimension; ${\mathfrak{F}}$_t(·), Fourier transform along the time dimension; |·|, magnitude; ∠, angle;. i(z, t), h(z, t) without noise; N, the number of spectral channels; M is the number of time samples; f_N, subscript indicating frequency domain with normal distribution; f_R, subscript indicating frequency domain with Rayleigh distribution; θ_vib, subscript indicating vibrational phase.

Download Full Size | PDF

A time series of spectral interferograms, H(k, t) is collected synchronously with sound output. We assume the DC term is subtracted from H(k, t) since H(k, t) is commonly preprocessed to remove it. Considering a single vibrating reflector and neglecting any autocorrelation,

The constant, b takes on the value 1 if both sides of the interferometer are collected and 2 if only one side is detected. In practice most swept-source based systems collect both sides using a balanced detector which has the added benefit of canceling the majority of the DC component and autocorrelation artifact. The remaining small DC component is typically subtracted in the preprocessing step noted above. Most spectrometer based systems only collect one side. In this case, the autocorrelation artifact remains and a large DC component must be subtracted in preprocessing.

Next we take the inverse Fourier transform of the interferometric signal, Eq. (1), to generate the complex signal, h(z, t), where z is the optical pathlength difference, related to the tissue depth (Δz) by z = 2nΔz. The magnitude, |h(z, t)|, is the time series of A-scans or line images which are typically plotted against Δz. We would commonly take the mean along time and display the resulting line image $|{\bar{h}(z )} |$ as the structural image.

The remaining processing steps are completed identically at each depth, hence for simplicity let us consider the interferometric signal at a single depth. Then, the L point discrete Fourier transform (DFT) yields,

Small scale vibrations (subresolution) of the reflector (δz) induce small changes in the phase of h(z) over time. The effect of the noise, n_z(z), on the phase of h(z), can be illustrated by complex domain analysis as shown in Fig. 2. It can be seen from Fig. 2 that n_z(z) distorts the ideal phase with the phase noise, ϕ_noise(z), given by

 

Fig. 2. Graphical phase noise analysis in the complex domain. |i(z)| and ϕ_i(z) are the magnitude and the phase of the ideal A-scan without noise. |n_z(z)| and ϕ_n(z) are the magnitude and the phase of the noise, n_z(z). ϕ_n-i(z) is the phase difference between the noise and ideal A-scan.

Download Full Size | PDF

In the limit that |h(z, t)| >> |n_z(z, t)| Eq. (3) simplifies to

The standard deviation of the phase noise in the time-domain is expressed as

As seen in Eq. (5), the phase noise is expressed as the weighted sum of real and imaginary components of n_z(z) which are n_z(z, t)cos(ϕ_n(z, t)) and n_z(z, t)sin(ϕ_n(z, t)), respectively. Because those real and imaginary components are stationary zero mean uncorrelated white Gaussian random processes, the probability density function (PDF) of the phase noise is also stationary zero mean uncorrelated white Gaussian random process since it is the weighted sum of those processes [21], having a variance of 1/(2SNR_z). These results are in agreement with those of Joo et al. [14] in terms of the distribution derived in the context of spectral domain optical coherence phase microscopy.

The next step is to compute the Fourier transform of the phase ∠h(z, t) = ϕ(z,t) along time to yield the magnitude and phase of the vibratory response as a function of depth and frequency. Again, since the same processing is done for each depth, we will explicitly consider the phase at only a single depth, (z_r) i.e. ϕ(z_r,t), without losing any generality. The time dependent phase at z_r is then

The time domain phase in Eq. (7) is described with a sinusoidal function to express an ideal vibration [22]. Generally, any sound can be represented by a sum of sinusoidal waves, hence this does not significantly limit the applicability of the derivation. As shown above, n_t(t) is a stationary white Gaussian noise that has a zero mean and standard deviation of (2SNR_z)^−1/2.

The next step is to take the Fourier transform of ϕ(t) to yield the complex vibratory response, ϕ(f) at z_r as shown in Fig. 1, i.e. ${\mathfrak{F}}$_t(ϕ(z,t)). Assuming M samples in the time series, the frequency domain vibrational signal is expressed as

The measured vibrational amplitude, A_mea, is the magnitude of Eq. (9). As can be seen in Eq. (9), the measured amplitude differs from true amplitude, A_vib, due to noise, i.e.

The statistical properties of A_noise depends on the magnitude of A_vib. When A_vib is equal to 0, A_mea is Rayleigh distributed [16], hence A_noise follows the Rayleigh distribution. The sensitivity of an OCT vibrometry system is typically measured in this limit by observing the noise of a frequency band where there is no vibration. In this case, the root mean square (RMS) value of A_noise is

When A_vib ≠ 0, A_mea is Rician distributed, however it can be approximated as a Gaussian distribution in the limit that the signal is significantly larger than the standard deviation of the noise. It has been shown previously that a s/n > 3 (i.e. ${\textrm{A}_{vib}}\, \ge \,3{\sigma _t}\,/\,{k_0}n\sqrt {2M}$) was sufficient to allow this approximation [16]. In OCT vibrometry, a measured displacement signal is typically in the range of few to one hundred nanometers [1,3,23,24] while noise standard deviations are commonly less than one nanometer [3,24]. Practically, the majority of measurements will have s/n that satisfy this criteria. Thus, the distribution of the detected magnitude (A_mea) in OCT vibrometry can be considered to be Gaussian whose mean and standard deviation are the square root of A_vib² + σ_t²/2Mk₀²n² and σ_t²/2Mk₀²n², respectively [16], leading to the result in Eq. (12).

3. Simulation and experimental results

3.1 Simulation results

Next we decided to numerically simulate interferograms of a vibrating reflector, that is, raw data acquired by M-scan that records the motion of a sample. This allowed us, in a very controlled setting where we know for certain that all of the noise is additive, to verify that the equations we have derived hold. The simulation was executed in MATLAB (Mathworks, Inc).

The set of interferograms were modeled using Eq. (1) such that ρS(k)√(R_RR_S)/b = 1, $\varDelta$z = 1.713 mm, k₀ = 2π/1310 rad. nm⁻¹, and n = 1. The wavelength bandwidth and the total number of spectral channels (N) were set to be 100 nm and 2000, respectively. The vibrational signal δz(t) was modeled using Eq. (7) with A_vib = 10 nm, θ_vib = 30 degree. And the vibrational frequency, f_vib, was defined to be the digital frequency, 0.1, i.e. f_vib = 0.1×A-scan rate. The number of points for the vibrational signal (M) was set to be 1000. The vibrational signal was embedded in the modeled interferograms to make one set of M-scan data as shown in Eq. (1). The k-domain noise n_k(k, t) in Eq. (1) was shaped to be additive Gaussian which has zero mean and standard deviation (σ_k) ranging from 0 to 0.1, and to be stationary so that each interferogram has the same statistical properties. The change in the standard deviation (σ_k) was intended to see the effects of additive noise on vibratory measurements.

Figure 3(a) shows one of the modeled interferograms, H(k, t), in the time-series with a k-domain noise standard deviation of 0.1. Figure 3(b) displays the magnitude response of the fast Fourier transform of that signal, |h(z, t)|, where a reflector is located at index 200 and has a SNR of 40.97 dB calculated from 20log₁₀(s/n). From the complex signal the time-domain phase, ϕ(200,t), can be extracted at the location of the reflector to yield the trace in Fig. 3(c). At this stage the sinusoidal vibration of the reflector is apparent. A fast Fourier transform yields the magnitude response of the vibration |ϕ(200,f)| shown in Fig. 3(d). The vibrational amplitude and phase are extracted at the normalized digital frequency of 0.1, ϕ(200,f = 0.1), showing 9.99 nm and 29.86°, respectively. Note, the displacement in Fig. 3(d) was placed on a log-scale in order to make the noise more visible. The k-domain noise was generated with the MATLAB function of randn multiplied by the desired standard deviation to generate Gaussian distributed random numbers with zero mean. For one M-scan data set, the noise was generated 1000 times with the same mean and standard deviation to be added to each modeled interferogram as seen in Eq. (1). These M-scan data were iteratively acquired to calculate statistical properties of experimental amplitude and phase detection noises. Then the theoretical mean and standard deviation of the amplitude, |ϕ(z, f)|, were calculated for comparison to simulation results.

 

Fig. 3. (a) Example of modeled interferograms with additive noise. (b) A-scan of (a). (c) Time domain vibrational signal extracted from a set of A-scans. (d) Frequency domain vibrational signal transformed from (c).

Download Full Size | PDF

First, the amplitude detection noise (Anoise) was extracted in the frequency region ranging from 0.25 to 0.5 where Avib was equal to 0 and its second moment was compared to Eq. (11) as shown in Fig. 4. In this figure, the top row shows the results from 300 M-scans while the bottom one displays from 3000 M-scans where each M-scan has 251 samples for A_noise. Hence 75,300 (251 × 300) and 753,000 (251 × 3000) samples were used for the top row and the bottom row, respectively. As seen in Figs. 4(a) and 4(c), the second moment of A_noise is on the rise as the k-domain noise standard deviation (σ_k) increases, matching the theoretical line clearly. This increase was caused by the rise of the standard deviation of the time domain phase noise (σ_t) and reduced SNR. Not only the second moment but also the distribution of A_noise was checked by comparing the histogram (blue box) to the fitted Rayleigh distribution (red line) obtained by the MATLAB function of fitdist as shown in Figs. 4(b) and 4(d). As would be expected, the histogram more closely matches the fit as the number of samples increases. The mean and standard deviation of the amplitude detection noise (A_noise) at ϕ(200,f = 0.1) was calculated from the modeled M-scans and plotted in Fig. 5 as a function of k-domain noise standard deviation, σ_k. The bottom x-axis on each plot is σ_k with the corresponding SNR plotted on the top x-axis except Figs. 5(c) and 5(f). Note that the SNR is always well above 9.5 dB (s/n = 3), so we are in the normal distribution limit. The y-axis is the noise mean in Figs. 5(a) and 5(d) and noise STD in Figs. 5(b) and 5(e). In Figs. 5(c) and 5(f), the x-axis and y-axis are A_noise and its frequency, respectively. The top and bottom rows in the figure differ by the number of M-scans used to generate the statistics, 300 top and 3000 bottom. A qualitative comparison of the top and bottom rows shows the expected trend: as we increase the number of samples used to generate the statistic, the results of the simulation more closely match the theory. The noise mean oscillates above and below the theoretical values. The standard deviation likewise oscillates above and below the theoretical line and clearly follows the same trend, showing the increase in a linear scale along σ_k. Also, it is shown that mean and standard deviation of amplitude detection noise increases as expected with an increase in the k-domain noise standard deviation. This is due to increase in the standard deviation (σ_t) of the time domain phase noise n_t(t) in Eq. (7) caused by decrease in SNR in z-domain. The histogram of A_noise was compared to the normal distribution fit in Figs. 5(c) and 5(f) calculated from the data of σ_k = 0.02, showing that the histogram matches the normal distribution as derived in Eq. (12).

 

Fig. 4. Simulation results for A_noise when A_vib = 0. Top row is for 300 simulated M-scans and bottom row for 3000 M-scans. (a, c) are the measured and theoretical second moment as a function of the STD of the noise in the spectral interferogram, σ_k. (b, d) are histograms (blue) with Rayleigh distribution fits (red line) of frequency domain noise from the data of σ_k = 0.02.

Download Full Size | PDF

 

Fig. 5. Simulation results for noise in the measure of A_vib. Top row is for 300 simulated M-scans and bottom row for 3000 M-scans. (a, d) are the measured and theoretical mean (µ_fR) as a function of the STD of the noise in the spectral interferogram, σ_k. (b, e) are the measured and theoretical STD (σ_θfR) as a function of the STD of the noise in the spectral interferogram, σ_k. (c, f) are histograms (blue) with normal distribution fits (red line) of A_noise from the data of σ_k = 0.02.

Download Full Size | PDF

The phase detection noise (θ_noise) was calculated by subtracting the true vibrational phase (θ_vib) from the measured one, θ_mea = ∠ϕ(200, f = 0.1). Its mean and standard deviation were plotted in Fig. 6 as a function of σ_k. As seen in the results of A_noise in Fig. 5, the mean and standard deviation of θ_noise match the theory more closely as the number of used samples increases. Both mean and standard deviation show the same trend while fluctuating slightly about the theory. However, the mean oscillates about 0 while the standard deviation fluctuates about the theoretical line. This shows the same trend as the amplitude detection noise because time domain phase noise directly affects the standard deviation while it has no effect on the mean as seen in Eq. (13). As seen in the amplitude detection noise plots, the increase in the standard deviation of θ_noise results from the decrease in SNR since this brings about rise up of σ_t. The distribution of θ_noise was verified by comparing the histogram acquired from the data of σ_k = 0.02 to the normal distribution fit as seen in Figs. 6(c) and 6(f). It is shown from these plots that the histogram of θ_noise follows the shape of the fitted distribution as derived in Eq. (13). Therefore, these MATLAB simulation results give us confidence that the derived equations are correct within the assumptions of additive noise and the statistical stationary state.

 

Fig. 6. Simulation results for noise in the measure of θ_vib. Top row is for 300 simulated M-scans and bottom row for 3000 M-scans. (a, d) are the measured and theoretical mean (µ_θvib) as a function of the STD of the noise in the spectral interferogram, σ_k. (b, e) are the measured and theoretical STD (σ_θvib) as a function of the STD of the noise in the spectral interferogram, σ_k. (c, f) are histograms (blue) with normal distribution fits (red line) of θ_noise from the data of σ_k = 0.02

Download Full Size | PDF

3.2 Experimental results

In order to verify the derived equations under normal experimental conditions, the movement of a piezo electric element was measured with a swept laser system as described in [25] where a stereo-microscope was used instead of the endoscope The piezo was driven sinusoidally at 4 kHz and measured with A-scan (laser sweep) rate of 128.04 kHz over 50 ms, acquiring 6403 A-scans per M-scan. Two neutral density (ND) filters were used to change the incident light power on the piezo, allowing us to acquire A-scans with 3 different SNRs such as 56.33 dB, 49.25 dB, and 18.09 dB. True vibrational amplitude and phase were estimated by averaging 100 sets of M-scan data with the highest A-scan SNR (56. 33 dB) which are 2.788 nm and −144.556°, respectively.

Figure 7 shows averaged A-scans and frequency domain vibrational responses from 100 M-scans measured with the piezo. And 1^st, 2^nd, and 3^rd columns of Fig. 7 display results from 3 different SNRs such as 56.33 dB, 49.25 dB, and 18.09 dB respectively. Figures 7(d)–7(f) show frequency domain vibrational responses calculated from A-scans seen in Figs. 7(a)–7(c). To obtain frequency domain vibrational responses, time domain vibrational signals were extracted at peak locations red-arrowed in Figs. 7(a)–7(c). The unwrapped phase was fit to a sixth-order polynomial and subtracted in order to reduce low frequency drift. Those subtracted signals were fast Fourier transformed after applying the Hann window for the frequency domain representation. As can be seen in Fig. 7, the more attenuation, the more the peak intensity of the A-scan reduces. This is accompanied by an increase in optical pathlength due to the glass thickness of the ND filters. Multiple reflections in the ND filter also led to a few extra peaks, particularly apparent in the lowest SNR signal. The second moment of the frequency domain noise was acquired for each case with magnitudes ranging from 6 kHz to 10 kHz extracted along 100 M-scan data set by calculating $\sqrt {{\mu ^2}\, + \,{\sigma ^2}}$ with their mean and standard deviation. As shown in Figs. 7(d)–7(f), the predicted RMS noise increases, 6.40 pm, 9.89 pm, and 288.44 pm, as the SNR decrease.

 

Fig. 7. Averaged A-scans and frequency domain vibrational responses measured from 100 M-scan with a piezo electric element driven by 4 kHz. (a, d) are A-scan and 100 frequency domain vibrational responses, respectively when no ND filter is used while (b, e) and (c, f) are from the use of 1 ND and 2 ND filters, respectively. In (d, e, f), thick black line shows the RMS value of the frequency domain noise.

Download Full Size | PDF

As a first step, frequency domain noise (or A_noise when A_vib = 0) was analyzed as shown in Fig. 8 to verify Eq. (11). In this figure, the second moments of A_noise were calculated for each SNR using data from different frequency regions such as 6∼10 kHz, 22∼32 kHz, and 50∼60 kHz. Those frequency regions were selected to see how A_noise changes with frequency since low frequency noise in Figs. 7(d) and 7(e), is clearly not Gaussian white. Above 10 kHz (not shown) the noise appears to converge toward Gaussian white.

 

Fig. 8. Second moments and histograms of amplitude detection noise extracted in the frequency region where A_vib = 0. (a, b) show theoretical and experimental second moments of A_noise and their percentage errors, respectively, where the second moments were measured in three frequency regions, 6∼10 kHz, 22∼32 kHz, and 50∼60 kHz. (c, d, e) display histograms (blue) with Rayleigh distribution fits (red line) of the data with a SNR of 56.33 dB acquired from 6∼10 kHz, 22∼32 kHz, and 50∼60 kHz, respectively.

Download Full Size | PDF

The experimental second moments were compared in Fig. 8(a) to the theoretical ones calculated from Eq. (11), following the trend of the theory. To investigate the difference between experimental and theoretical results in more detail, percent errors were computed using |Experimental second moment - Theoretical second moment|/Theoretical second moment × 100 and displayed in Fig. 8(b). In this plot, the second moments from the data with an SNR of 18.09 dB show similar errors between the three frequency ranges, as might be anticipated from the plot in Fig. 7(f). However, the two data sets with higher SNR shows larger error as the frequency decreases. At relatively low SNR or higher SNR and frequencies above ∼10 kHz the assumption of Gaussian white noise appears to hold. Below the 10 kHz threshold a ∼1/f noise component dominates leading to poorer agreement with the theoretical values. Nevertheless, as seen in Figs. 8(c)–8(e) the histograms for the highest SNR of 56.33 dB are well fit by a Rayleigh distribution. The theoretical distribution matches the experimental distribution at all frequencies and converges toward the theoretical second moment as characterized by the width of the distributions as the frequency increases.

Experimental mean and standard deviation of amplitude and phase detection noises are calculated from the 3 sets of 100 M-scan data and compared to theoretical ones as displayed in Fig. 9. In this figure, experimental results follow the trend of theoretical ones, showing slightly larger values than them. Also, histograms of measured A_noise and θ_noise from the data with a SNR of 56.33 dB are compared to their normal distribution fits in Figs. 9(c) and 9(f), showing again that they follow the correct distribution.

 

Fig. 9. Mean and standard deviation of amplitude and phase detection noises in accordance with Ascan SNR. (a, b) show in picometer scale mean and standard deviation of amplitude detection noise while (d, e) display in degree those of phase detection noise. Blue lines and red circles represent theoretical results calculated using Eqs. (12) and (13), and experimentally measured ones, respectively. (c, f) show histograms (blue box) with normal distribution fits (red line) of measured amplitude and phase detection noises from the data with a SNR of 56.33 dB.

Download Full Size | PDF

To investigate the difference between theoretical and experimental results in Fig. 9 in more detail, the percentage errors of A_noise and θ_noise were calculated using |Theoretical RMS – Experimental RMS|/Theoretical RMS × 100 and displayed in Fig. 10. In this figure, the second moments $(\sqrt {{\mu ^2}\, + \,{\sigma ^2}} )$ were used to consider the mean and the standard deviation simultaneously. It is shown in Fig. 10 that the percentage errors increase with increasing SNR, meaning that experimental results deviate more from theoretical ones. In this experiment, theoretical results were computed using the measured SNR in z-domain and Eqs. (6), (12), and (13). Since the SNR was determined with the magnitudes of the A-scans, any noise source that only impacts the phase would not be accounted for. Hence, other noise sources, which may not be additive white noise, were added to interferometric phases. For this reason, as seen in Figs. 7(d) and 7(e), clearly the total noise is not white, but rather more heavily weighted in the low frequency range. On the other hand, the measurement with the lowest SNR exhibits a noise profile that qualitatively looks white as can be seen in Fig. 7(f). Based on this observation, additive white noise in the k-domain is dominant when SNR is low whereas other non-white noise contributes significantly at higher SNRs. This results in greater deviation from the theoretical predictions for the higher SNR experiments and agrees with the finding from Fig. 8.

 

Fig. 10. RMS errors of amplitude and phase detection noises as function of A-scan SNR. (a) and (b) show percentage errors of amplitude detection noise and phase detection noise, respectively. RMS error was calculated by |Theoretical RMS – Experimental RMS|/Theoretical RMS × 100 where RMS was computed by the square root of mean squared plus standard deviation squared.

Download Full Size | PDF

4. Discussion

The sensitivity of OCT for vibrometry is clearly tied to the signal to noise ratio of the OCT system, i.e. σ_t.=(2SNR)^−1/2. The ultimate limit to the OCT system SNR is shot-noise, σ_k = (2ρR_rSB)^1/2, where B is the noise equivalent bandwidth and S is the source average power. Using the magnitude of Eq. (2) and recognizing that s(0) = NS for an ideal flat-top laser power profile, the SNR can be shown to be, (|i(z)|/E[|n(z)|²])² = NρSR_s/(4b²B). This is equivalent to the result in [10] for b = 1. Assuming an anti-aliasing filter has been used, then B is half of the sampling frequency and the SNR can be rewritten in terms of the sweep frequency (f_sweep) and duty cycle (D) or sweep time ($\varDelta$t),

While we have focused on revealing the statistical properties of the signals in OCT based vibrometry we have not directly addressed sensitivity. In our work [3,26,27] we typically use ${\mu _{{f_R}}} + 3{\sigma _{{f_R}}}$ as a threshold. In other words, any signal below this threshold is neglected. This metric has the advantage that it can readily be measured in each acquisition by considering the mean and standard deviation in a frequency region near the stimulus frequency where there is no stimulus and it corresponds to 98% confidence interval of the Rayleigh distribution. In terms of the relationships of Fig. 1 the sensitivity is,

OCT systems are commonly specified by the signal to noise ratio expected for a perfect reflector (R_s = 1) on a decibel scale. For instance, the maximum SNR of an OCT system might be reported as 107 dB, with a corresponding theoretical shot-noise limit of 110 dB. Developing a similar metric for vibrometry performance would be helpful in comparing systems with different architectures and newly developed algorithms for extracting the vibratory response. Assuming n = M = 1 and using Eqs. (15) and (16), this system has a potential sensitivity of 8.4 pm based on measured SNR and 5.9 pm in the shot-noise limit. Tissue reflectivity is typically in the range of −40 to −80 dB, hence in order to realize the 8.4 pm sensitivity M would need to be set in the range 10⁴–10⁸. For our 128.04 kHz system, that corresponds to an acquisition time of 78 ms – 13 min.

Small changes in the optical pathlength due to a wide range of sources, e.g. air-currents and building vibrations, introduce phase-noise to the system. These tend to be concentrated in the low frequency range, producing an ∼1/f dependence. Obviously, this deviates from our assumption of Gaussian white noise, however at some higher frequency, the mean and standard deviation will become approximately constant with increasing frequency. Beyond this limit our experimental measures match well with the derived equations because the noise is approximately Gaussian white. It is possible to suppress this low frequency phase noise using an interferometer where the reference and sample arm largely co-propagate so that the small scale vibrations are common to both, i.e. a common mode interferometer.

Larger scale changes in the optical pathlength can be introduced by patient/sample motion during data acquisition. Again, this will introduce phase-noise and could violate the condition of a stationary state utilized in the derivation of noise statistics. As we have noted elsewhere [9], under these circumstances, a broad band increase in phase-noise can be observed, however time-domain averaging of the interferometric phase substantially lowers the phase-noise since the motion artifact is incoherent.

A window function (w) is typically multiplied with the time-domain phase before taking the Fourier transform, i.e. w·ϕ(z, t), in order to suppress sidelobes in the frequency domain. Strong sidelobes can obscure weaker signals and produce errors in measured vibrational amplitude and phase. The window has an impact on statistical the properties of amplitude and phase detection noises since it changes intensities of signal and noise in the time domain. The equations in Fig. 1 were derived without using a window, which is equivalent to assuming a rectangular window. More generally we would want to choose a window function based on our tolerance for sidelobe intensities and need for frequency resolution. The frequency domain equations in Fig. 1 can be generalized by replacing M with W²(0)/E_w, where E_w is the total energy of the window function calculated from ${E_w} = \sum\limits_{i = 1}^M {{{({w({t_i})} )}^2}}$ and W(0) is the DC component of W(f) calculated from $W(0) = \sum\limits_{i = 1}^M {w({t_i})}$.

For example, consider Eq. (11) with an arbitrary window and then assuming the Hann window

It is shown from the derivation that the time domain phase noise (n_t(t)) in Eq. (7) is white Gaussian and the magnitude of the frequency domain noise (n_f(f)) is Rayleigh distributed when A_vib is equal to 0. Since the noise power is preserved between these domains, their second moments have the relationship of E[n_t²(t)] = 1/M*E[n_f²(f)], meaning that E[n_f²(f)] can be predicted by calculating E[n_t²(t)]. This relationship assumes that other noise sources that induce small changes in the optical pathlength are not included. Thus, it is important to process data to reduce the effects of other noise sources not only to have the consistent relationship between n_t(t) and n_f(f) but also to increase measurement accuracy.

The preprocessing (second step in Fig. 1) typically done is the subtraction of a background signal to remove the constant term, which will appear at z = 0. The side-lobes from a strong signal at z = 0 can sometimes overlap with signals of interest. We typically measure the background as the mean over a large number of acquisitions and bandpass filter. It is essentially noise free and therefore does not appreciably contribute to the noise. We also typically preprocess the time-domain phase, ϕ(z, t), before calculating the Fourier transform to get ϕ(z, f). At a minimum we calculate and then subtract the mean of ϕ(z, t), to suppress the signal at 0 Hz. This approach is fast and does not contribute to the noise above 0 Hz. In post-processing where computational time is not critical, we sometimes will instead fit ϕ(z, t) to a polynomial and subtract the fit from ϕ(z, t). This suppresses the signal at 0 Hz and removes some of the low frequency drift again without introducing noise. Some researchers have taken the derivative by computing $\varDelta$ϕ(z, t) = ϕ(z, t₂) - ϕ(z, t₁), similar to what is done for Doppler flow measurements. This also removes the 0 Hz portion of the signal, however it introduces an additional √2 noise so that σ_t = (4SNR)^−1/2.

5. Conclusions

In this article, additive noise is analyzed to investigate how it affects the measurement of vibratory amplitude and phase in the frequency domain. The statistical properties were derived in the limit of 0 vibratory signal and for vibratory signals with s/n > 3 and shown to follow a Rayleigh and Gaussian distribution, respectively. The derived equations were verified by MATLAB simulation and experimentally with a vibrating piezo electric element measured with a swept-source OCT system. The results show that derived equations are accurate within the assumption of additive Gaussian white noise, thus providing a theoretical tool to predict the performance of an OCT vibrometry system. These equations also provide a theoretical framework to investigate alternate signal processing algorithms to improve the measurement sensitivity. A measure of sensitivity is proposed that corresponds to the 98% confidence interval for the Rayleigh distribution. The adoption of this metric would make it easier to compare different systems and algorithms.