EURASIP Journal on Applied Signal Processing 2002:1, 67–72 © 2002 Hindawi Publishing Corporation Approximation of the Wigner Distribution for Dynamical Systems Governed by Differential Equations

A conceptually new approximation method to study the time-frequency properties of dynamical systems characterized by linear ordinary differential equations is presented. We bypass solving the differential equation governing the motion by writing the exact Wigner distribution corresponding to the solution of the differential equation. The resulting equation is a partial differential equation in time and frequency. We then show how it lends itself to effective approximation methods because in the time frequency plane there is a high degree of localization of the signal. Numerical examples are given and compared to exact solutions.


INTRODUCTION
Many dynamical systems are governed by an equation of motion that is an ordinary differential equation with a known driving function. In particular, by a n d n x dt n + a n−1 d n−1 x dt n−1 + · · · + a 1 dx dt where f (t) is the driving force and x(t) the state function.
Often, one wants to study the time-frequency properties of the solution. That would be done by solving (1) and putting the solution into a time-frequency distribution such as the Wigner distribution. However, only in rare cases one is able to solve (1) exactly, and even in those cases f (t) must be necessarily of a simple form (constant, sinusoid, polynomial, etc.). Alternatively, one can attempt to solve the differential equations approximately and substitute the solution into the Wigner distribution. However this is generally problematic because of the many possible regimes and we point out that even in the relatively simple case of the so-called gliding tone problem (to be discussed in Section 3) approximate solutions are quite involved.
Since the introduction of time-frequency methods, it has been realized that signals which may be complicated as a function of time or frequency are often simple in the timefrequency plane. We have developed an approach that takes advantage of this in a direct way. Our procedure is as follows.
In contrast to the standard methods where one solves the differential equation and then uses a time-frequency distribution, for example the Wigner distribution, to ascertain the time-frequency properties of the solution, we show that one can obtain a differential equation for the Wigner distribution of the solution and hence bypass the necessity for solving (1). That is, if the Wigner distribution is defined by we obtain an exact equation of motion for W (t, ω) directly and show that one can approximate the Wigner distribution effectively and directly.
In Section 2, we outline the method for obtaining an exact equation of motion for the Wigner distribution corresponding to the solution of (1). Then we show the exact solution to the gliding tone problem because its solution is crucial to our approximation method. Subsequently, we present our approximation method.

FROM EQUATIONS IN TIME TO EQUATIONS IN TIME-FREQUENCY
We have developed a method that allows us to rewrite the original problem, (1), in the Wigner domain [1,2]. We now describe briefly our solution. We rewrite (1) as a n D n + a n−1 D n−1 + · · · + a 1 D + a 0 x(t) = f (t) (3) and put it in the polynomial notation We have shown that the differential equation for the Wigner and the star sign stands for complex conjugation of the constants a 0 , . . . , a n . The distribution W f, f (t, ω) is the Wigner distribution of the driving force f (t). Using (5) we have been able to solve exactly the gliding tone problem which will be described in Section 3. The exact solution to the gliding tone problem is important not only for its own sake but because the approximation method makes use of the result obtained solving this problem. We point out that the gliding tone problem has a long history.
The possibility of writing a differential equation for the Wigner distribution was initiated by Wigner wherein he wrote the equation of motion for the Wigner distribution corresponding to the solution to the Schrödinger equation. From the very start of the field approximation methods were developed and in fact Wigner himself addressed that issue in his original paper where he calculated the quantum correction to the second virial coefficient of a gas [3,4,5]. However, these methods are based on expansion of the Wigner distribution in powers of Planck's constant, and the recent books by Schleich, Scully, and Zubairy and the references therein describe these methods [6,7]. However, our approximation goals and methods are different because we are not expanding in a small parameter but are trying to approximate the solution of (5) for arbitrary forcing functions.

THE EXACT SOLUTION TO THE GLIDING TONE PROBLEM
The gliding tone problem is the solution to The problem arises in many situations and was first studied by Barber and Ursell [8] and Hok [9]. No exact solution exists in time for this equation, but we have found the exact analytic solution of the Wigner distribution for the variable x(t) [10], that is, W x, x (t, ω). We first rewrite (7) as (8) and factor the differential operator acting on x(t) where The Wigner distribution equation of motion associated with the differential equation is [1] The explicit exact solution is and where u(t) is the step function. Explicit expressions for the underdamped, overdamped, and critically damped cases are given in [10].
In Figure 1, we plot the Wigner distribution for the underdamped case. Note the simplicity of the behavior in the time-frequency plane. The energy peak concentrated about f = 0.25 is due to the interaction of the driving force with the internal resonance of the harmonic oscillator. The response goes to zero as time goes to infinity, and this happens because the transfer function of the system goes to zero when the frequency goes to infinity. The oscillatory terms in the frequency region f = 0 − 0.15 are cross terms generated by the interference of the response for positive frequencies with the symmetric response for negative frequencies. Moreover, note that while the input linear chirp is a complicated function in time, its Wigner representation is a simple delta distribution. This relatively simple behavior of chirp signals in the time-frequency plane is the foundation of our approximation method.

APPROXIMATION METHOD
We now proceed to discuss the approximation method. In this paper, we restrict our attention to driving functions that are one-component. However, the method can be straightforwardly extended to the case of multicomponent signals [11]. We first highlight three important points regarding our approach: • As can be seen from (5) there are no derivatives with respect to frequency and therefore the frequencies in the solution of W x, x (t, ω) are independent of each other. This means that the solution for a certain frequency ω 1 can be evaluated independently of the solution for another frequency ω 2 = ω 1 . 1 • In the Wigner distribution solution of the gliding tone, if one changes the slope β of the input chirp, then also the Wigner distribution of the output changes simply as per (13). Therefore, changing β delays the response at each ω, and the new Wigner distribution can be thought of as a delayed version of the original one.
• The concentration in the time-frequency plane of a linear chirp is a delta function distribution along the instantaneous frequency. While the Wigner distribution of a signal of the form s(t) = e jϕ(t) is not a delta function along the instantaneous frequency, the delta function is nontheless a very good approximation to it. In time-frequency analysis this signal is said to be one-component [12], because its time-frequency spectrum shows an energy concentration (the component) located along the instantaneous frequency ω i (t) = ϕ (t). The case of a multicomponent driving force can be handled by writing it as a sum of delta functions each 1 The same thing happens when one transforms an equation in the class defined by (1) with the Fourier transform. In this case, in fact, the output at frequency ω 1 can be evaluated by multiplying the frequency component of the input at ω 1 by the transfer function of the system at the same frequency, thus been independent to what happens at ω 2 . one centered about the instantaneous frequency of that component. The details will be presented in a future publication.
The key idea of the method is to replace the Wigner distribution of the forcing function with a linear chirp centered at its instantaneous frequency. The independence of the different frequencies and the delay property discussed above assures that we are approximating the Wigner W x, x (t, ω) of the solution of the differential equation. We now formulate the method for obtaining the approximated Wigner distribution of the output of a linear system defined by an ordinary differential equation as in (1) when a one component signal is set as input. We start with a driving function with instantaneous frequency ϕ (t). The steps are: • First compute the exact Wigner W (t, ω) of the solution when a linear quadratic phase signal is the input. As can be seen from (13) the solution can be written as is evaluated by inversion of ϕ (t). This is simply obtained by setting ω = ϕ (t) and then solving for t, getting t = Φ(ω). If the analytic inversion of ϕ (t) is not possible, a polynomial function can be used to approximate Φ(ω).
This method gives the approximated Wigner distribution of x(t) up to a constant. The method will be valid as long as the instantaneous frequency ϕ (t) of the input driving force is slowly varying. When this condition is satisfied the delta function can be considered a good approximation to the Wigner distribution of the driving force f (t). The examples shown in Section 5 will confirm this statement.

EXAMPLES
As an example we take with the instantaneous frequency being quadratic in time and given by According to the method we approximate the Wigner distribution of f (t) by a delta function distribution centered at the instantaneous frequency ω f (t), that we write with respect to time according to the second step of the method presented in Section 4 where inversion yields and we have where A is the normalizing constant and where momentarily we only consider nonnegative values of time. We rewrite the approximated Wigner distribution of the driving force as where τ = t + β/γ − (β/γ) 2 + ω/γ. (This notation is chosen on purpose to be the same as in (13).) Applying the considerations outlined above the approximated output for the harmonic oscillator when the input is a quadratic component can be written, up to a constant and for t ≥ 0, as the Wigner distribution of the gliding tone problem, (13), by substituting τ = t + β/γ − (β/γ) 2 + ω/γ. For negative times the method can be repeated in the same way by substituting

Comparison with exact solution
We now compare the above approximate solution with the exact solution. By exact solution we mean that we first solved the differential equation numerically for x(t) and then numerically calculated the Wigner distribution of the solution. This is an involved procedure but we have done it for numerical comparison and indeed it is this procedure that our method avoids. We have done this for a number of different values of β and γ.
Case 1. β = 3 × 2π , γ = 1.5. In Figure 2, we plot the approximation obtained with the proposed method, and, in Figure 3, we plot the Wigner distribution of the solution computed by numerical integration as just described. The approximation is very good, and this behavior can be observed experimentally to be true in general for increasing values of β.
Case 2. β = 0.6 × 2π , γ = 1.5. In Figure 4, we plot the approximation, and, in Figure 5, we represent the Wigner distribution of the solution computed by numerical integration. The approximation becomes better as t and ω goes to infinity. Also we notice that there are some low energy frequency components before the instantaneous frequency of the input chirp, especially in the frequency interval f = 0.05 − 0.15.
Case 3. β = 0.01 × 2π , γ = 0.4. In Figure 6, we represent the approximation obtained by the proposed method, and, in Figure 7, we represent the Wigner distribution of the solution computed by numerical integration. Again we notice as in the previous experiment, that the quality of the approximation increases as t, ω → +∞. At low values of time and frequency, we notice that anti-causal terms in the Wigner dis- tribution of the solution cannot be neglected. By anti-causal terms we mean the energy located in the time-frequency region n = 20 − 100, f = 0.05 − 0.15. We call these oscillatory terms anti-causal because they arise well before the main interaction between the input driving force and the system, located on the instantaneous frequency of f (t). They represent the well-known cross terms, artifacts generated by the Wigner distribution due to its quadratic formulation (the signal is multiplied by itself).
Based on the numerical experiments presented, we conclude that generally the method works very well for large β, and becomes better as t, ω → ∞. The fact that for small β and for low values of time and frequency the approximation is poor can be explained considering that the Wigner distribution of the quadratic chirp is the Airy function centered at the instantaneous frequency of the chirp [12]. When the slope of the quadratic chirp is low (small β), the cross terms in the Wigner distribution become larger and, hence, they cannot be neglected in the approximation.

CONCLUSION
The approximation method presented takes advantage of the fact that, while solutions to differential equations may be involved and complicated the Wigner distribution of the solution may be relatively simple. In addition, the method takes advantage that in the time-frequency plane monocomponent forcing terms can be effectively approximated. Extension to multicomponent forcing terms are now being investigated. Also, we point out that of particular importance are partial differential equations such as wave equations with driving forces. We have recently presented a method for directly writing the equation for the Wigner distribution corresponding to the solution of a linear partial differential equation [13]. Our aim is to also develop approximation methods for partial differential equations along the same lines as we have developed here for ordinary differential equations.