Introduction

There are a variety of extensions of Optical Coherence Tomography (OCT) that take advantage of the interferometric phase to garner additional information. Some example applications include the measurement of blood flow with Doppler OCT [1–3] and tissue mechanical properties with OCT elastography [4, 5]. In particular, we are interested in using a technique called Phase-sensitive Optical Coherence Tomography (PhOCT) to accomplish vibrometry with picometer sensitivities. This exquisite sensitivity is enabling 3-D vibrometry in the intact cochlea of animal models of hearing, revealing new insights into cochlea function. In recent work [6], we have quantitatively measured differential motion in the mouse Organ of Corti using techniques collectively called Volumetric Optical Coherence Tomography and Vibrometry (VOCTV).

PhOCT has most often been prosecuted with spectrometer based or Spectral-Domain OCT (SDOCT) systems owing to their inherent phase stability [7–11]. Nevertheless, there are several advantages to the swept-laser source or Swept-Source OCT (SSOCT) architecture. For example, swept laser sources typically have very narrow instantaneous line widths that provide long coherence length and therefore low sensitivity rolloff as a function of depth. Likewise, they are compatible with architectures using some variant of the Mach-Zehdner interferometer and balanced detection. Balanced detection cancels the DC part of the signal and suppresses common mode noise. Unfortunately, swept laser sources used for OCT typically suffer sweep instabilities that lead to phase noise and therefore reduced sensitivity.

Specifically, the mechanical wavelength tuning of the laser source causes the wavenumber to be swept nonlinearly and to vary within subsequent sweeps which are referred to as intra- and inter-sweep variability, respectively. Several hardware and software approaches have been proposed to mitigate intra/inter-sweep variability. Postaid et al. and Choi et al. [12, 13] utilized a variable rate optical clock as a sampling clock for the Analog to Digital Converter (ADC) to remove intra-sweep variability. This method simplifies not only required hardware in the system but also the necessary post-processing. However, the ADC suffers from clock jitter generated by the variability of the optical clock. The clock jitter engenders a frequency dependent degradation of the signal-to-noise ratio (SNR) of the input signal [14]. Since the input signal is an interferometric signal, its frequency is a function of the optical pathlength difference. Therefore, the signal suffers more SNR degradation as the optical pathlength difference becomes larger, consequently deteriorating phase sensitivity as well as image intensity. There are software based approaches to deal with intra-sweep variability. Huber et al. [15] used a nearest-neighbor check algorithm to linearize the interferometric signal. This method is very efficient but not accurate. A more accurate algorithm was proposed by Gora et al. [16]. This method involves linearizing the interferometric signal using the phase of a reference interferometer. The phase of the reference interferometer is extracted by Fast Fourier Transform-Inverse Fast Fouirer Transform (FFT-IFFT) based Hilbert transform. This method is very accurate but requires substantial computation due to the FFT-IFFT process.

Several hardware based approaches have been used to remove inter-sweep variability. Vakoc et al. [17] and Baumann et al. [18] used a phase reference signal created by additional optics placed in the sample arm. The phase calibration signal is then used to measure and remove the inter-sweep variability. However, this approach requires careful adjustment so that the magnitude of the phase reference signal is sufficient for calibration but weak enough to avoid autocorrelation artifacts. A new type of hardware approach was proposed by Choi et al. [13]. The method places a Fiber Bragg Grating (FBG) in the optical system to act as a wavenumber reference signal. The inter-sweep variability is compensated by using the wavenumber reference signal. However, this method reflects some portion of the interferometric signal and needs careful selection and alignment of the FBG to work well. From the software point of view, the inter-sweep variability can be eliminated by a correlation method used by Braaf et al. [19]. In this method, two adjacent phases of k-clock signals are correlated with each other to calculate a shift magnitude. The shift magnitude is applied to the secondly measured k-clock signal to remove variations between sweeps. This method does not need additional optical hardware but has a considerable computational load due the calculation of correlation between two adjacent sweeps. In summary, hardware only based methods do not require significant post-processing but they increase hardware and alignment complexities and can be accompanied by unwanted artifacts or signal degradation. On the other hand, software based methods require considerable post-processing to guarantee accuracy and additional hardware to generate a reference interferometer signal, but they do not degrade the signal or introduce artifacts.

In our work which requires very high phase stability, we have adopted an approach similar to Gora et al. [16], where each sweep of the laser is calibrated by using a reference interferometer. The Hilbert phase of the reference interferometer is proportional to the instantaneous wavenumber of the laser and can be used to resample the interferometric signal such that δk is constant. However, the traditional Fast Fourier Transform-Inverse Fast Fourier Transform (FFT-IFFT) based Hilbert transform comes with considerable computational overhead. An alternative approach is to use a complex Finite Impulse Response (FIR) filter to extract the phase of the reference interferometer. The complex FIR filter is composed of two FIR filters used to generate an analytic signal [20]. The use of the FIR filter has the advantage of digital hardware efficiency compared to the FFT-IFFT method [21]. Other methods that have been utilized include the delay element and Hilbert FIR filter approach [22], and a bandpass and Hilbert FIR filter approach [23]. The complex FIR filter has been implemented in ultrasound imaging to detect the envelope of a beam formed signal and in OCT imaging to detect blood velocity [22, 23]. However, to the best our knowledge we present here the first application of the complex FIR filter to calibrate via a reference interferometer each sweep of a swept laser for SSOCT. We explicitly consider the accuracy of the phase of the different FIR based Hilbert transforms. In the context of PhOCT, we consider the impact of phase errors on phase noise. Furthermore we mitigate the computational load of the FIR filters by utilizing a field programmable gate array (FPGA) which can reduce processing time by working in parallel with the digitization [24]. Finally, we make use of the magnitude of the Hilbert transform to calculate an “ideal” window function, which allows us to carefully control the sidelobes of our axial point spread function.

Method

Proposed architecture

The proposed architecture is composed of two parts: phase extraction and resampling (as shown in Fig. 1). The phase extraction part uses a complex FIR filter to make an analytic signal, and coordinate rotation digital computer (CORDIC) to calculate the phase of the analytic signal. CORDIC allows a simple and efficient hardware implementation by using shift and add operations to calculate trigonometric functions [25]. The resampling part linearizes an interferometric signal henceforth called the OCT signal, by linearly interpolating it with the extracted reference interferometer phase and a linearized phase.

 

Fig. 1 The proposed architecture for continuous real-time spectral calibration and linearization.

Download Full Size | PDF

Generation and phase extraction of an analytic signal

The phase extraction from the reference interferometer is the most important task in the proposed architecture because it affects the accuracy of the spectral wavenumber calibration. Specifically, if the extracted phase has errors, those errors are propagated through the calculation of the depth resolved interferometric phase in the PhOCT image and reduces the system sensitivity. Therefore, it is important to extract the phase from reference interferometer accurately.

There are several types of FIR filter based methods for generating an analytic signal such as the delay element and Hilbert FIR approach, the bandpass and Hilbert FIR approach, and the complex FIR approach as shown in Fig. 2. Delay element and Hilbert FIR approach uses a delay element, Z^-d, for the in-phase signal to synchronize with the quadrature signal, and Hilbert FIR to make the quadrature signal. The bandpass and Hilbert FIR approach utilizes a bandpass FIR not only for synchronization but also for matching the bandwidth with the quadrature signal. The complex FIR approach consists of two filters: one from real coefficients, the other from imaginary coefficients. Ideally, these two filters exhibit an exact π/2 phase relationship as well as a matched frequency spectrum. The π/2 relationship is important for an analytic signal to be correctly defined as in Euler's formula. However, this property is not conserved in the other two FIR based methods, hence in this respect the complex FIR approach is the most accurate method among them.

 

Fig. 2 Block diagram illustrating several types of FIR filter approaches to generate an analytic signal (a): Delay element and Hilbert FIR approach (b): Bandpass and Hilbert FIR approach (c): complex FIR approach. In (a), Z^-d is a delay element to delay an input signal by d samples and Hilbert FIR is Hilbert transform implemented by FIR filter. The complex FIR filter has two kinds of filter coefficients: real and imaginary coefficients which are described as Complex FIR - Real and Complex FIR - Imag in (c), respectively.

Download Full Size | PDF

In order to analyze and compare the FIR approaches mathematically, a reference interferometer signal was modeled using the following equation,

Equation (2) shows that the reference interferometer has both a bandpass signal and a DC signal. In principle, the DC signal is cancelled when a balanced detector is used. However, the imperfect spectral response of the fiber coupler and difference in responsivity of photodetectors in a balanced detector keep the DC signal from being removed completely. The imperfect DC cancelation is the origin of the parameter α. The DC component should be removed in order to extract the phase from the reference interferometer accurately.

Assuming the Hilbert FIR has no delay, the impulse response for the delay and Hilbert FIR approach in Fig. 2(a) can be expressed as

with the corresponding phase

As shown in Eq. (4), the DC component, αA, included in the reference interferometer signal is not removed by the delay element and thus ends up with the cosine term in the calculation of the arctangent and corrupts the phase estimation. On the other hand, the αA term is eliminated by Hilbert FIR filter whose magnitude response is characteristic of a bandpass filter. The bandpass and Hilbert FIR approach in Fig. 2(b) can overcome this issue and remove the DC component in the reference interferometer signal. Nevertheless it has an inherent phase error which prevents it from reproducing the correct π/2 phase relationship between the bandpass and the Hilbert FIR filters. That is,

In contrast, the complex FIR approach in Fig. 2(c) generates real and imaginary filters dependently to ensure the exact π/2 phase relationship. In order to design the complex FIR filter, a low pass FIR filter (LPF) is designed first and then multiplied by a complex sinusoidal function. The resultant impulse response complex FIR filter can be expressed as

As shown in Eq. (6), the real and imaginary filters have the same envelope, h_LPF(n) and the exact π/2 phase relationship. The analytic signal and its phase from the complex FIR approach are obtained using the appropriate h_LPF(n) and f₀ as

with the corresponding phase

Design of complex FIR filter

The proper selection of h_LPF(n) and f₀ is important to the design of the complex FIR filter. There are two approaches; a wideband complex FIR filter approach [20] and a narrowband complex FIR filter approach [23]. The wideband approach sets the passband of h_LPF(n) to be close to the normalized digital frequency of 0.25 and f₀ as 0.25. In the narrowband approach the passband of h_LPF(n) and f₀ are determined by the bandwidth and the center frequency of an input signal, respectively. Although the narrowband approach has better SNR by rejecting the out of band noise, it requires optimization for a given swept laser and reference interferometer. On the other hand, the wideband approach can cover the total frequency range even though it does not reject the out of band noise. The SNR of the reference interferometer signal tends to be high, hence the improved noise rejection in the narrow band case is not particularly helpful. We have chosen to take the wideband approach for our spectral calibration algorithm since a single set of filter parameters could readily work for multiple systems and is therefore more general.

In order to design the complex FIR filter, the prevalent Park-McClellan algorithm (equiripple design) was used. The specifications such as the passband and stopband frequencies (f_p and f_s), the passband ripple (δ_p), and the stopband attenuation (δ_s) are first determined. From them, a specification for h_LPF(n) is developed by left-shifting the frequency response of the complex FIR by the normalized digital frequency of 0.25 as shown in Fig. 3.

 

Fig. 3 (a): The desired frequency response for complex FIR filter (b): The frequency response for the prototype LPF after left-shifting (a) by 0.25.

Download Full Size | PDF

In the wideband approach, two stopband frequencies are easily defined as 0 and 0.5. Two passband frequencies are set as close to 0 and 0.5 as possible and should have the relationship between them of f_p1 = 0.5-f_p2 in order to make h_LPF(n) a real filter. However, getting closer to 0 and 0.5 requires greater filter length, incurring greater computational load and therefore hardware resources, Therefore, f_p1 and f_p2 are heuristically set to 0.1 and 0.4, respectively. In order to set δ_s, the desired attenuation of the DC and the negative frequency component must be specified. Reilly et al. used the stopband attenuation of 0.17% below the passband [20]. By using this value, δ_s for both the DC and the negative frequency are reduced to 55 dB below the passband as depicted in Fig. 4. Figure 4(a) shows the simplified frequency response of the reference interferometer from Eq. (2) where α is assumed to be 1 (worst case) generating a DC component 6dB higher than the frequency component of interest. Thus, δ_s is set to 0.085% (−61dB) in order to decrease the DC and the negative frequency components to −55dB and −61dB, respectively, as shown in Fig. 4(b). Since the equiripple design method requires the same ripples in both the passband and the stopband, δ_p is set as 0.085%.

 

Fig. 4 (a): The frequency spectrum of reference interferogram (b): the desired frequency spectrum of reference interferogram after complex FIR filtering.

Download Full Size | PDF

Using the stated filter specification and Kaiser’s estimation method [26], the appropriate filter order was found to be 34.

Linearization of OCT signal

After extracting the phase from reference interferometer, the OCT signal is resampled such that δk is constant to enable the use of the fast Fourier transform. We have found that linear interpolation is sufficient for our purposes where we typically image near the zero pathlength difference. It is known that the performance of linear interpolation degrades as the depth approaches the Nyquist frequency, however the SNR is largely insensitive to the interpolation method [27]. The advantage of linear interpolation compared to polynomial interpolation methods is its computational simplicity. Assuming that the DC and auto-correlation terms of the signal are removed, the OCT signal for one reflector can be modeled as,

Intersweep variability in k(n) leads to an apparent sweep to sweep change in the phase of the interferogram, i.e. the argument of the cosine in Eq. (8). Assuming a fixed number of samples are collected after the sweep trigger, this variability can also be thought of as a sweep to sweep change in the center frequency of the laser source. In order to mitigate this issue the interpolant was fixed and encompassed a range of k that was reliably spanned by each sweep of the laser. This approach fixes the beginning, ending, and center wavenumber of every sweep at a predetermined value, thus eliminating the associated phase error.

Results and discussion

Simulation result

The performance of the FIR filter methods was first investigated by MATLAB simulation. For comparison, three FIR methods in Fig. 2 were designed with the same filter order of 34 and the same passband of 0.1~0.4. Figure 5 shows the magnitude responses in the frequencydomain of the three FIR approaches. Clearly, the DC component is not eliminated and the response is not flat in the negative frequency range for the delay element and Hilbert FIR approach, Fig. 5(a). Although the DC component is removed in the bandpass and Hilbert FIR approach, Fig. 5(b), the response is not flat in the negative frequency range. On the other hand, the complex FIR approach not only removes the DC component but also has flat magnitude response in the negative frequency range as depicted in Fig. 5(c). Errors in the estimation of the phase arise both from incomplete suppression of the DC component and residual negative frequency components.

 

Fig. 5 Frequency spetra of three FIR approaches. (a): delay element and Hilbert FIR filter (b): bandpass and Hilbert FIR filter (c): complex FIR filter.

Download Full Size | PDF

Next we compared the full calibration algorithm using the three different FIR approaches. Raw data was acquired from an OCT system with an identical architecture as [6] except the laser source had a 50 kHz sweep rate (SSOCT-1310, Axsun) and the digitizer was a 14-Bit, 250 MS/s Adapter Module (NI5761, National Instrument). The optical layout resembled a Mach-Zehdner interferometer with an optical circulator in the sample arm to preserve sample signal power. For additional details on the hardware setup and signal processing steps see [6]. The raw data consisted of a set of 300 reference interferograms and the corresponding OCT interferograms that comprised an OCT B-scan image of a mouse cochlea. Figure 6 shows the frequency spectra of reference interferograms filtered by delay element and Hilbert FIR, bandpass and Hilbert FIR, and complex FIR approaches as depicted in Fig. 6(c), 6(d) and 6(e), respectively. Each plot contains 300 reference spectral interferograms. The spectra processed with the delay element and Hilbert FIR approach, Fig. 6(c), and the bandpass and Hilbert FIR approach, Fig. 6(d), show a strong DC component and incomplete suppression of negative frequencies, respectively. In contrast, both DC and negative frequency components are eliminated in the complex FIR approach as seen in Fig. 6(e).

 

Fig. 6 (a) Example reference interferogram overlaid with its Hilbert phase. The signal had a center frequency of 19.36 MHz and 0.0775 in digital frequency. (b) Reference interferograms in the frequency domain (c): filtered by delay element and Hilbert FIR approach (d): filtered by bandpass and Hilbert FIR approach (e): filtered by complex FIR approach. These spectra are obtained from 300 reference interferograms.

Download Full Size | PDF

The B-scan images processed via the traditional FFT-IFFT approach and the three FIR approaches described here are shown in Fig. 7. The delay element and Hilbert FIR approach result in Fig. 7(b) shows strong ghost images repeated in depth which result from the high DC component as shown in Fig. 6(c). The bandpass and Hilbert FIR and complex FIR approaches (Fig, 7(c) and 7(d)) show similar results that compared well to the traditional FFT-IFFT Hilbert approach (Fig. 7(a)). That is because the negative component is suppressed enough not to adversely affect the image in bandpass and Hilbert FIR approach. However, in principle the spectrum of the reference signal can be located at any frequency range. If the frequency of the reference interferometer falls in a range where the negative components are not well suppressed, artifacts will arise.

 

Fig. 7 B-scan images of a mouse cochlea. (a): image from FFT-IFFT Hilbert transform (b): image from the delay element and Hilbert FIR approach (c): image from the bandpass and Hilbert FIR approach (d): image from the complex FIR approach. For reference the top of the images are near zero delay with digital frequency 0.024,and the bottom has a zero delay of 0.292.

Download Full Size | PDF

For illustration, the unfiltered reference signal was shifted down by the normalized digital frequency of 0.03, and A-scans recomputed using the down-shifted reference signal. Figure 8(a) shows the down shifted reference signal. Figures 8(b) and 8(c) display the filtered signals by the bandpass and Hilbert FIR, and complex FIR approaches, respectively. Compared to Fig. 6(d) and Fig. 6(e), the negative frequency component is increased by about 20 dB for the bandpass and Hilbert FIR approach while it is essentially unchanged for the complex FIR approach. A-lines from the FFT-IFFT Hilbert Transform, bandpass and Hilbert FIR and complex FIR approaches are shown in Figs. 8(d)-8(f), respectively. In the case of the bandpass and Hilbert FIR approach, strong artifacts arise at around 3.50 mm and 3.96 mm marked with an arrow in Fig. 8(e). These artifacts are not present in the result using the complex FIR approach which closely matches the result using the more traditional FFT-IFFT Hilbert transform.

 

Fig. 8 Effect of spectral down-shifting of the reference signal: frequency spectra of the down-shifted reference signal (a) before filtering, (b) bandpass and Hilbert FIR approach and (c) complex FIR approach. A-scans from (d) FFT-IFFT Hilbert transform, (e) bandpass and Hilbert FIR approach and (f) complex FIR approach. Red arrow indicates an artifact that arises due to incomplete suppression of the negative frequencies in the bandpass and Hilbert FIR approach.

Download Full Size | PDF

In order to assess the different approaches for making vibratory measurements with PhOCT, we utilized a piezo electric element driven at 8 kHz. An M-scan consisting of 10000 A-scans was recorded and the data analyzed using each of the four approaches. The results are shown in Fig. 9.

 

Fig. 9 Frequency responses of displacement of a piezo, from (a) FFT-IFFT Hilbert transform (b) delay element and Hilbert FIR approach (c) bandpass and Hilbert FIR approach (d) complex FIR approach.

Download Full Size | PDF

The phase noise of each approach is characterized by its mean ± standard deviation calculated in the range from 12.5 kHz to 25 kHz. The theoretical shot noise limited phase noise [28, 29] is 25.0 ± 13.1 pm. The result from the traditional FFT-IFFT Hilbert transform, Fig. 9(a), showed a peak at 8 kHz with a magnitude of 221 pm and 29.8 ± 15.4 pm phase noise. The result from the delay element and Hilbert FIR approach has no discernible peak and considerably worse phase noise, 112.5 ± 58.0 pm. The bandpass and Hilbert FIR and complex FIR approaches show similar results with a 221 pm magnitude peak at 8 kHz and 29.5 ± 15.3 pm and 29.6 ± 15.3 pm phase noise, respectively. As with the OCT B-scan image, the strong DC component in the delay element FIR approach severely compromises the results, however incomplete suppression of the negative frequencies in the bandpass filter FIR approach has no discernible impact. Based on the phase noise results both the bandpass and Hilbert FIR approach and complex FIR approach performed slightly better than the traditional FFT-IFFT Hilbert transform. The explanation lies in the difference in the number of data points of transient response between them. The FFT-IFFT Hilbert transform operation can be thought of as a circular convolution in the time domain with a FIR filter having the same length as the input signal [20]. Specifically, since we used 2000 data points per one reference interferogram, there were 1999 transient response data points for FFT-IFFT Hilbert transform while 34 points for bandpass FIR and complex FIR approaches. As a result, the FFT-IFFT Hilbert transform ends up with 1965 more unstable data points than the FIR approaches, which negatively affects the phase extraction of the reference signal. Presumably if we shifted the reference signal as in Fig. 8, the results from bandpass and FIR Hilbert approach would begin to deteriorate.

The bandpass and Hilbert FIR and complex FIR approaches provide comparable performance for both OCT imaging and OCT based vibrometry, at least to the limitations of the assessment presented here. Likewise, they have similar computational complexity. Nevertheless, going forward we have chosen to implement the complex FIR approach because it theoretically is more accurate, conserving the π/2 relationship between the complex components of the analytical signal, and completely suppresses the negative frequencies thereby mitigating potential problems related to the frequency of the reference interferometer signal.

FPGA implementation result

Our motivation for exploring the FIR approach to the Hilbert transform was to reduce the computational load from the traditional FFT-IFFT Hilbert transform. To that end, we have implemented the complex FIR approach on an FPGA and quantitatively compared the resource usage to our current FFT-IFFT method. The code was implemented on a Virtex-5 SX95T FPGA in an NI PXIe-7965R with 18-bit precision. Table 1 shows device utilization for the FPGA implementing post processing using the FFT-IFFT Hilbert and complex FIR approaches. The FFT-IFFT Hilbert approach is based on the structure in Fig. 2(a) where the Hilbert FIR filter is replaced by FFT-IFFT Hilbert transformation. The values in Table 1 include the operations in Fig. 2 as well as linear interpolation and DC subtraction.

Table 1. Comparison of device utilization of FPGA

View Table

Since both sets of code were compiled with the same clock frequency, they output data at the same rate however the latency for the FIR approach is 17 clock cycles while that of the FFT-IFFT approach is 12,000 clock cycles. This difference is further manifest in the resource usage. Since the complex FIR approach requires a smaller number of additions and multiplications than the FFT-IFFT Hilbert approach, slice registers, slice LUTs and DSP48s are reduced by around 20%. In addition, since the FFT-IFFT Hilbert approach requires block RAM in order to match the latency between real and imaginary signals as in Fig. 2(a), the usage of block RAM decreased by 11.5%. The results described here were for a system with a 200 kHz swept laser source.

In order to test the implementation, a piezo was driven at 12 kHz and its displacement was measured using a 200 kHz PhOCT system described in detail in [6]. The displacement was measured 10 times with 4000 A-scans for each M-scan. Figure 10 shows two frequency spectra of the measured displacement: one using the FFT-IFFT Hilbert Transform and the other from the complex FIR approach. The phase noise of each method was calculated as the square root of an averaged noise energy in the frequency domain along the range of 20 kHz to 100 kHz. The result from the traditional FFT-IFFT Hilbert transform displays a peak at 12 kHz with a magnitude of 3.78 ± 0.06 nm and phase noise of 40.2 ± 5.5 pm, where ± indicates the standard deviation among the 10 samples. The complex FIR approach showed similar results with a magnitude of 3.73 ± 0.06 nm and 32.6 ± 2.5 pm phase noise. In agreement with previous results, the complex FIR approach demonstrated slightly better phase noise than the FFT-IFFT Hilbert transform.

 

Fig. 10 Frequency responses of displacement of a piezo, (a) from FFT-IFFT Hilbert transform (b) from complex FIR approach

Download Full Size | PDF

Reducing the number of coefficients used in the FIR filter can reduce the required resources. Owing to the fairly efficient multiply and addition operations on the FPGA, there is a limited potential for further gains over what we have shown in Table 1. Nevertheless, small reductions in the resources can be the difference between an FPGA code successfully compiling and not. In order to explore the limits we have experimented with how few coefficients could be used before the image quality and/or the phase sensitivity were compromised. We repeated the processing of the mouse cochlea and piezo using the complex FIR approach systematically reducing the number of coefficients while retaining the passband and stopband frequencies in order to only affect passband ripple and stopband atteunuation. Figure 11(a) and Fig. 11(b)-11(e) show phase noise and B-scan images at various filter orders, respectively. The phase noise was measured from filter orders ranging from 4 to 40 and quantified by the root mean square (RMS) in order to simultaneously consider mean and standard deviation. In Fig. 11(a), the phase noise converges to around 33 pm at the filter order of 8. The plateau above the order of 8 in Fig. 11(a) comes from the fact that the added noise caused by the complex FIR filter is below the system noise. The B-scan images were processed using the same filter orders as the phase noise measurement. Images using filter orders, 4, 8, 12 and 18 are displayed in Fig. 11(b)-11(e), respectively. In these images there is clearly ghosting at the lower filter orders which essentially completely disappears with filter order of 18. These results taken together leads us to the conclusion that the filter order can be reduced to 18 without any noticeable impact. However, decreasing the filter order reduces the performance of the filter by increasing the passband ripple and the stopband attenuation if the passband and stopband frequencies are unchanged. Hence, careful selection of the minimum filter order is required, depending on the relative strength of the DC component to the frequency components in the reference interferogram.

 

Fig. 11 Phase noise and B-scan image in accordance with filter order. (a) shows phase noise calculated by Root Mean Square (RMS) according to the filter order ranging from 4 to 40. (b), (c), (d), (e) show B-scan images in accordance with the filter order of 4, 8, 12 and 18, respectively.

Download Full Size | PDF

We have also explored the accuracy of the magnitude of the Hilbert transform, |H|, as a measure of the instantaneous swept laser power. We typically use this to calculate a custom window to be applied to the OCT spectral interferogram before computing the IFFT. The custom window is meant to engender the spectral shape of a standard window function to the interferogram. For instance, if we desired the sideband structure of a Hanning window, we would calculate a custom window by W_cust = W_Hann/|H|. We do this to carefully control the sidelobe behavior of our axial point spread function. The phenomena of phase leakage [30] is well documented and known to cause errors in the estimation of depth resolved phase in OCT. Careful control of the sidebands can help mitigate this issue and reduce ambiguity in the interpretation of results of OCT vibratory measurements.

Our standard procedure is to place a mirror reflector in the sample arm of the OCT interferometer, calculate the Hilbert transform of the OCT spectral interferogram and use the measured |H| to calculate the W_cust. Since most swept laser sources have a relatively stable power and spectral shape at least after they have warmed-up we infrequently update this measurement. Therefore, our primary motivation was to determine if the complex FIR approach was accurate enough for this task in order to provide for simplicity and consistency in writing the code to run our system. Nevertheless, we recognized that there could be some advantage to recalculating the window for every sweep of the laser since any fluctuation in the laser power or spectral shape would automatically be compensated. Likewise, it would be a nearly trivial addition to the code already running on the FPGA if we could use the |H| from the reference interferometer instead of a mirror reflector.

In order to test this, the reflection from a mirror was measured using the 200 kHz PhOCT system. 12 M-scans were collected with 100 signal and reference interferograms. Among them, the first set of 100 A-scans was used to calculate dispersion compensation coefficients and a custom window using our standard procedure, as described above. Using this same set of data, coefficients to compensate for differences in chromatic losses between the signal and reference interferograms were calculated.

In the proposed new algorithm, |H| was calculated and then processed by low pass FIR filter (LPF) to remove additive noise before calculating the ideal window. The stopband of LPF was set to the normalized digital frequency of 0.05 because |H| is a baseband signal with narrow bandwidth. We compared the new algorithm with and without correcting for chromatic losses to our current algorithm on the remaining 11 sets of data. Figure 12 shows the results from one data set, in the depth range of the mirror reflection. The results were quantified in terms of Full width at half maximum (FWHM) resolution, SNR and Main to Side lobe Ratio (MSR) for the 1st side lobes on either side. We found no difference for the FWHM resolution or SNR among the three algorithms. Specifically, the mean ± standard deviation of the FWHM resolution over the 11 data sets was 8.56 ± 0.18, 8.56 ± 0.18 and 8.58 ± 0.19 pixels, for the current method, without chromatic correction and with chromatic correction, respectively. The FWHM of the perfect Hanning window, W_Hann, was 8.38 pixels which compares well with the results obtained here with the custom window designed to impart the shape of a Hanning window to the interferogram. Similarly, the SNR was 61.47 ± 1.76, 61.74 ± 1.67 and 61.43 ± 1.73dB, respectively. However, the MSR is degraded significantly if the chromatic correction is not included. From the data set in Fig. 12 the sidelobes increase by 1.85 dB and 3.44 dB. Likewise, it is fairly obvious in Fig. 12(a) that the sidelobe structure varies significantly when the chromatic correction is not applied. In contrast the sidelobe structure in Fig. 12(b) is faithfully reproduced even beyond the 1st sidelobes. These results indicate that the magnitude of the complex FIR Hilbert transform can be used to build a custom window and moreover it can be done using the reference interferometer as long as a chromatic correction is applied. We used a low-pass filter to suppress additive noise so that we could calculate a custom window on every sweep, however we found that using simple averaging of ~100 laser sweeps works equally as well if sweep to sweep variations are not a concern.

 

Fig. 12 Comparison of methods for calculating a custom window. (a): axial point spread functions from current and proposed procedures without chromatic correction (b): axial point spread functions from current and proposed procedures with chromatic correction. For comparison, the axial point spread function with no windowing is also shown.

Download Full Size | PDF

Conclusion

In this article we have evaluated several FIR based methods for spectral calibration of every sweep of a swept laser in real time. We found the complex FIR based approach was the most robust of the three considered. It provided as good or better results compared to the more traditional IFFT-FFT based Hilbert transform on OCT image quality and phase sensitivity. Our implementation on an FPGA also demonstrated that it used fewer hardware resources than the traditional IFFT-FFT approach. Our FPGA code used a filter order of 34 as recommended by Kaiser’s estimation method. We also explored using lower filter orders and found that we could reduce the order to as low as 18 without impacting image quality and phase sensitivity. We also explored using the Hilbert magnitude of the reference interferometer signal as an instantaneous measure of the spectral power of the swept laser. We used this information to calculate an ideal window function that carefully controlled the sidelobes on the axial point spread function. We compared results derived from imaging a mirror with the OCT interferometer and an IFFT-FFT based Hilbert transform to those obtained from the reference interferometer and the complex FIR filter. We found that after a chromatic correction that accounted for differences in the wavelength dependent attenuation between the reference interferometer and the OCT interferometer, we were able to generate essentially equivalent results using both methods. However, the reference interferometer based measurement can be done in real time and is independent of the sample being imaged with the OCT interferometer. In conclusion, the complex FIR approach to the Hilbert transform allows for highly accurate real time wavelength and magnitude calibration of a swept laser source.