Approximation and analysis of transient responses of a reverberation chamber by pulsed excitation

Reverberation chambers show transient behaviour when excited with a pulsed signal. The field intensities can in this case be significantly higher than in steady state, which implies that a transient field can exceed predefined limits and render test results uncertain. Effects of excessive field intensities of short duration may get covered and not be observable in a statistical analysis of the field characteristics. In order to ensure that the signal reaches steady state, the duration of the pulse used to excite the chamber needs to be longer than the time constant of the chamber. Initial computations have shown that the pulse width should be about twice as long as the time constant of the chamber to ensure that steady state is reached. The signal is sampled in the time domain with a sampling frequency according to the Nyquist theorem. The bandwidth of the input signal is determined using spectral analysis. For a fixed stirrer position, the reverberation chamber, wires, connectors, and antennas can jointly be considered as a linear time-invariant system. In this article, a procedure will be presented to extract characteristic signal properties such as rise-time, transient overshoot and the mean value in steady state from the system response. The signal properties are determined by first computing the envelope of the sampled data using a Hilbert transform. Subsequent noise reduction is achieved applying a Savitzky–Golay filter. The point where steady state is reached is then computed from the slope of the envelope by utilising a cumulative histogram. The spectral analysis is not suitable to examine the transient behaviour and determine the time constants of the system. These constants are computed applying the method of Prony, which is based on the estimation of a number of parameters in a sum of exponential functions.

An alternative to the Prony Method is the Time-Domain Vector-Fit method. In contrast to the first mentioned variant, it is now also possible to determine the transfer function of the overall RC system. Differences and advantages of the methods will be discussed.

1 Introduction

Reverberation Chambers are commonly used for immunity and emission tests, where an important class of immunity tests is performed with pulsed modulated sinusoidal signals. A reverberation chamber (RC) is characterised by its quality factor Q, which is dependent on the operating frequency f and the time constant τ of the chamber (Krauthäuser, 2007; Pfennig, 2015). The quality factor affects the time required until steady state is reached when the excitation is switched on. This characteristic places a limit on the minimum duration of the pulse with which the chamber can be excited if steady state is desired. Figure 1 shows a typical pulsed signal s with a sinusoidal carrier of frequency f_c and magnitude A.

Figure 1 Pulsed signal.

The duration of the pulse is denoted t_H and the time between two pulses t_L. The period T of the signal is, as a result, given by t_H+t_L. A typical value for the pulse duration is t_H≥2.5τ, i.e. a multiple of the time constant (Deutsche Kommission Elektrotechnik Elektronik Informationstechnik im DIN und VDE, 2011). Transient behaviour is observable in the response of the RC if a switching event takes place. The field intensities can in this case be significantly higher than in steady state and possibly exceed predefined limits in an immunity test which consequently can damage the device under test (DUT).

A large number of publications have already dealt with transient behaviour in   RCs (Artz, 2017; Manicke and Krauthäuser, 2012; Krauthäuser and Nitsch, 2002; Amador et al., 2010).

The possible overshoot can best be observed by measurements in time-domain. Artz has shown in his work that strong overshoot is accompanied by a low amplitude in the steady state and vice versa.

In contrast to the work of Artz (2017) it is not necessary to determine an ensemble mean value in advance of the evaluation. Each time series is evaluated individually.

The focus of this work is on signal analysis which depends on the stirrer position. First, in Sects. 2 and 3 a signal theoretical description of the overall RC system is given. The measurement setup and signal preprocessing is presented in Sect. 4.

An algorithm that extracts signal characteristics such as the amplitude of the overshoot when switching on and off as well as the amplitude at steady state is presented in Sect. 5.

The parameters thus obtained are required for further signal analysis.

The chamber response is recorded with an oscilloscope for each stirrer position and is therefore time-discrete. This time series of measurements is interpolated to a function that is continuous over time. The function depends on the location within the test volume, the time and the stirrer position.

In Sect. 6 the Method of Prony is used for this purpose. The transient behaviour of the RC exhibits exponential behaviour. By means of Prony's Method, which represents a sum of weighted exponential functions, this transient and decay behaviour can be represented intrinsically and is thus a well suited choice. For large data sets, the use of Prony's Method is not practical. This is due to the point wise nature of the method, where each data point is considered. A high time resolution and a long measurement duration result in large data sets, which leads to the problems mentioned above.

An alternative is the Time-Domain Vector-Fit method (Grivet-Talocia, 2004). Based on initial frequencies, a rational function is adapted to correspond to the measurement process in the time domain. An advantage over the Prony Method is the extraction of the transfer function of the entire system consisting of RC, antennas, cables and connectors. This procedure is presented in Sect. 7.

2 Reverberation chamber

A typical RC which can be regarded as rectangular cavity of width a, height b, and length c, that is equipped with a mechanical stirrer, a transmitting and a receiving antenna is shown in Fig. 2.

Figure 2 Reverberation chamber with stirrer and antennas.

Resonant cavities possess a linear time-invariant system and can therefore be considered as a linear time-invariant (LTI) system (Artz and Hirsch, 2014), as displayed in Fig. 3.

Figure 3 Reverberation chamber as linear time-invariant system.

The system response g of the chamber,

$$\begin{array}{}\text{(1)}& g\left(t\right)=s\left(t\right)\star h\left(t\right),\end{array}$$

is given by the convolution of the input signal s and the system impulse response h of the RC, where g, s and h are functions of time t. h can be determined by measurement (Artz and Hirsch, 2014). The system impulse response includes the impulse response of the receiving and transmitting antenna, cables and the feedthrough. LTI systems transform a time-harmonic input signal of a certain frequency to a time-harmonic output signal of the same frequency (Hoffmann, 1998). The resonant frequencies are determined by Hill (2009)

$$\begin{array}{}\text{(2)}& f_{mnp}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathit{\mu }\mathit{\epsilon }}}}\sqrt{{\left({\displaystyle \frac{m}{a}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{n}{b}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{p}{c}}\right)}^{\mathrm{2}}},\end{array}$$

where μ denotes the permeability and ε the permittivity of medium in the RC. The indices m, n, p determine the field mode. A TE- and a TM-mode exists for every combination of m, n, and p unless one of the indices is zero. In this case, only a TE- or a TM-mode exists. Any equipment in the chamber changes the boundary conditions for the field and consequently affects the resonant frequencies so that modes are not known exactly.

To alleviate this problem an approximation of the mode density

$$\begin{array}{}\text{(3)}& D_{\mathrm{s}}\left(f\right)=\mathrm{8}\mathit{\pi }abc{\displaystyle \frac{f^{\mathrm{2}}}{c_{\mathrm{0}}^{\mathrm{3}}}}-{\displaystyle \frac{a+b+c}{c_{\mathrm{0}}}}\end{array}$$

can be used as a performance indicator for the RC, where c₀ denotes the speed of light (Liu et al., 1983). D_s(f) provides information on the number of resonant frequencies in the neighbourhood of f. The mode density increases if the frequency increases. The chamber used for the measurements has a width a of 3.7 m, a height b of 3.0 m, a length c of 5.3 m, which amounts to a volume V of 58.83 m³. For a frequency of 400 MHz for example, the mode density D_s(f) in the given chamber equals 8.74 modes per MHz.

A key performance indicator of a RC is the quality factor Q, see Krauthäuser (2007).

A theoretical value for Q is given by

$$\begin{array}{}\text{(4)}& Q={\displaystyle \frac{\mathrm{3}{\mathit{\mu }}_{\mathrm{0}}V}{\mathrm{2}{\mathit{\mu }}_{\mathrm{m}}S{\mathit{\delta }}_{\mathrm{S}}}}\end{array}$$

with the chamber volume V, the skin depth δ_S and the walls permeability μ_m. The skin depth of the wall material can be given as

$$\begin{array}{}\text{(5)}& {\mathit{\delta }}_{\mathrm{S}}=\sqrt{{\displaystyle \frac{\mathrm{2}}{\mathit{\omega }{\mathit{\mu }}_{\mathrm{m}}\mathit{\sigma }}}},\end{array}$$

with the wall conductivity σ.

A RC can electrically be viewed as a network of RLC circuits, see Bernard Demoulin (2011). Each excited mode can be modelled with an RLC circuit.

3 Signal theory

The signal displayed in Fig. 1 can be written as

$$\begin{array}{}\text{(6)}& s\left(t\right)=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}\left(A\cos \left({\mathit{\omega }}_{\mathrm{c}}(t-kT)+\mathit{\phi }\right)\mathrm{rect}\left({\displaystyle \frac{t-kT}{t_{\mathrm{H}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\right),\end{array}$$

where ω_c=2πf_c and φ represents the phase shift of the carrier at t=0 s. The rectangular function vanishes for t<0 s and t>t_H, equals a half for t=0 s and t=t_H, and equals one for $\mathrm{0}<t<t_{\mathrm{H}}$. In order to generate a pulsed signal s, the single pulse is repeated periodically. The spectrum of the signal is determined by computing the Fourier coefficients

$$\begin{array}{}\text{(7)}& S_n={\displaystyle \frac{\mathrm{1}}{T}}\underset{T}{\int }s\left(t\right)e^{-jn{\mathit{\omega }}_{\mathrm{0}}t}\mathrm{d}t.\end{array}$$

The identity

$$\begin{array}{}\text{(8)}& \cos \left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(e^{jx}+e^{-jx}\right)\end{array}$$

enables to rewrite s as

$$\begin{array}{}\text{(9)}& s\left(t\right)={\displaystyle \frac{A}{\mathrm{2}}}\left({\stackrel{\mathrm{‾}}{s}}_{\mathrm{1}}\left(t\right)e^{j\mathit{\phi }}+{\stackrel{\mathrm{‾}}{s}}_{\mathrm{2}}\left(t\right)e^{-j\mathit{\phi }}\right)\end{array}$$

with

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{s}}_{\mathrm{1}}\left(t\right)=\mathrm{rect}\left({\displaystyle \frac{t}{t_{\mathrm{H}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)e^{j{\mathit{\omega }}_{\mathrm{c}}t}\text{(11)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{s}}_{\mathrm{2}}\left(t\right)=\mathrm{rect}\left({\displaystyle \frac{t}{t_{\mathrm{H}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)e^{-j{\mathit{\omega }}_{\mathrm{c}}t}.\end{array}$$

The frequency domain representation of ${\stackrel{\mathrm{‾}}{s}}_{\mathrm{1}}$ is given by

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{S}}_{\mathrm{1}n}& ={\displaystyle \frac{t_{\mathrm{H}}}{T}}\underset{\mathrm{0}}{\overset{\frac{T}{t_{\mathrm{H}}}}{\int }}{\stackrel{\mathrm{‾}}{s}}_{\mathrm{1}}\left(t^{\prime }\right)e^{-jn{\mathit{\omega }}_{\mathrm{0}}^{\prime }{t}^{\prime }}\mathrm{d}{t}^{\prime }\\ & ={\displaystyle \frac{t_{\mathrm{H}}}{T}}\mathrm{si}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(n{\mathit{\omega }}_{\mathrm{0}}^{\prime }-{\mathit{\omega }}_{\mathrm{c}}^{\prime }\right)\right)e^{-j\frac{\mathrm{1}}{\mathrm{2}}\left(n{\mathit{\omega }}_{\mathrm{0}}^{\prime }-{\mathit{\omega }}_{\mathrm{c}}^{\prime }\right)}\end{array}\end{array}$$

and similarly

$$\begin{array}{}\text{(13)}& {\stackrel{\mathrm{‾}}{S}}_{\mathrm{2}n}={\displaystyle \frac{t_{\mathrm{H}}}{T}}\mathrm{si}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(n{\mathit{\omega }}_{\mathrm{0}}^{\prime }+{\mathit{\omega }}_{\mathrm{c}}^{\prime }\right)\right)e^{-j\frac{\mathrm{1}}{\mathrm{2}}\left(n{\mathit{\omega }}_{\mathrm{0}}^{\prime }+{\mathit{\omega }}_{\mathrm{c}}^{\prime }\right)},\end{array}$$

where the primed symbols refer to normalised variables, i.e.

$$\begin{array}{}\text{(14)}& \begin{array}{c}{\mathit{\omega }}_{\mathrm{0}}^{\prime }={\mathit{\omega }}_{\mathrm{0}}{t}_{\mathrm{H}},\\ {\mathit{\omega }}_c^{\prime }={\mathit{\omega }}_{\mathrm{c}}{t}_{\mathrm{H}},\end{array},t^{\prime }={\displaystyle \frac{t}{t_{\mathrm{H}}}}\mathrm{and}\mathrm{d}t=t_{\mathrm{H}}\mathrm{d}{t}^{\prime }.\end{array}$$

The coefficients S_n can finally be written as

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}{S}_n& ={\displaystyle \frac{A}{\mathrm{2}}}\left({\stackrel{\mathrm{‾}}{S}}_{\mathrm{1}n}{e}^{j\mathit{\phi }}+{\stackrel{\mathrm{‾}}{S}}_{\mathrm{2}n}{e}^{-j\mathit{\phi }}\right)\\ & ={\displaystyle \frac{A}{\mathrm{2}}}{\displaystyle \frac{t_{\mathrm{H}}}{T}}\left[\mathrm{si}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(n{\mathit{\omega }}_{\mathrm{0}}^{\prime }-{\mathit{\omega }}_{\mathrm{c}}^{\prime }\right)\right)e^{-j\left(n{\mathit{\omega }}_{\mathrm{0}}^{\prime }-{\mathit{\omega }}_{\mathrm{c}}^{\prime }\right)/\mathrm{2}}{e}^{j\mathit{\phi }}\right.\\ & \left.+\mathrm{si}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(n{\mathit{\omega }}_{\mathrm{0}}^{\prime }+{\mathit{\omega }}_{\mathrm{c}}^{\prime }\right)\right)e^{-j\left(n{\mathit{\omega }}_{\mathrm{0}}^{\prime }+{\mathit{\omega }}_{\mathrm{c}}^{\prime }\right)/\mathrm{2}}{e}^{-j\mathit{\phi }}\right].\end{array}\end{array}$$

These coefficients determine the spectrum of the idealised pulsed signal with

$$\begin{array}{}\text{(16)}& S=\sum _{n=\mathrm{0}}^N{S}_n,\end{array}$$

and is shown in Fig. 4. The rise and fall times are finite in practice which limits the bandwidth (Pozar, 2005).

Figure 4 Spectrum s(t) with N=40.

To generalise the processing of the signals, it is beneficial to decompose the signal into its carrier and envelope parts. In a first step the analytic signal u_a(t) must be computed from the real signal u(t) with the Hilbert transform.

The Hilbert transform is defined as

$$\begin{array}{}\text{(17)}& v\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{u\left(\mathit{\tau }\right)}{t-\mathit{\tau }}}\mathrm{d}\mathit{\tau }\end{array}$$

and the inverse transform as

$$\begin{array}{}\text{(18)}& u\left(t\right)=-{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{v\left(\mathit{\tau }\right)}{t-\mathit{\tau }}}\mathrm{d}\mathit{\tau }.\end{array}$$

The integral representation of the Hilbert transform in Eq. (17) corresponds to the convolution integral, thus the transform can also be represented as

$$\begin{array}{}\text{(19)}& v\left(t\right)=u\left(t\right)\star {\displaystyle \frac{\mathrm{1}}{\mathit{\pi }t}},\end{array}$$

where ⋆ represents the convolution operation similar to Eq. (17). Using ⋆ the inverse Hilbert transform is

$$\begin{array}{}\text{(20)}& u\left(t\right)=-v\left(t\right)\star {\displaystyle \frac{\mathrm{1}}{\mathit{\pi }t}},\end{array}$$

with $\frac{\mathrm{1}}{\mathit{\pi }t}$ as the kernel of the Hilbert transformation.

In contrast to the Fourier Transform the Hilbert Transform does not change the domain, u(t), v(t)∈ℝ, see Eqs. (17) and (18). With the usage of Hilbert transform pairs the analytic signal u_a(t) becomes

$$\begin{array}{}\text{(21)}& {\displaystyle }{u}_{\mathrm{a}}\left(t\right)& {\displaystyle }=u\left(t\right)+jv\left(t\right)\text{(22)}& {\displaystyle }& {\displaystyle }=u\left(t\right)+j\mathcal{H}\mathit{\left\{}u\right(t\left)\mathit{\right\}}.\end{array}$$

The Fourier transform of the kernel $\frac{\mathrm{1}}{\mathit{\pi }t}$ is given by

$$\begin{array}{}\text{(23)}& \mathcal{F}\left\{{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }t}}\right\}=-j\mathrm{sgn}\left(f\right),\end{array}$$

with the signum function defined as

$$\begin{array}{}\text{(24)}& \mathrm{sgn}\left(x\right)=\left\{\begin{array}{rc}\mathrm{1}& x>\mathrm{0}\\ \mathrm{0}& x=\mathrm{0}\\ -\mathrm{1}& x<\mathrm{0}\end{array}\right..\end{array}$$

Expressing the spectrum of u_a(t) in terms of the Fourier Transform U(f) of u(t) leads to

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }{u}_{\mathrm{a}}\left(t\right)=u\left(t\right)+j{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }t}}\star u\left(t\right)\text{(26)}& {\displaystyle }& {\displaystyle }{U}_{\mathrm{a}}\left(f\right)=U\left(f\right)+\mathrm{sgn}\left(f\right)U\left(f\right).\end{array}$$

In Eqs. (25) and (26) use has been made of the linearity property of the Fourier transform, see Hoffmann (1998), and of Eq. (23).

Evaluating Eq. (26) leads to

$$\begin{array}{}\text{(27)}& U_{\mathrm{a}}\left(f\right)=\left\{\begin{array}{ll}\mathrm{2}U\left(f\right)& f>\mathrm{0}\\ U\left(\mathrm{0}\right)& f=\mathrm{0}\\ \mathrm{0}& f<\mathrm{0}\end{array}\right..\end{array}$$

The analytic signal u_a(t) in its polar form

$$\begin{array}{}\text{(28)}& {\displaystyle }{u}_{\mathrm{a}}\left(t\right)& {\displaystyle }={\widehat{u}}_{\mathrm{a}}\left(t\right)e^{j(\mathrm{2}\mathit{\pi }f+\mathit{\phi }(t\left)\right)}\text{(29)}& {\displaystyle }& {\displaystyle }={\widehat{u}}_{\mathrm{a}}\left(t\right)e^{j\mathit{\phi }\left(t\right)}{e}^{j\mathrm{2}\mathit{\pi }f}\text{(30)}& {\displaystyle }& {\displaystyle }={\underset{\mathrm{¯}}{\widehat{u}}}_{\mathrm{a}}\left(t\right)e^{j\mathrm{2}\mathit{\pi }f},\end{array}$$

with the envelope ${\widehat{u}}_{\mathrm{a}}\left(t\right)$, the complex envelope ${\underset{\mathrm{¯}}{\widehat{u}}}_{\mathrm{a}}\left(t\right)$, the instantaneous phase φ(t) and the carrier frequency f. The complex envelope is given by

$$\begin{array}{}\text{(31)}& {\underset{\mathrm{¯}}{\widehat{u}}}_{\mathrm{a}}\left(t\right)={\widehat{u}}_{\mathrm{a}}\left(t\right)e^{j\mathit{\phi }\left(t\right)},\end{array}$$

the envelope

$$\begin{array}{}\text{(32)}& {\widehat{u}}_{\mathrm{a}}\left(t\right)=\sqrt{u^{\mathrm{2}}\left(t\right)+v^{\mathrm{2}}\left(t\right)}\end{array}$$

and the instantaneous phase

$$\begin{array}{}\text{(33)}& \mathit{\phi }\left(t\right)=\arctan {\displaystyle \frac{v\left(t\right)}{u\left(t\right)}}.\end{array}$$

The transition from a time continuous signal u(t) to a time discrete signal u[k] is denoted by

$$\begin{array}{}\text{(34)}& {\displaystyle }& {\displaystyle }{t}_k=k\mathrm{\Delta }t\text{(35)}& {\displaystyle }& {\displaystyle }u\left[k\right]=u\left(t_k\right),\end{array}$$

with index k∈ℤ and the sample interval Δt.

As noted in Sect. 2 the RC is an oscillating system and can be modelled as a network of RLC circuits. Therefore the different excited modes can be modelled as a sum of exponential functions, see Sects. 6 and 7.

4 Measurements and signal processing

4.1 Measurement setup

Figure 5 shows the measurement setup.

Figure 5 Measurement setup and signal path.

A generator and an amplifier delivers the signal s to excite the RC employing the tx antenna. To record the chamber response g, a real-time oscilloscope is wired to the rx antenna of the same type as the tx antenna.

The carrier frequency f_c of the input signal is 400 MHz, with t_H=3 µs and t_L=7 µs. The power P of the input signal is set to P=50 dBm, which implies A=70.7 V in a 50 Ω system. The mode density is 8.74 modes per MHz, as calculated in Sect. 2. Figure 6 depicts the spectrum of the input signal in the frequency range from 395 to 405 MHz and the resonant frequencies, indicated by the vertical lines in the lower bar.

Figure 6 Spectrum of the input signal s(t) and resonant frequencies.

The input signal excites several modes simultaneously which can lead to high field intensities in the chamber.

4.2 Signal preprocessing

Figure 7 depicts the actual signal applied to excite the chamber when the carrier is switched on and Fig. 8 displays the signal when the carrier is switched off, where the curves clearly reflect the finite rise and fall times. The signal described in Eq. (6) has infinite rise and fall times it is therefore an idealisation of the actual signal.

Figure 7 Switching on the input signal, measured at the port of the signal generator.

Figure 8 Switching off the input signal, measured at the port of the signal generator.

The discrete Fourier Transform (DFT) of the actual input signal is computed employing the Fast Fourier Transform (FFT) algorithm. Figure 9 shows the spectrum of the signal sampled with the sample frequency f_s=25 GHz. At 400 MHz the peak of the carrier frequency is clearly visible. The dashed red box indicates the operational frequency range of the antennas according to the data sheet as their transfer functions are not known. The red arrows mark the harmonic distortions caused by elements with a non-linear transfer function, for instance the amplifier, other intermodulation products are also visible.

Figure 9 Spectrum of the input signal S(f). The red arrows mark the harmonic distortions. The dashed red box shows the frequency range of the antennas.

These frequency components will excite modes in the RC, because of the LTI nature of the system. To filter these unwanted frequency components from the chamber response a 4th-order Butterworth band-pass filter is used, with a bandwidth of 200 MHz.

The spectra of the idealised and of the actual signal are symmetrical with respect to the carrier frequency in displayed frequency range.

Figure 10 shows the response of the chamber in the time domain, sampled with f_s=2.5 GHz, where the dashes lines mark the trigger points when the carrier is switched on and off. After the trigger signal a delay of approximately 140 ns is observed, this signal propagation delay t_D is caused by the amplifier, cables and connectors. The propagation delay is frequency dependent because the transfer functions of the amplifier, cables and connectors are frequency dependent.

Figure 10 Chamber response to the actual input signal, dashed lines mark the trigger points.

The spectra of the actual input signal and the chamber response differ which implies that transfer function of the RC is frequency dependent. The peak value at the carrier frequency is decreased by 25 dB corresponding to the insertion loss of the chamber.

Figure 11 shows two chamber responses and their Fourier transforms for different stirrer positions. In time-domain the responses show different characteristics, such as settling time, overshoot and quasi-steady-state value. In frequency-domain the only observable difference is the derivation of the peak value at the carrier frequency, corresponding to the different quasi-steady-state values, indicated by the blue lines in the figure.

Figure 11 Comparison of chamber responses in time- and frequency-domain.

In Fig. 12 the black trace indicates the envelope of the chamber response, it is visible that the envelope is noisy. Figure 13 shows a chamber response for the same stirrer position as Fig. 12, where the acquisition mode is changed from single shot to average mode. In average mode the oscilloscope averages over N=150 successive traces. The envelope appears to be less noisy. This indicates that the envelope is superimposed by a mean-free white noise which does not contain any information.

Figure 12 Chamber response g[k] with envelope $\widehat{g}\left[k\right]$, acquisition mode single shot.

Figure 13 Chamber response g_avg[k] with envelope ${\widehat{g}}_{\mathrm{avg}}\left[k\right]$, acquisition mode average.

Figure 14 shows the spectra of the single shot and the averaged trace, the spectrum of the averaged signal has a sharper peak at the carrier frequency and the noise floor is significantly reduced and so the signal-to-noise ratio increased. The peak at 800 MHz is the second harmonic of the carrier frequency, generated by the signal generator or the amplifier.

Figure 14 Spectra of chamber response, recorded with single shot G(f) and average mode G_avg(f).

Figure 15 shows the spectrum of the filtered chamber response g[k] and its complex envelope $\widehat{g}\left[k\right]$, see Eq. (31). The spectra is band-limited and centred around f=0 Hz.

Figure 15 Spectra of the averaged chamber response g_avg and its complex envelope ${\widehat{g}}_{\mathrm{avg}}$.

It is now possible to downsample the complex envelope of the chamber response, because the upper frequency limit f_ge of the complex envelope is much smaller than the upper frequency limit f_g of g[k].

The signal g[k] has P samples with k∈I_k and

$$\begin{array}{}\text{(36)}& I_k:=[\mathrm{0},P-\mathrm{1}].\end{array}$$

To downsample g[k] it must first be multiplied with a D-periodic sequence of unit pulses

$$\begin{array}{}\text{(37)}& g_{\mathrm{D}\downarrow }=\widehat{g}\left[k\right]\cdot \sum _{m=-\mathrm{\infty }}^{\mathrm{\infty }}\mathit{\delta }[k-mD]\end{array}$$

with the decimation factor D

$$\begin{array}{}\text{(38)}& D=⌊{\displaystyle \frac{f_{\mathrm{g}}}{f_{\mathrm{ge}}}}⌋.\end{array}$$

The downsampled signal consists of every Dth value of g[k]

$$\begin{array}{}\text{(39)}& g_u\left[k\right]=g\left[kD\right]=g_{\mathrm{D}\downarrow }\left[kD\right],\end{array}$$

and the sample interval Δt_u

$$\begin{array}{}\text{(40)}& \mathrm{\Delta }{t}_u=D\cdot \mathrm{\Delta }t.\end{array}$$

The signal g_u[k] consists of $P^{\prime }=⌊\frac{P}{D}⌋$ samples with $k\in I_k^{\prime }$ and

$$\begin{array}{}\text{(41)}& I_k^{\prime }:=[\mathrm{0},P^{\prime }-\mathrm{1}].\end{array}$$

5 Signal analysis

In this section an algorithm is presented that extracts signal properties like rise time and overshoot amplitude. Table A1 lists all properties of interest. This algorithm works on a single chamber response, in contrast to Artz (2017) which needs to process a whole set of measurements first.

The averaged downsampled envelope of the input signal s_u,avg is designated as the input signal x and the averaged downsampled envelope of the chamber response ${\widehat{g}}_{u,\mathrm{avg}}$ is defined as the output signal y

$$\begin{array}{}\text{(42)}& {\displaystyle }& {\displaystyle }x\left[k\right]:={\widehat{s}}_{u,\mathrm{avg}}\left[k\right]\text{(43)}& {\displaystyle }& {\displaystyle }y\left[k\right]:={\widehat{g}}_{u,\mathrm{avg}}\left[k\right].\end{array}$$

5.1 Outline of the algorithm

In order to extract the signal properties of interest the signal is partitioned into four parts

1. rise phase,

2. quasi-steady state (QSS),

3. decay phase,

4. off phase.

First, the limits of the quasi-steady state are identified. Knowledge of the quasi-steady state allows to derive the start and end points of the remaining phases. A sequence of tasks is necessary to determine the start and end point of the quasi-steady state. These steps are

• i.

  smoothing of the envelope,

• ii.

  clipping,

• iii.

  computation time-derivative,

• iv.

  deriving threshold for detecting the quasi-steady state,

• v.

  identification of start and end points.

With the location of the quasi-steady state the amplitude is estimated. The algorithm concludes with the determination of the requested signal properties from the previously computed results.

5.2 Algorithm

In the following section the single steps for the determination of the quasi-steady state are presented.

5.2.1 Smoothing

Signal smoothing is accomplished by employing a Savitzky-Golay filter, see Savitzky and Golay (1964). In our example we used a window length of 83 and a polynomial degree of 3.

5.2.2 Clipping

The quasi-steady state has to be within the boundary of the trigger time and the switch off time, see Fig. 16, which mark the utmost limits of the QSS. This time-span covers the Rise phase and the QSS. All signal values outside of the interval are neglected. The index of the trigger time, the trigger point k_T, and the pulse width t_H is known. The start point $k_{\mathrm{s}}^{\prime }$ and the end point $k_{\mathrm{e}}^{\prime }$ of the clipped signal are

$$\begin{array}{}\text{(44)}& k_{\mathrm{S}}^{\prime }=k_{\mathrm{T}}{k}_{\mathrm{e}}^{\prime }=k_{\mathrm{T}}+⌊{\displaystyle \frac{t_{\mathrm{H}}}{\mathrm{\Delta }{t}_u}}⌋.\end{array}$$

In Sect. 4.2 a propagation delay was observed in the chamber response right after the trigger point. This delay is taken into account by shifting the former start and end points (Eq. 44) to

$$\begin{array}{}\text{(45)}& k_{\mathrm{s}}^{\prime }:=k_{\mathrm{s}}^{\prime }+⌊{\displaystyle \frac{t_{\mathrm{D}}}{\mathrm{\Delta }{t}_u}}⌋k_{\mathrm{e}}^{\prime }:=k_{\mathrm{e}}^{\prime }+⌊{\displaystyle \frac{t_{\mathrm{D}}}{\mathrm{\Delta }{t}_u}}⌋,\end{array}$$

where t_D=140 ns. This gives the clipped signal

$$\begin{array}{}\text{(46)}& y_{\mathrm{c}}\left[k\right]=y\left[k\right],\mathrm{with}k\in \left[k_{\mathrm{s}}^{\prime },k_{\mathrm{e}}^{\prime }\right].\end{array}$$

Figure 16 The envelope of the chamber response in grey, the smoothed signal y in black and the time-span of the clipped signal y_c is indicated by the vertical dashed lines.

5.2.3 Time-derivative

An ideal QSS is defined by a constant amplitude over a certain time period. The derivative for this time period is zero. In practice a deviation around zero is expected and an alteration limit needs to be defined, see Fig. 17.

Figure 17 Time-derivative $y_{\mathrm{c}}^{\prime }\left[k\right]$. The dotted lines mark the threshold ε.

The time-derivative is calculated with a forward-difference quotient

$$\begin{array}{}\text{(47)}& y_{\mathrm{c}}^{\prime }\left[k\right]={\displaystyle \frac{y_{\mathrm{c}}[k+\mathrm{1}]-y_{\mathrm{c}}\left[k\right]}{\mathrm{\Delta }{t}_u}}.\end{array}$$

5.2.4 Threshold determination

In order to determine the maximum deviation around zero for the QSS a cumulative histogram is used, see Fig. 18. The ordinate is restricted to the main interval I which is partitioned into N sub-intervals

$$\begin{array}{}\text{(48)}& I=\bigcup _{j=\mathrm{1}}^N{I}_j.\end{array}$$

The boundaries of the sub-intervals $I_j\in [k_{l,j}$, k_u,j] are the edges of the bin b_j, which is the cardinality of I_j

$$\begin{array}{}\text{(49)}& b_j=\left|I_j\right|\in \mathbb{N}.\end{array}$$

A threshold for the QSS value is obtained from the lower boundary of the interval $I_{n_{\mathrm{QSS}}}=[k_{l,n_{\mathrm{QSS}}},k_{u,n_{\mathrm{QSS}}}]$. The index of this interval n_QSS has to fulfil the condition

$$\begin{array}{}\text{(50)}& n_{\mathrm{QSS}}=\min \left\{n\left|s\right(n)\ge \mathrm{0.8},n=\mathrm{1},\mathrm{2},\mathrm{\dots },N\right\},\end{array}$$

where

$$\begin{array}{}\text{(51)}& s\left(n\right)={\displaystyle \frac{\mathrm{1}}{B}}\sum _{j=\mathrm{1}}^n{b}_{j}B=\sum _{j=\mathrm{1}}^N{b}_j.\end{array}$$

The threshold is

$$\begin{array}{}\text{(52)}& \left|\mathit{\epsilon }\right|\le y_{\mathrm{c}}^{\prime }\left[k_{l,n_{\mathrm{QSS}}}\right].\end{array}$$

Figure 18 Histogram with 100 bins of equal width, the cumulative percentage and the threshold ε. The x-axis is truncated at 0.06 V µs⁻¹ to increase the readability.

5.2.5 Start and End Point of the QSS

In the final step of the algorithm the start and end points of the QSS are determined. The QSS is defined by coherent regions where every value has to be smaller or equal than the threshold. These regions are

$$\begin{array}{}\text{(53)}& R_i:=\left\{(k,x^{\prime })|\left|y_{\mathrm{c}}^{\prime }\left[k\right]\right|\le \mathit{\epsilon },\forall k\in \left[k_i,k_i+L\right]\right\},\end{array}$$

with $\mathrm{1}\le i\le M$, the window length L and

$$\begin{array}{}\text{(54)}& M=⌊{\displaystyle \frac{k_{\mathrm{e}}^{\prime }-k_{\mathrm{s}}^{\prime }}{L}}⌋.\end{array}$$

If multiple regions satisfy the criterion and if these regions contain subsets of the previous regions they are merged into one coherent region. If the distance Δt_G between subsequent but disjoint regions is small enough the gap is neglected. In our case we used a maximum gap size of

$$\mathrm{\Delta }{t}_{\mathrm{G}}\le \mathit{\delta }\cdot \mathrm{\Delta }{t}_u=\mathrm{768}\mathrm{ns}$$

with δ=200. The merging process can have one of three possible outcomes. Either no region, one region or multiple disjoint regions are found. When no region is found no QSS exists.

The optimal outcome is that only one coherent segment results from the process. The identification of the QSS is obvious.

If more than one coherent segment have been found, the longest segment is chosen as the actual QSS.

Where y[k] is now divided into five segments

$$\begin{array}{}\text{(55)}& y\left[k\right]=\left\{\begin{array}{ll}\mathrm{Off}\mathrm{phase},& \mathrm{0}\le k<k_{\mathrm{s},\mathrm{rise}},\\ \mathrm{Rise}\mathrm{phase},& k_{\mathrm{s},\mathrm{rise}}\le k<k_{\mathrm{s},\mathrm{QSS}},\\ \mathrm{QSS},& k_{\mathrm{s},\mathrm{QSS}}\le k<k_{\mathrm{s},\mathrm{decay}},\\ \mathrm{Decay}\mathrm{phase},& k_{\mathrm{s},\mathrm{decay}}\le k<k_{\mathrm{s},\mathrm{off}},\\ \mathrm{Off}\mathrm{phase},& k_{\mathrm{s},\mathrm{off}}\le k<P^{\prime }-\mathrm{1},\end{array}\right.\end{array}$$

with

$$k_{\mathrm{s},\mathrm{rise}}=\min \left\{k\left|y\right[k]\ge {\displaystyle \frac{y_{\mathrm{QSS}}}{\mathrm{10}}},k=\mathrm{0},\mathrm{\dots },k_{\mathrm{s},\mathrm{QSS}}\right\},$$

$$k_{\mathrm{s},\mathrm{decay}}=k_{\mathrm{e},\mathrm{QSS}}$$

and

$$k_{\mathrm{s},\mathrm{off}}=\max \left\{k\left|y\right[k]\ge {\displaystyle \frac{y_{\mathrm{QSS}}}{\mathrm{10}}},k=k_{\mathrm{e},\mathrm{QSS}},\mathrm{\dots },P^{\prime }-\mathrm{1}\right\}.$$

The QSS mean value is

$$\begin{array}{}\text{(56)}& y_{\mathrm{QSS}}={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{e},\mathrm{QSS}}-k_{\mathrm{s},\mathrm{QSS}}}}\sum _{n=k_{\mathrm{s},\mathrm{QSS}}}^{k_{\mathrm{e},\mathrm{QSS}}}y\left[n\right].\end{array}$$

5.3 Results

The previously determined steady-state value y_QSS is now used as a normalisation parameter to simplify comparability between different chamber responses

$$\begin{array}{}\text{(57)}& \stackrel{\mathrm{‾}}{y}\left[k\right]={\displaystyle \frac{y\left[k\right]}{y_{\mathrm{QSS}}}}.\end{array}$$

An exemplary result is presented in Fig. 19 consisting of a measured noisy signal y[k] and its noise reduced part. The extracted boundaries for the four signal states are indicated by dashed vertical lines, further properties are listed in Table 1.

Figure 19 The envelope of the chamber response in black, the smoothed signal y in red. The vertical dashed lines indicate the limits of the signal parts. The Rise and the Decay phase are approximated with damped oscillations.

Due to the oscillating nature of the measured signals, the system responses are interpreted as oscillating systems. Based on this assumption it is possible to determine time constants for the rise and fall of each signal. A simple damped harmonic oscillation is

$$\begin{array}{}\text{(58)}& \begin{array}{rl}y\left(t\right)=A+e^{-T\left(t-t_{\mathrm{o}}\right)}& (B\cos \left(W\left(t-t_{\mathrm{o}}\right)\right)\\ & -C\sin \left(W\left(t-t_{\mathrm{o}}\right)\right)),T>\mathrm{0}\end{array}\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }B=\mathit{\varrho }\cos \mathit{\phi }\\ {\displaystyle }& {\displaystyle }C=\mathit{\varrho }\sin \mathit{\phi }\end{array}$$

where A, B, C, T and W are unknowns and t_o is a known shift along the t-axis. The sought-after time constant is then T⁻¹.

Computation of the unknowns in Eq. (58) is accomplished with the help of a Levenberg–Marquardt routine, an iterative trust-region algorithm (Marquardt, 1963). A necessity for convergence are well-posed starting values. Those values are obtained through a reduced model for A and T

$$\begin{array}{}\text{(59)}& y\left(t\right)=\stackrel{\mathrm{̃}}{A}+\stackrel{\mathrm{̃}}{B}{e}^{-\stackrel{\mathrm{̃}}{T}\left(t-t_{\mathrm{o}}\right)},\stackrel{\mathrm{̃}}{T}>\mathrm{0},\end{array}$$

described in Appendix A. The starting values for the fitting process are then

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{A}_{\mathrm{0}}=\stackrel{\mathrm{̃}}{A}\\ {\displaystyle }& {\displaystyle }{B}_{\mathrm{0}}=\mathrm{1}\\ {\displaystyle }& {\displaystyle }{C}_{\mathrm{0}}=\mathrm{1}\\ {\displaystyle }& {\displaystyle }{T}_{\mathrm{0}}=\stackrel{\mathrm{̃}}{T}\\ {\displaystyle }& {\displaystyle }{W}_{\mathrm{0}}={\displaystyle \frac{\mathit{\pi }}{x\left[\lfloor \frac{k_{\mathrm{s},\mathrm{QSS}}+k_{\mathrm{e},\mathrm{QSS}}}{\mathrm{2}}\rfloor \right]}}.\end{array}$$

The determined time constants for the chamber response depicted in Fig. 19 are for the the Rise and Decay phases

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{rise}}=\mathrm{674}\mathrm{ns}\\ {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{decay}}=\mathrm{2.24}\mathrm{\mu }\mathrm{s}.\end{array}$$

Table 1 Extracted signal properties.

The described algorithm is now applied to a series of measurements of L chamber responses. Measurements were taken over a full stirrer rotation at a stirrer step size of 1^∘, this leads to L=360 chamber responses. The measurement setup corresponds to that described in Sect. 4.1.

The distribution of incidences of selected properties are shown in Figs. 20–22. Those properties include the QSS value y_QSS, QSS mean value normalised to the pulse width

$$\mathrm{\Delta }{\stackrel{\mathrm{̃}}{t}}_{\mathrm{QSS}}={\displaystyle \frac{\mathrm{\Delta }{t}_{\mathrm{QSS}}}{t_{\mathrm{H}}}}$$

and the peak value in the Decay phase ${\widehat{y}}_{\mathrm{decay}}$. The bins are of equal width.

Figure 20 Histogram of the QSS mean value y_QSS with 50 bins. The dashed line indicates the mean value.

In Table 2 the population mean

$$\begin{array}{}\text{(60)}& \mathit{\mu }\left(y\right)={\displaystyle \frac{\mathrm{1}}{L}}\sum _{k=\mathrm{1}}^L{y}_k,\mathrm{with}y\in {\mathbb{R}}^L,\end{array}$$

and the population standard deviation

$$\begin{array}{}\text{(61)}& \mathit{\sigma }\left(y\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{L}}\sum _{k=\mathrm{1}}^L{\left(y_k-\mathit{\mu }\right)}^{\mathrm{2}}},\mathrm{with}y\in {\mathbb{R}}^L,\end{array}$$

as a measure of spread are computed for the three properties.

Table 2 Population mean μ and standard deviation σ computed for selected properties.

The standard deviation indicates a wide spread around the mean value of the y_QSS. This means that the field strength varies notably over the DUT for a full stirrer rotation.

Figure 21 shows the distribution of incidences of the duration of the quasi steady state in relation to the pulse width, it can be seen that the QSS value is present on average for only 57.1 % of the pulse width (Table 2).

Figure 21 Histogram of the QSS mean value $\mathrm{\Delta }{\stackrel{\mathrm{̃}}{t}}_{\mathrm{QSS}}$ normalised to the pulse width t_H, with 50 bins. The dashed line indicates the mean value.

The spread of the peak values in the Rise and Decay phase around their mean value is smaller than for the QSS mean value. In the histogram this can be seen by an accumulation of values around the mean value in Fig. 22.

Figure 22 Histogram of the peak values in the Decay phase ${\widehat{y}}_{\mathrm{decay}}$, with 50 bins. The dashed line indicates the mean value.

The relationship between the overshoot and the QSS mean value is shown in Fig. 23. Each QSS mean value y_QSS is the reference value to the corresponding overshoot ${\widehat{y}}_{\mathrm{rise}}$ of the data set. A corresponding example is given in Table 1 and the explanation of the symbols is given in Table A1. The overshoot decreases with increasing QSS mean value, this is in agreement with Artz (2017). This indicates that in this case, with the appropriate choice of test parameters, over-testing will not occur.

Figure 23 Plot of the signals overshoot ${\widehat{y}}_{\mathrm{rise}}$ normalised to the corresponding mean value in the quasi steady-state y_QSS over the mean value in the quasi steady-state.

The results of the algorithm are the foundation for further investigations.

6 Method of Prony

With Prony's Method a decomposition of a signal into a sum of weighted complex exponentials can be performed (Marple, 2019; Hauer et al., 1990). From this representation Amplitude, phase shift, attenuation and angular frequency can be determined.

Prony's Method is sensitive to noise (Marple, 2019; Peter, 2013), it is recommended to reduce noise in the input signal before applying the method. Noise reduction is achieved by the signal pre-processing described in Sect. 4.2.

6.1 Description of the method

The signal y[k] is approximated with

$$\begin{array}{}\text{(62)}& y\left[k\right]\approx \stackrel{\mathrm{̃}}{y}\left[k\right]=\sum _{n=\mathrm{0}}^{p-\mathrm{1}}{h}_n{z}_n^k,k\in [\mathrm{0},M-\mathrm{1}],\end{array}$$

where M≥2p and

$$\begin{array}{}\text{(63)}& {\displaystyle }& {\displaystyle }{z}_n=e^{\left({\mathit{\alpha }}_n+j\mathrm{2}\mathit{\pi }{f}_n\right)\mathrm{\Delta }t}\text{(64)}& {\displaystyle }& {\displaystyle }{h}_n={\widehat{y}}_n{e}^{j{\mathrm{\Omega }}_n}.\end{array}$$

${\widehat{y}}_n$ is an amplitude, Ω_n a phase shift, α_n a damping factor and f_n is a frequency.

The unknown factors appear linear as pre-factors h_n and non-linear as arguments of an exponential function z_n. These unknown parameters are determined in two steps. First the exponents are determined and then the coefficients.

The core of Prony's method is the assumption that the sum of the exponential functions satisfies a pth order difference equation

$$\begin{array}{}\text{(65)}& \left(\begin{array}{cccc}{y}_p& y_{p-\mathrm{1}}& \mathrm{\dots }& y_{\mathrm{0}}\\ y_{p+\mathrm{1}}& y_p& \mathrm{\dots }& y_{\mathrm{1}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ y_{\mathrm{2}p-\mathrm{1}}& y_{\mathrm{2}p-\mathrm{2}}& \mathrm{\dots }& y_{p-\mathrm{1}}\end{array}\right)\left(\begin{array}{c}{a}_{\mathrm{0}}\\ a_{\mathrm{1}}\\ \mathrm{⋮}\\ a_p\end{array}\right)=\left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{⋮}\\ \mathrm{0}\end{array}\right)\end{array}$$

(Hamming, 1987). The characteristic polynomial of this difference equation is

$$\begin{array}{}\text{(66)}& \mathrm{0}=a_{\mathrm{0}}\cdot z^p+a_{\mathrm{1}}\cdot z^{p-\mathrm{1}}+\mathrm{\dots }+a_p=\sum _{m=\mathrm{0}}^p{a}_m{z}^{p-m},\end{array}$$

where a_m∈ℝ. The z of Eq. (62) are the zeros of Eq. (66), which are either real or conjugate complex. To determine these zeros, the coefficients a₀, …, a_p must be known. If a₀=1, a linear system for the unknowns a₁, …, a_p is obtained from Eq. (65) as

$$\begin{array}{}\text{(67)}& \left(\begin{array}{cccc}{y}_{p-\mathrm{1}}& y_{p-\mathrm{2}}& \mathrm{\dots }& y_{\mathrm{0}}\\ y_p& y_p& \mathrm{\dots }& y_{\mathrm{1}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ y_{\mathrm{2}p}& y_{\mathrm{2}p-\mathrm{1}}& \mathrm{\dots }& y_{p-\mathrm{1}}\end{array}\right)\left(\begin{array}{c}{a}_{\mathrm{1}}\\ a_{\mathrm{2}}\\ \mathrm{⋮}\\ a_p\end{array}\right)=-\left(\begin{array}{c}{y}_p\\ y_{p+\mathrm{1}}\\ \mathrm{⋮}\\ y_{\mathrm{2}p-\mathrm{1}}\end{array}\right).\end{array}$$

With the exponents z determined from Eqs. (66) and (67), the coefficients h₁, …, h_p can be calculated from the data. Re-writing Eq. (62) in matrix notation gives

$$\begin{array}{}\text{(68)}& \left(\begin{array}{cccc}{z}_{\mathrm{0}}^{\mathrm{0}}& z_{\mathrm{1}}^{\mathrm{0}}& \mathrm{\dots }& z_{p-\mathrm{1}}^{\mathrm{0}}\\ z_{\mathrm{0}}^{\mathrm{1}}& z_{\mathrm{2}}^{\mathrm{1}}& \mathrm{\dots }& z_{p-\mathrm{1}}^{\mathrm{1}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ z_{\mathrm{0}}^{M-\mathrm{1}}& z_{\mathrm{2}}^{M-\mathrm{1}}& \mathrm{\dots }& z_{p-\mathrm{1}}^{M-\mathrm{1}}\end{array}\right)\left(\begin{array}{c}{h}_{\mathrm{0}}\\ h_{\mathrm{1}}\\ \mathrm{⋮}\\ h_{p-\mathrm{1}}\end{array}\right)=\left(\begin{array}{c}{y}_{\mathrm{0}}\\ y_{\mathrm{1}}\\ \mathrm{⋮}\\ y_{p-\mathrm{1}}\end{array}\right).\end{array}$$

From the roots z_n, n=1, …, p the searched for arguments of the exponential functions are

$$\begin{array}{}\text{(69)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_n={\displaystyle \frac{\ln \left|z_n\right|}{\mathrm{\Delta }t}}\text{(70)}& {\displaystyle }& {\displaystyle }{f}_n={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathrm{\Delta }t}}{\tan }^{-\mathrm{1}}\left({\displaystyle \frac{\mathrm{\Im }\left\{z_n\right\}}{\mathrm{\Re }\left\{z_n\right\}}}\right).\end{array}$$

The linear coefficients h_n give the parameters

$$\begin{array}{}\text{(71)}& {\displaystyle }& {\displaystyle }{\widehat{y}}_n=\left|h_n\right|\text{(72)}& {\displaystyle }& {\displaystyle }{\mathrm{\Omega }}_n={\tan }^{-\mathrm{1}}{\displaystyle \frac{\mathrm{\Im }\left\{h_n\right\}}{\mathrm{\Re }\left\{h_n\right\}}}.\end{array}$$

The approximated data set is written as

$$\begin{array}{}\text{(73)}& \stackrel{\mathrm{̃}}{y}\left[k\right]=\sum _{n=\mathrm{0}}^{p-\mathrm{1}}{\widehat{y}}_n{e}^{j{\mathrm{\Omega }}_n}{e}^{\left({\mathit{\alpha }}_n+j\mathrm{2}\mathit{\pi }{f}_n\right)k\mathrm{\Delta }t}.\end{array}$$

In Eqs. (67) and (68) the matrices are square and can get very large for big sets of data. But in practice the number of data points exceeds the number of unknowns in the exponential model (Eq. 62). This will lead to non-square matrices where the number of data points exceeds the number of unknowns. Various methods are available to solve these over-determined linear systems, for example least squares (Marple, 2019; Peter, 2013). In the present case the least squares approach is used. The method was tested on the transient and decay behaviour of the curve in Fig. 19. The necessary parameters such as the beginning and end of the corresponding phases were taken from Table 1. The data set for the transient oscillation comprises M=1133 data points and p=380 summands were chosen. The actual and approximate course is shown in Fig. 24. There is a very good agreement up to t=3 µs, after that the curves deviate noticeably. The approximation can be improved by increasing the summands in Eq. (62).

Figure 24 Rise phase of signal y[k] and Prony interpolation $\stackrel{\mathrm{̃}}{y}\left[k\right]$. The black line is the rise phase of signal y[k] and the dashed red line is the Prony interpolation $\stackrel{\mathrm{̃}}{y}\left[k\right]$.

The data set for the Decay phase includes M=1579 data points and p=450 parameters were used. The approximation and the actual course agree well up to about 13 µs (Fig. 25).

Figure 25 Decay phase of signal y[k] and Prony interpolation $\stackrel{\mathrm{̃}}{y}\left[k\right]$. The black line is the decay phase of signal y[k] and the dashed red line is the Prony interpolation $\stackrel{\mathrm{̃}}{y}\left[k\right]$.

Figure 26 shows the poles according to Eq. (66). A strong spread of the angular frequency is visible. However, since it is a smoothed and band-limited signal, some of the poles found lie outside the permissible frequency band. This means that the curves interpolated with the Prony method contain frequency components which are not physical. Band-limiting leads to a deterioration of the interpolation.

Figure 26 Conjugate complex poles of $\stackrel{\mathrm{̃}}{y}\left(t\right)$ with t in   to 5.35 µs.

An alternative interpolation method is the TD-VF where pole locations within a given frequency band are selected.

7 Time-domain vector-fit

Prony's method represents the measured chamber response as a sum of exponentials. This sum is the result of a convolution of a transfer-function and an excitation. Purely employing Prony's method makes it difficult to extract information about underlying characteristics of the signal transfer. An alternative method is Time-Domain Vector-Fit (TD-VF), which is based on the Vector-Fit (VF) algorithm (Gustavsen and Semlyen, 1999; Triverio, 2019).

Vector-Fit offers the possibility to approximate a transfer-function to a measured response. In the following sub-sections are VF and TD-VF briefly summarised. Afterwards an application to measured data is presented.

7.1 Fundamental vector-fit

The VF method approximates a rational function to frequency domain data using an iterative least squares technique. The transfer-function in frequency domain is

$$\begin{array}{}\text{(74)}& Y\left(s\right)=H\left(s\right)X\left(s\right),\end{array}$$

with Y(s) the result, H(s) the transfer function, X(s) the excitation and s a complex frequency parameter. The VF approximation of H(s) is

$$\begin{array}{}\text{(75)}& H\left(s\right)\approx d+\sum _{n=\mathrm{1}}^N{\displaystyle \frac{R_n}{s-p_n}},\end{array}$$

where d is a constant, R_n the nth residue to the pole p_n and N the number of real poles. The iterative least squares approximation is based on an initial set of poles, which are distributed either linearly or logarithmicaly over the investigated frequency range (Grivet-Talocia and Gustavsen, 2016). The main idea in VF is to weigh the input with a function σ

$$\begin{array}{}\text{(76)}& \mathit{\sigma }\left(s\right)Y\left(s\right)=\mathit{\sigma }\left(s\right)H\left(s\right)X\left(s\right)\end{array}$$

and enforce the product σ(s)H(s) to match the left-hand-side, with

$$\begin{array}{}\text{(77)}& \mathit{\sigma }\left(s\right)H\left(s\right)=d+\sum _{n=\mathrm{1}}^N{\displaystyle \frac{K_n}{s-p_n}},\end{array}$$

and

$$\begin{array}{}\text{(78)}& \mathit{\sigma }\left(s\right)=\mathrm{1}+\sum _{n=\mathrm{1}}^N{\displaystyle \frac{{\stackrel{\mathrm{̃}}{R}}_n}{s-p_n}}.\end{array}$$

Letting X(s) to be unity, then the unknowns d, K_n and ${\stackrel{\mathrm{̃}}{R}}_n$ form a linear least-squares problem for every data point Y(s_i)

$$\begin{array}{}\text{(79)}& \sum _{n=\mathrm{1}}^N{\displaystyle \frac{K_n}{s-p_n}}+d-\sum _{n=\mathrm{1}}^N{\displaystyle \frac{{\stackrel{\mathrm{̃}}{R}}_n}{s-p_n}}Y\left(s_i\right)=Y\left(s_i\right),\end{array}$$

i=1, 2, …, M_f, where M_f represents the number of frequency samples. When the least-squares problem has reached a solution, then the new poles for H(s) are the zeros of σ(s)

$$\begin{array}{}\text{(80)}& \mathrm{1}+\sum _{n=\mathrm{1}}^N{\displaystyle \frac{{\stackrel{\mathrm{̃}}{R}}_n}{s-p_n}}={\displaystyle \frac{\prod _{n=\mathrm{1}}^N(s-\stackrel{\mathrm{̃}}{z})}{\prod _{n=\mathrm{1}}^N\left(s-p_n\right)}}=\mathit{\sigma }\left(s\right)\end{array}$$

$$\begin{array}{}\text{(81)}& \begin{array}{rl}{\displaystyle \frac{\prod _{n=\mathrm{1}}^N(s-\stackrel{\mathrm{̃}}{z})}{\prod _{n=\mathrm{1}}^N\left(s-p_n\right)}}& Y\left(s\right){\displaystyle \frac{\prod _{n=\mathrm{1}}^N(s-z)}{\prod _{n=\mathrm{1}}^N\left(s-p_n\right)}}=\mathit{\sigma }\left(s\right)H\left(s\right)\\ & Y\left(s\right)={\displaystyle \frac{\prod _{n=\mathrm{1}}^N(s-z)}{\prod _{n=\mathrm{1}}^N(s-\stackrel{\mathrm{̃}}{z})}}.\end{array}\end{array}$$

With the relocated poles $p_n={\stackrel{\mathrm{̃}}{z}}_n$ the original problem (Eq. 75) is an over-determined linear system in the unknown residues R_n.

7.2 TD-VF

Similar to the VF in frequency domain the TD-VF solves a linearized least-squares problem in time domain based on a known excitation x(t) and measured response y(t) (Grivet-Talocia, 2004). Applying an inverse Laplace transformation on Eq. (74) leads to

$$\begin{array}{}\text{(82)}& y\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}h\left(\mathit{\tau }\right)x(t-\mathit{\tau })\mathrm{d}\mathit{\tau }.\end{array}$$

The approximated time domain transfer-function is

$$\begin{array}{}\text{(83)}& h\left(t\right)\approx d\mathit{\delta }\left(t\right)+\sum _{n=\mathrm{1}}^N{R}_n{e}^{p_{n}t}\mathbf{1}\left(t\right),\end{array}$$

with δ(t) the Dirac-function and 1(t) the unit step function. In order to determine the transfer function h(t) an estimate for initial poles has to be given first. Afterwards, a least-squares solution for the unknowns similar to Eq. (79) allows to determine relocated poles. Equation (79) in time domain is

$$\begin{array}{}\text{(84)}& \begin{array}{rl}dx\left(t\right)& +K_n\underset{\mathrm{0}}{\overset{t}{\int }}{e}^{p_n(t-\mathit{\tau })}x\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }\\ & -\stackrel{\mathrm{̃}}{R}\underset{\mathrm{0}}{\overset{t}{\int }}{e}^{p_n(t-\mathit{\tau })}y\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }=y\left(t\right).\end{array}\end{array}$$

The computationally demanding task of evaluating the convolution integrals simplifies with the usage of a recursive convolution algorithm (Semlyen and Dabuleanu, 1975). With the previously determined poles equation (Eq. 82) allows for the determination of the residues R_n in the least-squares sense.

In the case of conjugate complex poles Eq. (83) becomes

$$\begin{array}{}\text{(85)}& \begin{array}{rl}h\left(t\right)\approx d\mathit{\delta }\left(t\right)& +\sum _{n=\mathrm{1}}^{N_{\mathrm{r}}}{R}_n{e}^{p_{\mathrm{r},n}t}\mathbf{1}\left(t\right)\\ & +\mathrm{2}\sum _{n=\mathrm{1}}^{N_{\mathrm{c}}}\mathrm{\Re }\left\{C_n{e}^{p_{\mathrm{c},n}t}\right\}\mathbf{1}\left(t\right)\\ & -\mathrm{2}\sum _{n=\mathrm{1}}^{N_{\mathrm{c}}}\mathrm{\Im }\left\{C_n{e}^{p_{\mathrm{c},n}t}\right\}\mathbf{1}\left(t\right),\end{array}\end{array}$$

with N_r the number of real poles, N_c the number of complex poles, p_r,n are real poles and p_c,n are complex poles. In Eq. (85) the complex poles and residues appear in complex conjugate pairs, the total number of complex poles and likewise the total number of residues is 2N_c.

If the deviation between approximated response $\stackrel{\mathrm{̃}}{y}\left(t\right)$ and original response y(t) is smaller than a previously determined threshold the iteration stops. The deviation is evaluated in terms of a two-norm

$$\begin{array}{}\text{(86)}& \mathrm{RMS}={\displaystyle \frac{\mathrm{1}}{\sqrt{M}}}\left|\right|\mathit{y}-\stackrel{\mathrm{̃}}{\mathit{y}}|{|}_{\mathrm{2}},\mathit{y},\stackrel{\mathrm{̃}}{\mathit{y}}\in {\mathbb{R}}^M,\end{array}$$

which is called root mean square (RMS). M is the number of time samples. The vectors y and $\stackrel{\mathrm{̃}}{\mathit{y}}$ denote function values for each time step. A normalised error measure is

$$\begin{array}{}\text{(87)}& \mathrm{RMSE}={\displaystyle \frac{\left|\right|\mathit{y}-\stackrel{\mathrm{̃}}{\mathit{y}}|{|}_{\mathrm{2}}}{\left|\right|{\mathit{y}}_{\mathrm{2}}\left|\right|}}.\end{array}$$

7.3 Application

The TD-VF is applied to two measurements for two different but fixed stirrer positions and a fixed antenna position within the chamber. Before TD-VF is used a low-pass filter is applied to reduce noise and distortions, the procedure is described in Sect. 4.2.

A common difficulty is the prior knowledge of the number of required poles. So far there is no general procedure available that can solve this problem, apart from iterative solutions. It is possible to locate regions in frequency domain within which poles ought to occur, these regions are in between two local minima. This method is only suitable for frequency responses (Heindl et al., 2009).

The envelope of the excitation signal is depicted in Fig. 27. All following investigations relate to the envelope of the input signal. The reconstructed transfer-function (Eq. 85) is actually the transfer-function of the envelope of the measured response.

Figure 27 Excitation function.

In order to determine the transfer-function from Eq. (84) and afterwards from Eq. (82) the following abbreviations are used

$$\begin{array}{}\text{(88)}& {\displaystyle }& {\displaystyle }{t}_k=k\mathrm{\Delta }t\text{(89)}& {\displaystyle }& {\displaystyle }x\left[k\right]=x\left(t_k\right)\text{(90)}& {\displaystyle }& {\displaystyle }y\left[k\right]=y\left(t_k\right)\text{(91)}& {\displaystyle }& {\displaystyle }{x}_n\left(k\right)=\underset{\mathrm{0}}{\overset{t_k}{\int }}{e}^{p_n\left(t_k-\mathit{\tau }\right)}x\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }\text{(92)}& {\displaystyle }& {\displaystyle }{y}_n\left(k\right)=\underset{\mathrm{0}}{\overset{t_k}{\int }}{e}^{p_n\left(t_k-\mathit{\tau }\right)}y\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }.\end{array}$$

The linear system (Eq. 84) is of dimension $M\times (\mathrm{1}+\mathrm{2}{N}_{\mathrm{r}}+\mathrm{2}{N}_{\mathrm{c}})$ where one row corresponds to a time step t_k

$$\begin{array}{}\text{(93)}& \begin{array}{rl}dx\left[k\right]& +\sum _{n=\mathrm{0}}^{N_{\mathrm{r}}}{K}_{\mathrm{r},n}{x}_n\left(k\right)+\mathrm{2}\sum _{n=\mathrm{0}}^{N_{\mathrm{c}}}{K}_{\mathrm{c},n}^{\prime }\mathrm{\Re }\left\{x_n\left(k\right)\right\}\\ & -\mathrm{2}\sum _{n=\mathrm{0}}^{N_{\mathrm{c}}}{K}_{\mathrm{c},n}^{\prime \prime }\mathrm{\Im }\left\{x_n\left(k\right)\right\}-\sum _{n=\mathrm{0}}^{N_{\mathrm{r}}}{\stackrel{\mathrm{̃}}{R}}_{\mathrm{r},n}{y}_n\left(k\right)\\ & -\mathrm{2}\sum _{n=\mathrm{0}}^{N_{\mathrm{c}}}{\stackrel{\mathrm{̃}}{R}}_{\mathrm{c},n}^{\prime }\mathrm{\Re }\left\{y_n\left(k\right)\right\}+\mathrm{2}\sum _{n=\mathrm{0}}^{N_{\mathrm{c}}}{\stackrel{\mathrm{̃}}{R}}_{\mathrm{c},n}^{\prime \prime }\mathrm{\Im }\left\{y_n\left(k\right)\right\}\\ & =y\left[k\right].\end{array}\end{array}$$

Here K_r,n are real residues according to Eq. (84). K_c,n and ${\stackrel{\mathrm{̃}}{R}}_{\mathrm{c},n}$ denote complex residues

$$\begin{array}{}\text{(94)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{c},n}=K_{\mathrm{c},n}^{\prime }+j{K}_{\mathrm{c},n}^{\prime \prime }\text{(95)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{R}}_{\mathrm{c},n}={\stackrel{\mathrm{̃}}{R}}_{\mathrm{c},n}^{\prime }+j{\stackrel{\mathrm{̃}}{R}}_{\mathrm{c},n}^{\prime \prime }.\end{array}$$

The pole relocation process reduces to an eigenvalue problem (Gustavsen and Semlyen, 1999). The numerically determined eigenvalues have to fulfill the following requirements:

1. The poles have to be stable, which means that the real part has to be negative. If a pole has a positive real part the stability condition

   $$\begin{array}{}\text{(96)}& p_{\mathrm{c},n}=-p_{\mathrm{c},n}^{\prime }+j{p}_{\mathrm{c},n}^{\prime \prime },p_{\mathrm{c},n}^{\prime }>\mathrm{0},\end{array}$$

   shifts the pole to the left of the imaginary axis.

2. Poles outside of the frequency range (Eq. 99) are not valid poles and are discarded.

3. Each pole has to fulfill

   $$\begin{array}{}\text{(97)}& \underset{\left|s-p_n\right|=\mathit{\epsilon }}{\oint }H\left(s\right)\mathrm{d}s>\mathit{\delta },\mathit{\delta }\le {\mathrm{10}}^{-\mathrm{13}}.\end{array}$$

Equation (97) makes use of Eq. (75), where a linearization of the integral leads to a simplified evaluation in the complex s-plane. The problem simplifies further, when a linearized circle with small radius around p_n is used as integration domain. This test has been carried out in the AAA-Fit algorithm to detect and discard spurious poles (Nakatsukasa et al., 2018). Grivet-Talocia and Bandinu (2006) investigate further validity conditions also for data that includes noise.

Figure 28 Selection of a curve segment and resulting segment that is approximated with TD-VF.

Figure 29 Approximated function $\stackrel{\mathrm{̃}}{y}\left(t\right)$ in red and measured response y(t) in black. Below is the absolute error.

Both data sets consist of M=5000 sample points, with Δt=4 ns. A direct fit of the whole data set did not lead to an acceptable result and should be avoided. It is more convenient to partition the data and fit consecutive parts. An algorithm has been developed that locates minima in the data, which are used as partition boundaries. The overall error is reduced when the approximation is based on a delayed transfer-function (Charest et al., 2010). In Eq. (82) h(t) becomes

$$\begin{array}{}\text{(98)}& h\left(t\right)\approx \sum _{i=\mathrm{0}}^{N_{\mathrm{s}}}{h}_i\left(t-{\mathit{\tau }}_i\right)\left(\mathbf{1}\left(t-{\mathit{\tau }}_i\right)-\mathbf{1}\left(t-{\mathit{\tau }}_{i+\mathrm{1}}\right)\right),\end{array}$$

where each h_i(t−τ_i) corresponds to Eq. (85). N_s is the number of segments bounded by [τ_i, τ_i+1], i=1, 2, …, N_s.

In Fig. 28 is the segment selection process shown. A specific segment is chosen by multiplying the data set with a rectangle, the resulting curve segment is fitted with TD-VF. The start poles are chosen between

$$\begin{array}{}\text{(99)}& {\mathit{\omega }}_{\min }={\displaystyle \frac{\mathit{\pi }}{{\mathit{\tau }}_{i+\mathrm{1}}-{\mathit{\tau }}_i}}{\mathit{\omega }}_{\max }={\displaystyle \frac{\mathit{\pi }}{\mathrm{\Delta }t}},\end{array}$$

where τ_i+1 and τ_i are the boundaries for the segment. Since the curve segment is the product of the original data and a rectangle function, a near infinite frequency width would be contained in this signal section. Frequencies that exceed the band limit, which is fixed due to the signal preprocessing (Sect. 4.2), are excluded and the fitting procedure is repeated until a good curve approximation is obtained or until the iteration limit is reached. An approximation is shown in Fig. 29 together with the deviation of the smoothed data set. The RMS is 0.138 mV with N_r=179 and N_c=560 poles. A second example is depicted in Fig. 30 with an RMS of 0.163 mV.

Figure 30 Piece-wise approximated $\stackrel{\mathrm{̃}}{y}\left(t\right)$ and absolute error.

Both approximations fit well to the original curve. A majority of the deviations appear at segment boundaries. Charest et al. (2010) make use of a wavelet decomposition for the response signal, which allows for an algorithmic search for ideal segmentation points. This could not be applied yet and is left for future work.

7.4 Outlook

The TD-VF has been applied on a single measured signal. The signal belongs to a fixed position within the RC-chamber where the stirrer is at rest. The reconstructed signal $\stackrel{\mathrm{̃}}{y}\left(t\right)$ is based on a transfer-function which is defined through residues and poles according to Eqs. (98) and (85). This transfer-function is not the transfer-function of the RC-chamber rather of the whole system consisting of cables, antennas and the chamber. Knowledge of the transfer-functions of the sub-systems – cables and antennas – give a more detailed description of the chamber, like its eigenmodes.

A multi-resolution in frequency and time gives the possibility to computationally partition the response signals. An automated TD-VF evaluation of signals, each measured at a fixed but different stirrer position, could then construct the transfer-functions systematically.

8 Conclusion

A system theoretical approach has been presented to explain the transient behaviour of a reverberation chamber excited with a pulsed signal. A spectral analysis on a pulsed signal and on the chamber response has been employed.

The spectrum of the chamber response shows that resonant frequencies in the neighbourhood of the carrier frequency are excited. Chamber responses of different stirrer positions provide nearly identical spectra. Differences between the chamber responses are more visible in the time domain.

An algorithm for the extraction of defined parameters was developed. Distinction of the individual responses is automatically possible on the basis of these parameters. The parameters include, for example, the duration and amplitude of the steady state as well as the time constants for the rise and decay. The results of the algorithm constitute the basis for further investigations, in particular the identification of the steady state is essential.

For further analysis of the chamber the data were approximated by a continuous function. This approximation is based on a superposition of damped oscillations for which the Prony interpolation method was used. The practical implementation is realised by a piecewise approximation of the data whose segmentation requires the knowledge of the steady state.

A similar approach allows to construct the transfer function of the overall RC system from the measurements with known excitation. This method is called Time-Domain Vector-Fit. The transfer function also represents an oscillating system. Similar to the Prony Method, the data set is segmented and approximated piecewise.

However, the approximations are based on linear time-invariant systems, i.e. the stirrer position must be kept constant for a measurement. Further investigations include the simulation of the chamber behaviour based on the approximated transfer functions as a function of the stirrer position. The transfer function of the chamber can be improved by knowing the transfer functions of the individual system components such as antennas, cables and connectors.

Furthermore, a comparison between the two interpolation methods is particularly useful with regard to the parameters obtained, such as the pole points.

Appendix A: Start values for time-constant estimates

To estimate the time constants for the Rise and the Decay phase from a damped oscillation (Eq. 58) start values for the unknown parameters have to be determined first. This is accomplished if a reduced model of a simple exponential process (Eq. 59) is used

$$\begin{array}{}\text{(A1)}& y\left(t\right)=\stackrel{\mathrm{̃}}{A}+\stackrel{\mathrm{̃}}{B}{e}^{-\stackrel{\mathrm{̃}}{T}\left(t-t_{\mathrm{o}}\right)},\stackrel{\mathrm{̃}}{T}>\mathrm{0}.\end{array}$$

One of the searched for quantities $\stackrel{\mathrm{̃}}{T}$ still appears in the argument of the exponential function, which would again require a start value for $\stackrel{\mathrm{̃}}{T}$. This can be avoided if Eq. (A1) is rewritten as an integral equation. As a consequence appears the former nonlinear quantity as a linear unknown and the application of classic least-squares is possible (Jacquelin, 2014). Integration of Eq. (A1) with regard to t leads to

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}I\left(y\right(t\left)\right)& :=\underset{t_{\mathrm{0}}}{\overset{t}{\int }}y\left(t\right)\mathrm{d}t\\ & =\stackrel{\mathrm{̃}}{A}\left(t-t_{\mathrm{0}}\right)-{\displaystyle \frac{\stackrel{\mathrm{̃}}{B}}{\stackrel{\mathrm{̃}}{T}}}{e}^{-\stackrel{\mathrm{̃}}{T}\left(t-t_{\mathrm{o}}\right)}+{\displaystyle \frac{\stackrel{\mathrm{̃}}{B}}{\stackrel{\mathrm{̃}}{T}}}{e}^{-\stackrel{\mathrm{̃}}{T}\left(t_{\mathrm{0}}-t_{\mathrm{o}}\right)}.\end{array}\end{array}$$

The exponential term is expressed with Eq. (A1) and the integral equation

$$\begin{array}{}\text{(A3)}& y\left(t\right)-y_{\mathrm{0}}=\stackrel{\mathrm{̃}}{C}\left(t-t_{\mathrm{0}}\right)-\stackrel{\mathrm{̃}}{T}I\left(y\right(t\left)\right),\end{array}$$

is established, where

$$\begin{array}{}\text{(A4)}& {\displaystyle }& {\displaystyle }{y}_{\mathrm{0}}=y\left(t_{\mathrm{0}}\right)=\stackrel{\mathrm{̃}}{A}+\stackrel{\mathrm{̃}}{B}{e}^{-\stackrel{\mathrm{̃}}{T}\left(t_{\mathrm{0}}-t_{\mathrm{o}}\right)}\text{(A5)}& {\displaystyle }& {\displaystyle }{y}_k=y\left(t_k\right)\text{(A6)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{C}:=\stackrel{\mathrm{̃}}{T}\stackrel{\mathrm{̃}}{A}\end{array}$$

and t₀ is the starting time value. The sum of the squared error for each time step is

$$\begin{array}{}\text{(A7)}& \sum _{k=\mathrm{0}}^{M-\mathrm{1}}{\mathit{\epsilon }}_k^{\mathrm{2}}=\sum _{k=\mathrm{0}}^{M-\mathrm{1}}{\left(\stackrel{\mathrm{̃}}{C}\left(t_k-t_{\mathrm{0}}\right)-\stackrel{\mathrm{̃}}{T}I\left(y_k\right)-\left(y_k-y_{\mathrm{0}}\right)\right)}^{\mathrm{2}}.\end{array}$$

Deriving Eq. (A7) with regard to $\stackrel{\mathrm{̃}}{C}$ and $\stackrel{\mathrm{̃}}{T}$ results in

$$\begin{array}{l}{\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial {\mathit{\epsilon }}^{\mathrm{2}}}{\partial \stackrel{\mathrm{̃}}{C}}}& =\mathrm{2}\sum _{k=\mathrm{0}}^{M-\mathrm{1}}\left(t_k-t_{\mathrm{0}}\right)\left(\stackrel{\mathrm{̃}}{C}\left(t_k-t_{\mathrm{0}}\right)\right.\\ & \left.\mathrm{2}\sum _{k=\mathrm{0}}^{M-\mathrm{1}}\left(t_k-t_{\mathrm{0}}\right)-\stackrel{\mathrm{̃}}{T}{I}_k-\left(y_k-y_{\mathrm{0}}\right)\right)\stackrel{!}{=}\mathrm{0}\end{array}\\ {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial {\mathit{\epsilon }}^{\mathrm{2}}}{\partial \stackrel{\mathrm{̃}}{T}}}& =\mathrm{2}\sum _{k=\mathrm{0}}^{M-\mathrm{1}}\left(-I_k\right)\left(\stackrel{\mathrm{̃}}{C}\left(t_k-t_{\mathrm{0}}\right)\right.\\ & \left.\mathrm{2}\sum _{k=\mathrm{0}}^{M-\mathrm{1}}\left(-I_k\right)-\stackrel{\mathrm{̃}}{T}{I}_k-\left(y_k-y_{\mathrm{0}}\right)\right)\stackrel{!}{=}\mathrm{0}\end{array}\end{array}$$

which is rewritten as

$$\begin{array}{}\text{(A8)}& \left(\begin{array}{cc}\sum {\mathit{\tau }}_k^{\mathrm{2}}& -\sum {\mathit{\tau }}_k{I}_k\\ -\sum {\mathit{\tau }}_k{I}_k& \sum I_k^{\mathrm{2}}\end{array}\right)\left(\begin{array}{c}\stackrel{\mathrm{̃}}{C}\\ \stackrel{\mathrm{̃}}{T}\end{array}\right)=\left(\begin{array}{c}\sum {\mathit{\tau }}_k{\mathit{\gamma }}_k\\ -\sum I_k{\mathit{\gamma }}_k\end{array}\right),\end{array}$$

with

$$\begin{array}{}\text{(A9)}& {\displaystyle }& {\displaystyle }\sum :=\sum _{k=\mathrm{0}}^{M_{\mathrm{1}}}\text{(A10)}& {\displaystyle }& {\displaystyle }{I}_k=I\left(y_k\right)\text{(A11)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_k=t_k-t_{\mathrm{0}}\text{(A12)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_k=y_k-y_{\mathrm{0}}.\end{array}$$

The quantity $\stackrel{\mathrm{̃}}{C}$ is discarded in favor of a second minimization based on the determined value of $\stackrel{\mathrm{̃}}{T}$. Simple least squares for $\stackrel{\mathrm{̃}}{A}$ and $\stackrel{\mathrm{̃}}{B}$ on the basis of Eq. (A1) leads to the desired first order approximation.

Table A1 Properties.