## Abstract

The high acquisition speed of state-of-the-art optical coherence tomography (OCT) enables massive signal-to-noise ratio (SNR) improvements by signal averaging. Here, we investigate the performance of two commonly used approaches for OCT signal averaging. We present the theoretical SNR performance of (a) computing the average of OCT magnitude data and (b) averaging the complex phasors, and substantiate our findings with simulations and experimentally acquired OCT data. We show that the achieved SNR performance strongly depends on both the SNR of the input signals and the number of averaged signals when the signal bias caused by the noise floor is not accounted for. Therefore we also explore the SNR for the two averaging approaches after correcting for the noise bias and, provided that the phases of the phasors are accurately aligned prior to averaging, then find that complex phasor averaging always leads to higher SNR than magnitude averaging.

Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

## 1. Introduction

Optical coherence tomography (OCT) [1] performs biomedical imaging at high speed and with high sensitivity. State-of-the art OCT systems provide data rates of $\sim$100 million pixels per second and frame rates of $\sim$100-200 frames per second and beyond [2,3]. OCT provides a distinctively high sensitivity of typically $\sim$100 dB, meaning that backscatter signals as weak as $10^{-10}$ of a mirror reflection can still be detected [4–6]. At the same time, OCT enables imaging with a very high dynamic range spanning several tens of decibels between the strongest and the weakest signal in the image.

The strong performance of OCT in detecting weak signals can be improved even more by image processing. Noise in OCT images limits the detection capabilities for weakly scattering structures, and thus a variety of different approaches for reducing OCT image noise have been proposed [7–15]. The high acquisition speeds of modern OCT systems enable the improvement of the signal-to-noise ratio by averaging multiple signals, for instance by fusing OCT frames or even volumes quickly repeated at the same sample position. In recent years, several methods were proposed for improving the detection sensitivity of OCT by averaging the complex-valued OCT signals rather than just averaging their magnitudes [16–20]. In 2013, Szkulmowski and Wojtkowski published a thorough analysis of signals and noise subject to different averaging approaches [21]. In their analysis, the authors found a much stronger reduction of the noise floor by complex averaging as compared to magnitude averaging but also observed a heterogeneous outcome in terms of signal-to-noise performance for different imaging scenarios.

In this article, we set out to answer the question: Which signal averaging approach works best for improving signal-to-noise in OCT images, and when? We first introduce signals and noise in OCT (section 2.1) based on the analysis in Ref. [21]. We then analyze the signals and noise after magnitude and complex averaging, respectively, for a given pixel in an OCT image (section 2.2) and present simple expressions for the resulting signal-to-noise ratios. Next, we compare the somewhat surprising theoretical performance of the averaging schemes for different input signal levels and for different numbers of averaged signals and introduce a noise bias corrected signal-to-noise ratio (SNR) analysis (section 2.3). After substantiating our theoretical analysis with data from simulations (section 3.1) and experimentally acquired OCT data (section 3.2), we conclude the paper with a brief summary of the findings and their implications on actual OCT image processing (section 4). An overview of the terminology as well as a brief section on phase correction required for complex phasor averaging is provided in the appendix (Table 1).

## 2. Analysis of signals and noise in OCT

#### 2.1 OCT signals and noise

OCT signals $S_{OCT}$ can be described by complex phasors of the form

where $A(\mathbf {x},t)$ denotes the amplitude and $\phi (\mathbf {x},t)$ the phase of the OCT signal at location $\mathbf {x}$ and time $t$. Here the amplitude $A(\mathbf {x},t)$ represents the length of the phasor and the phase $\phi (\mathbf {x},t)$ corresponds to the polar angle in the complex plane. The noise $S_{noise}$ in OCT images can be specified by phasors too, namely by random phasors. In the complex plane, the real and imaginary parts of random phasors, $r_{noise}$ and $i_{noise}$, can be described by normal distributions with zero mean and $\sigma ^2$ variance [22]. Hence, the complex noise signal is characterized by a binormal (or Beckmann) distibution $p_{ri}$ in the complex plane [22,23]:In every OCT image pixel, the detected signal is a combination of the actual OCT signal $S_{OCT}$ and the noise $S_{noise}$. By choosing the OCT signal to have $\phi = 0$ and thus to point into the direction of the positive real axis (Fig. 1(e)) similar to Goodman [22], the PDF of the measured amplitudes incorporating both signal and noise is represented by a Rice distribution

OCT amplitude data are usually squared in order to display a quantity proportional to sample reflectivity, $I = A^2$, henceforth called intensity. Therefore it makes sense to investigate the behavior of the PDFs of the squared signal and noise amplitude data. For this purpose, a variable transform of $A = \sqrt []{I}$ is performed such that the PDFs are rendered into $p_I(I) = p_A(A = \sqrt []{I})\vert \frac {dA}{dI}\vert$. The resulting PDFs for noise ($p_I(I_{noise})$) and signal intensities ($p_I(I,I_{noise})$) are

#### 2.2 Averaging magnitudes and phasors

For averaging OCT data, usually the absolute values (magnitudes) of the signals are used. For some applications, the complex signals have been exploited for adding or increasing image contrast in one way or another [31–36]. Here, we systematically analyze the effect of averaging magnitude and complex signals in OCT images.

### 2.2.1 Magnitude averaging

Magnitude averaging uses the absolute values of $N$ spatially and/or temporally separated signals (e.g., $N$ pixels in a 2D or 3D kernel or $N$ repeated measurements at the same spatial pixel location but at different time points):

The PDF representing the average of $N$ data points can be computed numerically by convolving the PDF of a single data point $(N-1)$-times [21]. Unfortunately, no closed form expressions are available for the distribution describing an $N$-fold average of the above PDFs and only some approximations and bounds have been derived [37]. Nevertheless, the mean intensity values and intensity variances can be calculated for Rayleigh distributed noise and Rician signals affected by random noise, respectively, as### 2.2.2 Complex averaging

In contrast to magnitude averaging, complex averaging also includes the phase information into the averaging process. This inclusion exploits the fact that the noise phasors (in absence of a signal) randomly fluctuate about the origin according to a 2D Gaussian PDF with a standard deviation of $\sigma$, see Eq. (2). By averaging $N$ spatially and/or temporally independent complex noise phasors, the standard deviation $\sigma$ of the Beckmann distribution in Eq. (2) is reduced by a factor $1/\sqrt []{N}$ [22]. This reduction of the noise amplitude variance translates into a reduction of both the mean noise intensity as well as of its variance:

Recalling that we initially assumed that the signal vector $S_{OCT}$ was pointing in the direction of the positive real axis (i.e. Fig. 1(e)), the precondition for averaging signals in a complex fashion is that the phase of $S_{OCT}$ remains constant at $\phi = 0$. (Likewise, a different signal phase $\phi = \phi _0$ could be chosen, as long as the phases of all signals to be averaged are aligned at the same angle $\phi _0$. The signal phase ($\phi = 0$) is chosen here out of convenience to ensure a purely real signal which facilitates the calculations described above.) For such perfectly aligned signals in the presence of random noise, the following mean signal intensity and intensity variance will be observed:

#### 2.3 Signal-to-noise performance

The signal-to-noise ratio is the measure of choice to describe the influence of noise on signal measurements. An overview of methods for determining the SNR and in particular the sensitivity of an OCT system has recently been published by Agrawal et al. [24]. The SNR in OCT is typically defined as the ratio of the average signal intensity to the standard deviation of the noise intensity, $\textrm {SNR} = \overline {I}/\sigma _{I_{noise}}$ [24–28] (or alternatively $\overline { \langle I \rangle }/\langle \sigma _{I_{noise}}\rangle$ for averaged signals). As shown in the previous sections and visualized in Fig. 2, the measured signal intensity $\overline {I}$ is the sum of the pure signal intensity $I$ and the noise floor $\overline {I_{noise}}$. Historically, the bias caused by the noise floor is not specifically factored out from the SNR analysis, however some investigations of OCT image data did exclude the background [29,30]. Because of the multiple approaches for OCT data processing, it makes sense to ask, "how important is noise bias correction when measuring SNR?". In the following section, we will investigate the SNR for the total measured signal (i.e. including the contribution of the noise floor) with and without averaging. Subsequently, we will account for the signal offset caused by the noise floor and evaluate the SNR with noise bias correction.

### 2.3.1 Signal-to-noise performance without noise bias correction

Using the signal and noise pairs in Eqs. (13) and (12), (18) and (17), and (22) and (21), and the relation $\overline {I_{noise}} = \sigma _{I_{noise}} = 2\sigma ^2$ in Eqs. (11) and (12), we find the following SNRs for single signals, $N$ magnitude averaged signals, and $N$ complex averaged signals, respectively:

The relative SNR improvements SNR$_{MAG,CPX}$/SNR$_1$ and SNR$_{CPX}$/SNR$_{MAG}$ are plotted for $N$ up to 100 and for typical OCT image dynamics with SNR$_1$ between 5 dB and 20 dB in Fig. 3. In particular for small $N$, complex averaging provides a considerable SNR improvement over magnitude averaging. At the same time, the SNR$_1$ dependence of complex averaging has a considerable impact on its performance: While for strong signals with SNR$_1 \gg 1$, complex averaging outperforms magnitude averaging by a factor $\sqrt []{N}$, the signal-to-noise improvement becomes much less for weaker input signals.

In the next section, we are going to explore the SNR enhancement of the two averaging approaches particularly for very weak signals on the order of the noise background, however still without correcting for the signal bias caused by the noise floor. Will magnitude averaging or complex averaging be superior in terms of recovering such small signals?

### 2.3.2 Averaging domains and borderline SNR

When the signal bias caused by the noise is not accounted for, the SNR performance averaging depends on the approach taken. Equations (25) and (26) revealed a dependence on the number of averaged signals, $N$, and – for complex averaging – on the SNR of the input signals, SNR$_1$. In this section, we are investigating this dependence for different numbers of averaged signals, $N$, and for a great range of signals – from much smaller than the noise variance to several orders of magnitude greater. Recalling the identities for the noise variance and the observed signal intensity in Eqs. (12) and (13), we choose the quantity $I/2\sigma ^2$, which is identical to SNR$_1 - 1$, as the benchmark for the input signal.

Figure 4(a) shows three plots of the signal-to-noise ratios calculated for magnitude averaging (red) and complex averaging (blue) where relatively small signals $I/2\sigma ^2= 1$, 0.5 and 0.1 were chosen as the respective inputs. Note that for $I = 2\sigma ^2$ (left plot), complex averaging always yields greater SNR than magnitude averaging for $N > 1$. However, for $I/2\sigma ^2 < 1$ and a small number of signals $N$, magnitude averaging provides a more effective SNR improvement. Figure 4(b) plots three pairs of SNR profiles for a fixed number of averaged signals ($N = 2$, 10 and 100), this time as a function of the relative signal intensity $I/2\sigma ^2$. Again, for weak input signals, the superior performance of magnitude averaging can be observed, while complex averaging excels beyond a crossover point of $I/2\sigma ^2 = 1/\sqrt []{N}$. This relative borderline signal is plotted for $N$ up to 100 in Fig. 4(c) alongside the borderline input SNR$_1$

which is the input SNR$_1$ for which magnitude averaging and complex averaging perform equally well.An overview of the relative SNR performance SNR$_{CPX}/$SNR$_{MAG}$ at various inputs $I/2\sigma ^2$ and $N$ is shown as a heat map in Fig. 5. At first glance, two domains can be differentiated: For stronger signals $I/2\sigma ^2$ and larger $N$, complex averaging provides a greater SNR improvement than magnitude averaging (see the area in blue in Fig. 5). For signals comparable to or smaller than the noise level, magnitude averaging yields a better SNR enhancement (see the area in red). In between these two domains, the borderline with equal SNR performance of the two approaches is indicated in white, SNR$_{CPX}/$SNR$_{MAG}$ = 1.

### 2.3.3 SNR in absence of a signal (noise floor only)

An interesting scenario is posed by the calculation of the SNR when the signal intensity $I$ is zero. Then, using the relation $\overline {I_{noise}} = \sigma _{I_{noise}} = 2\sigma ^2$ from Eqs. (11) and (12), the SNRs for a single signal, after magnitude averaging and after complex averaging, respectively, become:

### 2.3.4 SNR analysis with noise bias correction

In the limit of a very weak OCT signal, the small meaningful signal contribution sits on top of a comparatively large noise floor. This noise bias is evident as the term related to $\overline {I_{noise}}$ in Eqs. (24–26) and is also visualized in Fig. 2. Thus, it makes sense to actually consider the impact of this noise bias in the SNR analysis – in particular for weak signals. In this section, we introduce and evaluate an SNR analysis with noise bias correction. A similar approach has for instance also been used by Makita et al. [39] and recently been modified by our group [40] to correct for noise contributions in PS-OCT images and improve the computation of the degree of polarization uniformity for weak signals. To some extent, this SNR with noise bias correction resembles the SNR definitions in Refs. [29,30] and has formal similarities with previous definitions of the contrast-to-noise ratios in Refs. [21,27].

The noise bias of the average OCT signals can be removed by subtracting the average noise intensity ($\overline {I_{noise}}$, $\overline {\langle I_{noise}\rangle _{MAG}}$, and $\overline {\langle I_{noise}\rangle _{CPX}}$ in Eqs. (11), (16) and (20), respectively) from the average signal intensity ($\overline {I}$, $\overline {\langle I\rangle _{MAG}}$, and $\overline {\langle I\rangle _{CPX}}$ in Eqs. (13), (18) and (22), respectively). The noise bias corrected average signal intensities then read

## 3. Experimental validation

#### 3.1 Numerical simulation

OCT signals and noise were simulated according to the probability density functions described in Eqs. (3) and (6), respectively. For the noise contribution, binormally distributed complex signals with similar standard deviations $\sigma$ along both the real axis and the imaginary axis were generated. For the signal contribution, the binormal distribution was generated around a real-valued signal $A$ such that the phasor distribution was shifted from the origin of the complex plane to $(A,0)$ (see also Fig. 1(e)). A total of $M=100$ averaged signals were computed for every simulation run. For complex averaging, in every run $N=1$ to $N=100$ complex-valued phasors were averaged. For magnitude averaging, the absolute values were computed for every single noise and signal phasor prior to averaging. Finally, the respective SNRs and SNR$'$s were calculated from the average signal intensity and the variance of the noise signals for every $N$.

Results of the simulations performed in MATLAB (R2014a, MathWorks) are shown in Fig. 7. To showcase the effect of averaging for strong and weak signals, SNR simulations for two distinct relative signal levels of $I/2\sigma ^2=10$ and $I/2\sigma ^2=0.1$ are presented in Fig. 7(a) and (b), respectively. Unlike magnitude averaging, complex averaging markedly decreased the signal intensity but at the same time also reduced noise much more by averaging. For the standard deviation of the noise fluctuations, a decrease by $1/\sqrt []{N}$ and $1/N$ was observed for magnitude and complex averaging, respectively. The resulting SNRs are shown in the four panels in Fig. 7. Here, for averaging up to 100-fold, SNR$_{CPX}$ reveals a better SNR improvement for the strong input signal (Fig. 7(a)) when the noise bias is not corrected for, while magnitude averaging appears more effective for the weaker input signal in Fig. 7(b). After noise bias correction, however, no more dependence upon the input signal level can be observed, and complex averaging provides an $N$-fold increase in SNR$'$ whereas the SNR$'$ is only enhanced by a factor $\sqrt {N}$ after magnitude averaging.

#### 3.2 Averaging experimental OCT data

For demonstrating the effect of the different averaging approaches on OCT images, we used a spectral domain (SD) OCT system in our lab [38]. This polarization sensitive SD-OCT system was used for imaging a stationary eye phantom. This system is based on a multiplexed superluminescent diode ($\lambda$ = 840 nm, $\Delta \lambda$ = 100 nm) as a light source, a free-space Michelson-type interferometer, and a polarization-sensitive detection unit including two spectrometers. For the investigations presented here, only the co-polarized detection channel was analyzed; thus the PS-OCT system was reduced to what would be considered a standard SD-OCT with high-resolution imaging capabilities - 3.6 $\mu$m axial resolution (assuming a refractive index of 1.35), 83 kHz A-scan rate - similar to state-of-the-art commercial SD-OCT scanners for retinal imaging.

The eye phantom, composed of several layers of transparent nail polish on a glass bead to produce a laminar reflectance pattern [38], was imaged using B-scan and M-scan protocols at different levels of attenuation (from 0 dB to -40 dB) to observe the effects of averaging in low-signal conditions. The B-scan protocol scanned a 1 mm lateral range with 100 repeats to generate cross-sectional images of the phantom (Fig. 8(a-d)). Representative depth profiles are shown in Fig. 8(e).

The M-scan protocol was used to acquire 1000 repeated A-lines from a fixed position within the phantom for a precise comparison of averaging methods. The M-scan protocol was chosen to minimize phase differences between consecutive A-lines and ensure that averaging represents a best-case real-world scenario. The 1000 A-lines were split into 10 temporally independent bins of 100 consecutive A-lines each. Complex and magnitude averaging for $N$ up to 100 was performed on each of the 10 bins and the resulting signal curves were averaged together before calculating SNR and SNR$'$ in order to yield more stable SNR curves, particularly for lower numbers of averages. Here, the noise variances were computed from the noise pixel data of 100 consecutive A-lines, similar to the simulation above. Three characteristic SNR and SNR$'$ curves from the magnitude-averaging dominant, borderline, and complex-averaging dominant domains are shown in Fig. 8(f) and (g). Theoretical SNR profiles based on SNR$_{1}$ as described in Eqs. (25) and (26) demonstrate good agreement with the experimental data in Fig. 8(f). These results show that the theoretical magnitude-averaging dominant domain is reproducible with real world OCT data when the noise floor induced bias is not accounted for. Similarly, the SNR$'$ plots from the experimental data with noise bias correction follow the expected slopes of $\sqrt {N}$ and $N$ (Fig. 8(g)).

## 4. Discussion and conclusion

Signal averaging is one of the most commonly used procedures for OCT image processing. Here, we investigated the SNR performance of magnitude and complex averaging in theory and backed our findings with experimental results from simulations and actual OCT image data. We also studied the effect of the background noise bias on the measured SNR and calculated a noise bias corrected SNR. In the following paragraphs, the main observations are summarized and discussed in order to provide a better understanding of the strengths and weaknesses of the two approaches.

#### 4.1 Impact of averaging on noise and signals

Complex averaging reduces the noise variance by $1/N$, and thus provides a much stronger noise reduction than magnitude averaging, which only reduces the noise by $1/\sqrt []{N}$. At the same time, complex averaging also reduces the signal – in contrast to magnitude averaging, which maintains the signal level and only reduces noise. The main cause of the average signal decrease by complex averaging is the reduction of the average noise level. Both signal and noise characteristics have to be considered to understand the impact of averaging on OCT data.

#### 4.2 Leaving or removing the noise floor bias for the SNR calculation

The signal-to-noise ratio is a commonly used measure to assess the image quality and to describe the system performance of OCT machines. Usually, the ratio of the average intensity of a signal peak and the standard deviation of the noise intensity is used to calculate the SNR [24,28]. While it is rather clear how to estimate the noise variance (e.g., by considering the noise intensity in a signal-free image area), the definition of the average intensity is not that obvious. As visualized in Fig. 2, taking the entire measured signal intensity from the zero line to the (average) peak intensity $\overline {I}$ includes two portions, namely the actual signal intensity term $I$ and the background noise level $\overline {I_{noise}}$.

For relatively strong signals, the measured signal is dominated by the intensity term $I$ and $\overline {I} \approx I \gg \overline {I_{noise}}$. In contrast, when the actual signal contribution to the measured signal is small (i.e. $I \ll \overline {I}$), the noise floor bias $\overline {I_{noise}}$ prevails. We hence investigated the signal-to-noise performance of magnitude and complex averaging for (a) leaving the noise floor bias as part of the measured signal and (b) correcting for the noise bias. As discussed in detail in the following sections, similar results were observed for both SNR analyses when strong signals were investigated. However, for rather weak signals, the impact of $\overline {I_{noise}}$ manifested in the observation of dissimilar SNR characteristics.

#### 4.3 Noise-afflicted signal-to-noise ratios after averaging strong and weak signals

For sufficiently large signals $\overline {I} = I + \overline {I_{noise}}$, the $N\cdot {\textrm{SNR}}_1$ term in Eq. (26) dominates such that complex averaging converges on an $N$-fold SNR improvement, while magnitude averaging only enhances the SNR $\sqrt []{N}$-fold. In contrast, for recovering weak signals comparable to or smaller than the noise, magnitude averaging appeared to outperform complex averaging – in particular for small $N$. When the input SNR$_1$ was just slightly larger than unity (i.e., $I/2\sigma ^2$ just slightly greater than zero), the right-hand side of Eq. (26) was essentially neutralized such that the SNR improvement was far less than the $\sqrt []{N}$-fold enhancement yielded by magnitude averaging, especially in the limit of small $N$. Hence, magnitude averaging could be considered more effective for boosting small signals as they may be found in the outer nuclear layer in the retina when the noise floor induced signal bias is not accounted for. However, by taking a look at the SNR in absence of an actual signal ($I = 0$), odd SNR characteristics were observed (see section 2.3.3) which underscored the importance of a noise bias correction.

#### 4.4 SNR calculations with noise bias correction

When SNRs are calculated in the classical way described in section 2.3.1, the average noise floor level $\overline {I_{noise}}$ contributes to the intensity measured as "signal". This noise bias impacts the measured SNR, in particular for weak signals. Akin to noise offset removal in PS-OCT processing [39,40], we performed SNR calculations with noise bias correction in sections 2.3.4 and 3. By using this modified approach, a superior SNR performance was always observed for complex averaging, regardless of the input signal strength. Compared to the noise bias corrected SNR of a single signal prior to averaging, magnitude and complex averaging improved the SNR (after noise bias correction) by a factor of $\sqrt {N}$ and $N$, respectively. The theoretically predicted performance agreed well with the performance observed in simulated and experimental OCT data (Figs. 7 and 8). This finding suggests that noise bias correction definitely is an important processing step in SNR analyses – especially when it comes to averaging weak signals.

#### 4.5 Implications for practical OCT image averaging

State-of-the-art swept source and spectral domain OCT devices deliver complex-valued signals (see Eq. (1)) right out of the box. Hence both magnitude and complex averaging can be easily implemented using the image data. For effective complex averaging, it is imperative that the phases of the signals are accurately aligned before averaging as described in more detail in Appendix B. While this means additional computational steps, it may be well worth the effort when the SNR is to be increased in images of scattering structures such as the retina where OCT has been frequently applied. Retinal OCT images may easily span a dynamic range of 40 dB with hyperscattering structures such as the nerve fiber layer and the pigment epithelium that can serve as landmarks for phasor alignment [18]. Other retinal tissues such as the ganglion cell layer and the outer nuclear layer [41], and also structures in the vitreous, which are usually much less reflective [42], may then be visualized by averaging multiple OCT images together. Complex averaging has also been found promising to detect signals in settings with multiple scattering [20].

Magnitude averaging on the other hand can be implemented at less computational expense and may be the averaging approach of choice for images containing mostly weak signals of interest, which would render phase alignment for complex averaging difficult or even impossible. Also in scenarios where images include a lot of varying motion, e.g. when visualizing weakly scattering structures such as flow in lymphatic vessels [43], albeit theoretically less effective by $\sqrt {N}$ in terms of SNR improvement, magnitude averaging may likely trump complex averaging.

#### 4.6 Future perspectives

Finally, we would like to point out that, while the analysis presented here particularly focused on the SNR improvement for different averaging approaches, it may be interesting to also investigate other image metrics such as their contrast-to-noise characteristics and/or their efficiency in terms of speckle reduction. For this purpose and to explore scenarios with tissue-like scattering, recently described OCT signal models could be particularly interesting candidates [28,44–46]. Additionally, OCT images are often displayed on a logarithmic scale in order to compress signals covering a large dynamic range. The logarithm pulls up weak signals but also the noise floor [21], such that the SNR performance for averaged logarithmic amplitude data will deviate from that for uncompressed OCT data discussed here. An in-depth analysis similar to the one performed for linear data in section 2.2 may provide more insight in the advantages and disadvantages of magnitude- and phasor-based averaging approaches for enhancing OCT image quality.

## Appendix A – Glossary

The table below provides an overview of the variables and symbols used in this article.

## Appendix B – Complex signal averaging in case of axial motion: Effect of axial motion on complex averaging

Keeping the signal phase $\phi$ constant is essential for complex signal averaging. If $\textrm {d}^2\phi (\mathbf {x},t)/\textrm {d}\mathbf {x}\textrm {d}t \neq 0$, the complex sum of the signal vectors will have a length smaller than the sum of their absolute values. In other words, if the phases of the complex signals are not matched prior to averaging, the length of the average signal vector will be reduced. The condition $\Delta \phi = 0$ is fulfilled for simultaneously acquired signals at the same $z$-position within the same speckle. When averaging signals among different speckles acquired at the same time point, first the phase offset between different speckles has to be eliminated [48]. Similarly, when signals acquired at the same location $\mathbf {x}$ but at different time points $t$ are averaged (e.g., signals from repeatedly acquired OCT images), the influence of axial motion has to be removed prior to averaging [18].

Axial motion introduces a proportional phase shift in the complex signals (Fig. 9(a)). This phase shift has been extensively exploited for velocity measurements in Doppler OCT, for displacement measurements in OCT elastography and, recently most prominently, for some OCT angiography approaches [32–36]. The phase shift of an OCT signal is proportional to $2k\Delta z$ where $k =2\pi /\lambda$ is the central wavenumber and $\Delta z$ is the axial displacement between successive measurements. When $N$ complex signals with displacements $\Delta z_j$ relative to the first signal are averaged in a complex fashion, the resulting relative signal reduction corresponds to the penalty factor $P(N)$:

where $W_j = 2\Delta z_j / \lambda$ denotes the displacement of the $j$-th signal in terms of wavelengths. In case of constant axial motion, the displacement will increase by $(j-1)\Delta z$ for the $j$-th signal compared to the first signal. The penalty factor $P(N)$ when averaging $N$ signals with relative separation $W_j = W$ is## Compensation of axial motion

In order to remove the detrimental effect of axial motion and thus to fully exploit the SNR advantage of complex averaging, the phases $\phi _j$ of the signals to be averaged have to be aligned. Ju et al. proposed two approaches: one compensating phase offsets A-line-wise and another one compensating phase offsets in a pixel-wise fashion [18]. Assuming a constant "bulk" displacement of the simultaneously acquired signals within one A-scan, the bulk phase shift $\Delta \phi _{bulk}$ between the $m$-th A-line in a reference frame ($ref$) and the corresponding $m$-th A-scan in other frames to be averaged can be calculated by computing the angle of the weighted complex signal average along every A-scan:

If now, in a third step, the locally averaged phase difference map (Eq. (42)) is used to correct the phasors of the $j$-th frame with respect to the reference frame,

## Funding

Austrian Science Fund (P25823-B24, P26553-N20); H2020 European Research Council (ERC Starting Grant 640396 OPTIMALZ).

## Acknowledgements

The authors would like to acknowledge fruitful discussions with Antonia Lichtenegger, Danielle J. Harper, Pablo Eugui, Laurin Ginner, and Florian Beer. We would like to express our particular gratitude to Stanislava Fialova for her contributions to the imaging system and eye phantom. We also acknowledge the thoughtful feedback and inputs provided by the anonymous reviewers.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. Fujimoto, “Optical coherence
tomography,” Science **254**(5035),
1178–1181
(1991). [CrossRef]

**2. **W. Drexler, M. Liu, A. Kumar, T. Kamali, A. Unterhuber, and R. A. Leitgeb, “Optical coherence
tomography today: speed, contrast, and
multimodality,” J. Biomed.
Opt. **19**(7),
071412 (2014). [CrossRef]

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

**4. **R. Leitgeb, C. K. Hitzenberger, and A. F. Fercher, “Performance of fourier
domain vs. time domain optical coherence
tomography,” Opt. Express **11**(8),
889–894
(2003). [CrossRef]

**5. **M. A. Choma, M. V. Sarunic, C. Yang, and J. A. Izatt, “Sensitivity advantage
of swept source and Fourier domain optical coherence
tomography,” Opt. Express **11**(18),
2183–2189
(2003). [CrossRef]

**6. **J. F. de Boer, B. Cense, B. H. Park, M. C. Pierce, G. J. Tearney, and B. E. Bouma, “Improved
signal-to-noise ratio in spectral-domain compared with time-domain
optical coherence tomography,” Opt.
Lett. **28**(21),
2067–2069
(2003). [CrossRef]

**7. **A. Sakamoto, M. Hangai, and N. Yoshimura, “Spectral-Domain Optical
Coherence Tomography with Multiple B-Scan Averaging for Enhanced
Imaging of Retinal Diseases,”
Ophthalmology **115**(6),
1071–1078.e7
(2008). [CrossRef]

**8. **P. Puvanathasan and K. Bizheva, “Interval type-II fuzzy
anisotropic diffusion algorithm for speckle noise reduction in optical
coherence tomography images,” Opt.
Express **17**(2),
733–746
(2009). [CrossRef]

**9. **W. Wu, O. Tan, R. R. Pappuru, H. Duan, and D. Huang, “Assessment of
frame-averaging algorithms in OCT image
analysis,” Ophthalmic Surg. Lasers
Imag. Retin. **44**(2),
168–175
(2013). [CrossRef]

**10. **C.-L. Chen, H. Ishikawa, G. Wollstein, R. A. Bilonick, L. Kagemann, and J. S. Schuman, “Virtual Averaging
Making Nonframe-Averaged Optical Coherence Tomography Images
Comparable to Frame-Averaged Images,”
Trans. Vis. Sci. Tech. **5**(1), 1
(2016). [CrossRef]

**11. **D. Xu, N. Vaswani, Y. Huang, and J. U. Kang, “Modified compressive
sensing optical coherence tomography with noise
reduction,” Opt. Lett. **37**(20),
4209–4211
(2012). [CrossRef]

**12. **M. Sugita, S. Zotter, M. Pircher, T. Makihira, K. Saito, N. Tomatsu, M. Sato, P. Roberts, U. Schmidt-Erfurth, and C. K. Hitzenberger, “Motion artifact and
speckle noise reduction in polarization sensitive optical coherence
tomography by retinal tracking,”
Biomed. Opt. Express **5**(1),
106–122
(2014). [CrossRef]

**13. **H. Zhang, Z. Li, X. Wang, and X. Zhang, “Speckle reduction in
optical coherence tomography by two-step image
registration,” J. Biomed. Opt. **20**(3), 036013
(2015). [CrossRef]

**14. **A. C. Chan, K. Kurokawa, S. Makita, M. Miura, and Y. Yasuno, “Maximum a posteriori
estimator for high-contrast image composition of optical coherence
tomography,” Opt. Lett. **41**(2),
321–324
(2016). [CrossRef]

**15. **B. Tan, A. Wong, and K. Bizheva, “Enhancement of
morphological and vascular features in OCT images using a modified
Bayesian residual transform,” Biomed.
Opt. Express **9**(5),
2394–2406
(2018). [CrossRef]

**16. **P. H. Tomlins and R. K. Wang, “Digital phase
stabilization to improve detection sensitivity for optical coherence
tomography,” Meas. Sci.
Technol. **18**(11),
3365–3372
(2007). [CrossRef]

**17. **J. W. Jacobs and S. J. Matcher, “Digital phase
stabilization for improving sensitivity and degree of polarization
accuracy in polarization sensitive optical coherence
tomography,” Proc. SPIE **7889**, 788938
(2011). [CrossRef]

**18. **M. J. Ju, Y.-J. Hong, S. Makita, Y. Lim, K. Kurokawa, L. Duan, M. Miura, S. Tang, and Y. Yasuno, “Advanced multi-contrast
Jones matrix optical coherence tomography for Doppler and polarization
sensitive imaging,” Opt.
Express **21**(16),
19412–19436
(2013). [CrossRef]

**19. **T. Pfeiffer, W. Wieser, T. Klein, M. Petermann, J. P. Kolb, M. Eibl, and R. Huber, “Flexible A-scan rate
MHz OCT: computational downscaling by coherent
averaging,” Proc. SPIE **9697**, 96970S
(2016). [CrossRef]

**20. **L. Thrane, S. Gu, B. J. Blackburn, K. V. Damodaran, A. M. Rollins, and M. W. Jenkins, “Complex decorrelation
averaging in optical coherence tomography: a way to reduce the effect
of multiple scattering and improve image contrast in a dynamic
scattering medium,” Opt. Lett. **42**(14),
2738–2741
(2017). [CrossRef]

**21. **M. Szkulmowski and M. Wojtkowski, “Averaging techniques
for OCT imaging,” Opt. Express **21**(8),
9757–9773
(2013). [CrossRef]

**22. **J. W. Goodman, * Statistical
Optics* (John Wiley &
Sons, 2000).

**23. **P. Beckmann, “Statistical
distribution of the amplitude and phase of a multiply scattered
field,” J. Res. Natl. Bur. Stand. (U.
S.) **66D**(3),
231–240
(1962). [CrossRef]

**24. **A. Agrawal, T. J. Pfefer, P. D. Woolliams, P. H. Tomlins, and G. Nehmetallah, “Methods to assess
sensitivity of optical coherence tomography
systems,” Biomed. Opt. Express **8**(2),
902–917
(2017). [CrossRef]

**25. **X. Zhu, Y. Liang, Y. Mao, Y. Jia, Y. Liu, and G. Mu, “Analyses and
calculations of noise in optical coherence tomography
systems,” Front. Optoelectron.
China **1**(3-4),
247–257
(2008). [CrossRef]

**26. **J. Izatt and M. Choma, “Theory of optical
coherence tomography,” in Optical
Coherence Tomography - Technology and Applications,
(2008), pp.
47–72.

**27. **M. Szkulmowski, I. Gorczynska, D. Szlag, M. Sylwestrzak, A. Kowalczyk, and M. Wojtkowski, “Efficient reduction of
speckle noise in Optical Coherence Tomography,”
Opt. Express **20**(2),
1337–1359
(2012). [CrossRef]

**28. **M. Sugita, A. Weatherbee, K. Bizheva, I. Popov, and A. Vitkin, “Analysis of scattering
statistics and governing distribution functions in optical coherence
tomography,” Biomed. Opt.
Express **7**(7),
2551–2564
(2016). [CrossRef]

**29. **D. M. Stein, H. Ishikawa, R. Hariprasad, G. Wollstein, R. J. Noecker, J. G. Fujimoto, and J. S. Schuman, “A new quality
assessment parameter for optical coherence
tomography,” Br. J.
Ophthalmol. **90**(2),
186–190
(2006). [CrossRef]

**30. **A. Lozzi, A. Agrawal, A. Boretsky, C. G. Welle, and D. X. Hammer, “Image quality metrics
for optical coherence angiography,”
Biomed. Opt. Express **6**(7),
2435–2447
(2015). [CrossRef]

**31. **M. Wojtkowski, A. Kowalczyk, R. Leitgeb, and A. F. Fercher, “Full range complex
spectral optical coherence tomography technique in eye
imaging,” Opt. Lett. **27**(16),
1415–1417
(2002). [CrossRef]

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

**33. **S. Makita, K. Kurokawa, Y.-J. Hong, M. Miura, and Y. Yasuno, “Noise-immune complex
correlation for optical coherence angiography based on standard and
Jones matrix optical coherence tomography,”
Biomed. Opt. Express **7**(4),
1525–1548
(2016). [CrossRef]

**34. **K. V. Larin and D. D. Sampson, “Optical coherence
elastography - OCT at work in tissue biomechanics
[Invited],” Biomed. Opt.
Express **8**(2),
1172–1202
(2017). [CrossRef]

**35. **J. Xu, S. Song, Y. Li, and R. K Wang, “Complex-based OCT
angiography algorithm recovers microvascular information better than
amplitude- or phase-based algorithms in phase-stable
systems,” Phys. Med. Biol. **63**(1), 015023
(2017). [CrossRef]

**36. **B. Braaf, S. Donner, A. S. Nam, B. E. Bouma, and B. J. Vakoc, “Complex differential
variance angiography with noise-bias correction for optical coherence
tomography of the retina,” Biomed. Opt.
Express **9**(2),
486–506
(2018). [CrossRef]

**37. **S. Nadarajah, “A review of results on
sums of random variables,” Acta Appl.
Math. **103**(2),
131–140
(2008). [CrossRef]

**38. **S. Fialová, M. Augustin, M. Glösmann, T. Himmel, S. Rauscher, M. Gröger, M. Pircher, C. K. Hitzenberger, and B. Baumann, “Polarization properties
of single layers in the posterior eyes of mice and rats investigated
using high resolution polarization sensitive optical coherence
tomography,” Biomed. Opt.
Express **7**(4),
1479–1495
(2016). [CrossRef]

**39. **S. Makita, Y.-J. Hong, M. Miura, and Y. Yasuno, “Degree of polarization
uniformity with high noise immunity using polarization-sensitive
optical coherence tomography,” Opt.
Lett. **39**(24),
6783–6786
(2014). [CrossRef]

**40. **B. Baumann, M. Augustin, A. Lichtenegger, D. J. Harper, M. Muck, P. Eugui, A. Wartak, M. Pircher, and C. K. Hitzenberger, “Polarization-sensitive
optical coherence tomography imaging of the anterior mouse
eye,” J. Biomed. Opt. **23**(08), 1
(2018). [CrossRef]

**41. **Z. Liu, K. Kurokawa, F. Zhang, J. J. Lee, and D. T. Miller, “Imaging ganglion cells
in the living human retina,” Proc. Nat.
Acad. Sci. **114**(48),
12803–12808
(2017). [CrossRef]

**42. **J. J. Liu, A. J. Witkin, M. Adhi, I. Grulkowski, M. F. Kraus, A.-H. Dhalla, C. D. Lu, J. Hornegger, J. S. Duker, and J. G. Fujimoto, “Enhanced Vitreous
Imaging in Healthy Eyes Using Swept Source Optical Coherence
Tomography,” PLoS One **9**(7), e102950
(2014). [CrossRef]

**43. **C. Blatter, E. F. J. Meijer, A. S. Nam, D. Jones, B. E. Bouma, T. P. Padera, and B. J. Vakoc, “In vivo label-free
measurement of lymph flow velocity and volumetric flow rates using
Doppler optical coherence tomography,”
Sci. Rep. **6**(1),
29035 (2016). [CrossRef]

**44. **M. Almasian, T. G. van Leeuwen, and D. J. Faber, “OCT amplitude and
speckle statistics of discrete random media,”
Sci. Rep. **7**(1),
14873 (2017). [CrossRef]

**45. **T. B. Dubose, D. Cunefare, E. Cole, P. Milanfar, J. A. Izatt, and S. Farsiu, “Statistical models of
signal and noise and fundamental limits of segmentation accuracy in
retinal optical coherence tomography,”
IEEE Trans. Med. Imaging **37**(9),
1978–1988
(2018). [CrossRef]

**46. **M. Sugita, R. A. Brown, I. Popov, and A. Vitkin, “K-distribution
three-dimensional mapping of biological tissues in optical coherence
tomography,” J. Biophotonics **11**(3), e201700055
(2018). [CrossRef]

**47. **W. C. van Etten, * Introduction to Random Signals
and Noise* (John Wiley &
Sons, 2005).

**48. **Y. Lim, M. Yamanari, S. Fukuda, Y. Kaji, T. Kiuchi, M. Miura, T. Oshika, and Y. Yasuno, “Birefringence
measurement of cornea and anterior segment by office-based
polarization-sensitive optical coherence
tomography,” Biomed. Opt.
Express **2**(8),
2392–2402
(2011). [CrossRef]