1. Introduction

Measurement and characterization of ultrashort laser pulses has experienced major breakthroughs in the past decade, enabling retrieval of amplitude envelope and phase of the electric field of a short pulse with methods such as frequency-resolved optical gating (FROG [1]) and spectral phase interferometry for direct electric-field reconstruction (SPIDER [2, 3]). While the former method iteratively reconstructs the pulse shape from measured two-dimensional spectrograms of autocorrelations, SPIDER allows for direct retrieval of the spectral phase from one-dimensional interferograms. This non-iterative retrieval can be accomplished using only Fourier transforms and elementary arithmetic operations, which makes SPIDER the method of choice for video-rate characterization of femtosecond laser pulses [4]. In fact, data acquisition with a posteriori phase processing has been demonstrated at rates ranging up to 1 kHz [5]. For most applications, however, acquisition with concurrent phase retrieval at about 20 Hz is completely sufficient for online monitoring of dispersion compensation, e. g., for aligning amplifier laser chains.

SPIDER is based on acquisition of a spectral interferogram S(ω), arising from the interference between two mutually temporally delayed and spectrally shifted replicas of the pulse under test [2, 3]. The spectral envelope of this interferogram is identical to the spectrum of an individual pulse, with the appearance of additional interference fringes when both pulses enter the spectrograph. Without the spectral shear, the fringe pattern is absolutely equidistant at a period Δω₀ = 2π/Δt ₀. Introducing a spectral shear between the replicas gives rise to deviations from equidistance. Differences between the spectrally varying Δω(ω) and the fringe period Dw0 prior to application of the spectral shear then gives access to the spectral phase of the pulse under test. For determining the spectrally varying modulation period Δω(ω) = 2π/Δt(ω), one typically uses the equation

first suggested by Takeda et al. for applications in holography [6]. Equation (1) processes sampled data analogous to frequency demodulation in radio electronics. This equation can be efficiently computed using two fast Fourier transforms and also effectively allows for noise rejection by suitable choice of the filter time τ _f. In the absence of noise a maximum passband is obtained with τ _f = Δt ₀/2. As soon as experimental data has to be processed, narrowing of the passband filters out noise in the retrieved phase. However, as filtering must not affect the encoded signal, the frequency modulation depth of the fringe pattern imposes a lower limit for τ _f, i.e., it is impossible to reject noise directly in the signal band. One further problem of Eq. (1) is the use of the complex argument function, which is restricted to a range of 2π and requires phase unwrapping prior to further processing. Moreover, the Fourier formalism yields the phase ϕ _mod rather than Δt directly, which requires numerical differentiation and makes the Takeda algorithm susceptible to experimental and numerical noise. As we will demonstrate in the subsequent section, this becomes a particular problem at low signal-to-noise ratios, i.e., in the spectral wings of the interferograms and for signals with strongly modulated spectral amplitude. Eventually, the algorithm may not be able to track the modulation signal anymore, a scenario that will be discussed as carrier collapse in the following. If a carrier collapse occurs in between two decodable segments of the spectrum, this may give rise to arbitrary multiple 2π phase jumps in the retrieved spectral phase. Such a cycle slip of the phase then results in an erroneous kink of the reconstructed Δt(ω).

Recently, the use of wavelets has been suggested as an alternative for phase retrieval in SPIDER [7]. Deng et al. showed in numerical simulations that wavelet based frequency analysis generally outperforms the Takeda algorithm in terms of precision and also demonstrated their method experimentally [8]. Despite its proven advantage, however, wavelet-based phase retrieval is not commonly used for phase retrieval in interferometry. One of the major drawbacks of this method is its relatively low computational efficiency, requiring computation times comparable to iterative pulse retrieval as used for FROG. Even though widely similar acceleration strategies as the bit-reverse shuffling used in fast Fourier transforms have also been suggested for wavelet transformation, their use is restricted to particular choices of wavelets [9] or gives access to only a limited subset of frequency coordinates in dyadic steps [10]. These restrictions render fast wavelet transforms unattractive for frequency analysis applications, which require complex wavelets with real and imaginary parts in quadrature [11].

Here we propose a numerically efficient strategy for application of a fast wavelet transform for phase retrieval. Our strategy is based on the Gabor wavelet [11, 12] that allows for equally good localization in frequency and time, and it is not restricted to particular coordinates in either domain. Our method can compute the transform in fractions of a second, allowing for direct integration of such an algorithm into a video-rate SPIDER apparatus. We will further show which type of wavelet delivers the best performance, and we will demonstrate by experimental example that wavelet-based phase retrieval significantly reduces the susceptibility of the method towards experimental noise, therefore combining the virtues of computational speed and an increased precision.

2. Performance of Takeda phase retrieval for supercontinuum characterization

For a demonstration of the Takeda algorithm, we recorded 60 SPIDER interferograms measured at the output of a single-stage hollow fiber compressor [13, 14, 15]. Using chirped mirrors the pulses of the supercontinuum were compressed to about 9.5 fs pulse duration. For acquisition of the interferograms, an uncooled line-scan CCD camera with 2048 pixels and a dynamic range of 8 bit was used. Characterizing supercontinua with the SPIDER method is a particularly challenging task as the octave-spanning spectra from such a compressor may exhibit spectral drop-outs, i.e., small spectral regions with power densities of 10% of the peak value [16]. In these regions of low fringe contrast, carrier collapse may occur, and the coherence between different spectral sections of the pulse will then be lost. We therefore chose experimental conditions where this problem is imminent, as shown in Fig. 1.

 

Fig. 1. a) Example of a SPIDER interferogram of the hollow-fiber continuum out of a series of 60 measurements. (b) Average signal-to-noise ratio ∑ observed in the series. (c) Standard deviation of phases retrieved by the Takeda algorithm. Cycle slips have been eliminated in the evaluation. (d) Relative occurrence of cycle slips in Takeda-type phase retrieval.

Download Full Size | PDF

Figure 1(a) shows a measured SPIDER interferogram out of a series of 60 individual measurements that are typical for characterization measurements of ultrabroad continua. The spectrum shows a well-defined modulation over most of the central part, however, with some spectral drop-outs in between. For a deeper analysis of this trace, we isolate both, the fringe modulation A _mod and the dc part A _dc of the trace via Fourier filtering, i.e.,

We use these quantities to calculate the signal-to-noise ratio ∑ = A _mod/(A _dc-A _mod) depicted in Fig. 1(b). In this representation, ∑ = 0dB corresponds to a perfect 100% modulation of the fringes, whereas a value of ∑ = -24 dB indicates the detection limit for an 8-bit camera system. The central spectral drop-outs are also clearly discernible in Fig. 1(b). Using Eq. (1) for an analysis of the measured data, we retrieved the phases ϕ _mod from all 60 interferograms and evaluated their standard deviation, as shown in Fig. 1(c). We removed cycle slips prior to calculation of the standard deviation and plotted their occurence vs. frequency separately in Fig. 1(d). For ∑ > -10dB, we do not observe any cycle slips in our retrieved phases. In this regime, however, the spectral drop-outs become clearly visible giving rise to an increased standard deviation of the phase. This is coupling of amplitude noise into phase noise and has been discussed before, see, e.g. [17]. For -15dB < ∑ < -10dB, we observe occasional drop outs on the order of about 1%; for even lower quality of the signal, the Takeda algorithm starts to completely fail, producing 1 or several cycle slips within each modulation cycle. This behavior is accompanied by a dramatic increase of the phase error in Fig. 1(c), i.e., the algorithm is not able to retrieve any meaningful data out of the noisy data below ∑ = -15dB. This may serve as an illustration for the limitations of the Takeda algorithm. In the following, we will discuss a wavelet-based formalism that avoids cycle slips in the first place and is also less susceptible to noise.

3. Fast wavelet transforms for phase retrieval

Avoidance of cycle slip artifacts can be guarenteed using a wavelet-based approach that directly retrieves the modulation period Δω(ω). This is accomplished by choice of the wavelet

with suitable definition of a gating function g(ξ) that allows for localization in the ξ domain. If g̃, i.e. , the Fourier transform of g, is also localized, the parameter δ enables an analysis in the Fourier domain as $\tilde{\Psi }$ (κ) = g̃(κ-δ), with ξ and κ being a Fourier pair. Localization in both domains is ensured when g is a bell-shaped single maximum function. ξ is a dimensionless variable, which, for our application of the analysis of spectral interferograms S(ω), spells out as

where ω enables translation (analysis of different parts of the measured spectrum) and Δω dilation, i.e., adjustment for detection of particular fringe periods Δω in the interferogram. In the following, we will use ω and Δω as the fundamental two coordinates for our wavelet-based analysis of SPIDER interferograms.

Real and imaginary part of the wavelet modulation function in Eq. (6) appear in quadrature, and a straightforward choice for the gating function is a Gaussian, which leads to the Gabor wavelet

with σ = 2log2 and the Heisenberg scaling parameter h, which determines the localization in both domains [8]. Optimization of this parameter will be discussed in Sect. 5. Using the above definition, we define the Gabor wavelet transform (GWT) of a given interferogram S(ω) as

For numerical evaluation of the GWT, input data S_m = S(ω_m) is sampled at N discrete equidistant spectral positions ω_m = ω ₀ + mδω with m = 0…N -1. Here δω is equivalent to the spectral resolution imposed by the used CCD camera and spectrograph, and ω ₀ defines the grating position of the spectrograph. Similarly, we define Δω_n = nδω with n = 1…N/2. In this representation the Gabor wavelet is then evaluated as

and the GWT computes as

From this representation, it is evident that a complete evaluation of the GWT is cubic in N, i.e., requires evaluation of an N-term sum for N possible values of ω_m and N/2 sensible values of Δω_n. Moreover, as each individual calculation requires computation of two real-valued trigonometric functions and one exponential, on the order of 10¹² floating-point operations are necessary at N = 2048 for a full evaluation of the GWT, resulting in several minutes computation time on standard PCs.

On the other hand, however, two optimizations of the GWT are quite obvious for application in spectral interferometry: first, Ψ can only be called with N ² different independent input values as ξ = (j - m)/n in the discrete representation of the wavelet. This allows for usage of a look-up table, i.e., Ψ_mnj is computed for every possible combination of indices m, n and j at initialization of the program and permanently stored in the computer’s random access memory. For example, considering the grid size of N = 2048 used below in our measurements, this requires 64 MByte of memory at most, which is readily available on modern PCs. This optimization step renders calculation of Eq. (9) into the evaluation of a simple scalar product, roughly gaining a factor 100 over repeated calculation of the transcendental functions, which allows for a complete evaluation of the GWT in a few seconds on a modern PC.

As a further optimization step specific for application in spectral interferometry, it is important to see that the modulation period Δω typically stays within a narrow ±50% range around Δω ₀. Moreover, SPIDER typically works best in terms of noise rejection when a rather small modulation period is chosen [13]. With 8-bit data this allows truncation of the wavelet at 1/256 and constriction of the sum in Eq. (9) to a few ten terms. Moreover, the small variation of fringe spacing may be used for alleviation of storage demands for the look-up table, but more importantly, this finding may also serve to render the phase extraction via the GWT a linear algorithm, i.e., an algorithm that requires ∝= N evaluations. This is accomplished with the ridge tracking algorithm discussed below.

Let us assume some prior knowledge on the rough location of the strongest relevant modulation signal W(Δω_s,ω_t), which is easily obtainable from inspection of a single measured interferogram. Initially setting ω_m = ω_t, we scan the Δω_n coordinate and determine the center of gravity

with n ₁ < s < n ₂. Typically, we find it sufficient to include some ten terms in the search for the ridge position, i.e., |n ₁ - n ₂| ≈ 10. It is important to note that computation of the center of gravity is essential for reasonable precision in ridge tracking as simple determination of the maximum value along the Δω coordinate may otherwise repeatedly yield the same integer coordinate n along the ridge and show no variation with ω. The procedure in Eq. (10) is therefore repeated along the ω_m coordinate, retrieving the complete spectral dependence Δω _max(ω_m) from the interferogram under test. For this purpose, however, it is typically not necessary to determine the center of gravity for each possible value of ω_m. We find it sufficient to sparsely track the ridge, scanning the ω_m coordinate in steps of the period Δω _max(ω_m) or even slightly coarser. In the SPIDER method, variations of the spectral phase on scales finer than the spectral shear are meaniningless. Therefore rather coarse scanning can be used. We use the setup described in [13] with a spectral shear Ω_shear = 81 × 10¹² rad/s= 12.9 THz, which is to be compared to a camera pixel spacing corresponding to only 0.15THz. Therefore, even sparsely tracking the ridge at a rate of 10 pixels ≈ 1.5 THz still oversamples the data by a factor 10.

With the described optimizations, phase retrieval along the ω_m coordinate becomes essentially a linear algorithm, requiring computation of two multiplications and two additions per point carried into the sum of Eq. (9), which then has to be evaluated about N times for a complete phase retrieval. Depending on the chosen truncation and sparsity, this allows reduction of the computational time to 10 to 100 ms, which is to be compared with about 1000 s for a complete unoptimized GWT on the same computer.

4. Comparison of the fast GWT with the Takeda algorithm

For a comparison of the fast GWT with the Takeda algorithm, we again use the 60 data sets already discussed above. Other than before, we now apply the concatenation procedure to retrieve the spectral phase of the pulse, see Chapter 3.2 of Ref. [3]. This means that the phases shown in Fig. 2 are spectrally integrated as compared to the modulation phase shown in Fig. 1.

 

Fig. 2. a) Average phases retrieved with the GWT (red curve) and the Takeda algorithm (black curve) and spectral power density (blue curve) of the fundamental. Zero phase was arbitrarily chosen at 390 THz. Error bars indicate the standard deviation in 60 individual measurements. (b) Spectral dependence of the standard deviation for the GWT (red curve) and the Takeda algorithm (black curve). The data from Fig. 1(a) is shown for comparison as the blue curve.

Download Full Size | PDF

We first analyzed the data again with the Takeda algorithm and then redid the analysis with the GWT algorithm discussed above. Average spectral phases reconstructed with either method are depicted in Fig. 2(a). Quite clearly, the retrieved phases agree very well in the center between 370 and 440 THz. Below 370 THz there is already a marked deviation between the retrieved phases, in the extreme wings of the spectrum both methods roll off into opposite directions. Apart from this systematic deviation, we also analyzed statistical fluctuations of both algorithms, see Fig. 2(b). Starting from the fixed phase at pulse center, phase standard deviation grows much faster for the Takeda algorithm, with fluctuations being roughly twice as high as with the GWT. Moreover, fluctuations first tend to grow linearly with distance from spectral center up to the collapse of the method. Again this collapse appears to be delayed for the GWT method, gaining more than 10 THz on the Stokes side of the spectrum and about 5 THz on the other, which corresponds to an increase of 20% in detection bandwidth.

The higher immunity of the GWT against amplitude-noise conversion into phase noise has to be attributed to the better filtering capacities of the GWT. While the Takeda algorithm can only globally reject noise by adjusting τ _f in Eq. (1), the wavelet is localized in ω, which allows for adaptive noise rejection along this coordinate. The GWT therefore allows for partial rejection of noise in the signal band. For the case of an equidistant fringe pattern, one could certainly set the limits in Eq. (1) reasonably tight without suppressing any information, yielding the same noise suppression as with the GWT. However, as soon as information is encoded in the fringe spacing, one has to use less selective filtering with the Takeda algorithm. We therefore expect the relative noise immunity of the GWT to be most pronounced for a SPIDER apparatus with strongly modulated fringe spacings, i.e., in particular for broadband set-ups and those using relatively large amounts of dispersion for generating the reference pulse.

It should be emphasized that the observed 2:1 ratio of the phase standard deviations for experimental data is still quite a dramatic improvement, as real experimental data should also show inherent phase jitter which may easily blur the improvement gained by switching to the GWT. The very pronounced improvement in the experimental data, however, quite clearly confirms the trends predicted by numerical simulation in Ref. [7].

5. Optimum Heisenberg box

As an example, the structure of a complete GWT of a SPIDER interferogram together with its tracked ridge is shown in Fig. 3(c) and (d) for Heisenberg parameters h = 0.7 and h = 3.0, respectively. In this representation, the central drop-outs and their consequences on ridge tracking are also clearly visible, and ultimately the ridge tracking algorithm also collapses in the far wings of the spectrum at above 440 THz and below 335 THz.

The Heisenberg parameter h determines the width of the gating functions g and g̃ in either Fourier domain, i.e., σ _rms and $\tilde{\sigma }$ _rms, respectively. For our choice of unchirped Gaussian gating functions, the product of the two is fixed via

which is symbolically represented by the rectangular Heisenberg boxes in Fig. 3(a).

Now the question on the optimum aspect ratio of the Heisenberg box arises. It may appear intuitive to choose much higher resolution along the Δω axis because the signal is strongly oversampled in ω, and ultimately SPIDER cannot resolve features of the spectral phase that are finer than the spectral shear. To investigate the influence of the aspect ratio of the Heisenberg box, we reconstructed the spectral phase with different h parameters, yielding the results shown in Fig. 3(a). For comparison, the phase error is also shown as a function of the width of the Heisenberg box along the ω coordinate in Fig. 3(b). In both representations, a minimum phase error is observed for particular values of h ≈ 2 and σ _rms ≈ 0.4THz. While the phase error changes only by little in a broad plateau region at 1 < h < 3, it does exhibit a sharp rise at σ _rms < 0.3THz. This behavior is caused by the discretization of the signal and by the determination of the ridge via center of gravity formation. Once the localization in the Δω direction becomes smaller than the discretization along this coordinate, only one significant term enters into computation of the center of gravity [Eq.(10)], and noise enters everywhere else. The other way around, when the signal is too strongly localized in ω, the signal extends too far along the Δω domain, with little difference between the individual elements of the sum and poor rejection of spurious signals at other fringe periods Δω. This makes it quite understandable that there is an optimum value of h ≈ 2, which corresponds to equal localization in both domains at about the discretization in these domains.

 

Fig. 3. a) Dependence of the average phase error on the Heisenberg parameter h. The shape of the Heisenberg box is illustrated for three different values of h. b) same as a function of the width of the Heisenberg box in ω direction. c) Complete GWT shown as logarithmic hues, with the ridge position indicated by a white line for h = 0.7. d) same but with h = 3.0.

Download Full Size | PDF

6. Conclusion

For phase retrieval in spectral interferometry and SPIDER, the Gabor wavelet transform can be utilized as a one-to-one replacement for the Takeda algorithm. Using an optimized algorithm, evaluation times in the 10 to 100 ms range can be readily obtained on standard PCs, making video-rate retrieval possible. Our experiments clearly indicate that ridge tracking via the GWT is significantly more robust towards carrier collapse than the traditional Takeda algorithm. The increased robustness immediately gains some 10% to 20% in bandwidth under otherwise identical conditions. Moreover, the GWT also appears to be less susceptible to conversion of amplitude noise into phase noise under less fatal conditions, as indicated by the low noise levels in the spectral center of our interferogram with its drop-outs.

With all these advantages and execution times only slightly inferior to Fourier filtering techniques, we expect the GWT to replace Takeda-based retrieval for applications in spectral inter-ferometry. Apart from the slightly higher complexity of the algorithm, we cannot see any real disadvantages of the method. In all our experiments, we never encountered a situation where the Takeda algorithm would perform better than the GWT. Finally, we also believe that the very same optimization of the algorithm may find equal application in digital holography and other applications of shearing interferometry [19, 20].