---
abstract: 'For many years, the Simplified Refined Instrumental Variable method for Continuous-time systems (SRIVC) has been widely used for identification. The intersample behaviour of the input plays an important role in this method, and it has been shown recently that the SRIVC estimator is not consistent if an incorrect assumption on the intersample behaviour is considered. In this paper, we present an extension of the SRIVC algorithm that is able to deal with input signals that cannot be interpolated exactly through hold reconstructions. The proposed estimator is generically consistent for any input reconstructed through zero or first-order-hold devices, and we show that it is generically consistent for continuous-time multisine inputs as well. The statistical performance of the proposed estimator is compared to the standard SRIVC estimator through extensive simulations.'
address:
- 'Division of Decision and Control Systems, KTH Royal Institute of Technology, 10044 Stockholm, Sweden'
- 'School of Electrical Engineering and Computing, University of Newcastle, Callaghan, 2308 NSW, Australia'
author:
- 'Rodrigo A. González'
- 'Cristian R. Rojas'
- Siqi Pan
- 'James S. Welsh'
bibliography:
- 'References.bib'
title: 'Consistent identification of continuous-time systems under multisine input signal excitation'
---

,

,

,

System identification; Continuous-time systems; Instrumental variable method; Consistency

Introduction
============

System identification involves using measured input and output data for building mathematical models that characterise a system’s behaviour. Different approaches to system identification have been developed depending on whether a discrete-time (DT) or continuous-time (CT) model is needed. Continuous-time system identification has applications in many areas of science and engineering such as economics, biology, physics and control, with comprehensive literature written on the subject [@rao2006identification; @Garnier2008book; @young2012recursive]. Although the system identification community has focused mainly in DT setups, as it has been investigated during a predominantly digital era, there are many reasons why CT system identification has had a renewed interest during the last decades [@garnier2015direct]. For example, model coefficients are directly linked to physical parameters, and more parsimonious models can be obtained as knowledge of the relative degree of the CT system can be accommodated. Also, contrary to the DT system identification using the forward shift operator, irregular and fast sampling can be easily handled, since the associated parameters remain invariant with respect to the varying sampling period and the model poles do not become statistically ill-defined as the sampling period decays to zero.

One of the main difficulties in CT system identification is the treatment of time derivatives. Since the goal is to obtain an estimate of a CT system, knowledge of the derivatives of the input and output are, either explicitly or implicitly, required. However, these derivatives are not exactly computable when only sampled input-output data is obtained. To overcome this problem, many algorithms have been suggested (see, e.g., [@sinha1991identification; @rao2006identification] and the references therein). One of the most popular algorithms is the Simplified Refined Instrumental Variable method for Continuous-time systems (SRIVC), which was first presented in [@young1980refined]. This method has been suggested for general use due to its robustness and accuracy in practical applications [@garnier2014advantages]. Many further extensions of this method also exist in the literature, for example, to handle non-uniformly sampled data [@huselstein2002approach] or multi-input systems [@garnier2007optimal]. Extensions to output error (OE) and Box-Jenkins (BJ) models [@chen2013refined], unification of DT and CT transfer function estimation [@young2015refined], and a comprehensive consistency analysis [@pan2020consistency] have also been presented.

The SRIVC algorithm uses interpolation of the input and output data in order to compute filtered regressor and instrument vectors in an iterative estimation procedure. This reconstruction of the CT input and output signals is usually implemented through simple interpolation schemes like zero-order hold (ZOH) or first-order hold (FOH) devices, independently of the nature of the true signals [@garnier2004time]. For inputs that can be described exactly with these reconstruction schemes, the SRIVC estimator has recently been shown to be generically consistent [@pan2020consistency]. However, when the intersample behaviour assumption on the model input does not match that of the system input, continuous-time estimation methods can deliver large estimation errors if the sampling period is large [@schoukens1994identification], and in particular, the SRIVC estimator is known to be generically inconsistent in this case. Important input signals for identification that cannot be described by holds are band limited signals such as multisines. These signals are advantageous due to their flexibility regarding power spectrum design, time domain averaging possibilities, simplification of the model validation step and finite sample estimation performance [@schoukens1994advantages]. For these input signals, the complete CT input signal is known to the practitioner, but the SRIVC procedure only performs simple interpolations of the input, which impact its consistency regardless of the sampling period.

In summary, in this paper,

-   we present a refinement of the SRIVC method that is shown to yield generic consistency of the estimated model parameters for CT multisine input signal excitations;

-   we prove that, given knowledge of the CT multisine input signal and measured output samples, the exact computation of the input regressors is necessary and sufficient for a generically consistent estimate of the CT system;

-   we propose a computationally efficient algorithm for computing the regressors under the multisine case, and introduce an extension of the SRIVC algorithm for arbitrary CT inputs, which are not necessarily constructed through hold devices; and

-   we exemplify the consistency of the proposed estimator through extensive Monte Carlo simulations.

The remainder of this paper is organised as follows. The identification problem is formulated in Section \[sec2\]. Section \[sec3\] provides a description of the SRIVC estimator and its consistency properties. The proposed SRIVC-type method is presented and analysed in Section \[sec4\], and Section \[sec5\] illustrates this method with extensive numerical examples. Finally, conclusions are drawn in Section \[sec6\].

Problem formulation {#sec2}
===================

Consider a linear and time-invariant (LTI), causal, stable, proper, single-input single-output, CT system $$\label{eq1}
x(t) = \frac{B^*(p)}{A^*(p)} u(t), \notag$$ where $p$ is the Heaviside operator, i.e., $p g(t) := \text{d}g(t)/\text{d}t$, and the numerator and denominator polynomials are coprime and given by $$\begin{aligned}
B^*(p) &= b_0^* p^{m^*} + b_1^* p^{{m^*}-1} + \dots + b_{m^*}^*, \notag \\
A^*(p) &= a_1^* p^{n^*} + a_2^* p^{{n^*}-1} + \dots + a_{n^*}^* p+1. \notag\end{aligned}$$ Suppose that the CT input $u(t)$ is known from $t=t_1$ to $t=t_N$, and that $N$ noisy measurements of the output $x(t)$ are obtained at the instants $\{t_k\}_{k=1}^N$. In other words, the output observations are given by $$\label{outputdescription}
y(t_k) =x(t_k) + v(t_k), \quad k=1,\dots,N,$$ where it is assumed that the sampled noise sequence $\{v(t_k)\}$ can be described as a zero-mean and finite variance random process. Due to the nature of the sampled signals and the difficulty of computing the time-derivative of CT white noise, which does not have finite variance [@aastrom2012introduction], we only consider DT noise in this paper.

To identify the system, we propose the model structure $$G(p)= \frac{b_0 p^{m} + b_1 p^{m-1} + \dots + b_m}{a_1 p^{n} + a_2 p^{n-1} + \dots + a_n p+1},  \notag$$ where the parameter vector $$\theta := \begin{bmatrix}
a_1, & a_2, &\dots, & a_n, & b_0, & b_1, &\dots, & b_m 
\end{bmatrix}^\top \notag$$ needs to be estimated. The goal is to obtain an accurate model of the CT system $G^*(p):=B^*(p)/A^*(p)$ given the knowledge of $N$ samples of the output measurements and the CT input signal. Note that in this framework the input signal is not limited to hold reconstructions. Hence, the description includes the standard framework where $u(t)$ is assumed to be obtained through a ZOH or FOH and extends to more general inputs, such as multisines and band limited signals [@schoukens1994advantages].

The identification of the system $G^*(p)$ can be done by obtaining the data points $\{u(t_k),y(t_k)\}$ and applying a method for CT system identification, such as in [@chen2013refined], or as in [@young1980refined; @maruta2013projection; @gonzalez2018asymptotically] for regular sampling schemes. In most of these algorithms, however, the hold reconstructions of the input and output are assumed, and they are independent of the exact nature of the signals. In this work, we show that the knowledge of the exact intersample behaviour of the input can provide further insights for a better design of the identification procedure.

The Simplified Refined Instrumental Variable method for Continuous-time systems (SRIVC) {#sec3}
=======================================================================================

The SRIVC estimator is an adaptive instrumental variable algorithm where parameter-dependent CT filters are updated iteratively. In each step, the instruments are computed using the parameter estimate obtained in the previous iteration until the model parameters have converged. The iterative procedure of the SRIVC algorithm is designed so that the sum of squares of the residuals (also called the generalised equation errors or GEEs) $\varepsilon(t_k)$, is minimised. The residuals are written as $$\begin{aligned}
    \varepsilon(t_k) :&= y(t_k)-\frac{B(p)}{A(p)}u(t_k) \notag \\
    \label{gee}
    &= A(p)y_f(t_k)-B(p)u_f(t_k),  
    \end{aligned}$$ where $$y_f(t_k) = \frac{1}{A(p)}y(t_k) \textnormal{, and } u_f(t_k) = \frac{1}{A(p)}u(t_k). \label{filt1}$$ Note that in and we have adopted a mixed notation of CT operators and DT data. Since this dichotomy is repeatedly encountered in this paper, we formalise it in the following remark.

**.** \[remark1\] In this paper, $G(p)x(t_k)$ means that the DT signal $x(t_k)$ is interpolated in some manner, e.g., using a ZOH or FOH, and the resultant output through the CT filter $G(p)$ is sampled at $t=t_k$. On the other hand, $\{G(p)x(t)\}_{t_k}$ (or $[G(p)x(t)]_{t_k}$ in the vector-valued case) means that the CT signal $x(t)$ is filtered through $G(p)$, and later sampled at $t=t_k$.

The SRIVC method is described in Algorithm \[algorithm1\], where we denote $\varphi_f(t_k)$ as the filtered regressor vector, $\hat{\varphi}_f(t_k)$ as the filtered instrument vector, and $y_f(t_k)$ as the filtered output. Note that line 8 of Algorithm \[algorithm1\] requires the DT signals to be prefiltered by CT transfer functions. This is usually done by assuming a ZOH or FOH reconstruction for the input and output signals and then simulating the response by using, for example, the `lsim` command in MATLAB. Although this approach has provided a quick procedure to compute the filtered regressor and instrument vectors, it is prone to approximation errors that can jeopardise the statistical properties of the method.

\[algorithm1\]

Input: $\{(u(t_k),y(t_k))\}_{k=1}^N$, model order $(n,m)$, initial vector estimate $\theta_1$, tolerance $\epsilon$ and maximum number of iterations $M$ Using $\theta_1$, form the estimated system polynomials $A_1(p)$ and $B_1(p)$ $j\gets 1$, $\textnormal{flag}\gets 1$ Reflect the unstable poles of $1/A_j(s)$ into the stable region of the complex $s$-plane Prefilter the (DT) input $\{u(t_k)\}_{k=1}^N$ and output $\{y(t_k)\}_{k=1}^N$ by CT filters to form

\_f(t\_k) & \^,&\
\_f(t\_k) & \^, &\
\[filteredoutput\] y\_f(t\_k) &y(t\_k)&

Compute the parameter estimate $$\label{iterations}
        \theta_{j+1} \gets \hspace{-0.05cm}\left[\sum_{k=1}^N  \hat{\varphi}_f(t_k) \varphi_f^\top(t_k)\right]^{-1}\hspace{-0.1cm}\left[ \sum_{k=1}^N \hat{\varphi}_f(t_k) y_f(t_k)  \right]$$ $\textnormal{flag} \gets 0$ $j \gets j+1$ Output: $\theta_{j}$ and its associated model $B_j(p)/A_j(p)$.

\[algorithm1\]

**.** In the SRIVC method, the user has several choices regarding the intersample behaviour assumptions. In particular, the intersample behaviour of the input in both and can be chosen, as well as the reconstruction of the output signal for the filtering steps in  and . Usually the output is selected to have a FOH behaviour, since it is argued that it typically gives rise to a satisfactory approximation if the sampling period is small [@chen2017robust].

Consistency Analysis of the SRIVC estimator
-------------------------------------------

Previous works [@young2008refined; @young2015refined] have suggested that the SRIVC estimator uses the optimal instrumental variable terms, and that it minimises the prediction error and maximises the likelihood function, but they lack rigorous theoretical analysis regarding the influence of the interpolation of the input and output for the prefiltering step. Only recently [@pan2020consistency] has the intersample behaviour of the signals been taken into account for the consistency analysis. In [@pan2020consistency Theorem 1], the generic consistency of the SRIVC estimator was proven for inputs that can be exactly interpolated by FOH or ZOH devices. More precisely, under mild assumptions regarding the sampling period and persistence of excitation of the input, the following statements are true for an input that is exactly reconstructible with FOH or ZOH interpolation:

1.  The matrix $E\{\hat{\varphi}(t_k)\varphi^\top(t_k)\}$ is generically non-singular.[^1]

2.  The true parameter $\theta^*$ is the unique converging point.

3.  As the sample size $N$ approaches infinity, $\theta_{j+1}$ in  converges to $\theta^*$ for $j\geq 1$.

Also, the effect of choosing a different intersample behaviour than that of the system input was also analysed in [@pan2020consistency Corollary 3]. In the following, we say that a *correct intersample behaviour in the input signal is assumed* whenever the intersample behaviour of such signal in the SRIVC algorithm matches that of the system input. In [@pan2020consistency] it was shown that the SRIVC estimator

1.  remains generically consistent if an incorrect assumption on the intersample behaviour is used for generating the filtered signals in the instrument vector $\hat{\varphi}_f(t_k)$, and

2.  is generically inconsistent if an incorrect assumption on the intersample behaviour is used for filtering the input signal in the regressor vector $\varphi_f(t_k)$.

This result indicates that the intersample behaviour of the input signal needs to be correctly taken into account for the consistency of the SRIVC estimator. In particular, it implies that if the system input is a signal that is not produced by a hold mechanism, the estimator will be generically inconsistent. This argument holds regardless of whether the additive noise $v(t_k)$ is white or coloured.

Consistent SRIVC-type method {#sec4}
============================

As mentioned in the previous section, a correct assumption of the intersample behaviour of the input (ZOH of FOH) in the regressor vector $\varphi(t_k)$ guarantees generic consistency under mild conditions. The extension of this principle constitutes our main contribution. In this work, we propose an extension of the SRIVC method that computes the filtered regressors exactly for a wide class of input signals whose intersample behaviour is known and prove its consistency for multisine input excitations.

The generalised equation error for the proposed approach is $$\label{gee2}
    \varepsilon(t_k) = A(p)y_f(t_k)-\{B(p)u_f(t)\}_{t_k},$$ where $u_f(t) = \frac{1}{A(p)}u(t)$. In , the predicted output measurement is explicitly calculated by first computing the underlying CT signal, and later evaluating it at $t=t_k$. The proposed estimator follows the procedure described in Algorithm \[algorithm1\], but the filtered regressor and instrument vectors in Equations and now become

&\_f(t\_k) = \^,&

and

&\_f(t\_k) = \^\_[t\_k]{}.&

Note that the $t_k$ in follows the notation in Remark \[remark1\].

**.** The proposed estimator is an extension of the standard SRIVC estimator that uses the complete CT input signal for identification. For input signals that are reconstructed exactly through a ZOH or FOH, this estimator is equivalent to the SRIVC estimator. Thus, the SRIVC-type estimator with prefiltering stage given by  and is generically consistent under the same assumptions as in [@pan2020consistency].

In order to further analyse the asymptotic properties of the proposed estimator, we first study its consistency for multisine inputs. Later, an extension for arbitrary input excitations is presented.

Multisine input signal
----------------------

We consider multisine input signals of the form $$\label{inputmultisine}
u(t)=\alpha_0+\sum_{l=1}^{m_u} \alpha_{l} \cos(\omega_l t + \psi_l),$$ where $m_u,\{\alpha_l\}_{l=0}^{m_u},\{\omega_l\}_{l=0}^{m_u}$ and $\{\psi_l\}_{l=0}^{m_u}$ are input parameters. The frequencies $\omega_l$ are assumed to be positive and distinct, and without loss of generality we assume that the weights $\alpha_l$ are positive as well. It is well known that the output in stationary regime of an asymptotically stable LTI filter $H(s)$ when $u(t)$ is applied is also a multisine, given by $$\label{eval3}
y(t)=H(0)\alpha_0+\sum_{l=1}^{m_u} \alpha_l|H(i\omega_l)| \cos(\omega_l t + \psi_l+\angle H(i\omega_l)).$$ This property of LTI systems provides a natural way to obtain exact values for the signal evaluations in and , and it is of low computational cost, since the prefiltering is directly obtained by evaluating with the corresponding filter. Another advantage of this approach is that it extends naturally to non-uniformly sampled data. For such type of sampling, the proposed method is not as computationally intensive as the standard SRIVC method, since the algorithm only obtains approximations of the filtered output $p^i A_j^{-1}(p)y(t_k), i=0,\dots,n$, instead of computing approximations of the filtered values of both $u(t_k)$ and $y(t_k)$. The filtered output computations can be carried out by, e.g., an adaptive Runge-Kutta method (as in [@chen2013refined]), or by any oversampling technique with intersample behaviour assumptions.

We now prove the consistency of the proposed estimator for the multisine input. The assumptions we use during the analysis are the following:

1.  \[assumption1\] The true system $B^*(p)/A^*(p)$ is proper ($n^* \geq m^*$) and asymptotically stable with $A^*(p)$ and $B^*(p)$ being coprime.

2.  \[assumption2\] The disturbance sequence $\{v(t_k)\}$ is a zero-mean stationary random process.

3.  \[assumption3\] The number of sinusoids of the input, $m_u$, satisfies $m_u\geq (n+m)/2$, and the input offset, $\alpha_0$, is different from zero.

4.  \[assumption4\] All the zeros of $A_j(p)$ have strictly negative real parts, $n\geq m$, with $A_j(p)$ and $B_j(p)$ being coprime.

5.  \[assumption5\] The degrees of the polynomials in the model satisfy $\min(n-n^*,m-m^*)=0$.

Assumptions \[assumption1\] and \[assumption2\] are standard. The condition in Assumption \[assumption3\] is a persistence of excitation requirement, where $\alpha_0 \neq 0$ is set only for simplicity in our derivations and can be removed.[^2] Given that the poles of unstable models are reflected in line 6 of Algorithm \[algorithm1\], Assumption \[assumption4\] is met in practice. Assumption \[assumption5\] takes into account the model structure, as it ensures a unique solution of the model parameters to be obtained.

Since deterministic inputs will be considered in conjunction with stochastic noise processes, our analysis uses the standard definition of expectation for quasi-stationary signals [@ljung1998system pp. 34], which is $$\overline{E}\{g(t)\} := \lim_{N\to \infty} \frac{1}{N} \sum_{t=1}^N E\{g(t)\}. \notag$$

**.** \[theorem1\] Consider the SRIVC-type estimator with a fixed sampling period $h$ and filtered regressor and instrument vectors given by and respectively, and suppose that Assumptions \[assumption1\] to \[assumption5\] hold. Then, the following statements are true:

1.  There exists a maximum sampling period $h^*>0$ such that, if $h\leq h^*$, the matrix $\overline{E}\{\hat{\varphi}_f(t_k)\varphi_f^\top(t_k)\}$ is generically non-singular.

2.  If $h\leq h^*$, the true parameter $\theta^*$ is the unique converging point.

3.  As the sample size $N$ approaches infinity, $\theta_{j+1}$ converges to $\theta^*$ for $j\geq 1$.

[Proof.]{}

*Proof of Statement 1*. By substituting $$y(t_k) = \left\{\frac{B^*(p)}{A^*(p)}u(t)\right\}_{t_k} + v(t_k) \notag$$ into , we find that $\varphi_f(t_k)=\varphi_{f1}(t_k)-V(t_k)$, where $$\begin{aligned}
    \hspace{-0.21cm}\varphi_{f1}(t_k)\hspace{-0.05cm}:= \hspace{-0.12cm}\Bigg[\dfrac{-p^n}{A_j(p)}&\left\{ \dfrac{B^*(p)}{A^*(p)}u(t) \right\}_{t_k} \hspace{-0.2cm},\hspace{0.1cm} \dots, \hspace{0.1cm}\dfrac{-p}{A_j(p)}\left\{ \dfrac{B^*(p)}{A^*(p)}u(t) \right\}_{t_k} \notag \\
    &\hspace{0.1cm} \left\{\dfrac{p^m u(t)}{A_j(p)} \right\}_{t_k}  \hspace{-0.1cm}, \hspace{0.1cm} \dots,\hspace{0.1cm} \left\{\dfrac{u(t)}{A_j(p)} \right\}_{t_k}
    \Bigg]^\top,
    \end{aligned}$$ and $$\label{Vtk}
V(t_k):= \begin{bmatrix}
\dfrac{p^n}{A_j(p)} v(t_k), & \hspace{0.05cm}\dots, & \hspace{0.1cm} \dfrac{p}{A_j(p)} v(t_k), & 0, & \hspace{0.05cm}\dots, & 0  
\end{bmatrix}^\top.$$ On the other hand, we also have $$\hat{\varphi}_f(t_k)=S(-B_j,A_j) \left[ \dfrac{U_{du}(t)}{A_j^2(p)} \right]_{t_k}, \notag$$ where $$\label{Udu}
    U_{du}(t) := \begin{bmatrix}
    \dfrac{\textnormal{d}^{n+m}}{\textnormal{d}t^{n+m}}u(t), & \dfrac{\textnormal{d}^{n+m-1}}{\textnormal{d}t^{n+m-1}}u(t), & \hspace{0.1cm}\dots, & \hspace{0.1cm} u(t)
    \end{bmatrix}^\top \hspace{-0.1cm},$$ and $S(-B_j,A_j)$ is the Sylvester matrix associated with the polynomials $-B_j(p)$ and $A_j(p)$, whose non-singularity follows from the same analysis done in [@pan2020consistency], where Assumption \[assumption4\] is used. With this, we compute $$\begin{aligned}
    \overline{E}\{\hat{\varphi}_f(t_k)\varphi_f^\top(t_k)\} &= S(-B_j,A_j)\underbrace{\overline{E}\left\{\left[ \dfrac{U_{du}(t)}{A_j^2(p)} \right]_{t_k}\varphi_{f1}^\top(t_k) \right\}}_{=:\Phi} \notag \\
    \label{Psidef}
    &\hspace{-1.5cm}-S(-B_j,A_j)\underbrace{\overline{E}\left\{\left[ \dfrac{U_{du}(t)}{A_j^2(p)} \right]_{t_k}V^\top(t_k)\right\}}_{=:\Psi}.
    \end{aligned}$$ Thus, for showing that $\overline{E}\{\hat{\varphi}_f(t_k)\varphi_f^\top(t_k)\}$ is generically non-singular for a small enough sampling period $h$, it is sufficient to show that $\Psi=0$ and $\Phi$ is generically non-singular for a small enough $h$. The difference between the analysis in [@pan2020consistency Theorem 1] and the proof in the current paper is that the signals of interest are hybrid in nature: some are evaluations of CT signals, whereas others are DT signals interpolated with a reconstruction device, such as a FOH.

The proof of $\Psi=0$ can be found in Lemma \[lemma1\] in the Appendix. Regarding the invertibility of $\Phi$, we will conveniently write $\varphi_{f1}(t_k)$ as $\varphi_{f2}(t_k)+\Delta(t_k)$, where $\Delta(t_k) \in \mathbb{R}^{n+m+1}$ has entries $$\label{deltai}
    \Delta_i(t_k)= \begin{cases}
    \left\{\frac{p^{n+1-i}}{A_j(p)}x(t)\right\}_{t_k}-\frac{p^{n+1-i}}{A_j(p)} x(t_k) &\hspace{-0.2cm},\hspace{0.07cm} i=1,\dots,n \\
    0 &\hspace{-2.2cm},\hspace{0.07cm} i=n+1,\dots,n+m+1,
    \end{cases}$$ and $$\begin{aligned}
    \varphi_{f2}(t_k) &= \Bigg[
    \dfrac{-p^nB^*(p)}{A_j(p) A^*(p)}u(t), \hspace{0.1cm}\dots, \hspace{0.1cm}\dfrac{-pB^*(p)}{A_j(p)A^*(p)}u(t), \notag \\
    &\hspace{2cm} \quad \dfrac{p^m}{A_j(p)}u(t),  \hspace{0.1cm}\dots, \hspace{0.1cm} \dfrac{1}{A_j(p)}u(t)
    \Bigg]^\top_{t_k} \notag \\
    &= S(-B^*,A^*) \left[ \frac{U_{du}(t)}{A_j(p)A^*(p)} \right]_{t_k}, \notag 
    \end{aligned}$$ with $S(-B^*,A^*)$ being the Sylvester matrix associated with the polynomials $-B^*(p)$ and $A^*(p)$, which is non-singular since $A^*(p)$ and $B^*(p)$ are coprime. Hence, we can write the expected value of interest as $$\begin{aligned}
    \overline{E}\{\hat{\varphi}_f(t_k)\varphi_f^\top(t_k)\} &= S(-B_j,A_j)\Phi_1 S^\top(-B^*,A^*) \notag \\
    \label{righthandside}
    &\hspace{-1.9cm}+S(-B_j,A_j)\overline{E}\left\{\left[ \dfrac{U_{du}(t)}{A_j^2(p)} \right]_{t_k}\Delta^\top(t_k)\right\},
    \end{aligned}$$ where $$\Phi_1 := \overline{E}\left\{ \left[\frac{U_{du}(t)}{A_j^2(p)}\right]_{t_k} \left[\frac{U_{du}(t)}{A_j(p)A^*(p)} \right]_{t_k}^\top \right\}. \notag$$ It is shown in Lemma \[lemma2\] that $\Phi_1$ is generically non-singular, which means that the first summand of the right hand side of is generically non-singular.

Finally, as $h$ tends to zero, the infinity norm of the difference between the direct evaluation of a CT signal and its interpolated counterpart also tends to zero. Thus, $\Delta(t_k)\to 0$ as $h\to 0$. This, together with the fact that (generic) non-singularity of a matrix is preserved under small-enough matrix perturbations [@Horn2012 Chap. 6], leads to the first statement of the theorem.

*Statement 2*. Suppose that $\bar{\theta}$ is a limiting point of the iteration in , where $\varphi_f(t_k)$ and $\hat{\varphi}_f(t_k)$ are defined as in  and respectively, and the corresponding polynomials are denoted by $\bar{A}(p)$ and $\bar{B}(p)$. These polynomials are coprime by Assumption \[assumption4\]. The ergodic lemmas in [@soderstrom1975ergodicity Lemma 3.1] and [@soderstrom1983instrumental Lemma A4.3] permit us to write the iteration equation , at the converging point and as $N$ tends to infinity, as $$\label{eqlimit}
    \overline{E}\{\hat{\varphi}_f(t_k,\bar{\theta})\varphi_f^\top(t_k,\bar{\theta})\}^{-1} \overline{E}\{\hat{\varphi}_f(t_k,\bar{\theta})\varepsilon(t_k,\bar{\theta})\}=0,$$ where $\varepsilon(t_k,\bar{\theta})$ is the GEE evaluated at the converging point. Since the matrix inverse in is assumed to be non-singular, the second expectation in must be zero, i.e., $$\label{eqcon2}
    \overline{E}\{\hat{\varphi}_f(t_k,\bar{\theta})\varepsilon(t_k,\bar{\theta})\}=0.$$

Let $\bar{A}(p)B^*(p)-\bar{B}(p)A^*(p)=h_0p^r+h_1 p^{r-1}+\dots + h_r$, where $r = \max(n+m^*,n^*+m)=n+m$. Then, the GEE in can be rearranged as $$\begin{aligned}
    \varepsilon(t_k,\bar{\theta})&=\left\{\frac{\bar{A}(p)B^*(p)-\bar{B}(p)A^*(p)}{\bar{A}(p)A^*(p)}u(t)\right\}_{t_k}+v(t_k) \notag \\
    &= \left\{\frac{U_{du}^\top(t)}{\bar{A}(p)A^*(p)}H \right\}_{t_k} + v(t_k), \notag
    \end{aligned}$$ where $H = \begin{bmatrix}
    h_0, & h_1, & \hspace{0.1cm}\dots, & h_{n+m}\end{bmatrix}^\top$. Now, note that the instrument vector $\hat{\varphi}_f(t_k)$ can be written as $$\hat{\varphi}_f(t_k) = S(-\bar{B},\bar{A})\left[\frac{U_{du}(t)}{\bar{A}^2(p)}\right]_{t_k}, \notag$$ where $S(-\bar{B},\bar{A})$ is a Sylvester matrix associated with the polynomials $\bar{B}(p)$ and $\bar{A}(p)$, which again is non-singular. So, we can express as $$\begin{aligned}
    &0 = S(-\bar{B},\bar{A}) \underbrace{\overline{E}\left\{\left[\frac{U_{du}(t)}{\bar{A}^2(p)}\right]_{t_k} \left[\frac{U_{du}(t)}{\bar{A}(p)A^*(p)}\right]_{t_k}^\top\right\}}_{:=\bar{\Phi}} H \notag \\
    \label{eqcon3}
    &\hspace{0.2cm}+ S(-\bar{B},\bar{A})\underbrace{\overline{E}\left\{ \left[\frac{U_{du}(t)}{\bar{A}^2(p)}\right]_{t_k}v(t_k) \right\}}_{:=\bar{\Psi}}.
    \end{aligned}$$ Following a similar approach as in Lemma \[lemma1\], we conclude that $\bar{\Psi}=0$, and by Lemma \[lemma2\], $\bar{\Phi}$ is generically non-singular. Thus, for to hold we need $H=0$, which implies that $$\begin{aligned}
    \bar{A}(p)B^*(p)&-\bar{B}(p)A^*(p)=0 \notag \\
    \frac{\bar{B}(p)}{\bar{A}(p)}&= \frac{B^*(p)}{A^*(p)}, \notag 
    \end{aligned}$$ i.e., $\theta^*$ is the unique limiting point.

*Statement 3*. The proof follows from the analysis made for proving Statement 3 of Theorem 1 in [@pan2020consistency].

Note that if the commonly used FOH (or ZOH) were chosen as the intersample behaviour of the signals when discretising the prefilters, the reconstruction of $u(t)$ would suffer from high frequency distortion, which usually leads to inaccuracies in the computation of $\varphi_f(t_k)$ and $\hat{\varphi}_f(t_k)$. As stated next, only an inaccurate computation of the regressor vector $\varphi_f(t_k)$ causes generic inconsistency of the proposed method under CT multisine input excitation.

**.** Assume that an incorrect intersample behaviour in the input signal is assumed, but nevertheless satisfies $G(p)u(t_k) = \{G(p)u(t)\}_{t_k}$ as $h\to 0$. The SRIVC-type estimator with filtered regressor and instrument vectors given by and respectively

1.  remains generically consistent if an incorrect assumption on the intersample behaviour is used for generating the filtered signals in the instrument vector $\hat{\varphi}_f(t_k)$, and

2.  is generically inconsistent if an incorrect assumption on the intersample behaviour is used for filtering the input signal in the regressor vector $\varphi_f(t_k)$.

[Proof.]{}

*Statement 1*: The result follows from the same logic as in the proof in [@pan2020consistency Corollary 3, Statement 1].

*Statement 2*: Statement 1 of Theorem \[theorem1\] still holds by following the same steps as before, but this time the vector $\Delta(t_k)$ in will also have non-zero elements in its bottom $m+1$ entries. Namely, the $i$-th component of $\Delta(t_k)$, with $i=n+1,\dots,n+m+1$, is now $$\Delta_i(t_k) = \frac{p^{m+n+1-i}}{A_j(p)}u(t_k)-\left\{\frac{p^{m+n+1-i}}{A_j(p)}u(t)\right\}_{t_k}, \notag$$ which still satisfies $\Delta_i(t_k)\to 0$ as $h\to 0$. Thus, Theorem \[theorem1\] is valid for this case as well. However, Statement 2 of Theorem \[theorem1\] does not yet hold. This fact follows from the same analysis done in the proof in [@pan2020consistency Corollary 3, Statement 2].

**.** A similar procedure to and could be proposed for the computation of $y_f(t_k)$ by exploiting the fact that the noiseless output also corresponds to a multisine (thus, a more adequate reconstruction scheme could be designed). However, Remark 5 of [@pan2020consistency] suggests that, as the number of iterations tends to infinity, the GEE at the converging point does not depend on the intersample behavior of the output. Thus, if the iterations converge, a more precise filtering of the output is not needed.

Extension to arbitrary input signals {#arbitrary}
------------------------------------

The previous method is naturally suited for multisine inputs due to the simplicity of the filtered outputs at stationary regime, since they are also multisine signals. By introducing a Delta transform description, the computations in and can be generalised for an arbitrary input signal with arbitrary accuracy. For this procedure, only regular sampling is considered, although extensions to irregular sampling are also possible.

The proposed algorithm for computing $\varphi_f(t_k)$ and $\hat{\varphi}_f(t_k)$ at the $j$-th iteration of the SRIVC-type method goes as follows:

1.  Given the sampling period $h$ of $y(t)$, (over)sample $u(t)$ with sampling period $\delta h$, where $\delta\ll 1$.

2.  From $\theta_j$, form the prefilters of interest, namely $p^i/A_j(p)$ and $B_j(p)p^l/A_j^2(p)$ for $i=0,\dots,m$; $l=1,\dots,n$.

3.  Compute the Delta equivalent [@middleton1990digital] of the prefilters, and calculate the response at instants $t_k$ of each filter to the fast-sampled version of $u(t)$ in the Delta domain.

The Delta domain description leads to an exact simulation of the underlying CT system when $\delta h\to 0$. Thus, the SRIVC-type method with prefilters computed as in steps 1, 2 and 3 above calculates the filtered regressor and instrument vectors accurately if the oversampling period $\delta h$ is chosen small as compared to the sampling period. Due to potentially high sampling rates, the use of the Delta operator is needed for ameliorating rounding errors and ill-conditioning problems regarding the sensitivity of the coefficients of the prefilters. Note that step 3 of the algorithm can be easily performed by using the Delta domain Toolbox in MATLAB.

**.** The precision of this procedure will depend on the over-sampling period $\delta h$. Via extensive simulations, we have found that sampling at least 100 times faster is usually enough to provide reliable estimates.

Simulation examples {#sec5}
===================

Via numerical simulations we compare the performance of the standard SRIVC method with the proposed SRIVC-type method, which is labeled SRIVC-c. For a multisine input, we examine the consistency of both methods for different regular sampling period and also for irregular sampling. We also study the consistency of each method under an arbitrary input excitation for different regular sampling periods. For the following tests, we consider the system $$\label{systemexample}
G^*(p) = \frac{1.25}{0.25p^2+0.7p+1},$$ where the parameters of interest are $a_1^*=0.25$, , and $b_0^*=1.25$. Regarding the implementation of the standard SRIVC method, we have used the `srivc` command from the CONTSID toolbox version 7.3 for MATLAB [@garnier2018contsid], under default initialisation and tolerance settings. It was set to estimate the best model among the correct model structure with a FOH as the intersample behaviour.

Multisine input: Regular sampling
---------------------------------

We first test if the algorithms provide consistent estimates of the parameter vector $[a_1, \hspace{0.1cm} a_2, \hspace{0.10cm} b_0]^\top$. The system in is excited with the CT input $$u(t) = \sin(0.714t)+ \sin(1.428t)+\sin(2.142t). \notag$$ The noiseless output is computed analytically by assuming that it corresponds to the output of the system at the stationary regime, i.e., $$x(t) = \sum_{k=1}^3|G^*(i\omega_k)|\sin(\omega_k t+\angle G^*(i\omega_k)), \notag$$ where $(\omega_1,\omega_2,\omega_3) = (0.714, 1.428,2.142)[\textnormal{rad}/\textnormal{s}]$. This output is sampled at $h = 0.3[\textnormal{s}]$ and is contaminated with additive noise, which is set as an i.i.d. Gaussian sequence of variance 0.1. Sixty different sample sizes are considered, ranging logarithmically from $N=100$ to $N=25500$, and 300 Monte Carlo runs are performed for each value of $N$.

Figures \[fig1\] and \[fig2\] show the sample means and sample mean square errors (MSEs) of each estimated parameter. The SRIVC-c estimator accurately identifies all parameters while the standard SRIVC method fails to recover the true parameter vector as $N$ increases. The MSEs for the SRIVC-c estimator decrease to zero as expected for consistency, while at least two out of the three estimated parameters given by the SRIVC method are biased, which empirically indicates that the SRIVC estimator is not consistent in this example.

Multisine input: Different sampling periods
-------------------------------------------

We now study the effect of the intersample behaviour on the SRIVC-type estimates. Under the same input and noise variance as the previous simulation, we test the performance of each algorithm for a fixed number of output measurements ($N=2000$) with different regular sampling periods. Since the rise time of the system is approximately 2 seconds, a good choice for the sampling period should be between 0.2 and 0.5 seconds according to the criterion suggested in [@astroem1984computer]. In order to cover fast, normal and slow sampling, we test with sampling periods $h=0.06,0.2$ and $0.6$\[s\].

The sample mean and mean square error of each parameter over $300$ Monte Carlo runs for each sampling period are shown in Table \[table1\]. On average, the SRIVC-c estimator delivers the true values of every parameter for all sampling periods in this study, whereas the SRIVC estimator only performs well when the sampling period is small. For $h=0.6[\textnormal{s}]$, the large sampling period exaggerates the interpolation error of the input signal in the standard SRIVC estimator, which severely degrades its performance. This is confirmed by the order of magnitude of difference in MSE of the parameters given by the two estimators.

\[table1\]

[|@m[0.92cm]{}|m[0.96cm]{}|@m[0.63cm]{}@|m[1.13cm]{}|m[1.13cm]{}|m[1.13cm]{}|]{} Method &

  ---------
   Param.
   (Value)
  ---------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

& Stats. & & $h=0.2$ & $h=0.6$\
& $a_1 (0.25)$ &

  ------
   Mean
   MSE
  ------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.253$
   $1.8\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.251$
   $1.1\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.248$
   $1.6\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

\
SRIVC & $a_2 (0.7)$ &

  ------
   Mean
   MSE
  ------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.697$
   $7.0\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.694$
   $9.3\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.668$
   $1.1\cdot 10^{-3}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

\
& $b_0 (1.25)$ &

  ------
   Mean
   MSE
  ------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $1.244$
   $1.6\cdot 10^{-4}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $1.251$
   $1.2\cdot 10^{-4}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $1.286$
   $1.5\cdot 10^{-3}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

\
& $a_1 (0.25)$ &

  ------
   Mean
   MSE
  ------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.250$
   $1.1\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.250$
   $1.1\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.250$
   $1.2\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

\
& &

  ------
   Mean
   MSE
  ------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.700$
   $6.0\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.699$
   $6.3\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $0.700$
   $6.7\cdot 10^{-5}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

\
& $b_0 (1.25)$ &

  ------
   Mean
   MSE
  ------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $1.249$
   $1.3\cdot 10^{-4}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $1.250$
   $1.3\cdot 10^{-4}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

&

  --------------------
        $1.250$
   $1.3\cdot 10^{-4}$
  --------------------

  : Sample mean and MSE of each parameter, for SRIVC and SRIVC-c, when $h=0.06, 0.2$ and $0.6[\textnormal{s}]$.

\

Multisine input: Irregular sampling
-----------------------------------

We consider the same system described before, with the same input and noise variance. In this simulation study, 2000 irregularly sampled output measurements are obtained. The sampling interval is distributed uniformly between $h_{lb}$ and $h_{hb}$, where the lower bound is fixed at $h_{lb}=0.05$, while the upper bound is varied from $0.1$ to $0.6$. A total of 6 Monte Carlo simulations are performed with each simulation containing 300 runs.

Figure \[fig3\] shows the mean value of each parameter, with their standard deviation around this value. As expected, the SRIVC-c estimator provides accurate estimates for all sampling period ranges in this study. On the other hand, the SRIVC estimator has a degrading performance as the sampling range increases, which could be attributed to the approximation errors in the prefilter calculations due to incorrect assumptions on the intersampling behaviour.

Arbitrary input application: Chirp signals
------------------------------------------

The next goal is to check whether the algorithms can provide accurate estimates for arbitrary input signals. Now, the system in is excited with an up-chirp signal, which is a CT signal that increases in frequency with time. These signals are widely used in signal processing applications such as radar systems and seismology and have been used for system identification [@xia1997system; @muller2001transfer]. The chirp signal used in this example is $$u(t)=\cos\left(f_0 \left[\frac{f_1}{f_0}\right]^{t/T_f} 2\pi t\right), \notag$$ where $f_0=0.1$\[Hz\] and $f_1=0.6$\[Hz\], and $T_f=500$\[s\] is the length of one period of the chirp signal. In this case, we determine the true system output using the explicit Runge-Kutta formulae RK5(4) [@dormand1980family]. The output is sampled every $h=0.5$\[s\] with $\delta=0.001$, and the measurement noise has variance $0.05$, which corresponds to approximately $10\%$ of the energy of the noiseless output. For the computation of the SRIVC-c estimate, we follow the algorithm described in Section \[arbitrary\]. The number of periods of the input signal vary from $1$ to $10$, which is equivalent to sample sizes ranging from 1000 to 10000, and $300$ Monte Carlo runs are performed for each case. The empirical evidence in Figure \[fig4\] suggests that the SRIVC-c estimator can also lead to accurate estimates for signals that are not described exactly by hold reconstructions nor multisines.

Conclusions {#sec6}
===========

In this paper, we have derived an algorithm for continuous-time system identification that is consistent for a wide class of input signals that have a known intersample behaviour. This estimator extends the applicability of the standard SRIVC method to inputs that are not exactly described by hold devices. We put forward a comprehensive analysis of the generic consistency of this estimator for multisine inputs, and extensive simulations have confirmed the advantages of this estimator over the widely popular SRIVC method. Further research on this topic concerns variance analyses of these estimators, and theoretical guarantees for irregular sampling schemes.

This work was partially supported by the Swedish Research Council under contract number 2016-06079 (NewLEADS) and by the Australian government Research Training Program (RTP) scholarship.

Appendix
========

**.** \[lemma1\] Consider $u(t)$ as in , and $V(t_k)$ and $U_{du}(t)$ as defined in and respectively. Under Assumption \[assumption2\], the matrix $\Psi$ defined in is equal to zero.

[Proof.]{} From the definition of $V(t_k)$, we directly obtain that all the entries $\Psi_{il}$ of $\Psi$, with , are equal to zero. For the other entries, we see that an arbitrary entry of this matrix is of the form $$\label{expec1}
\Psi_{il} = \overline{E}\left\{ \left\{\frac{p^{n+m+1-i}u(t)}{A_j^2(p)}\right\}_{t_k} \frac{p^{n+1-l}}{A_j(p)}v(t_k) \right\}.$$ If we define $g_i(t)$ as the inverse Laplace transform of $s^{n+m+1-i}A_j^{-2}(s)$, the first term in the expectation in  can be written as $$\left\{\frac{p^{n+m+1-i}u(t)}{A_j^2(p)}\right\}_{t_k} = \int_0^{t_k} g_i(t_k-\tau) u(\tau)\textnormal{d}\tau. \notag$$ Note that this is a DT signal, as a function of the time measurements $\{t_k\}$. On the other hand, the second term  can be described by $$\frac{p^{n+1-l}}{A_j(p)}v(t_k) = \sum_{r=1}^k \beta_{k-r,l} v(t_r). \notag$$ where $\{\beta_{j,l}\}_{j=0}^{k-1}$ are the first $k$ values of the impulse response of the FOH DT equivalent of $p^{n+1-l}A_j^{-1}(p)$. So, we compute $\Psi_{il}$ as $$\begin{aligned}
\Psi_{il} &= \overline{E}\left\{ \int_0^{t_k} g_i(t_k-\tau) u(\tau)\textnormal{d}\tau \sum_{r=1}^k \beta_{k-r,l} v(t_r) \right\} \notag \\
&=\lim_{N\to \infty}\frac{1}{N} \sum_{k=1}^N    \sum_{r=1}^k \int_0^{t_k} g_i(t_k-\tau)u(\tau) \beta_{k-r,l} E\{v(t_r)\} \textnormal{d}\tau \notag \\
&= 0, \notag\end{aligned}$$ where we have used the fact that the disturbance signal has zero mean.

**.** \[lemma2\] Under Assumptions \[assumption1\] to \[assumption5\], with $u(t)$ described as in , the following matrix is generically non-singular: $$\bar{\Phi}:= \overline{E}\left\{\left[\frac{U_{du}(t)}{\bar{A}^2(p)}\right]_{t_k} \left[\frac{U_{du}(t)}{\bar{A}(p)A^*(p)}\right]_{t_k}^\top\right\}. \notag$$

[Proof.]{} Similar to [@pan2020consistency], we follow an analyticity argument. We must first prove that $$\bar{\Phi}^*:= \overline{E}\left\{\left[\frac{U_{du}(t)}{{A^*}^2(p)}\right]_{t_k} \left[\frac{U_{du}(t)}{{A^*}^2(p)}\right]_{t_k}^\top\right\} \notag$$ is positive definite. For this, let $\textbf{z}\in \mathbb{R}^{n+m+1}$. We write $$\textbf{z}^\top \bar{\Phi}^* \textbf{z} = \overline{E}\left\{\left(\left\{\frac{B_{\textbf{z}}(p)}{{A^*}^2(p)}u(t)\right\}_{t_k}\right)^2 \right\} \geq 0. \notag$$ Since $u(t)$ is a multisine of the form , we have in stationary regime $$\frac{B_{\textbf{z}}(p)}{{A^*}^2(p)}u(t) = \tilde{\alpha}_0 + \sum_{l=1}^{m_u} \tilde{\alpha}_l \cos(\omega_l t + \tilde{\phi}_l), \notag$$ where $\tilde{\alpha}_0 = \alpha_0 B_{\textbf{z}}(0)/{A^*}^2(0)$, $\tilde{\alpha}_l = \alpha_l|B_{\textbf{z}}(i\omega_l)/{A^*}^2(i\omega_l)|$, and $\tilde{\phi}_l = \phi_l + \angle B_{\textbf{z}}(i\omega_l)/{A^*}^2(i\omega_l)$. Therefore,

$$\begin{aligned}
    \textbf{z}^\top \bar{\Phi}^* \textbf{z} &= \lim_{N\to \infty} \frac{1}{N} \sum_{k=1}^N \left(\tilde{\alpha}_0 + \sum_{l=1}^{m_u} \tilde{\alpha}_l \cos(\omega_l kh + \tilde{\phi}_l)\right)^2 \notag \\
    \label{cosinesums1}
    &\hspace{-1cm}= \lim_{N\to \infty} \frac{1}{N} \sum_{k=1}^N \Bigg(\tilde{\alpha}_0^2 + 2\tilde{\alpha}_0\sum_{l=1}^{m_u} \tilde{\alpha}_l \cos(\omega_l kh+ \tilde{\phi}_l) \\
    \label{cosinesums2}
    &\hspace{-0.7cm}+\sum_{j,l=1}^{m_u} \tilde{\alpha}_j\tilde{\alpha}_l \cos(\omega_j kh + \tilde{\phi}_j) \cos(\omega_l kh + \tilde{\phi}_l)\Bigg).
    \end{aligned}$$

Recall the formula for a geometric series $$\begin{aligned}
    \lim_{N\to \infty} \frac{1}{N}\sum_{k=1}^{N}\cos(\omega k + \phi) &= \textnormal{Re}\left\{\lim_{N\to \infty} \frac{e^{i\phi}}{N} \sum_{k=1}^N e^{i\omega k} \right\} \notag \\
        &= 0. \notag
    \end{aligned}$$ Using this result, and the identity $\cos(\alpha)\cos(\beta) = [\cos(\alpha+\beta)+\cos(\alpha-\beta)]/2$, the second term in the sum in is zero. Moreover, in the term for $j\neq l$ is a sum of sinusoids whose sum tends to zero as $N$ tends to infinity, while for $j= l$ constants appear. Thus, $$\label{expectationcomputation}
    \textbf{z}^\top \bar{\Phi}^* \textbf{z} = \tilde{\alpha}_0^2 + \frac{1}{2}\sum_{j=1}^{m_u} \tilde{\alpha}_j^2.$$ This computation leads to stating that $\textbf{z}^\top \bar{\Phi}^* \textbf{z}=0$ occurs if and only if $\tilde{\alpha}_0=\tilde{\alpha}_1 = \cdots = \tilde{\alpha}_{m_u}=0$, which in turn is equivalent to imposing $$\frac{B_{\textbf{z}}(0)}{{A^*}^2(0)} = 0, \hspace{0.15cm} \frac{B_{\textbf{z}}(i\omega_l)}{{A^*}^2(i\omega_l)} = \frac{B_{\textbf{z}}(-i\omega_l)}{{A^*}^2(-i\omega_l)}= 0, \hspace{0.15cm} l = 1,\dots, m_u. \notag$$ Since $m_u \geq (n+m)/2$, the only rational function that satisfies all of these restrictions is the null transfer function. Thus, $B_{\textbf{z}}(p)=0$ and $\textbf{z}=0$. With this, we have shown that $\bar{\Phi}^*$ is positive definite.

We now show that the entries of the matrix $\bar{\Phi}$ are real analytic functions of the (real-valued) parameters $(\bar{a}_1,\dots,\bar{a}_n)$ in the domain where $\bar{A}(p)$ is a stable polynomial. We denote this domain as . The entries of the matrix $\bar{\Phi}$ are given by $$\bar{\Phi}_{jl}:= \overline{E}\left\{\left\{\frac{p^{n+m+1-j}u(t)}{\bar{A}^2(p)}\right\}_{t_k}\hspace{-0.15cm} \left\{\frac{p^{n+m+1-l}u(t)}{\bar{A}(p)A^*(p)}\right\}_{t_k}\right\}, \notag$$ where $j,l=1,2,\dots,n+m+1$. By computing the expectation similarly to the derivation of , we find that $$\bar{\Phi}_{jl} =\tilde{\alpha}_0^j \tilde{\alpha}_0^l + \frac{1}{2} \sum_{r=1}^{m_u} \tilde{\alpha}_r^j \tilde{\alpha}_r^l \cos(\tilde{\phi}_r^j-\tilde{\phi}_r^l), \notag$$ where $$\begin{aligned}
    \tilde{\alpha}_0^j &= \begin{cases} 0, & \hspace{-0.2cm} j<n+m+1 \\
\alpha_0, & \hspace{-0.2cm}j= n+m+1\end{cases}, \hspace{0.15cm}
\tilde{\alpha}_0^l = \begin{cases} 0, & \hspace{-0.2cm}l<n+m+1 \\
\alpha_0, & \hspace{-0.2cm}l= n+m+1\end{cases}, \notag \\
\tilde{\alpha}_r^j &= \alpha_r\left|\frac{\omega_r^{n+m+1-j}}{\bar{A}^2(i\omega_r)}\right|, \hspace{0.85cm} \tilde{\phi}_r^j = \phi_r + \angle\left[ \frac{(i\omega_r)^{n+m+1-j}}{\bar{A}^2(i\omega_r)}\right] \notag \\
\tilde{\alpha}_r^l &= \alpha_r \left|\frac{\omega_r^{n+m+1-l}}{\bar{A}(i\omega_r)A^*(i\omega_r)}\right|, \hspace{0.24cm} \tilde{\phi}_r^l = \phi_r + \angle\left[\frac{(i\omega_r)^{n+m+1-l}}{\bar{A}(i\omega_r)A^*(i\omega_r)}\right]. \notag
    \end{aligned}$$

The coefficient $\tilde{\alpha}_r^j$ can be equivalently expressed as $$\begin{aligned}
    \tilde{\alpha}_r^j &= \frac{\alpha_r \omega_r^{n+m+1-j}}{\textnormal{Re}\{\bar{A}(i\omega_r)\}^2+\textnormal{Im}\{\bar{A}(i\omega_r)\}^2} \notag \\
    \label{ch5:alphar}
    &= \frac{\alpha_r\omega_r^{n+m+1-j}}{\left(1+ \sum\limits_{\substack{1\leq k \leq n \\ k\textnormal{ even}}}
         \hspace{-0.2cm}\bar{a}_k \omega_r^k (-1)^{\frac{k}{2}}\right)^{\hspace{-0.1cm}2}+\left(\sum\limits_{\substack{1\leq k \leq n \\ k\textnormal{ odd}}}\hspace{-0.2cm}\bar{a}_k \omega_r^k (-1)^{\frac{k-1}{2}}\right)^{\hspace{-0.1cm}2}}.
    \end{aligned}$$

From , we see that the denominator of $\tilde{\alpha}_r^j$ is a multivariate polynomial in the variables $(\bar{a}_1,\hspace{0.1cm} \dots,\hspace{0.1cm}\bar{a}_n)$, which is strictly positive in $\Omega$, since we know that $\bar{A}(p)$ is a stable polynomial for any $(\bar{a}_1,\hspace{0.1cm} \dots,\hspace{0.1cm}\bar{a}_n)\in \Omega$. This shows that the denominator of $\tilde{\alpha}_r^j$ is real analytic in $\Omega$, and since the quotient of real analytic functions is real analytic as long as the denominator does not vanish [@krantz2002primer Proposition 2.2.2], we have that $\tilde{\alpha}_r^j$ is real analytic in $\Omega$.

Similarly, the coefficient $\tilde{\alpha}_r^l$ can be written as $$\tilde{\alpha}_r^l = \frac{\alpha_r\omega_r^{n+m+1-l}}{|A^*(i\omega_r)|\sqrt{\textnormal{Re}\{\bar{A}(i\omega_r)\}^2+\textnormal{Im}\{\bar{A}(i\omega_r)\}^2}}. \notag$$ Here, the argument of the square root is a real analytic function, and the function $x\rightarrow 1/\sqrt{x}$ is real analytic for $x\in(0,\infty)$. Thus, by [@krantz2002primer Proposition 2.2.8], which states that the composition of real analytic functions is real analytic, we conclude that $\tilde{\alpha}_r^l$ is real analytic in $\Omega$.[^3]

Finally, note that $$\phi_r^j-\phi_r^l = \frac{\pi}{2}(l-j)+\angle A^*(i\omega_r)-\angle \bar{A}(i\omega_r), \notag$$ which leads to $$\begin{aligned}
    \cos(\phi_r^j-\phi_r^l) &= \frac{\textnormal{Re}\{\bar{A}(i\omega_r)\}}{|\bar{A}(i\omega_r)|} \cos\left(\frac{\pi}{2}(l-j)+\angle A^*(i\omega_r)\right) \notag \\
    &\hspace{-0.9cm}+ \frac{\textnormal{Im}\{\bar{A}(i\omega_r)\}}{|\bar{A}(i\omega_r)|} \sin\left(\frac{\pi}{2}(l-j)+\angle A^*(i\omega_r)\right). \notag 
    \end{aligned}$$ By the same justification above, $\textnormal{Re}\{\bar{A}(i\omega_r)\}/|\bar{A}(i\omega_r)|$ and $\textnormal{Im}\{\bar{A}(i\omega_r)\}/|\bar{A}(i\omega_r)|$ are real analytic functions for any . Therefore, $\cos(\phi_r^j-\phi_r^l)$ is real analytic for any $(\bar{a}_1, \dots,\bar{a}_n)\in \Omega$. Since any function defined by multiplication and addition of real analytic functions is real analytic [@krantz2002primer], we conclude that $\bar{\Phi}_{jl}$ is real analytic in the variables $(\bar{a}_1, \dots,\bar{a}_n)$, in the domain $\Omega$.

Hence, by Lemma A2.3 of [@soderstrom1983instrumental] and its corollary, the matrix $\bar{\Phi}$ is generically non-singular.

[^1]: In this context, generically non-singular means that the set $M=\{\theta_j\in \mathbb{R}^n\colon A_j(p)$ is a stable polynomial, $E\{\hat{\varphi}(t_k)\varphi^\top(t_k)\}$ is singular$\}$ has Lebesgue measure zero in $\mathbb{R}^n$.

[^2]: If no offset is considered, then at least $(n+m+1)/2$ sinusoids are required for our results.

[^3]: Note that the coefficients in $A^*(i\omega_r)$ do not play a role in the analyticity of $\tilde{\alpha}_r^l$, since $\tilde{\alpha}_r^l$ is viewed as a function of $(\bar{a}_1,\hspace{0.1cm} \dots,\hspace{0.1cm}\bar{a}_n)$ only.
