Optimization and selection of time-frequency transforms for wave-packet analysis in ultrafast spectroscopy

The analysis of quantum beats in time-resolved spectroscopic signals is becoming a task of primary importance because it is now clear that they bring crucial information about chemical reactivity, transport, and relaxation processes. Here we describe how to exploit the wide family of time-frequency transform methodologies to obtain information not only about the frequency but also about the dynamics of the oscillating components contributing to the overall beating signal. Several linear and bilinear transforms have been considered, and a general and easy procedure to judge in a non-arbitrary way the performances of different transforms has been outlined. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
In the last decades, spectroscopists have learned how to manipulate laser pulses with increasingly shorter time durations and how to use them in advanced photonics applications to capture transient physical phenomena with an increasing degree of subtleness and detail. A typical feature of the spectroscopic response in the ultrafast regime is the presence of oscillations of the signal amplitude as a function of delay times between exciting pulses. Such oscillations, typically referred to as quantum beats, are associated with the time evolution of coherent superpositions of states prepared by the ultrafast pulses, i.e. a wave-packet. The presence of quantum beats in ultrafast signals is a well-known textbook phenomenon, typically associated with coherent nuclear or electronic motions [1][2][3]. Quantum beats may appear in ultrafast spectroscopy when two levels are excited coherently and emit to a final common level; in these conditions, the detected spectrum oscillates with a frequency corresponding to the energy gap between the two levels.
The analysis of the quantum beats in different time regimes has proved to be particularly informative to extract information on various aspects of the energetics and the dynamics of the systems. For example, quantum beats in nanosecond spectroscopy have been exploited for the determination of molecular structure parameters, such as spin-orbital coupling constants, interaction constants in the Zeeman effect, dipole moments in excited states, rotation and anharmonic constants [4,5] and singlet fission mechanisms [6]. Moving to pico-and femtosecond resolution allows monitoring the dynamics of chemical reactions, transition states, energy transfer and vibrational dynamics in real time [7][8][9]. Making a further step towards the attosecond domain, information about the rapid motion of electron density inside a molecule, sub-cycle evolution of coherent electronic wavepacket and photoionization can be obtained [10][11][12][13].
Whatever the nature of the studied material or the particular technique and time resolution used, the analysis of the frequency and time behavior of the beatings is a key ingredient to extract information about the system dynamics. To this aim, methodologies based on linear and bilinear time-frequency transforms (TFT) have been already proposed and applied to ultrafast 1D [14][15][16][17][18][19] and 2D spectroscopy measurements [20][21][22][23][24][25][26]. The advantage of TFTs is that they analyze a non-stationary oscillating signal providing both time and frequency resolution, spreading out the information content along the two dimensions of the time-frequency plane [27]. TFT methods are model-independent. Thus, while other parametric methods independently developed to estimate wave-packet dynamics require a careful a priori definition of a suitable fitting function and initial parameters [28][29][30], TFT methods can be directly applied without any preliminary assumption. Therefore, they are ideal to achieve a quick and immediate visualization of the frequency and time behavior of the relevant modes contributing to the beating pattern and, for example, identifying families of long-or short-leaved components [31].
One should anyway take into account that TFTs, like Fourier transforms, are subjected to the time-frequency uncertainty principle. This implies that the time and frequency resolution are strongly correlated, and it is not possible to reach simultaneously ideal resolution along both dimensions. In order to maintain the joint information across the time-frequency plane, one has thus necessarily to compromise and, among the variety of strategies available in the time-frequency analysis, select from time to time the TFT and the parameters most suitable for the specific input signal to be analyzed and for the specific expected output. In other words, the TFT approach is necessarily application-dependent. In addition, these methods can be particularly prone to artifacts, both along time and frequency dimensions. In this context, in [32] we recently overviewed the most known TFT methods and qualitatively compared their performances when applied to typical signals obtained by coherent femtosecond spectroscopy. The comparison is particularly challenging because of the arbitrary choice of TFTs parameters, which greatly influence the outcome of the transforms. To make this approach more general and ease the exploitation of this powerful methodology, also in different time regimes, we developed an easy and general procedure to quantitatively optimize and evaluate the performances of various TFTs when applied to wave-packet signals. The procedure is exemplified for the analysis of typical femtosecond coherent responses.

Generalities about time-frequency transforms
In general, the time-domain representation s(t) and the frequency-domain representation S(ν) of a signal are related through the Fourier transform (FT) and inverse Fourier transform (IFT): In standard Fourier analysis, the square of the signal in time (|s(t)| 2 ) and in frequency domain (|S(ν)| 2 ), are called instantaneous energy and power spectrum, respectively. The power spectrum can also be expressed as the FT of the signal correlation function R(t) [33]: where τ is the time lag between the signal and its complex conjugate. The critical point in a conventional FT analysis is that in s(t) and S(ν) the time and frequency variables are treated as mutually exclusive: one can interpret the frequency representation as averaged over the entire time domain and vice versa. The representations in terms of time-frequency distributions (TFTs) overcome this limitation since they are joint functions of time and frequency, able to represent the energy or intensity per unit time per unit frequency. Table 1 summarizes the TFTs considered in this work.
To understand the basic properties of a TFT, it is useful to consider the best-known example of this class, the short-time Fourier transform (STFT) [34]. It examines the frequency content of the signal as a window h(t) is moved in time, providing the Fourier transform only of the portion of the signal contained in the window. This operation differs from the conventional FT only by the presence of the window, typically described as a Gaussian function centered at the time t. The addition of a window function allows retaining the time information performing a localization of the FT. A similar approach is adopted in the continuous wavelet transform (CWT) [35,36]. Here, instead than a fixed window h(t), a short time oscillating function φ(ν, t) called wavelet is used. The wavelet φ has a central frequency ν 0 and it is stretched and translated along time, leading to a variable frequency resolution. Often the squared modulus of STFT and CWT are employed. They are called spectrogram (SP=|STFT| 2 ) and scalogram (SCA=|CWT| 2 ), respectively. SP and SCA are used, rather than their linear analogs, to be compared with bilinear transforms, to maintain the quadratic dependence on the signal. The Margenau-Hill spectrogram (MHS) can be built similarly, starting from a product of two STFTs with independently defined windows.
Most bilinear TFTs belong to a more general class of transforms, the so-called Cohen's class, and they can be specified as [37]: This general form is reminiscent of Eq. (2) with the addition of a kernel function F(t , τ), whose definition characterizes different transforms, as summarized in Table 1. These bilinear TFTs are generally able to distribute in a more sophisticated way the energy of the signal over the time-frequency plane, than linear TFTs. The definition of the kernel function is the crucial point in the interpretation of TFT operations and performances, as more deeply discussed in Section 4. When F(t , τ) = 1, the Wigner-Ville transform (WV) is retrieved. WV is the basic function acting as a precursor for all the other functions of the Cohen's class family. WV achieves the best performances in terms of uncertainty principle and has the highest resolution among all the time-frequency distributions of the class. However, it presents strong cross-term interference as a result of the structure of Eq. (3) that basically compares the signal with a translated replica. For this reason it is very rarely employed for practical purposes. The issue of interference terms attenuation is a central topic in bilinear TFTs and results particularly critical in the context of ultrafast spectroscopy beating analysis. The cross-term interference is manifested as oscillating features arising in the time-frequency plane halfway between each pair of components [32]. One of the possible approaches to efficiently attenuate interference is to define a suitable kernel function F(t , τ). For example, in the smoothed-pseudo-Wigner-Ville distribution (SPWV), two time windows g and h are introduced and translated along the time dimension, introducing time and frequency smoothing, respectively. As a result, interference contributions are strongly mitigated. These two windows can be arbitrarily and independently tuned in order to obtain the best performances for each particular type of signal [32]. Another of the best-known transforms of the Cohen's class is the Choi-Williams (CW) distribution [38]. To suppress cross-term interferences, CW adopts an exponential kernel where the width parameter α (Table 1) tunes the amount of interference filtering. This strategy is very effective for signals where the different components are well separated in time and frequency, but strong interference cross-terms cannot be avoided for synchronized and frequency overlapping oscillations. For the sake of completeness, in addition to the most famous transforms previously mentioned, two other widely used time-frequency representations have been considered: the Butterworth distribution (BUD) [39] which has an exponential kernel similar to CW, and the Zhao-Athlas-Marks distribution (ZAM) that has a cone-shaped kernel [40]. In analogy with the definition of the SPWV, two additional window functions h and g have been included in the definition of CW, BUD and ZAM in Table 1. This is a standard procedure to mitigate the known limitations of exponential and cone-shaped kernels. All the considered transforms, summarized in Table 1, have been chosen based on their spread and frequency of utilization and on the availability of the associated algorithm in the most used data analysis software platforms. Our calculations are performed using MATLAB and the Time-Frequency Toolbox [41]. Table 1. Linear and bilinear transforms considered in this work.

Name
Acronym Definition

Optimization and selection procedure
A proper time-frequency signal processing procedure requires (i) the selection of the best TFT and (ii) the optimization of its parameters, based on the properties of the particular experimental signal to analyze. Given the differences in the definitions and in the parameters on which they depend, a meaningful comparison of different TFTs requires a preliminary step of optimization of the associated parameters. A different choice of parameters and windows (or wavelets) leads indeed to significantly different results, hindering the comparison of the outcomes. To select the best TFT in a non-arbitrary way the following procedure has been set up: (i) generation of a synthetic signal with known frequency and time properties. Working with a signal with known information content will help in the evaluation of the performances of different transforms as well as in the identification of possible artifacts arising during the analysis; (ii) definition of a theoretical TFT based on the synthetic signal and complying with basic properties of time-frequency distributions; (iii) optimization of the TFTs based on the best matching with the theoretical transform; (iv) application of the TFTs with optimized windows parameters and comparison of their results. For simplicity, following the standard procedure in the literature, Gaussian windows are considered in all the transforms and their width σ is the parameter to be optimized. Similarly, for |CWT| 2 a complex Morlet wavelet is used, defined as a product of a complex exponential and a Gaussian window [36]. The choice of a Gaussian shape is motivated by its peculiar properties. Indeed, it guarantees the optimal joint time-frequency resolution from the uncertainty point of view. Moreover, the balance between the time and the frequency resolution can be controlled by the width parameter σ. The performances of other window functions were already explored in a previous publication [32]. In any case, the optimization procedure described here can be easily generalized to other window functions.

Generation of a synthetic signal
First of all, in order to emulate the typical oscillating features appearing in time-resolved data, a discrete wave-packet signal has been considered: where the n−th component y n (t) of the wave-packet is characterized by the scalar amplitude A n , the damping time τ n and the oscillation frequency ν n . The signal is assumed to exist only for positive times and the associated FT is denoted with the capital letter, i.e. Y (ν) and Y n (ν). The choice of modeling the signal as a collection of exponentially damped oscillatory components is general and fully justified by the typical dephasing dynamics of the wave-packets observed in any time-resolved experiment [1,2,7]. The specific values of frequencies and time parameters can then be selected based on the particular frequency and time regime of the experimental signal we want to reproduce. As an example, here we have chosen a set of parameters (Table 2) and an investigated time window and sampling (from 0 to 1000 fs with a step of 5 fs, as in Fig. 1), that mimic the typical conditions encountered in a femtosecond experiment.  The amplitudes in Eq. (4) are defined as A n = 1/ √ τ n in order to guarantee that the total energy of each component is the same, i.e., each component has a peak with the same area in the power spectrum. Indeed, it is possible to analytically demonstrate that for each component y n (t): where This definition of amplitudes does not diminish the generality of the approach, but it represents an elegant choice to guarantee that in the ensuing optimization procedure all the wave-packet components have similar weight.

Definition of a theoretical transform
A theoretical bilinear time-frequency representation THEO (ν, t) is introduced. This function is specifically tailored for exponentially damped signals. The most important property imposed in the definition of this distribution is the fulfillment of the marginal conditions [27]: If marginals are fulfilled, then, we are in the ideal condition that (i) the sum of the energy distribution for all frequencies at a particular time gives the instantaneous energy, and (ii) the sum over all times at a particular frequency would give the energy density spectrum; with energy distribution and energy density spectrum defined as in Section 2.
The function THEO is defined as: note that this equation incorporates the analytical expressions reported in Eq. (6). The marginal properties are rigorously respected for a signal with N = 1, and they are still satisfactorily valid if the components of the signal do not strongly overlap in the time-frequency plane, as in our case. For the n−th component, it is possible to demonstrate that which implies that the total energy of the signal is conserved when distributed on the timefrequency plane through the function THEO. The idea behind the definition of this theoretical time-frequency distribution is straightforward; a Lorentzian shape is used to describe the frequency profile of the features and an exponential decay is used for the time profile. A factor (A 2 n τ n /2) −1 is introduced to guarantee the conservation of the energy. This theoretical TFT thus represents the hypothetical limit for any time-frequency transform: it is characterized by the highest possible frequency resolution of the Fourier transform and the exact time evolution of the originating signal, without the presence of any cross-term interference. Being defined as the hypothetical limit of a TFT, THEO perfectly incorporates all the characteristics of the synthetic signal without any artifact, but it clearly violates the time-frequency indetermination principle.

Optimization of TFTs
In order to find the optimal performance conditions, an unconstrained minimization problem over the set of parameters p = (p 1 , p 2 , ..., p M ) that characterize each investigated TFT can be solved: where M is the number of parameters. The minimum of the problem represents the TFT realization that gives the smallest difference with respect to the theoretical transform and thus delivers an optimal compromise between interference attenuation and time-frequency concentration of features. Consequently, the solution of the minimization problem in Eq. (10) provides the best parameters that optimize the different TFTs on the basis of a common performance criterion, allowing a meaningful comparison. As an example, in the case of SPWV, the expression inserted in the minimization problem in Eq. (10) is where p 1 is a scale factor to account for the different normalization conditions of the timefrequency representations, and p 2 = σ h and p 3 = σ g are the width parameters of the Gaussian window functions h(t) and g(t), respectively, as defined in Section 2 and Table 1. The procedure has been repeated for all the bilinear TFTs reported in Table 1. The optimized parameters for the various TFTs are summarized in Table 3. It is useful to clarify that these optimized parameters strongly depend on the specific signal we considered, thus they should be refined according to the experimental signals requirements. Moreover, it is important to stress that the minimization procedure outputs the parameters that simultaneously optimize the auto-terms for the best match with the theoretical representation, and optimally reduce interference fringes, accordingly to the TFT capability. That is because THEO does not contain any interference by definition, thus every interference contribution in the TFT increases the norm in the minimization problem.

Results
The output of the minimization problem for all the considered TFTs are reported in Fig. 2. The histogram in Fig. 2(a) clearly shows that SPWV is the TFT that guarantees the best match with the theoretical one. In the optimized SPWV the cross-term interference is mostly eliminated maintaining at the same time a high resolution along both time and frequency axes. This result can be rationalized giving a closer look to the optimized kernel functions F(t , τ) introduced in Eq. (3). In particular, the mathematical operations behind the Cohen's class distributions and the effect of kernel filtering can be studied exploiting the narrow-band ambiguity function (AM), which is directly connected to the base function of the Cohen's class (WV):   The representation of the various transforms in the (η, τ) space allows reinterpreting how the kernel function smooths the interference artifacts, as already quickly outlined in Section 2. In the (η, τ) space, indeed, the time and frequency smoothing operations can be performed simply multiplying the AM by the kernel function. In this picture, the optimization problem in Eq. (10) is acting by shaping the bidimensional kernel function, within the constraints of its parametrization, in order to remove the interference contributions as much as possible. Some examples of optimized kernel functions are reported in Fig. 3(e).
Both SP and SPVW use a simple Gaussian smoothing kernel. The difference is that in SPVW the bidimensional Gaussian function has independent widths along τ and η, which leads to a reliable and stable interference terms attenuation. In SP, instead, the widths along the two dimensions are controlled by a single parameter and therefore they cannot be optimized independently. Exponential kernels like in CW (and BUD) are in principle more efficient and can guarantee a higher time-frequency concentration of the features. Nevertheless, they do not attenuate interference terms efficiently in the specific case of wave-packet signals. In fact, the kernel has a shape elongated along the τ = 0 line, where interference of signals superimposed in time appears. Similar considerations can be done about the performance of cone-shaped kernels, like in ZAM. We can conclude that the SPWV distribution, among all the considered TFTs, is the most suited transform in the analysis of a sum of exponentially damped signals with frequency and damping times typically encountered in femtosecond optical experiments. It is the most stable and gives the most predictable interference attenuation, reducing the probability of feature deformations and distortions in noisy experimental signals. [32]  The qualitative prediction reported in [32] was thus confirmed and rationalized, also after the numerical optimization step. Note that, while this specific result applies to femtosecond spectroscopy responses, the proposed procedure is general and it can be used for wave-packet oscillating signals in any time domain, from attoseconds to nanoseconds and beyond, simply varying the expression and/or the parameters of the synthetic signal used for the definition of the THEO transform. The approach can be easily furtherer generalized considering different shapes for the window functions, and accounting for signals with different properties (for example including components overlapped in time or frequency and characterized by various levels of noise or chirp).
Whatever the properties of the beating signal to be analyzed, this procedure allows to identify the best TFT and optimize its parameters to achieve the best analysis of frequency and time behavior of the oscillating signal.

Conclusions
In the latest years, time-resolved experiments in different time-domains are producing a growing amount of evidence, proving the crucial role of wave-packet dynamics in chemical reactivity, transport and relaxation processes, manifested as multicomponent and non-stationary oscillating signals. Owing to the complexity of such oscillations, it is often challenging to extract information about their origin and the associated dynamics using only conventional spectral analysis. To overcome this difficulty, we described here how to exploit analysis methods borrowed from the signal-processing field, known as time-frequency transforms, to obtain information not only about the frequency but also about the dynamics of the oscillating components contributing to the overall beating signal. Several linear and bilinear transforms have been considered and a general and easy procedure to judge in a non-arbitrary way the performances of different transforms has been outlined. This numerical procedure relies on the definition of a synthetic signal with properties able to reproduce the most important features of the experimental signal to be analyzed. This synthetic signal is then used to define a theoretical transform retaining ideal resolution along both time and frequency axes. The minimization of the difference between this ideal theoretical transform and a generic TFT allows ascertaining the best TFT and optimizing its parameters to obtain the best analysis. As an example, the application of this methodology to synthetic signals mimicking the results of a typical femtosecond time-resolved experiment lead to the identification of the bilinear smoothed-pseudo Wigner-Ville (SPWV) transform as particularly advantageous for the analysis of signals in this time and frequency domain. The different performances of the considered TFTs have been justified and interpreted in the framework of the narrow-band ambiguity function, providing a clear picture of how different transforms, and the associated kernel functions, operate to minimize interference artifacts. This is an important information for the generalization of the proposed methodology to different time-resolved beating signals.