Multiple-order singularity expansion method

The first singularity expansions were developed in the aforementioned SEM [27, 28] based on the observation of a system's response to a sinusoidal input in the transient regime. It was shown to be a combination of damped sinusoidal functions associated with complex frequencies which were the singularities of the transfer function of the system. While Baum first derived the singularity expansion by taking into account the order of these singularities, most applications made use of (and later demonstrated) this expansion in the case of simple poles, i.e isolated singularities of order 1 [34, 35, 37]. The resulting expansion is referred to as the simple pole expansion (SPE), and it has the following expression:

$$\begin{align} H\left(\omega\right) = H_{\mathrm{bgrd}}\left(\omega\right) + \sum_{p} \frac{\mathrm{Res}\left(H,p\right)}{w-p} \end{align} \tag{ 1 }$$

where p denotes the simple poles of H. The SPE is composed of a background term $H_{\mathrm{bgrd}}(\omega)$ which does not possess poles [41], and a sum of resonant terms $\mathrm{Res}(H,p)/(\omega - p)$. This translates the idea that the shape of $H(\omega$) in a specific frequency range is chiefly affected by the nearby singularities, with the background term compensating for the offset introduced by the singularities with respect to the average (or background) value of H in that range. When all the poles p of H are known, the background term reduces to a constant. The contribution of the poles to the shape of H is illustrated in figure 2 for H defined as the reflection coefficient of a thin layer of silver illuminated from one side at normal incidence:

$$\begin{align} H\left(\omega\right) = r\left(\omega\right) - \frac{t\left(\omega\right)t^{\prime}\left(\omega\right)r\left(\omega\right)\mathrm{e}^{2\mathrm{i}\omega \frac{n\left(\omega\right)d}{c}}}{1 - r\left(\omega\right)^2\mathrm{e}^{2\mathrm{i}\omega \frac{n\left(\omega\right)d}{c}}} \end{align} \tag{ 2 }$$

where d = 70 nm is the thickness of the silver layer, c is the speed of light in the air, $n(\omega)$ is the refractive index of silver (see the supporting information (SI) for the detailed expression), $r(\omega) = (n(\omega)-1)/(n(\omega)-1)$ is the Fresnel reflection coefficient at the air/silver interface, and $t(\omega) = 2/(n(\omega)+1)$ and $t^{\prime}(\omega) = 2n(\omega)/(n(\omega)+1)$ are the Fresnel transmission coefficients at the air/silver and silver/air interfaces respectively. The zeros, the poles and the residues of H were determined numerically. Each pole p in the complex frequency plane, in figure 2(b) is associated with a resonant term $\mathrm{Res}(H,p)(\omega - p)$, whose moduli are plotted in figure 2(a). The sum of the resonant terms associated with the poles in the plotted complex frequency window gives the red curve in figure 2(a), which matches the local shape of the exact response (dashed line which corresponds to equation (2), i.e. the approximate position of the local minimum and maximum frequencies. The poor match between this sum of resonant terms and the exact expression is due to the missing background term $H_{\mathrm{bgrd}}(\omega)$ which should correct for the influence of the poles outside of that frequency window.

Zoom In Zoom Out Reset image size

Since most physical systems are described by poles of order 1, the SPE is well suited for the description of the transfer function of such systems, or their response to a temporal or spatial sinusoidal input [35, 37, 38, 48, 49]. Other expressions must be used when the order of the poles is increased. One commonly used alternative in electrical and electric engineering is what we will refer to as the PZF.

The PZF is obtained by making the assumption that the harmonic-domain function H is the ratio of two polynomials of the complex frequency variable ω [1, 2, 24], and that any system responds instantaneously to an input signal:

$$\begin{align} H\left(\omega\right) = \frac{N\left(\mathrm{i}\omega\right)}{D\left(\mathrm{i}\omega\right)} = \frac{N_0\prod_{\ell = 1}^{\mathrm{deg}\left(N\right)}\left(\mathrm{i}\omega - iz_\ell\right)}{D_0\prod_{\ell = 1}^{\mathrm{deg}\left(D\right)}\left(\mathrm{i}\omega - ip_\ell\right)} \end{align} \tag{ 3 }$$

where $iz_\ell$ and $ip_\ell$ are the zeros of N and D respectively, with potentially identical zeros and poles (in which case their order is higher). Let us point out that we consider, in this case, the variable $\mathrm{i}\omega$ instead of ω for N and D in order to remain consistent with the usual definition of the Laplace transform formalism. Using a partial fraction expansion, it is possible to write an expression similar to the SPE:

$$\begin{align} \begin{aligned} H\left(\omega\right) & = \alpha_{1,0} + \frac{\alpha_{1,1}}{\omega - p_1} + \ldots + \frac{\alpha_{1,2}}{\left(\omega - p_1\right)^{\nu_1}} \\ & \quad + \alpha_{2,0} + \frac{\alpha_{2,1}}{\omega - p_2} + \ldots + \frac{\alpha_{2,2}}{\left(\omega - p_2\right)^{\nu_2}}. \\ & \quad + \ldots \\ & = H_{\mathrm{bgrd}} + \sum_{\ell} \sum_{m = 1}^{\nu_\ell} \frac{\alpha_{\ell,m}}{\left(\omega - p_\ell\right)^{m}} \end{aligned} \end{align} \tag{ 4 }$$

where $H_{\mathrm{bgrd}} = \sum_\ell \alpha_{\ell,0}$ and the poles $p_\ell$ are considered with their order or multiplicity $\nu_\ell$. The constants $\alpha_{\ell,1}$ can be identified as the residues, and they can be calculated with the zeros and poles using equation (3) [27, 35]:

$$\begin{align} \begin{aligned} \alpha_{k,1} & = \mathrm{Res}\left(H, p_k\right) \\ & = \lim_{\omega \rightarrow p_k} \left(\omega - p_k\right)H\left(\omega\right) \\ & = -\mathrm{i}\frac{N_0\prod_{\ell = 1}^{\mathrm{deg}\left(N\right)}\left(\mathrm{i}p_k - iz_\ell\right)}{D_0\prod_{\ell \neq k}^{\mathrm{deg}\left(D\right)}\left(\mathrm{i}p_k - ip_\ell\right)} \end{aligned} \end{align} \tag{ 5 }$$

As stated in the introduction, the PZF and the resulting pole expansion are widely used in electronics and system control where they provide highly accurate results when studying systems at time scales larger than the transient time scale. However, this is not the case for other fields such as wave physics, in which case the PZF is missing a complex exponential factor depending on the frequency ω which can be interpreted as the result of a time delay required by a system to respond to different frequencies (i.e. a dispersive system) [35, 37]. This term, which is particularly important for systems which interact with signals with respect to both time and space cannot be obtained from the restrictive hypothesis that H is a ratio of two polynomials. In the next section, we present a singularity expansion taking into account the potentially infinite order of the singularities. From this expression, a more general SZF will then be derived.

We now consider a meromorphic function H, i.e. holomorphic everywhere on $\mathbb{C}$ except for a set of points $\mathcal{P}$ which is the set of its singularities which are all assumed to be poles or isolated essential singularities. Using the Cauchy integration theorem and the residue theorem, the function H can be expressed as the sum of an integral term and an expansion on its set of singularities (see equations (S1)–(S5) in the SI):

$$\begin{align} \begin{aligned} H\left(\omega\right) & = \frac{1}{2\mathrm{i}\pi}\int_{\gamma} \frac{H\left(z\right)}{z-\omega} \mathrm{d}z - \sum_{p} \mathrm{Res}\left(F,p\right) \\ F\left(z\right) & = \frac{H\left(z\right)}{z-\omega} \end{aligned} \end{align} \tag{ 6 }$$

where γ is a closed curve around the singularities $\{p\}$ of H, and $\mathrm{Res}(F, p)$ is the residue of F at the singularity p. The residues of F can be analytically determined from the Laurent series coefficients $\alpha(H, p, n)$ of H around the poles p, which provides local information regarding the behavior of H in the vicinity of its singularities. This leads to the following expression (details of the calculations can be found in equations (S6)–(S17) of the SI):

$$\begin{align} \begin{aligned} H\left(\omega\right) = \frac{1}{2\mathrm{i}\pi}\int_{\gamma} \frac{H\left(z\right)}{z-\omega} \mathrm{d}z + \sum_{p} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(H, p, -m\right)}{\left(\omega-p\right)^m} \end{aligned} \end{align} \tag{ 7 }$$

where ν_p is the potentially infinite order of the singularity p. Equation (7) shows two contributions to the behavior of H at the frequency ω: (i) the nearby singularities $\{p\}$, within the closed curve γ; (ii) the set of values of H on γ, represented by the integral term, which accounts for the contribution of all the singularities outside of γ. The bigger the closed curve gets, the less the integral affects the value of H since the singularities outside the closed curve become too far from ω, provided that H does not grow faster than $|\omega|$ in the complex plane (which is usually the case for physically realistic signals). When this hypothesis holds, we show that the integral term can be replaced by a constant by replacing the curve by a circle of infinite radius (equations (S18)–(S24) of the SI):

$$\begin{align} H\left(\omega\right) = H_{\mathrm{NR}} + H_{\mathrm{R}}\left(\omega\right) \end{align} \tag{ 8 }$$

where $H_{\mathrm{NR}}$ is a constant non-resonant term, and $H_{\mathrm{R}}(\omega)$ is a resonant shaping term which depends on the frequency ω. The non-resonant term $H_{\mathrm{NR}}$ is defined as:

$$\begin{align} H_{\mathrm{NR}} = H\left(a\right) - \sum_{p} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(H, p, -m\right)}{\left(a-p\right)^m} \end{align} \tag{ 9 }$$

where a is an arbitrary complex frequency which is not a pole. Let us stress that the choice of a has no influence over the value of the constant term $H_{\mathrm{NR}}$. The resonant term $H_{\mathrm{R}}(\omega)$ possesses the following expression:

$$\begin{align} H_{\mathrm{R}}\left(\omega\right) = \sum_{p} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(H, p, -m\right)}{\left(w-p\right)^m}. \end{align} \tag{ 10 }$$

If 0 is not a pole, we usually set a = 0 and the non-resonant term $H_{\mathrm{NR}}$ can be expressed using the static response H(0):

$$\begin{align} H_{\mathrm{NR}} = H\left(0\right) - \sum_{p} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(H, p, -m\right)}{\left(-p\right)^m}. \end{align} \tag{ 11 }$$

We refer to equations (8) and (9) as the multiple-order SEM (MOSEM). The accuracy of MOSEM is shown in the case of poles of order 2 in a purely theoretical example in the SI. When the order ν_p of all the poles is 1, the expression becomes:

$$\begin{align} H\left(\omega\right) = H_{\mathrm{NR}} + \sum_{p} \frac{\alpha\left(H,p,-1\right)}{\omega - p}. \end{align} \tag{ 12 }$$

Let us stress that the residue of H associated with the pole p is defined as $\mathrm{Res}(H,p) = \alpha(H,p,-1)$. Therefore by identifying $H_{\mathrm{NR}}$ as a constant background term $H_{\mathrm{bgrd}}$ when all the singularities are taken into account in the resonant terms, we obtain the SPE expression defined in equation (1):

$$\begin{align} H\left(\omega\right) = H_{\mathrm{bgrd}} + \sum_{p} \frac{\mathrm{Res}\left(H,p\right)}{\omega - p}. \end{align} \tag{ 13 }$$

If H is reconstructed using only a finite set of singularities, two strategies can be adopted to still obtain a good match between the resulting truncated SPE or MOSEM and the exact expression; (i) as mentioned before, $H_{\mathrm{bgrd}}$ can be considered as a holomorphic function, i.e. with no singularities, which models the influence of the singularities missing from the resonant term $H_{\mathrm{R}}(\omega)$; (ii) the influence of the poles outside of the region of interest can be accounted for by considering an additional virtual singularity within the resonant term as was done in [38].

The Weierstrass and Hadamard factorization theorems state that any meromorphic function can be written as the ratio of two complex polynomials multiplied by a complex exponential. Starting from this consideration, expressions such as the PZF (in which the complex exponential is usually removed) can be obtained (although not straightforwardly) to study the phase and amplitude variations at real frequencies. The poles and zeros can be defined in the resulting rational fractions as the zeros of the denominator and numerator respectively. Alternatively, a factorization involving the singularities and zeros can be obtained from MOSEM, providing a clearer link between the expanded and factorized forms as was done in [15, 35] in the case of poles of order 1.

Here, we aim at deriving a generalized factorized expression from the MOSEM expression of the function H in equations (8) and (9). H can always be written as:

$$\begin{align} H\left(\omega\right) = \omega^m G\left(\omega\right) \end{align} \tag{ 14 }$$

with $m\unicode{x2A7E}0$ the order of the zero $\omega=0$ of H, and G a meromorphic function which does not possess 0 as a zero. If m = 0, we have G = H. Let us consider F the log-derivative of G, with Gʹ the complex derivative of G:

$$\begin{align} \begin{aligned} &F = \frac{G^{\prime}}{G} \\ &G^{\prime} = \partial_\omega\left[G\right]. \end{aligned} \end{align} \tag{ 15 }$$

Gʹ is calculated by taking the derivative of MOSEM applied to G (equations (8) and (9)). The poles of F are the zeros $z_\ell$ and the singularities $p_\ell$ of G, and they are all of order 1 (see equations (S33)–(S36) in the SI):

$$\begin{align} \forall z_\ell \in \mathbb{C}, G\left(z_\ell\right) & = 0 \iff F\left(z\right) = \infty \nonumber\\ \forall p_\ell \in \mathbb{C}, G\left(p_\ell\right) & = \infty \iff F\left(z\right) = \infty. \end{align} \tag{ 16 }$$

In addition, their associated residues are the order ν_z and ν_p of the zeros and poles of G respectively:

$$\begin{align} \forall \left(z_\ell, \nu_z\right), ~ \mathrm{Res}\left(F, z\right) & = +\nu_z \nonumber\\ \forall \left(p_\ell, \nu_p\right), ~ \mathrm{Res}\left(F, p\right) & = -\nu_p. \end{align} \tag{ 17 }$$

We apply MOSEM to F taking these values into account:

$$\begin{align} F\left(\omega\right) : & = F_{\mathrm{NR}} + F_{\mathrm{R}}\left(\omega\right) \nonumber\\ F_{\mathrm{NR}} : & = F\left(a\right) - \sum_{p} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(F, p, -m\right)}{\left(a-p\right)^m} \nonumber\\ F_{\mathrm{R}}\left(\omega\right) : & = \sum_{p} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(F, p, -m\right)}{\left(w-p\right)^m} \end{align} \tag{ 18 }$$

with $p\in\{z_\ell, p_\ell\}$, $\nu_p = 1$, and $\alpha(F, p, -1) = 1$ if $p = z_\ell$ is a zero of H and $\alpha(F, p, -1) = -1$ if $p = p_\ell$ is a pole of H:

$$\begin{align} F\left(\omega\right) & = F_{\mathrm{NR}} + F_{\mathrm{R}}\left(\omega\right) \nonumber\\ F_{\mathrm{NR}} & = F\left(a\right) + \sum_{p_\ell} \frac{\nu_{p_\ell}}{a-p_\ell} - \sum_{z_\ell} \frac{\nu_{z_\ell}}{a-z_\ell} \nonumber\\ F_{\mathrm{R}}\left(\omega\right) & = -\sum_{p_\ell} \frac{\nu_{p_\ell}}{\omega-p_\ell} + \sum_{z_\ell} \frac{\nu_{z_\ell}}{\omega-z_\ell}. \end{align} \tag{ 19 }$$

Finally, we obtain the following expression of F:

$$\begin{align} F\left(\omega\right) & = \frac{G^{\prime}\left(a\right)}{G\left(a\right)} + \sum_{z_\ell} \left( \frac{\nu_{z_\ell}}{\omega-z_\ell} - \frac{\nu_{z_\ell}}{a-z_\ell} \right) - \sum_{p_\ell} \left( \frac{\nu_{p_\ell}}{\omega-p_\ell} - \frac{\nu_{p_\ell}}{a-p_\ell} \right). \end{align} \tag{ 20 }$$

By integrating F on a curve from the arbitrary complex frequency a to the frequency of interest ω, and applying the exponential function to the result, we derive the following expression (equations (S37)–(S41) in the SI):

$$\begin{align} H\left(\omega\right) & = \omega^m ~ G\left(a\right) ~ \frac{\prod_{z_\ell}\left(1 - \frac{\omega-a}{z_\ell-a}\right)^{\nu_{z_\ell}}}{\prod_{p_\ell}\left(1 - \frac{\omega-a}{p_\ell-a}\right)^{\nu_{p_\ell}}} ~ \mathrm{e}^{\mathrm{i}\tau\left(\omega-a\right)} \nonumber\\ \tau & = -\mathrm{i}\left(\frac{G^{\prime}\left(a\right)}{G\left(a\right)} - \sum_{z_\ell} \frac{\nu_{z_\ell}}{a-z_\ell} + \sum_{p_\ell} \frac{\nu_{p_\ell}}{a-p_\ell} \right) \end{align} \tag{ 21 }$$

where G(a) is defined as:

$$\begin{align} G\left(a\right) & = \frac{H\left(a\right)}{a^m} ~ \qquad \qquad \textrm{if } a \neq 0, \nonumber\\ G\left(a\right) & = \frac{\partial_\omega^m\left[H\right]\left(0\right)}{m!} ~\qquad \textrm{otherwise}. \end{align} \tag{ 22 }$$

and $G^{\prime}(a)$ is defined as:

$$\begin{align} G^{\prime}\left(a\right) & = \frac{\partial_\omega\left[H\right]\left(a\right)}{a^m}-m\frac{H\left(a\right)}{a^{m+1}} \quad \textrm{if}~ a\neq0,\nonumber\\ G^{\prime}\left(a\right) & = \frac{\partial_\omega^{\left(m+1\right)}\left[H\right]\left(0\right)}{\left(m+1\right)!} \qquad \qquad \textrm{otherwise}. \end{align} \tag{ 23 }$$

We refer to equation (21) as the SZF. Let us point out the presence of the aforementioned complex exponential $\mathrm{e}^{-\mathrm{i}\tau\omega}$ missing from the PZF, but also the presence of the known response G(a) (the static response in the case of a = 0). Let us also stress that in the case of an LTIS, the phase introduced by the complex exponential can be set to an arbitrary position by changing the time origin or the space-origin (for a system with coupled space and time variables). Shifting the time origin by τ₀ in the temporal domain results in a multiplication by $\mathrm{e}^{\mathrm{i}\omega\tau_0}$ of the input $X(\omega)$ in the harmonic domain. This is tantamount to considering the transfer function $H(\omega)\mathrm{e}^{\mathrm{i}\omega\tau_0}$ with the same input $X(\omega)$ as depicted in figure 3 for H defined as:

$$\begin{align} H\left(\omega\right) & = 5 + \frac{2 + 0.1\mathrm{i}}{\omega - \left(2 + 0.3\mathrm{i}\right)} + \frac{4 + 0.1\mathrm{i}}{\omega - \left(5 + 1.4\mathrm{i}\right)} + \frac{6 + 0.1\mathrm{i}}{\omega - \left(6.5 + 0.7\mathrm{i}\right)}. \end{align} \tag{ 24 }$$

In figure 3(a), the phase of the transfer function $H(\omega)$ is shown in the complex plane. In figure 3(b), it is plotted for $H(\omega)$ multiplied by $\mathrm{e}^{\mathrm{i}\omega\tau_0}$ with τ₀ arbitrarily set to $\pi/6$. We show that the position of the zeros and singularities remains the same after switching from $H(\omega)$ to $H(\omega)\mathrm{e}^{\mathrm{i}\omega\tau_0}$, but a non-constant phase-shift is introduced in the complete complex frequency plane by the phasor $\mathrm{e}^{\mathrm{i}\omega\tau_0}$. It is thus possible to set τ₀ in such a way that $\tau + \tau_0 = 0$, where τ is the time constant naturally appearing in the SZF in equation (21).

Zoom In Zoom Out Reset image size

By setting τ = 0, a = 0, and defining two constants N₀ and D₀ appropriately, we recover the PZF described in the previous section by considering only N_s poles and N_z zeros:

$$\begin{align} H\left(\omega\right) & = \frac{N_0\prod_{\ell = 1}^{N_z}\left(\mathrm{i}\omega - iz_\ell\right)}{D_0\prod_{\ell = 1}^{N_s}\left(\mathrm{i}\omega - i p_\ell\right)} \nonumber\\ N_0 & = \left[\prod_{\ell = 1}^{N_z}\frac{i}{z_\ell}\right] ~ G\left(0\right)\nonumber\\ D_0 & = \prod_{\ell = 1}^{N_s}\frac{i}{p_\ell}. \end{align} \tag{ 25 }$$

The MOSEM and SZF expressions, which are reminded in figure 4, are two equivalent means to characterize a function: (i) MOSEM expression, which relies on the behavior in the vicinity of the singularities only and thus depend on the singularities and Laurent series coefficients associated with them, is useful to get a fast and accurate approximation of a function in the harmonic-domain; (ii) the SZF is more convenient to look at the phase and better understand the behavior on the real frequency axis. It is easy to obtain the Laurent series coefficients, and thus the MOSEM expression, from the SZF using the definition of the Laurent series coefficients:

$$\begin{align} \alpha\left(H,p_0,-m\right) & = \lim_{\omega \rightarrow p_0} \frac{1}{m!}\partial_\omega^{m-1} \left[\eta\right]\left(\omega\right) \nonumber\\ \eta\left(\omega\right) &:= \left(\omega-p_0\right)^{\nu_{p_0}}H\left(\omega\right) \end{align} \tag{ 26 }$$

where p₀ is a pole of order $\nu_{p_0}$ of H, and $H(\omega)$ is given by the SZF in equation (21). It is however more difficult to obtain the zeros from the MOSEM expression, although they can be approximated in a specific frequency range by writing the MOSEM expression as a rational function involving only the singularities in the selected complex frequency window and solving for the zeros of the numerator.

Zoom In Zoom Out Reset image size

Let us consider a real-valued function h(t) in the temporal domain, associated with a complex function $H(\omega)$. Since h is real, H possesses a hermitian symmetry in the complex plane [35, 50]:

$$\begin{align} H\left(-\omega^*\right)^* = H\left(\omega\right). \end{align} \tag{ 27 }$$

This property leads to constraints on the distribution of the poles and singularities in the complex frequency plane. By evaluating the complex conjugate of the singularity expansion in equation (8) evaluated at $-\omega^*$, it can be shown that for any pole p of order ν_p , $-p^*$ is also a pole of order ν_p .

We can also obtain the Laurent coefficients associated with $-p^*$ to those of p via:

$$\begin{align} \forall m\in\mathbb{Z},~\alpha\left(H, -p^*, -m\right) = \left(-1\right)^m\alpha\left(H, p, -m\right)^* \end{align} \tag{ 28 }$$

and we show similarly that the non-resonant term $H_{\mathrm{NR}}$ must be real. The Laurent series coefficients of H at $-p^*$ are the opposite of the complex conjugate of those at p. This leads to the following MOSEM expression:

$$\begin{align} H\left(\omega\right) & = H_{\mathrm{NR}} + \sum_{p, ~ p\in \mathrm{i}\mathbb{R}} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(H, p, -m\right)}{\left(w-p\right)^m} + \sum_{p, ~ \mathrm{Re}\left[p\right]\gt 0} \sum_{m = 1}^{\nu_p} \left[\frac{\alpha\left(H, p, -m\right)}{\left(w-p\right)^m} + \left(-1\right)^m \frac{\alpha\left(H, p, -m\right)^*}{\left(w+p^*\right)^m} \right]. \end{align} \tag{ 29 }$$

Let us point out that equation (28) restricts $\alpha(H, p, -m)$ to $\mathrm{i}\mathbb{R}$ if $p\in \mathrm{i}\mathbb{R}$. The poles and Laurent series coefficients thus always come in pairs, as shown in figure 5 for the transfer function H defined in equation (2). In (a) and (b), the log-amplitude in a complex frequency window as well as its symmetric window is plotted, highlighting the symmetry of the amplitude and thus of the singularities. In (c) and (d), the same frequency windows were chosen for the phase plots. They show the antisymmetry of the phase relative to the imaginary axis in the complex frequency plane, which is tantamount to an antisymmetry of the Laurent series coefficients. Let us stress that the hermitian symmetry only arises from the fact that we have real-valued signals in the temporal domain. Considering a causal plane wave $\mathrm{e}^{-\mathrm{i}\omega_0 t}u(t)$, with u the Heaviside step function and $\omega_0 \gt 0$ is not equivalent to considering that time flows backward. It only means that the variations of the phase are opposed to those of $\mathrm{e}^{\mathrm{i}\omega_0 t}u(t)$.

Zoom In Zoom Out Reset image size

If we now look at the complex conjugate of the SZF (equation (21)) evaluated at $-\omega^*$, we can show that if z is a zero of order ν_z of H, then $-z^*$ is also a zero of order ν_z . In addition, let us show that the Hermitian symmetry forces the time constant τ to be real-valued. If τ were complex-valued, it could be written $\tau = \tau_\mathrm{R} + \mathrm{i}\tau_\mathrm{I}$, with $\tau_\mathrm{R}$ and $\tau_\mathrm{I}$ the real and imaginary parts of τ respectively. Since the opposite of the complex conjugate of zeros and poles are also zeros and poles of the same order, the Hermitian symmetry would thus lead to:

$$\begin{align} \omega^m~\mathrm{e}^{\mathrm{i}\tau_\mathrm{R}\omega}~\mathrm{e}^{-\tau_\mathrm{I}\omega} = \left(-1\right)^m~\omega^m~\mathrm{e}^{\mathrm{i}\tau_\mathrm{R}\omega}~\mathrm{e}^{\tau_\mathrm{I}\omega}. \end{align} \tag{ 30 }$$

This condition would not be satisfied as ω tends towards $+\infty$ unless $\tau_I = 0$, therefore, $\tau \in \mathbb{R}$. If we set a = 0, the SZF thus yields the following expression of τ:

$$\begin{align} \tau = &2\mathrm{Im}\left(\sum_{z_\ell, \mathrm{Re}\left[z_\ell\right]\gt 0} \frac{\nu_{z_\ell}}{z_\ell} - \sum_{p_\ell, \mathrm{Re}\left[p_\ell\right]\gt 0} \frac{\nu_{p_\ell}}{p_\ell} \right) -\mathrm{i} \left(\frac{G^{\prime}\left(0\right)}{G\left(0\right)} - \sum_{p_\ell\in \mathrm{i}\mathbb{R}} \frac{\nu_{p_\ell}}{p_\ell} \right)\end{align} \tag{ 31 }$$

with $G^{\prime}(0)/G(0) \in \mathrm{i}\mathbb{R}$. The time constant τ is thus the sum of two contributions: (1) the phase-shift introduced by the imaginary part of the singularities and the non-null zeros, (2) a constant term depending on the static response of the derivatives of H. In physical systems for which the space and time variables are coupled, the phase shift can be modified by moving the spatial or temporal origin. In this case, the zeros and poles are identical, and only the constant term linked to the static response is changed.

Stability and causality are linked but distinct principles which can both be expressed in terms of the position of the singularities in the complex plane depending on the convention used to perform a Fourier transform [50]. Causality states that any signal h must be generated at a certain time t_h , and that it cannot depend on its future values. If h is a causal signal, it can therefore be written, using the Heaviside step function u, as:

$$\begin{align} h\left(t\right) = h\left(t\right)u\left(t-t_h\right). \end{align} \tag{ 32 }$$

As long as h does not diverge faster than an exponential function for a long time t, it can thus always be regularized using a function h_γ , γ > 0, which converges:

$$\begin{align} \begin{aligned} &h\left(t\right) = \mathrm{e}^{\gamma t} h\left(t\right) \mathrm{e}^{-\gamma t} u\left(t\right) \\ &h\left(t\right) = \mathrm{e}^{\gamma t} h_\gamma\left(t\right). \end{aligned} \end{align} \tag{ 33 }$$

Let us point out that we set $t_h = 0$ in u, and we can do so without losing in generality. By construction, h_γ possesses a Fourier transform H_γ from which the harmonic domain function H associated with h can be determined:

$$\begin{align} H\left(\omega\right) = H_\gamma\left(\omega - \mathrm{i}\gamma\right). \end{align} \tag{ 34 }$$

h can be retrieved by integrating H (multiplied by a complex exponential) over a horizontal line within the region of convergence of H, lower-bounded by the amplitude of the smallest diverging exponential function diverging faster than h (see figure 6). Using the residue theorem, it can thus be shown that any pole possessing a positive imaginary part in the complex frequency plane is associated with a causal diverging exponential function in the temporal domain, i.e. an unstable signal. Therefore, the poles of the harmonic-domain function associated with any stable signal must have a negative imaginary part. The stability of the signals can be interpreted with energy considerations. As an input signal interacts with a system, it exchanges energy with it. This leads to a modification of the input signal which results in the output signal. For a passive system, the energy is transferred from the input to the system. The output thus has a lower energy and cannot diverge if the input is stable. For the output to diverge or become unstable, it must result from a sufficient energy transfer from the system to the input, or the interaction with an already diverging input with the system. Therefore, the only way to obtain an unstable output is through a high-energy, unstable input, or an active system. In terms of singularities, this means that the singularities of the transfer function of a passive system always possess a negative imaginary part (using our Fourier transform convention). If we inject energy into the system, i.e. the system is active, we move the singularities closer and closer to the real axis until the system is unstable and at least one singularity possesses a positive imaginary part.

Zoom In Zoom Out Reset image size

In many problems, the main focus is the global temporal response of the LTIS at a specific frequency. In these situations, the characterization of the system in the harmonic domain is sufficient since it is almost equivalent to a study of the temporal permanent-regime. We expect the response of a stable system to an excitation to be composed of that same input, scaled and phase-shifted in the permanent-regime, as well as exponentially decaying functions in the transient-regime [38]. It is worth demonstrating this known result by deriving the response y(t) to a sinusoidal input $x(t) = \mathrm{Re[e}^{-\mathrm{i}\omega_0t}u(t - t_\mathrm{e})]$ of an LTIS of impulse response h(t). Since y and x are real, h must be real, which leads to:

$$\begin{align} \begin{aligned} y\left(t\right) & = \left(h \ast \mathrm{Re}\left[x\right]\right)\left(t\right) = \mathrm{Re}\left[\left(h \ast x\right)\left(t\right)\right] \\ & = \mathrm{Re}\left[\mathcal{L}^{-1}\left[H.X\right]\right] = \mathrm{Re}\left[\mathcal{L}^{-1}\left[Y\right]\right]. \\ \end{aligned} \end{align} \tag{ 38 }$$

The response to a sinusoidal input is the real part of the response to a causal plane wave $t \mapsto \mathrm{e}^{-\mathrm{i}\omega_0t}u(t)$. In the harmonic domain, Y and X can be written as:

$$\begin{align} \begin{aligned} &Y\left(\omega\right) = X\left(\omega\right)H\left(\omega\right) \\ &X\left(\omega\right) = \frac{\mathrm{i}}{\omega - \omega_0}. \end{aligned} \end{align} \tag{ 39 }$$

Therefore, the set of isolated singularities of the harmonic response includes the singularities of the transfer function with the same order, as well as ω₀ of order 1 introduced by the harmonic input signal X. Applying MOSEM to Y gives:

$$\begin{align} Y\left(\omega\right) = Y_{\mathrm{NR}} + Y_{\mathrm{R}}\left(\omega\right) \end{align} \tag{ 40 }$$

with

$$\begin{align} Y_{\mathrm{NR}} & = Y\left(0\right) - \sum_{p} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(Y, p, -m\right)}{\left(-p\right)^m} + \frac{\alpha\left(Y,\omega_0,-1\right)}{\omega_0}, \nonumber\\ Y_{\mathrm{R}}\left(\omega\right) & = \sum_{p} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(Y, p, -m\right)}{\left(w-p\right)^m} + \frac{\alpha\left(Y,\omega_0,-1\right)}{\omega - \omega_0}. \end{align} \tag{ 41 }$$

The calculation of the Laurent series coefficients of Y leads to:

$$\begin{align} \alpha\left(Y, p, n\right) = -\sum_{\ell = -\nu_p}^{n} \frac{\mathrm{i}\alpha\left(H,p,\ell\right)}{\left(\omega_0 - p\right)^{n-\ell+1}} \end{align} \tag{ 42 }$$

for the singularity p of order ν_p . In addition, the residue of Y at ω₀ is proportional to the transfer function evaluated at the plane wave frequency ω₀:

$$\begin{align} \begin{aligned} &\alpha\left(Y, \omega_0, -1\right) = \lim_{\omega \rightarrow \omega_0} \left(\omega - \omega_0\right)X\left(\omega\right)H\left(\omega\right), \\ &\alpha\left(Y, \omega_0, -1\right) = \mathrm{i}H\left(\omega_0\right). \end{aligned} \end{align} \tag{ 43 }$$

By replacing the corresponding terms in equations (40) and (41), we obtain:

$$\begin{align} \begin{aligned} Y\left(\omega\right) & = Y_{\mathrm{NR}} + \mathrm{i}\frac{H\left(\omega_0\right)}{\omega - \omega_0} + \sum_{p} \sum_{m = 1}^{\nu_p} \frac{\alpha\left(Y, p, -m\right)}{\left(w-p\right)^m}. \end{aligned} \end{align} \tag{ 44 }$$

It follows that the inverse Laplace transform of Y can be written as the sum of a Dirac's delta function multiplied by $Y_{\mathrm{NR}}$, an expansion on the dynamic states similarly to the resonant term in equation (37) after replacing $\alpha(H, p, n)$ by $\alpha(Y, p, n)$, and the scaled and phase-shifted sinusoidal function $t \mapsto \mathrm{e}^{-\mathrm{i}\omega_0t}$:

$$\begin{align} \mathcal{L}^{-1}\left[Y\right]\left(t\right) & = \mathrm{Re} \left[Y_{\mathrm{NR}}\delta\left(t\right) + H\left(\omega_0\right)\mathrm{e}^{-\mathrm{i}\omega_0t} u\left(t\right) +\sum_{p} \sum_{m = 1}^{\nu_p} \alpha\left(Y, p, -m\right) f_{m, p}\left(t\right) \right]\end{align} \tag{ 45 }$$

with

$$\begin{align} f_{m,p}\left(t\right) := \frac{\left(-\mathrm{i}\right)^m}{\left(m-1\right)!}t^{m-1}\mathrm{e}^{-\mathrm{i}pt}~u\left(t\right). \end{align} \tag{ 46 }$$

The initial value theorem imposes that $Y_{\mathrm{NR}}$ must be null if x is a continuous input. Since $\lim_{t \rightarrow 0^+} y(t) = 0$, i.e. the system cannot immediately respond to a physical signal, the initial value theorem applied to y gives:

$$\begin{align} 0 = \lim_{t \rightarrow 0^+} y\left(t\right) = \lim_{\omega \rightarrow +\infty} \mathrm{i}\omega Y\left(\omega\right). \end{align} \tag{ 47 }$$

It follows that the constant, non-resonant term must be equal to 0, and that the sum of the residues of Y is null:

$$\begin{align} \begin{aligned} &\left( \omega Y_{\mathrm{NR}} \rightarrow 0 \right) \Rightarrow \left(Y_{\mathrm{NR}} = 0\right) \\ &0 = \alpha\left(Y,\omega_0,-1\right) + \sum_p \alpha\left(Y, p, -1\right). \end{aligned} \end{align} \tag{ 48 }$$

We obtain the aforementioned expected result: the response y of a stable system to a causal sinusoidal input is only composed of the scaled and phase-shifted input in the permanent regime. In the transient regime, y must be expanded on the set of dynamic states $f_{m,p}$ of the system:

$$\begin{align} \mathcal{L}^{-1}\left[Y\right]\left(t\right) & = \left[ H\left(\omega_0\right)\mathrm{e}^{-\mathrm{i}\omega_0t} u\left(t\right) + \sum_{p} \sum_{m = 1}^{\nu_p} \alpha\left(Y, p, -m\right) f_{m, p}\left(t\right). \right]\end{align} \tag{ 49 }$$

This expansion remains valid in the case of an unstable system, although the concept of transient and permanent regime would no longer hold. In this case, the dynamic states would hold information regarding the diverging speed of y.