# A frequency-domain sequential Bayesian filter for sparse and broadband force estimation problems

Mathieu Aucejo, Olivier de Smet

► To cite this version:

Mathieu Aucejo, Olivier de Smet. A frequency-domain sequential Bayesian filter for sparse and broadband force estimation problems. Mechanical Systems and Signal Processing, 2025, 232, pp.112729. 10.1016/j.ymssp.2025.112729. hal-05038580v2

HAL Id: hal-05038580

https://hal.science/hal-05038580v2

Submitted on 5 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

<figure><img src="image_3.png" /></figure>

Full length article

# A frequency-domain sequential Bayesian filter for sparse and broadband force estimation problems

M. Aucejo<sup>\*ID</sup>, O. De Smet

Laboratoire de Mécanique des Structures et des Systèmes Couplés, Conservatoire National des Arts et Métiers, F 75003 Paris, France

<figure><img src="image_4.png" /></figure>

Check for
updates

## ARTICLE INFO

Communicated by J. Awrejcewicz

**Keywords:**

Linear inverse problem

Force identification

Space-frequency approach

Bayesian filter

## ABSTRACT

This paper presents a novel method for estimating the external sources acting on a mechanical structure in the frequency domain. Under the assumption of spatially sparse and broadband sources, a sequential Bayesian filter is derived. Its general structure follows that of a sequential Kalman-like filter, which is commonly used for input-state estimation problems in the time domain. This paper also includes an original Bayesian method for computing the noise variances of each measurement channel, which is a key element for the proper tuning of the proposed filtering algorithm. The proposed method is validated by a numerical experiment and an experimental application. The numerical experiment considers a simply supported beam subjected to a broadband point force under different operating conditions, while the experimental application deals with the identification of a point force acting on a simply supported plate. The comparison made with approaches available in the literature shows that the proposed strategy is able to estimate the external forces acting on a mechanical structure with the best trade-off between computational time/resources and accuracy.

## 1. Introduction

In structural dynamics, force estimation is a class of inverse problems that attempts to reconstruct an external excitation field applied to a mechanical structure from a set of vibration measurements. From a broad perspective, there are many methods available today to solve such a problem in both time and frequency domains, ranging from least squares strategies [1–4] to machine learning techniques [5–10]. However, if we limit our scope to the frequency domain, force identification methods have been developed to deal with various applications such as Transfer Path Analysis [11–13], material properties characterization [14], stochastic source identification [15,16], virtual sensing [17], uncorrelated sources identification [6]. Despite the existence of a wide variety of strategies, the landscape is dominated by two main families of methods. The first one is based on the filtering of the (analytical) equation of motion of the structure. Among the methods belonging to this family, one can mention the Force Analysis Technique and its derivatives [18–20] or the strategies derived from the virtual work principle, such as the virtual fields method [21–23] or the selective source identification [24]. The second class of methods are the regularization-based strategies, which are certainly the most widely used in structural dynamics because of their ability to constrain the space of admissible solutions by incorporating some prior knowledge about the sources to be identified. This family includes additive and multiplicative regularizations, such as the $\ell_2$-regularization (a.k.a. Tikhonov regularization) [12,25], the $\ell_1$-regularization [26] or the $\ell_p$-regularization [27–29], as well as the Bayesian regularization [30–34].

<footnote>\* Corresponding author.</footnote>

<footnote>E-mail address: mathieu.aucejo@lecnam.net (M. Aucejo).</footnote>
Although the above-mentioned force reconstruction techniques have proven their efficiency, some drawbacks remain when dealing with the identification of the space-frequency characteristics of the excitation sources, since, the inverse problem is generally solved at each frequency separately. However, such a solution process ignores the fact that the excitation spectrum may have some continuity. Consequently, large reconstruction errors can be observed at the resonance frequencies of lightly-damped structures due to the near singularity of the linear system to be solved (non-uniqueness of the solution). A possible solution to this problem can be found in mixed-norm regularization strategies [35–38], which are able to introduce prior knowledge about the space-frequency characteristics of the excitation sources. While mixed-norm regularization provides better results than standard regularization strategies, this approach is not without limitations. In particular, the higher accuracy of the reconstructed excitation fields comes at the price of a higher computational cost, which can be prohibitive for large-scale problems since the space-frequency problem is solved all at once.

To overcome the previous shortcomings, it may be of great interest to approach the problem recursively, i.e. frequency by frequency while keeping track of the information obtained at the previous frequencies. In the time domain, Kalman-like filtering offers such a possibility with good reconstruction accuracy [39–42]. Therefore, considering the adaptation of Kalman filtering in the frequency domain seems to be a promising approach. In fact, such a strategy has already been applied with some success in the field of audio signal processing for echo cancellation [43,44] or system identification [45,46].

This paper proposes to revisit the idea of frequency-domain Kalman filtering by introducing an original Bayesian filter to identify broadband and spatially sparse excitation sources. To this end, the paper is organized as follows. Section 2 presents the problem definition and the state space representation of the system. Section 3 introduces the proposed frequency-domain Bayesian filter and its detailed implementation. This section also describes a novel Bayesian strategy for estimating the measurement noise covariance matrix, which plays a crucial role in the tuning of the proposed filtering algorithm. Section 4 presents a numerical experiment to validate the proposed method, consisting of a simply supported beam subjected to a broadband point force under different operating conditions. Finally, Section 5 presents an experimental application of the proposed frequency-domain Bayesian filter to the identification of a point force acting on a simply supported plate. All the proposed validations are supported by comparisons with some existing methods to show the advantages or the limitations of the proposed strategy in terms of accuracy and computational time/resources.

## 2. Frequency-domain state space representation

Kalman filtering, and more generally Bayesian filtering, requires that the system be represented in a state space form. To construct the state space representation of the system, let us consider a linear time-invariant mechanical structure. In the frequency domain, the noiseless vibration field $x(f)$ is related to the external excitation field $u(f)$ by the transfer function matrix $H(f)$, which completely characterizes the dynamic behavior of the considered structure. This relation is given explicitly by:

$$x(f) = H(f)u(f), \quad (1)$$

where $f$ is the frequency.

Basically, the state space representation consists of a state equation, describing the evolution of the state vector $x(f)$ between two frequency steps, and an observation equation, relating the measured data to the system state and the input.

In the frequency domain, the state equation can be obtained by estimating the noiseless vibration field at a frequency $f + \delta f$, where $\delta f$ is a small frequency increment. By applying a first-order Taylor series, one has:

$$x(f + \delta f) = x(f) + \frac{\partial x(f)}{\partial f} \delta f + w(f), \quad (2)$$

where $w(f)$ represents the modeling error related to unmodeled higher order terms. In the following, the modeling error is regarded as a process noise, which is assumed to be a complex Gaussian noise with zero mean and covariance matrix $Q(f)$.

The derivative of the noiseless vibration field with respect to the frequency is simply expressed as:

$$\frac{\partial x(f)}{\partial f} = \frac{\partial H(f)}{\partial f} u(f) + H(f) \frac{\partial u(f)}{\partial f}. \quad (3)$$

Applying a first-order forward finite difference scheme to the frequency derivatives of the input vector and transfer function matrix, it comes that:

$$\frac{\partial x(f)}{\partial f} \delta f \approx [H(f + \delta f) - 2H(f)]u(f) + H(f)u(f + \delta f). \quad (4)$$

Consequently, from Eqs. (2) and (4), the state equation can be written as:

$$x_{k+1} = x_k + \mathbf{B}_k^{-} u_k + \mathbf{B}_{k+1}^{+} u_{k+1} + \mathbf{w}_k, \quad (5)$$

where $x_{k+1} = x(f_k + \delta f)$, $x_k = x(f_k)$, $\mathbf{B}_k^{-} = \mathbf{H}_{k+1}^{-} - 2\mathbf{H}_k$ and $\mathbf{B}_{k+1}^{+} = \mathbf{H}_k$.

Regarding the observation equation, it is obtained by assuming that the noiseless vibration field is corrupted by an additive complex Gaussian noise $v(f)$ with zero mean and covariance matrix $R(f)$. Thus, the observation equation is such that:

$$y_k = x_k + v_k, \quad (6)$$

where $y_k$ is the measured vibration field at frequency $f_k$.
Putting it all together, the frequency-domain state space representation of the system is given by:

$$ egin{cases} \mathbf{x}_{k+1} = \mathbf{x}_k + \mathbf{B}_k^{-} \mathbf{u}_k + \mathbf{B}_{k+1}^{+} \mathbf{u}_{k+1} + \mathbf{w}_k \\ \mathbf{y}_k = \mathbf{x}_k + \mathbf{v}_k \end{cases} \qquad (7) $$

At this point, it is worth noting that the previous formulation is non-standard, since the state vector at frequency step $k+1$ requires the knowledge of the input (excitation) vector at frequencies $k$ and $k+1$. To reduce the state space representation to its standard form, a reduced state $z_{k+1}$ is commonly introduced [47,48]. It is mathematically expressed as:

$$ \mathbf{z}_{k+1} = \mathbf{x}_k + \mathbf{B}_k^{-} \mathbf{u}_k + \mathbf{w}_k \qquad (8) $$

In doing so, the state space representation considered in the rest of the paper is finally given by:

$$ egin{cases} \mathbf{z}_{k+1} = \mathbf{z}_k + \mathbf{B}_k \mathbf{u}_k + \mathbf{w}_k \\ \mathbf{y}_k = \mathbf{z}_k + \mathbf{D}_k \mathbf{u}_k + \mathbf{v}_k, \end{cases} \qquad (9) $$

where $\mathbf{B}_k = \mathbf{B}_k^+ + \mathbf{B}_k^-$ and $\mathbf{D}_k = \mathbf{B}_k^+$.

### 3. Frequency-domain Bayesian filter

This section introduces the proposed frequency-domain Bayesian filter. Its structure follows that of a sequential Kalman-like filter, which is commonly used for input-state estimation problems in the time domain. The section describes the key elements of the proposed filtering algorithm, initial from a general overview to a detailed description of some specific technical aspects.

#### 3.1. Bayesian formulation of the filtering algorithm

The derivation of the input-state estimation algorithm is based on the Bayesian formulation proposed by Sedehi et al. in Ref. [41] and further unified and adapted by the authors in Refs. [42,49] to develop a family of sequential Bayesian filters in the time domain. From a general point of view, the considered Bayesian formulation is divided into 6 steps:

**0. Model definition** – It is essentially a rewrite of the state space representation given by Eq. (9) from a Bayesian perspective, that is:

$$ egin{cases} \mathbf{z}_{k+1} \sim p(\mathbf{z}_{k+1} | \mathbf{z}_k, \mathbf{u}_k) = \mathcal{N}_c(\mathbf{z}_{k+1} | \mathbf{z}_k + \mathbf{B}_k \mathbf{u}_k, \mathbf{Q}_k) \\ \mathbf{y}_k \sim p(\mathbf{y}_k | \mathbf{z}_k, \mathbf{u}_k) = \mathcal{N}_c(\mathbf{y}_k | \mathbf{z}_k + \mathbf{D}_k \mathbf{u}_k, \mathbf{R}_k), \end{cases} \qquad (10) $$

where $\mathcal{N}_c(\mathbf{x} | \mu, \Sigma)$ is the multivariate complex Gaussian distribution over the variable $\mathbf{x}$ with mean $\mu$ and covariance matrix $\Sigma$, while the notation $p(\mathbf{y} | \mathbf{x})$ stands for the conditional probability distribution of $\mathbf{y}$ given $\mathbf{x}$.

**1. Initialization** – This step requires the prior knowledge of the initial reduced state and input vectors $\mathbf{z}_0$ and $\mathbf{u}_0$. This information is encoded in the prior probability distributions $p(\mathbf{z}_0)$ and $p(\mathbf{u}_0)$, which are assumed to be multivariate complex Gaussian distributions, namely:

$$ \mathbf{z}_0 \sim p(\mathbf{z}_0) = \mathcal{N}_c(\mathbf{z}_0 | \hat{\mathbf{z}}_0, \mathbf{P}_0^z), \qquad (11a) $$

$$ \mathbf{u}_0 \sim p(\mathbf{u}_0) = \mathcal{N}_c(\mathbf{u}_0 | \hat{\mathbf{u}}_0, \mathbf{P}_0^u), \qquad (11b) $$

where $\hat{\mathbf{z}}_0$ and $\hat{\mathbf{u}}_0$ are the known mean state and input vectors, while $\mathbf{P}_0^z$ and $\mathbf{P}_0^u$ are the corresponding covariance matrices. Practically, only the knowledge of the initial prior distribution of the input vector $\mathbf{u}_0$ is needed, since the mean and the covariance matrix of the prior distribution of the reduced state vector $\mathbf{z}_0$ can be computed from the mean and the covariance matrix of the prior distribution of the input vector. This is done by combining Eqs. (1), (5) and (8) to obtain:

$$ \hat{\mathbf{z}}_0 = \overline{\mathbf{H}}_0 \hat{\mathbf{u}}_0 \quad \text{and} \quad \mathbf{P}_0^z = \overline{\mathbf{H}}_0 \mathbf{P}_0^u \overline{\mathbf{H}}_0^H, \qquad (12) $$

where $\overline{\mathbf{H}}_0 = \mathbf{H}_0 - \mathbf{B}_0^+$ and $\mathbf{x}^H$ is the Hermitian adjoint of $\mathbf{x}$.

Finally, the initialization is completed by the prediction of the reduced state given the initial measurement vector. This information is encoded in the predictive probability distribution $p(\mathbf{z}_1 | \mathbf{y}_0)$, whose expression is:

$$ \begin{aligned} p(\mathbf{z}_1 | \mathbf{y}_0) &= \iint_{\mathbf{z}_0} p(\mathbf{z}_1 | \mathbf{z}_0, \mathbf{u}_0) p(\mathbf{z}_0) p(\mathbf{u}_0) d\mathbf{z}_0 d\mathbf{u}_0 \\ &= \mathcal{N}_c(\mathbf{z}_1 | \tilde{\mathbf{z}}_1, \tilde{\mathbf{P}}_1^z), \end{aligned} \qquad (13) $$

where $\tilde{\mathbf{P}}_1^z$ is the covariance matrix associated to the predicted mean reduced state vector $\tilde{\mathbf{z}}_1$. These quantities are written as:

$$ \tilde{\mathbf{z}}_1 = \hat{\mathbf{z}}_0 + \mathbf{B}_0 \hat{\mathbf{u}}_0, \qquad (14a) $$

$$ \tilde{\mathbf{P}}_1^z = [\mathbf{I}_{n_z} \quad \mathbf{B}_0] \begin{bmatrix} \mathbf{P}_0^z & \mathbf{P}_0^{zu} \\ \mathbf{P}_0^{uH} & \mathbf{P}_0^u \end{bmatrix} [\mathbf{I}_{n_z} & \mathbf{B}_0^{uH}] + \mathbf{Q}_0, \qquad (14b) $$

with $\mathbf{P}_0^{zu} = \overline{\mathbf{H}}_0 \mathbf{P}_0^u$, the cross-covariance matrix between the reduced state and input vectors. $\mathbf{I}_{n_z}$ is the identity matrix of size $n_z$ ($n_z$: number of states).
2. **Input prediction** – This step requires the user to define the predictive distribution of the input vector given all the measurements collected up to the previous frequency step $k-1$, denoted by $\mathbf{y}_{1:k-1} = \{\mathbf{y}_1, \dots, \mathbf{y}_{k-1}\}$. Since our goal is to promote the spatial sparsity of the excitation field, the chosen predictive distribution must reflect this behavior.

In practice, several choices can be made at this point. For input estimation problems, two options are generally retained. The first one is to assume that the components of the input vector are independent and identically distributed real generalized Gaussian random variables [32,49]. The second option is to assume that the components of the input vector are independent and identically distributed complex Gaussian random variables [50–52]. In this paper, the second option is chosen, because it provides more flexibility in tuning the corresponding hyperparameters and a simpler computational implementation [42].

It follows from the preceding that the predictive distribution of the input vector given all the measurements collected up to the previous frequency step $k-1$ is written as:

$$p(\mathbf{u}_k | \mathbf{y}_{1:k-1}) = \prod_{i=1}^{n_u} \mathcal{N}_c(u_{ik} | 0, \tau_{ik}^{-1}), \qquad (15)$$

where $u_{ik}$ is the $i$th component of the input vector at the frequency step $k$, $\tau_{ik}$ is the corresponding precision parameter, $n_u$ is the number of input identification points.

Finally, to simplify the implementation in the subsequent steps, it should be noted that Eq. (15) is rewritten in a more compact form as:

$$p(\mathbf{u}_k | \mathbf{y}_{1:k-1}) = \mathcal{N}_c(\mathbf{u}_k | \mathbf{0}, \tilde{\mathbf{P}}_k^u), \qquad (16)$$

where $\tilde{\mathbf{P}}_k^u = \mathbf{T}_k^{-1}$ with $\mathbf{T}_k = \text{diag}(\tau_{1k}, \dots, \tau_{n_u k})$.

3. **Input estimation** – This step consists of applying the Bayes' rule to update our prediction of the input vector, given the information provided by the measurement at the current frequency step $k$. All this information is summarized in the filtering probability distribution $p(\mathbf{u}_k | \mathbf{y}_{1:k})$, which is written as:

$$ \begin{aligned} p(\mathbf{u}_k | \mathbf{y}_{1:k}) &\approx \int_{\mathbf{z}_k} p(\mathbf{u}_k | \mathbf{z}_k, \mathbf{y}_{1:k}) p(\mathbf{z}_k | \mathbf{y}_{1:k-1}) d\mathbf{z}_k \\ &\propto \int_{\mathbf{z}_k} p(\mathbf{y}_k | \mathbf{z}_k, \mathbf{u}_k) p(\mathbf{u}_k | \mathbf{y}_{1:k-1}) p(\mathbf{z}_k | \mathbf{y}_{1:k-1}) d\mathbf{z}_k \\ &= \mathcal{N}_c(\mathbf{u}_k | \hat{\mathbf{u}}_k, \mathbf{P}_k^u). \end{aligned} \qquad (17) $$

where $\hat{\mathbf{u}}_k$ and $\mathbf{P}_k^u$ are the mean and covariance matrix of the filtering distribution of the input vector at the frequency step $k$. They are expressed as:

$$\hat{\mathbf{u}}_k = \mathbf{K}_k^u \mathbf{i}_k, \quad \mathbf{i}_k = \mathbf{y}_k - \tilde{\mathbf{z}}_k \qquad (18a)$$

$$\mathbf{P}_k^u = (\mathbf{I}_{n_u} - \mathbf{K}_k^u \mathbf{D}_k) \tilde{\mathbf{P}}_k^u + \mathbf{K}_k^u \tilde{\mathbf{P}}_k^z \mathbf{K}_k^{uH}, \qquad (18b)$$

$$\mathbf{K}_k^u = \tilde{\mathbf{P}}_k^u \mathbf{D}_k^{\mathrm{H}} (\mathbf{D}_k \tilde{\mathbf{P}}_k^u \mathbf{D}_k^{\mathrm{H}} + \mathbf{R}_k)^{-1}. \qquad (18c)$$

4. **Reduced state estimation** – This step consists of updating the prediction of the reduced state vector given the information provided by the measurement at the current frequency step $k$. This step follows a similar procedure to that of the input estimation, meaning that the filtering probability distribution $p(\mathbf{z}_k | \mathbf{y}_{1:k})$ is written as:

$$ \begin{aligned} p(\mathbf{z}_k | \mathbf{y}_{1:k}) &= \int_{\mathbf{u}_k} p(\mathbf{z}_k | \mathbf{u}_k, \mathbf{y}_{1:k}) p(\mathbf{u}_k | \mathbf{y}_{1:k}) d\mathbf{u}_k \\ &\propto \int_{\mathbf{u}_k} p(\mathbf{y}_k | \mathbf{z}_k, \mathbf{u}_k) p(\mathbf{z}_k | \mathbf{y}_{1:k-1}) p(\mathbf{u}_k | \mathbf{y}_{1:k}) d\mathbf{u}_k \\ &= \mathcal{N}_c(\mathbf{z}_k | \hat{\mathbf{z}}_k, \mathbf{P}_k^z), \end{aligned} \qquad (19) $$

where $\hat{\mathbf{z}}_k$ and $\mathbf{P}_k^z$ are the mean and covariance matrix of the filtering distribution of the reduced state vector at the frequency step $k$. They are expressed as:

$$\hat{\mathbf{z}}_k = \tilde{\mathbf{z}}_k + \mathbf{K}_k^z (\mathbf{y}_k - \tilde{\mathbf{z}}_k - \mathbf{D}_k \hat{\mathbf{u}}_k), \qquad (20a)$$

$$\mathbf{P}_k^z = (\mathbf{I}_{n_z} - \mathbf{K}_k^z) \tilde{\mathbf{P}}_k^z + \mathbf{K}_k^z \mathbf{D}_k \mathbf{P}_k^u \mathbf{D}_k^{\mathrm{H}} \mathbf{K}_k^{uH}, \qquad (20b)$$

$$\mathbf{K}_k^z = \tilde{\mathbf{P}}_k^z (\tilde{\mathbf{P}}_k^z + \mathbf{R}_k)^{-1}. \qquad (20c)$$

At this point of the derivation, the cross-covariance matrix $\mathbf{P}_k^{zu}$, which is required to continue the filtering process, can be computed as:

$$ \begin{aligned} \mathbf{P}_k^{zu} &= \iint_{\mathbf{z}_k \mathbf{u}_k} (\mathbf{z}_k - \hat{\mathbf{z}}_k)(\mathbf{u}_k - \hat{\mathbf{u}}_k)^H p(\mathbf{z}_k, \mathbf{u}_k | \mathbf{y}_{1:k}) d\mathbf{z}_k d\mathbf{u}_k \\ &= -\mathbf{K}_k^z \mathbf{D}_k \mathbf{P}_k^u. \end{aligned} \qquad (21) $$
5. **Reduced state prediction** – This step consists of predicting the state vector at the next frequency step $k+1$ given the estimated reduced state and input vectors at step $k$, i.e:

$$ 
\begin{aligned}
p(\mathbf{z}_{k+1} | \mathbf{y}_{1:k}) &= \iint_{\mathbf{z}_k \mathbf{u}_k} p(\mathbf{z}_{k+1} | \mathbf{z}_k, \mathbf{u}_k) p(\mathbf{z}_k, \mathbf{u}_k | \mathbf{y}_{1:k}) d\mathbf{z}_k d\mathbf{u}_k \\
&= \mathcal{N}_c(\mathbf{z}_{k+1} | \tilde{\mathbf{z}}_{k+1}, \tilde{\mathbf{P}}_{k+1}^z),
\end{aligned} 
\quad (22) $$

where $\tilde{\mathbf{z}}_{k+1}$ and $\tilde{\mathbf{P}}_{k+1}^z$ are the mean and covariance matrix of the predictive distribution of the reduced state vector at the frequency step $k+1$. Their expressions are given by:

$$ \tilde{\mathbf{z}}_{k+1} = \hat{\mathbf{z}}_k + \mathbf{B}_k \hat{\mathbf{u}}_k, \quad (23a) $$

$$ \tilde{\mathbf{P}}_{k+1}^z = [\mathbf{I}_{n_z} \quad \mathbf{B}_k] \begin{bmatrix} \mathbf{P}_k^z & \mathbf{P}_k^{zu} \\ \mathbf{P}_k^{zuH} & \mathbf{P}_k^u \end{bmatrix} \begin{bmatrix} \mathbf{I}_{n_z} \\ \mathbf{B}_k^{H} \end{bmatrix} + \mathbf{Q}_k. \quad (23b) $$

In this paper, only the input vector is of primary interest. However, if the true state vector statistics are also needed, it is necessary to combine Eqs. (5) and (8) to finally infer that:

$$ \hat{\mathbf{x}}_k = \hat{\mathbf{z}}_k + \mathbf{B}_k^+ \hat{\mathbf{u}}_k, \quad (24a) $$

$$ \mathbf{P}_k^x = [\mathbf{I}_{n_z} \quad \mathbf{B}_k^+] \begin{bmatrix} \mathbf{P}_k^z & \mathbf{P}_k^{zu} \\ \mathbf{P}_k^{zuH} & \mathbf{P}_k^u \end{bmatrix} \begin{bmatrix} \mathbf{I}_{n_z} \\ \mathbf{B}_k^{H} \end{bmatrix}. \quad (24b) $$

## 3.2. Practical implementation

The practical implementation of the proposed frequency-domain Bayesian filter (FDBF) deriving from the Bayesian formulation presented in the previous section is summarized in Alg. 1. The Julia code implementing the algorithm is made publicly available on Zenodo [53] or can be directly downloaded on the Github repository of the first author<fnref n="1" />.

<caption>**Algorithm 1:** Frequency-domain Bayesian Filter – FDBF</caption>

Input: $\mathbf{y}_k, \mathbf{B}_k^+, \mathbf{B}_k, \mathbf{D}_k, \hat{\mathbf{u}}_0, \mathbf{P}_0^u, \mathbf{Q}_k, \mathbf{R}_k$,

Output: $\hat{\mathbf{u}}_k, \mathbf{P}_k^u, \hat{\mathbf{z}}_k, \mathbf{P}_k^x$

1. Initialization

$$ \tilde{\mathbf{z}}_0 = \bar{\mathbf{H}}_0 \hat{\mathbf{u}}_0, \quad \mathbf{P}_0^z = \bar{\mathbf{H}}_0 \mathbf{P}_0^u \bar{\mathbf{H}}_0^H, \quad \bar{\mathbf{H}}_0 = \mathbf{B}_k^+ - \mathbf{B}_0^+ $$

$$ \tilde{\mathbf{z}}_1 = \hat{\mathbf{z}}_0 + \mathbf{B}_0 \hat{\mathbf{u}}_0, \quad \tilde{\mathbf{P}}_1^z = [\mathbf{I}_{n_z} \quad \mathbf{B}_0] \begin{bmatrix} \mathbf{P}_0^z & \mathbf{P}_0^{zu} \\ \mathbf{P}_0^{zuH} & \mathbf{P}_0^u \end{bmatrix} \begin{bmatrix} \mathbf{I}_{n_z} \\ \mathbf{B}_0^{H} \end{bmatrix} + \mathbf{Q}_0 $$

for each frequency step $k > 0$ do

2. Input estimation

$$ \hat{\mathbf{i}}_k = \mathbf{y}_k - \tilde{\mathbf{z}}_k $$

$$ \tilde{\mathbf{P}}_k^u = \mathbf{T}_k^{-1} \quad \text{with} \quad \mathbf{T}_k = \text{diag}(\tau_{1k}, \dots, \tau_{n_u k}) $$

$$ \mathbf{K}_k^u = \tilde{\mathbf{P}}_k^u \mathbf{D}_k^H (\mathbf{D}_k \tilde{\mathbf{P}}_k^u \mathbf{D}_k^H + \mathbf{R}_k)^{-1} $$

$$ \hat{\mathbf{u}}_k = \mathbf{K}_k^u \hat{\mathbf{i}}_k, \quad \mathbf{P}_k^u = (\mathbf{I}_{n_u} - \mathbf{K}_k^u \mathbf{D}_k) \tilde{\mathbf{P}}_k^u + \mathbf{K}_k^u \tilde{\mathbf{P}}_k^z \mathbf{K}_k^u $$

3. Reduced state estimation

$$ \mathbf{K}_k^z = \tilde{\mathbf{P}}_k^z (\tilde{\mathbf{P}}_k^z + \mathbf{R}_k)^{-1}, \quad \hat{\mathbf{z}}_k = \tilde{\mathbf{z}}_k + \mathbf{K}_k^z (\hat{\mathbf{i}}_k - \mathbf{D}_k \hat{\mathbf{u}}_k), \quad \mathbf{P}_k^z = (\mathbf{I}_{n_z} - \mathbf{K}_k^z) \tilde{\mathbf{P}}_k^z + \mathbf{K}_k^z \mathbf{D}_k \mathbf{P}_k^u \mathbf{D}_k^H \mathbf{K}_k^z $$

4. True state estimation

$$ \hat{\mathbf{s}}_k = \hat{\mathbf{z}}_k + \mathbf{B}_k^+ \hat{\mathbf{u}}_k, \quad \mathbf{P}_k^x = [\mathbf{I}_{n_z} \quad \mathbf{B}_k^+] \begin{bmatrix} \mathbf{P}_k^z & \mathbf{P}_k^{zu} \\ \mathbf{P}_k^{zuH} & \mathbf{P}_k^u \end{bmatrix} \begin{bmatrix} \mathbf{I}_{n_z} \\ \mathbf{B}_k^{H} \end{bmatrix} $$

5. Reduced state prediction

$$ \tilde{\mathbf{z}}_{k+1} = \hat{\mathbf{z}}_k + \mathbf{B}_k \hat{\mathbf{u}}_k, \quad \tilde{\mathbf{P}}_{k+1}^z = [\mathbf{I}_{n_z} \quad \mathbf{B}_k] \begin{bmatrix} \mathbf{P}_k^z & \mathbf{P}_k^{zu} \\ \mathbf{P}_k^{zuH} & \mathbf{P}_k^u \end{bmatrix} \begin{bmatrix} \mathbf{I}_{n_z} \\ \mathbf{B}_k^{H} \end{bmatrix} + \mathbf{Q}_k $$

end

At this stage, the FDBF algorithm makes four main assumptions that can greatly affect the quality of the filtering process. First, Alg. 1 implicitly suggests to run the FDBF in a forward mode, i.e., from the lowest frequency of interest to the highest. However, the recursive nature of the FDBF algorithm suggests that it is possible to run the filter in a reverse mode, i.e., from the highest to the lowest frequency of interest. This particular mode will be called R-FDBF (for Reverse FDBF), while the forward mode will be simply referred to as FDBF in the rest of the paper. Then, it is assumed that the precision parameters $\tau_{ki}$ used in the input estimation step are known at each frequency. However, setting these parameters a priori is not straightforward and requires special attention. The third assumption concerns the existence of an initial input vector $\hat{\mathbf{u}}_0$ and its corresponding covariance matrix $\mathbf{P}_0^u$. In contrast to the time domain, the statistics of the initial input vector cannot be easily set to reasonable values (e.g., zero mean vector and isotropic covariance matrix) without prior knowledge. Finally, the last assumption is that the process and measurement noise covariance

<footnote n="1">¹ https://github.com/maucejo/KalmanFreq.jl</footnote>
matrices $\mathbf{Q}_k$ and $\mathbf{R}_k$ are known at each step. Although this assumption is reasonable, since the measurement noise covariance matrix can be estimated from the measurements taken at the beginning of the experiment, it requires some user expertise to set it correctly. For the process noise covariance matrix, the tuning is less obvious because it is related to unknown and/or unmodeled dynamics and depends on the assumptions made to derive the state equation. In the remainder of this section, automatic procedures are proposed to address the last three main assumptions made to derive the FDBF.

### 3.2.1. Computation of the estimated input vector

As mentioned in the preamble of Section 3.2, Alg. 1 assumes that the precision parameters $\tau_{ik}$ used in the input estimation step are known at each frequency. However, choosing an appropriate set of precision parameters a priori is far from an easy task. To address this problem, a simple and effective procedure is used in this paper. It is based on the observation that the estimated input vector $\hat{\mathbf{u}}_k$ computed using Eq. (18a) is the solution of the following Bayesian optimization problem [49]:

$$ \hat{\mathbf{u}}_k = \underset{\mathbf{u}_k}{\operatorname{argmax}} p(\mathbf{u}_k | \mathbf{i}_k) = \underset{\mathbf{u}_k}{\operatorname{argmax}} p(\mathbf{i}_k | \mathbf{u}_k) p(\mathbf{u}_k), \qquad (25) $$

where $p(\mathbf{i}_k | \mathbf{u}_k)$ is the likelihood function and $p(\mathbf{u}_k)$ is the prior distribution of the input vector. In the context of the FDBF, these two probability distributions are assumed to be multivariate complex Gaussian distributions, namely:

$$ p(\mathbf{i}_k | \mathbf{u}_k) = \mathcal{N}_c(\mathbf{i}_k | \mathbf{D}_k \mathbf{u}_k, \mathbf{R}_k) \quad \text{and} \quad p(\mathbf{u}_k) = \mathcal{N}_c(\mathbf{u}_k | \mathbf{0}, \hat{\mathbf{P}}_k^u). \qquad (26) $$

Based on this observation, the idea here is to extend the previous Bayesian optimization problem to also estimate the precision parameters $\tau_{ik}$. This is done in a straightforward manner by finding the optimal parameters $\mathbf{u}_k$ and $\tau_{ik}$ as the solution of the following optimization problem:

$$ (\hat{\mathbf{u}}_k, \hat{\tau}_{ik}) = \underset{(\mathbf{u}_k, \tau_{ik})}{\operatorname{argmax}} p(\mathbf{i}_k | \mathbf{u}_k) \prod_{i=1}^{n_u} p(u_{ik} | \tau_{ik}) p(\tau_{ik}), \qquad (27) $$

where $p(u_{ik} | \tau_{ik}) = \mathcal{N}_c(u_{ik} | 0, \tau_{ik}^{-1})$ is the prior probability distribution over the components of the input vector, while $p(\tau_{ik})$ is the prior probability distribution over the precision parameters $\tau_{ik}$. With conjugacy in mind, and given that the precision parameters are strictly positive, a classical and convenient choice for $p(\tau_{ik})$ is the gamma distribution, which is written as:

$$ G(\tau_{ik} | \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} \tau_{ik}^{\alpha-1} \exp(-\beta \tau_{ik}), \qquad (28) $$

where $\alpha$ and $\beta$ are the strictly positive shape and rate parameters, respectively, and $\Gamma(x)$ is the gamma function. Because the order of magnitude of $\tau_{ik}$ is usually unknown, $\alpha$ and $\beta$ are chosen so as to make the prior distribution weakly informative. For this reason, the shape and rate parameters are set to $\alpha = 1$ and $\beta = 10^{-18}$ in this paper.

The solution of the optimization problem defined by Eq. (27) is obtained by maximizing the full conditional probability distribution associated to each parameter, that is:

$$ \hat{\tau}_{ik} = \underset{\tau_{ik}}{\operatorname{argmax}} p(u_{ik} | \tau_{ik}) p(\tau_{ik}) \qquad (29a) $$

$$ \hat{\mathbf{u}}_k = \underset{\mathbf{u}_k}{\operatorname{argmax}} p(\mathbf{i}_k | \mathbf{u}_k) \prod_{i=1}^{n_u} p(u_{ik} | \tau_{ik}). \qquad (29b) $$

Since Eqs. (29a) and (29b) have analytical expressions, the resolution process consists of iteratively solving these two subproblems until convergence. At a given iteration $j$, the precision parameter $\hat{\tau}_{ik}^{(j)}$ and the input vector $\hat{\mathbf{u}}_k^{(j)}$ are updated as follows:

$$ \hat{\tau}_{ik}^{(j)} = \frac{\alpha}{\beta + |\hat{\mathbf{u}}_{ik}^{(j-1)}|^2} \quad \text{and} \quad \hat{\mathbf{u}}_k^{(j)} = \mathbf{K}_k^{\mathbf{u}(j)} \mathbf{i}_k, \qquad (30) $$

where:

$$ \mathbf{K}_k^{\mathbf{u}(j)} = \tilde{\mathbf{P}}_k^{\mathbf{u}(j)} \mathbf{D}_k^{\mathbf{H}} (\mathbf{D}_k \tilde{\mathbf{P}}_k^{\mathbf{u}(j)} \mathbf{D}_k^{\mathbf{H}} + \mathbf{R}_k)^{-1} \qquad (31a) $$

$$ \tilde{\mathbf{P}}_k^{\mathbf{u}(j)} = \operatorname{diag}(\hat{\tau}_{ik}^{(j-1)}, \dots, \hat{\tau}_{n_u k}^{(j-1)}). \qquad (31b) $$

Finally, since the process is iterative, it is necessary to define an initial input vector $\hat{\mathbf{u}}_k^{(0)}$ and to monitor its convergence. In this paper, the initial input vector is the estimated input vector at the previous frequency, i.e., $\hat{\mathbf{u}}_k^{(0)} = \hat{\mathbf{u}}_{k-1}$. Regarding the monitoring of the convergence, the iterative process is stopped when the relative error between two iterations $\delta$, defined such that:

$$ \delta = \frac{\| \hat{\mathbf{u}}_k^{(j)} - \hat{\mathbf{u}}_k^{(j-1)} \|_1}{\| \hat{\mathbf{u}}_k^{(j-1)} \|_1}, \qquad (32) $$

is less than $10^{-6}$.
### 3.2.2. Computation of the initial input vector

The FDBF algorithm being recursive, it is necessary to define the initial input vector $\hat{u}_0$ and its corresponding covariance matrix $\mathbf{p}_0$. However, this choice cannot be arbitrary, since the quality of the filtering process depends on it. A good initial point can usually be obtained from a Bayesian optimization similar to the one implemented in the previous section. Here, the optimization problem is written as:

$$ (\hat{\mathbf{u}}_0, \hat{\tau}_{i0}) = \underset{(\mathbf{u}_0, \tau_{i0})}{\operatorname{argmax}} p(\mathbf{y}_0 | \mathbf{u}_0) \prod_{i=1}^{n_u} p(u_{i0} | \tau_{i0}) p(\tau_{i0}). \quad (33) $$

where $p(u_{i0} | \tau_{i0})$ and $p(\tau_{i0})$ are the same as those used defined in Eqs. (27) and (28), while the likelihood function $p(\mathbf{y}_0 | \mathbf{u}_0) = \mathcal{N}_c(\mathbf{y}_0 | \mathbf{H}_0 \mathbf{u}_0, \mathbf{R}_0)$.

This problem is solved by exactly the same iterative process as the one described in the previous section. The only difference is that the initial input vector $\hat{\mathbf{u}}_0^{(0)}$ is computed by a standard Tikhonov regularization, namely:

$$ \hat{\mathbf{u}}_0^{(0)} = \underset{\mathbf{u}_0}{\operatorname{argmin}} \|\mathbf{y}_0 - \mathbf{H}_0 \mathbf{u}_0\|_{\mathbf{R}_0}^2 + \lambda \|\mathbf{u}_0\|_2^2, \quad (34) $$

where $\|\mathbf{x}\|_{\mathbf{R}}^2 = \mathbf{x}^{\mathrm{H}}\mathbf{R}^{-1}\mathbf{x}$ is the squared Mahalanobis distance and $\lambda$ is the regularization parameter, whose optimal value is automatically chosen from the L-curve criterion [54].

### 3.2.3. Estimation of the measurement noise covariance matrix

The last practical point of the FDBF algorithm concerns the estimation of the measurement noise covariance matrix $\mathbf{R}_k$ which is assumed to be known. However, over- or underestimating the components of this matrix can significantly affect the filtering process. In practice, a large area of research in the literature related to Kalman-like filtering is devoted to the online estimation of the measurement noise (but also process) noise covariance matrices [55–58]. In this paper, the estimation of the measurement noise covariance matrix is performed a priori from a simple offline Bayesian optimization procedure inspired by the works of Pereira et al. [59] and Garcia [60].

Let us assume that the measurement $\mathbf{y}_i \in \mathbb{C}^{n_f}$ ($n_f$: number of frequency steps) made by a sensor at location $i$ is statistically independent from those made by the other sensors. Let us also assume that the measurement noise variance $\sigma_i^2 \in \mathbb{R}$ for the $i$th sensor is constant over the frequency range. Under these assumptions, the measurement matrix is diagonal and such that:

$$ \mathbf{R}_k = \mathbf{R} = \operatorname{diag}(\sigma_1^2, \dots, \sigma_{n_y}^2), \quad (35) $$

where $n_y$ is the number of sensors.

For each measurement channel $i$, the observed signal $\mathbf{y}_i$ is related to the noiseless signal $\mathbf{x}_i$ by the following equation:

$$ \mathbf{y}_i = \mathbf{x}_i + \mathbf{v}_i, \quad (36) $$

where $\mathbf{v}_i$ is the measurement noise vector on the sensor $i$.

Theoretically, the noiseless signal $\mathbf{x}_i$ can be estimated from the observed signal $\mathbf{y}_i$ by maximizing the following posterior probability distribution:

$$ p(\mathbf{x}_i | \mathbf{y}_i) \propto p(\mathbf{y}_i | \mathbf{x}_i) p(\mathbf{x}_i), \quad (37) $$

where $p(\mathbf{y}_i | \mathbf{x}_i)$ is the likelihood function and $p(\mathbf{x}_i)$ is the prior distribution of the noiseless signal. In the context of signal denoising, these two probability distributions can be chosen as multivariate complex Gaussian distributions, namely:

$$ p(\mathbf{y}_i | \mathbf{x}_i) = \mathcal{N}_c(\mathbf{y}_i | \mathbf{x}_i, \sigma_i^2 \mathbf{I}_{n_f}) \quad \text{and} \quad p(\mathbf{x}_i) = \mathcal{N}_c(\mathbf{x}_i | \mathbf{0}, \alpha_i^2 (\mathbf{D}^\mathrm{T}\mathbf{D})^{-1}), \quad (38) $$

where $\sigma_i^2$ and $\alpha_i^2$ are variance parameters, while $\mathbf{D} \in \mathbb{R}^{n_f \times n_f}$ is the second difference matrix, defined such that:

$$ \mathbf{D} = \begin{bmatrix} 1 & -1 & & & \\ -1 & 2 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 2 & -1 \\ & & & -1 & 1 \end{bmatrix}, \quad (39) $$

which allows to impose a smoothness constraint on the noiseless signal (here $\mathcal{C}^2$-continuity).

The goal here is to estimate the variance parameters $\sigma_i^2$ and $\alpha_i^2$ from the data only. This can be done by maximizing the following posterior probability distribution:

$$ p(\sigma_i^2, \alpha_i^2 | \mathbf{y}_i) \propto p(\mathbf{y}_i | \sigma_i^2, \alpha_i^2) p(\sigma_i^2, \alpha_i^2). \quad (40) $$

Assuming that $p(\sigma_i^2, \alpha_i^2) \propto 1$, i.e. that any set of values $(\sigma_i^2, \alpha_i^2)$ is equiprobable, it can be shown that [59]:

$$ \begin{align} p(\sigma_i^2, \alpha_i^2 | \mathbf{y}_i) &\propto p(\mathbf{y}_i | \sigma_i^2, \alpha_i^2) = \int_{\mathbf{x}_i} p(\mathbf{y}_i | \mathbf{x}_i, \sigma_i^2) p(\mathbf{x}_i | \alpha_i^2) d\mathbf{x}_i \tag{41} \\ &= \mathcal{N}_c(\mathbf{y}_i | \mathbf{0}, \alpha_i^2 [\mathbf{D}^\mathrm{T}\mathbf{D}]^{-1} + \lambda \mathbf{I}_{n_f}), \nonumber \end{align} $$
where $\lambda = \sigma_i^2/\alpha_i^2$.

Following Strang [61] and Yueh [62], $p(y_i | \sigma_i^2, \alpha_i^2)$ is explicitly given by:

$$p(y_i | \sigma_i^2, \alpha_i^2) = \frac{1}{\pi^{n_f} (\alpha_i^2)^{n_f} \prod_{k=1}^{n_f} g_k(\lambda)} \exp \left( -\frac{1}{\alpha_i^2} \sum_{k=1}^{n_f} \frac{|p_{ik}|^2}{g_k(\lambda)} \right), \quad (42)$$

where:

* $p_i$ is the discrete cosine transform of $y_i$;

* $g_k(\lambda) = \frac{1+\lambda s_k^2}{s_k^2}$;

* $s_k = 2 \left[ \cos \left( \frac{(k-1)\pi}{n_f} \right) - 1 \right]$ is the $k$th eigenvalue of D.

At this stage, maximizing $p(\sigma_i^2, \alpha_i^2 | y_i)$ w.r.t. the parameter $\alpha_i^2$ allows obtaining its optimal value $\hat{\alpha}_i^2$ as a function of $\lambda$ that is:

$$\hat{\alpha}_i^2(\lambda) = \frac{1}{n_f} \sum_{k=1}^{n_f} \frac{|p_{ik}|^2}{g_k(\lambda)}. \quad (43)$$

It follows that the process can only be completed by finding the optimal value of $\lambda$ given the measured data $y_i$, i.e. the value that maximizes the conditional probability distribution $p(\lambda | y_i)$. Its expression is obtained by marginalizing the joint probability distribution $p(\lambda, \alpha_i^2 | y_i)$ with respect to the parameter $\alpha_i^2$, i.e:

$$\begin{aligned} p(\lambda | y_i) &= \int_{\alpha_i^2} p(\lambda, \alpha_i^2 | y_i) d\alpha_i^2 \\ &= \int_{\alpha_i^2} p(\sigma_i^2, \alpha_i^2 | y_i) \left| \frac{\partial \sigma_i^2}{\partial \lambda} \right| d\alpha_i^2 \\ &= \int_{\alpha_i^2} p(\sigma_i^2, \alpha_i^2 | y_i) \alpha_i^2 d\alpha_i^2 \end{aligned} \quad (44)$$

After some calculations, which are not detailed here for the sake of brevity, it can be shown that:

$$p(\lambda | y_i) \propto \left[ \hat{\alpha}_i^2(\lambda)^{n_f-2} \prod_{k=1}^{n_f} g_k(\lambda) \right]^{-1}. \quad (45)$$

It follows that the optimal value of $\lambda$, denoted $\hat{\lambda}$, is computed by solving the following minimization problem:

$$\hat{\lambda} = \underset{\lambda}{\operatorname{argmin}} (n_f - 2) \log \hat{\alpha}_i^2(\lambda) + \sum_{k=1}^{n_f} \log g_k(\lambda). \quad (46)$$

The previous minimization problem is solved using a standard univariate optimization algorithm on a fixed interval, such as the Brent's method [63].

Finally, the optimal value of the noise variance $\sigma_i^2$, noted $\hat{\sigma}_i^2$ is obtained by:

$$\hat{\sigma}_i^2 = \hat{\lambda} \hat{\alpha}_i^2(\hat{\lambda}). \quad (47)$$

The proposed procedure is summarized in Alg. 2.

## Algorithm 2: Estimation of R

Input: y

Output: R

1. Compute $s_k = 2 \left[ \cos \left( \frac{(k-1)\pi}{n_f} \right) - 1 \right]$, $\forall k = \{1, \dots, n_f\}$

for $i = 1$ to $n_y$ do

2. Compute the discrete cosine transform of $y_i - p_i = \text{DCT}(y_i)$

3. Solve $\hat{\lambda} = \underset{\lambda}{\operatorname{argmin}} (n_f - 2) \log \hat{\alpha}_i^2(\lambda) + \sum_{k=1}^{n_f} \log g_k(\lambda)$

with $\hat{\alpha}_i^2(\lambda) = \frac{1}{n_f} \sum_{k=1}^{n_f} \frac{|p_{ik}|^2}{g_k(\lambda)}$ and $g_k(\lambda) = \frac{1+\lambda s_k^2}{s_k^2}$

4. Compute the optimal noise variance $-\hat{\sigma}_i^2 = \hat{\lambda} \hat{\alpha}_i^2(\hat{\lambda})$

end

5. Construct the measurement noise covariance matrix

$R = \text{diag}(\hat{\sigma}_1^2, \dots, \hat{\sigma}_{n_y}^2)$
<figure><img src="image_5.png" /><figcaption>Fig. 1. Graphical summary of the implementation of the FDBF algorithm.</figcaption></figure>

### 3.2.4. Tuning of the process noise covariance matrix

As classically done in the literature, the process noise covariance matrix $Q_k$ can be assumed to be isotropic and constant over the frequency range of interest. This means that it is defined by a single variance parameter, denoted $\sigma_x^2$. A heuristic value can be found by analyzing the process used to derive the state equation.

From a theoretical point of view, the state equation is derived from a first-order Taylor series applied to Eq. (1) for a frequency increment $\delta f$. This implies that for any function $h(f)$ and a “sufficiently small” increment $\delta f$:

$$h(f + \delta f) = h(f) + \frac{\partial h(f)}{\partial f} \delta f + \frac{1}{2} \frac{\partial^2 h(f)}{\partial f^2} \delta f^2 + \dots \quad (48)$$

In the present case, this means that the process noise $w(f)$ should be theoretically at least proportional to $\delta f^2$. However, because other sources of uncertainty can arise in the model, such as those related to the construction of the transfer function matrix or to the continuity of the quantities involved in the formulation of the state space model, it seems reasonable to consider that the process noise is proportional to $\delta f$.

Under this assumption, the variance parameters can be expressed as:

$$\sigma_x^2 = (\gamma \cdot \delta f)^2, \quad (49)$$

where $\gamma$ can be thought as an estimate of the amplitude of the unmodeled terms per Hz.

The value of the parameter $\gamma$ can be calibrated using automatic offline methods such as the L-curve principle [39] or innovation-based strategies [64]. However, some of the computational efficiency of the proposed approach will be lost if such strategies are implemented. Therefore, this paper proposes a heuristic value to achieve a trade-off between solution accuracy and computational efficiency. After a series of experiments conducted on a set of numerical and experimental data, we have found that $\gamma = 2.5 \times 10^{-3}$ provides a sufficiently good compromise for a relatively wide range of frequency resolution.

Finally, it should be noted that such approximation of the process noise variance is consistent with the intuition that our confidence in the state evolution model should increase as the frequency resolution decreases.

### 3.2.5. Graphical summary of the implementation

Fig. 1 summarizes the implementation of the FDBF algorithm. It is divided into three main steps:
1. **Data acquisition** – The first step is to acquire the vibration data of the structure, namely to measure the acceleration spectra $y_i(\omega)$ at each sensor $i$ and to measure/model the transfer function matrices $H_{ij}(\omega)$ between a sensor $i$ and a potential excitation point $j$ over the frequency range of interest.

2. **Data preprocessing** – The second step is the preprocessing of the data. This includes the construction of the matrices of the state space representation of the system, $B_k^+$, $B_k$, and $D_k$, calculating the initial solution $\hat{u}_0$ and its covariance matrix $P_0^u$, and estimating the process and measurement noise covariance matrix $Q$ and $R$.

3. **Filtering** – The last step is to run the FDBF algorithm given in Alg. 1.

# 4. Numerical experiment

This section aims at validating the proposed FDBF algorithm on a numerical example using an academic case study. In particular, the proposed numerical experiment is designed to assess the performance of the FDBF algorithm in various operating conditions by comparing it with some regularization strategies available in the literature.

## 4.1. Problem definition

The purpose of this numerical validation is to test the performance of the previously described method under different operating conditions. This section first presents the general information common to all the test cases considered, before detailing them in a second part.

### 4.1.1. General information

The structure under consideration is a simply supported stainless steel beam. The beam is 1 m long, 3 cm wide and 1 cm thick. The material properties are $E = 210$ GPa for the Young's modulus and $\rho = 7850$ kg m⁻³ for the density. A modal damping factor is considered and set to 1%. It is also assumed that the structure is excited at $x_0 = 62$ cm from the left end of the beam by a broadband point force between 100 Hz and 1 kHz. In this frequency range, the beam has 4 vibration modes (211 Hz, 375 Hz, 586 Hz and 844 Hz).

Regarding the data generation, the noiseless acceleration data were synthesized from a finite element model consisting of 100 beam elements. Selection matrices were further applied to fit the measurement mesh defined in Fig. 2. The data were then corrupted with an additive Gaussian white noise with a controlled signal-to-noise ratio (SNR). On the other hand, the transfer function matrices needed to construct the state space representation were calculated analytically from a classical mode expansion using the first 9 modes of the structure (i.e. up to 2 kHz).

To run the proposed filter, it remains to define the noise covariance matrices $R$ and $Q$. Regarding the measurement noise covariance $R$, it is computed a priori from the data using the procedure described in Section 3.2.3. As for the process noise covariance $Q$, which is supposed to be isotropic, the value of the variance parameters $\sigma_x^2$ is computed according to Eq. (49).

In this paper, three metrics are considered to measure the reconstruction accuracy of the proposed method:

* **Global Relative Error (GRE)**: it provides a measure of how close the estimated input vector is to the true input vector over the entire structure. It is defined as:

$$ \text{GRE} = \frac{\|\hat{\mathbf{u}} - \mathbf{u}_{\text{ref}}\|_1}{\|\mathbf{u}_{\text{ref}}\|_1}. \qquad (50) $$

$\mathbf{u}_{\text{ref}}$ and $\hat{\mathbf{u}}$ represent respectively the reference and estimated input vectors of dimensions $n_u \times n_f$ defined over the whole frequency range.

* **Frequency Relative Error (FRE)**: it provides a measure of how close the input spectrum estimated at the point force location, $\hat{\mathbf{u}}^{\text{exc}} \in \mathbb{C}^{n_f}$, is to the true input spectrum $\mathbf{u}_{\text{ref}}^{\text{exc}} \in \mathbb{C}^{n_f}$. It is defined as:

$$ \text{FRE} = \frac{\|\hat{\mathbf{u}}^{\text{exc}} - \mathbf{u}_{\text{ref}}^{\text{exc}}\|_1}{\|\mathbf{u}_{\text{ref}}^{\text{exc}}\|_1}. \qquad (51) $$

* **Correlation Coefficient (CC)**: it measures the space-frequency correspondence between the amplitudes of the reference and estimated input vectors. It is defined as:

$$ \text{CC} = \frac{|\text{vec}(\mathbf{u}_{\text{ref}})|^T |\text{vec}(\hat{\mathbf{u}})|}{\|\text{vec}(\mathbf{u}_{\text{ref}})\|_2 \|\text{vec}(\hat{\mathbf{u}})\|_2}, \qquad (52) $$

where $\text{vec}(\mathbf{x})$ is the vectorization of the matrix $\mathbf{x}$.

In addition to these accuracy metrics, it is also interesting to evaluate how long an algorithm takes to run and how much memory it uses. This information is obtained by monitoring the computation time $t_e$ and the memory footprint $m_f$ of the algorithm. These two additional metrics allow to estimate the trade-off between computational time/resources and accuracy.

Finally, to provide the reader with some points of comparison, the FDBF, run in forward and reverse modes, is compared with three regularization strategies available in the literature. The first one is a variant of the Relevant Vector Regression (RVR) proposed by Tipping in Ref. [50], which consists mainly in solving the Bayesian optimization problem defined by Eq. (33) at each frequency.
<table><caption>**Table 1**
Summary of the different test cases considered in the numerical experiment.</caption><thead><tr><th>Case ID</th><th>Excitation type</th><th>SNR (dB)</th><th>δf (Hz)</th><th>Sensor configuration</th></tr></thead><tbody><tr><td>1</td><td>Unitary</td><td>25</td><td>0.5</td><td>Fully collocated</td></tr><tr><td>2</td><td>Random</td><td>25</td><td>0.5</td><td>Fully collocated</td></tr><tr><td>3</td><td>Unitary</td><td>15</td><td>0.5</td><td>Fully collocated</td></tr><tr><td>4</td><td>Unitary</td><td>25</td><td>[0.5, 5]</td><td>Fully collocated</td></tr><tr><td>5</td><td>Unitary</td><td>25</td><td>0.5</td><td>Under-determined, non-collocated</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Geometry of the beam and reference sensors configuration – (●) Accelerometers and (■) Excitation.</figcaption></figure>

The second method is the multiplicative $\ell_q$-regularization, which consists in solving the following minimization problem at each frequency:

$$ \hat{\mathbf{u}}_k = \underset{\mathbf{u}_k \in \{0\}}{\operatorname{argmin}} \|\mathbf{y}_k - \mathbf{H}_k \mathbf{u}_k\|_{\mathbf{R}}^2 \cdot \|\mathbf{u}_k\|_q^q. \qquad (53) $$

Here, the value of $q$ is set to 0.5 in order to promote the spatial sparsity of the input vector. From a practical point of view, the problem is solved using a fixed-point algorithm (see Ref. [29] for the implementation details). The last strategy is the multiplicative (mixed-norm) $\ell_{p,q}$-regularization, which consists in solving the following minimization problem for all frequencies at once:

$$ \hat{\mathbf{U}} = \underset{\mathbf{U} \in \{0\}}{\operatorname{argmin}} \|\mathbf{Y} - \bar{\mathbf{H}}\mathbf{U}\|_{\bar{\mathbf{R}}}^2 \cdot \|\mathbf{U}\|_{p,q}^q, \qquad (54) $$

where $\mathbf{Y} = \operatorname{vec}(\mathbf{y})$ and $\mathbf{U} = \operatorname{vec}(\mathbf{u})$, while $\bar{\mathbf{H}} = \operatorname{diag}(\mathbf{H}_1, \dots, \mathbf{H}_{n_f})$ and $\bar{\mathbf{R}} = \mathbf{I}_{n_f} \otimes \mathbf{R}$ ($\otimes$ stands for the Kronecker product). Here, the value of $p$ is set to 2 and $q$ to 0.5 in order to promote the continuity of the force spectrum as well as its spatial sparsity. The problem is solved using the fixed-point algorithm described in Ref. [36].

### 4.1.2. Test cases description

A test case is defined by four main parameters, namely the type of excitation force, the sensor configuration, the measurement noise level, and the frequency resolution. This section describes the different test cases considered in the present numerical experiment. They are mainly variations of a reference configuration described below, where only one parameter differs. A summary of the different test cases is proposed in Table 1.

#### Case 1 – Reference configuration

In the reference configuration, a set of 20 accelerometer is mounted along the beam and one of these sensors if collocated with the excitation, which is supposed to be a unit point force (See Fig. 2).

In this configuration, the same grid is used for both measurement and identification, which corresponds to a fully collocated configuration. Finally, the SNR is set to 25 dB and the frequency resolution to 0.5 Hz.

#### Case 2 – Random excitation

The second test case aims to evaluate the influence of the type of excitation signal on the performance of the FDBF algorithm. For this purpose, the unit excitation is replaced by a white noise excitation applied at point $x_0$ and defined such that

$$ \mathbf{u}_{\text{ref}}^{\text{exc}} = \mu_u + \sigma_u \mathbf{r}, \qquad (55) $$

where $\mu_u = 2$ N, $\sigma_u = 0.1$ N and $\mathbf{r} \sim \mathcal{N}(\mathbf{r} | \mathbf{0}, \mathbf{I}_{n_f})$.

In this test case, the sensor configuration, the measurement noise level and the frequency resolution are those prescribed in the reference configuration.

#### Case 3 – Measurement noise level

The level of measurement noise can significantly affect the performance of estimation algorithms. In general, the estimation accuracy decreases as the noise level increases. To evaluate the robustness of the proposed algorithm with respect to the measurement
<figure><img src="image_7.png" /><figcaption>Fig. 3. Under-determined and non-collocated sensor configuration.</figcaption></figure>

<table><caption>**Table 2**
Performance metrics – Reference configuration.</caption><thead><tr><th>Method</th><th>GRE (%)</th><th>FRE (%)</th><th>CC (%)</th><th>t<sub>e</sub> (s)</th><th>m<sub>f</sub> (MiB)</th></tr></thead><tbody><tr><td>FDBF</td><td>3.22</td><td>3.22</td><td>99.96</td><td>0.71</td><td>533</td></tr><tr><td>R-FDBF</td><td>3.28</td><td>3.28</td><td>99.96</td><td>0.56</td><td>533</td></tr><tr><td>RVR</td><td>6.63</td><td>4.19</td><td>98.32</td><td>0.60</td><td>403</td></tr><tr><td>ℓ<sub>q</sub>-reg.</td><td>16.41</td><td>7.02</td><td>97.94</td><td>1.53</td><td>1070</td></tr><tr><td>ℓ<sub>p,q</sub>-reg.</td><td>2.75</td><td>2.71</td><td>99.97</td><td>1.07</td><td>646</td></tr></tbody></table>

noise level, the signal-to-noise ratio (SNR) of the vibration data is reduced to 15 dB, while the excitation, the sensor configuration and the frequency resolution are those prescribed in the reference configuration.

## Case 4 – Frequency resolution

A first-order Taylor expansion and a first-order forward finite difference scheme are used to derive the state equation. Consequently, the frequency resolution is a key parameter in the FDBF algorithm and it is legitimate to wonder to which extent the approximations made limit the frequency resolution. For this reason, this test case assesses the performance of the proposed approach by varying frequency resolution from 0.5 Hz to 5 Hz. The excitation, the measurement noise level and the sensor configuration are those prescribed in the reference configuration.

## Case 5 – Sensor configuration

An important parameter when implementing a force estimation method is the choice of the sensor configuration. The ideal case consists in the well-determined configuration, where the identification and measurement points are collocated and equal in number. However, in real applications such a configuration is rarely met (accessibility, limited number of available sensors, ...). Therefore, it is worthwhile to evaluate the performance of the proposed strategy in a more realistic setting. Here, it is proposed to apply the different identification strategies in the under-determined and non-collocated configuration presented in Fig. 3, for which $n_z = 5$ and $n_u = 20$.

In this test case, the excitation, the measurement noise level and the frequency resolution are those defined in the reference configuration.

## 4.2. Results

This section collects the results and analyses corresponding to the five test cases described in the previous section. Special attention is given to explaining the differences observed between the methodologies compared in this paper.

### 4.2.1. Case 1 – Reference configuration

As a first application, the FDBF and the R-FDBF are compared to the three regularization strategies described in the previous section on the reference configuration. As shown in Fig. 4, all the methods are globally successful in estimating the location of the force as well as its amplitude. However, the visual inspection of the results shows however that the FDBF/R-FDBF and the $ℓ_{p,q}$-regularization perform better than the RVR and the $ℓ_q$-regularization. A closer look at the identified force spectrum at the point force location presented in Fig. 5 shows that the FDBF/R-FDBF and the $ℓ_{p,q}$-regularization behave rather similarly, while the RVR and the $ℓ_q$-regularization exhibit large discrepancies near the resonance frequencies of the beam.

The quantitative analysis proposed in Table 2 shows that the $ℓ_{p,q}$-regularization is slightly more accurate than the FDBF/R-FDBF, which provide very similar results. However, the FDBF and the R-FDBF are slightly faster and have a smaller memory footprint than the $ℓ_{p,q}$-regularization. This is due to the repeated inversion of a large matrix of dimensions $(n_z \cdot n_f) \times (n_u \cdot n_f)$ for the $ℓ_{p,q}$-regularization, which is computationally more expensive than the repeated inversion of a $n_z \times n_u$ matrix at each frequency step for the FDBF/R-FDBF.

The fastest and lightest approach is the RVR, but its accuracy is worse than that of the FDBF and the $ℓ_{p,q}$-regularization due to large reconstruction errors at the resonance frequencies of the structure. Finally, the worst performance is obtained with the $ℓ_q$-regularization, for which the computation time and memory footprint reveal convergence difficulties of the resolution algorithm, mainly noticeable around the resonance frequencies of the beam.

To complete the analysis of the reference configuration, it is interesting to evaluate the ability of the different methods to estimate the state vector $x_k$, even though its estimation is not of primary interest in the present study. Such information is particularly relevant in the context of virtual sensing or vibration level monitoring and may paves the way for future applications of the proposed method. Here, it is important to note that the state vector is a direct output the filtering algorithm for the FDBF and its reverse version, while it is recovered from Eq. (1) for the other estimation methods.
<figure><img src="image_8.png" /><figcaption>Fig. 4. Estimated force distribution for the reference configuration.</figcaption></figure>

The performance metrics computed from the estimated state vector and collected in Table 3 show that all methods are able to reconstruct the noiseless acceleration field with high accuracy in the reference configuration, despite inaccurate estimations of the input vector for the RVR and the $ℓ_q$ regularization. At first sight, this result is surprising, but it may be related to the fact that the estimation of the input vector is the solution of an optimization problem, which aims to provide the best possible compromise between the fit to the measured data (data fidelity) and the prior information about the excitation field used to constrain the estimation.

For the sake of the completeness and to have a visual interpretation of the values of the performance metrics, the acceleration at the driving point computed with the FDBF is proposed in Fig. 6.
<figure><img src="image_9.png" /><figcaption>Fig. 5. Estimated force spectrum from the reference configuration.</figcaption></figure>

<table><caption>**Table 3**
Performance metrics – State estimation in the reference configuration.</caption><thead><tr><th>Method</th><th>GRE (%)</th><th>FRE (%)</th><th>CC (%)</th></tr></thead><tbody><tr><td>FDBF</td><td>2.08</td><td>2.68</td><td>99.99</td></tr><tr><td>R-FDBF</td><td>2.41</td><td>3.40</td><td>99.99</td></tr><tr><td>RVR</td><td>2.18</td><td>3.52</td><td>99.99</td></tr><tr><td>ℓ<sub>q</sub>-reg.</td><td>3.33</td><td>4.36</td><td>99.98</td></tr><tr><td>ℓ<sub>p,q</sub>-reg.</td><td>2.07</td><td>3.31</td><td>99.99</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 6. Acceleration estimated at the driving point by the FDBF in the reference configuration.</figcaption></figure>

<table><caption>Table 4
Performance metrics – White noise excitation.</caption><thead><tr><th>Method</th><th>GRE (%)</th><th>FRE (%)</th><th>CC (%)</th><th>t<sub>e</sub> (s)</th><th>m<sub>f</sub> (MiB)</th></tr></thead><tbody><tr><td>FDBF</td><td>3.26</td><td>3.26</td><td>99.96</td><td>0.59</td><td>543</td></tr><tr><td>R-FDBF</td><td>3.29</td><td>3.29</td><td>99.96</td><td>0.75</td><td>546</td></tr><tr><td>RVR</td><td>8.56</td><td>5.28</td><td>97.12</td><td>0.47</td><td>312</td></tr><tr><td>ℓ<sub>q</sub>-reg.</td><td>17.06</td><td>8.34</td><td>97.16</td><td>1.67</td><td>1035</td></tr><tr><td>ℓ<sub>p,q</sub>-reg.</td><td>2.73</td><td>2.70</td><td>99.97</td><td>1.10</td><td>676</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Fig. 7. Input force estimated by the FDBF in the case of a white noise excitation.</figcaption></figure>

### 4.2.2. Case 2 – Random excitation

Random excitations are generally challenging for recursive filtering methods because the filter must quickly adapt to the variation in excitation amplitude from one frequency step to another. The results presented in Table 4 are consistent with those obtained in the previous section, since the performance metrics are globally similar. Consequently, the type of excitation signal is not a discriminating factor in the present numerical experiment. However, it should be noted that the accuracy metrics of the RVR and the ℓ<sub>q</sub>-regularization are slightly worse than those obtained for the reference configuration. For the sake of the completeness, the excitation field and excitation spectrum estimated by the FDBF algorithm are shown in Fig. 7.
<figure><img src="image_12.png" /><figcaption>Fig. 8. Space-frequency distribution of the input force estimated by the FDBF and R-FDBF for a SNR of 15 dB.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 9. Force spectrum estimated by the FDBF and the R-FDBF for a SNR of 15 dB.</figcaption></figure>

### 4.2.3. Case 3 – Measurement noise level

Figs. 8 and 9 show that the FDBF fails to correctly estimate the excitation force for such an SNR, while the R-FDBF provides an accurate estimation of the input force.

A closer analysis reveals that these results are strongly correlated with the accuracy of the initial estimated input vector. As shown in Fig. 10. This result highlights the influence of the initial solution on the performance of the filtering process. This observation suggests that more accurate solutions can be obtained by choosing an initial frequency for which the solution is easier to compute. In fact, accurate solutions are difficult to compute at frequencies close to the resonances of the system, because the transfer function matrix becomes ill-conditioned in these regions, especially for lightly damped structures. In the present case, the initial frequency of 100 Hz is close to the second resonance frequency of the beam at 94 Hz. This is not the case at 1 kHz, as this frequency is sufficiently far away from the two adjacent resonance frequencies of the beam, namely 844 Hz and 1149 Hz. Further analysis also shows that the condition number of the transfer function matrix at 100 Hz is twice as large at 1 kHz.

To make the previous statement more concrete, the FDBF is run with the initial solution computed at 1 kHz. This fictitious situation allows the pivotal role played by the accuracy of the initial solution in the performance of the FDBF to be demonstrated. In this case, Fig. 11 shows that the FDBF is able to accurately estimate the input force when an accurate initial solution is available.

Finally, to give the reader a fair overview of the results obtained, the performance metrics computed for all the methods considered in this paper are given in Table 5. As expected, the results deteriorate with decreasing SNR. Again, the accuracy metrics of the RVR and the $ℓ_q$-regularization are significantly worse than those of the R-FDBF and the $ℓ_{p,q}$-regularization, which are on par in terms of accuracy. However, the R-FDBF achieves the best trade-off between the computation time/resources and accuracy.

### 4.2.4. Case 4 – Frequency resolution

In this section, only the influence of the frequency resolution on the performance of the FDBF, the R-FDBF and the $ℓ_{p,q}$-regularization (a.k.a. mixed norm regularization) is examined, since the performance of the RVR and the $ℓ_q$-regularization should be insensitive to this parameter.
<figure><img src="image_14.png" /><figcaption>Fig. 10. Comparison of the initial solution of the FDBF and the R-FDBF computed for a SNR of 15 dB.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 11. Input force estimated by the FDBF with an accurate initial solution for a SNR of 15 dB.</figcaption></figure>

<table><caption>Table 5
Performance metrics – SNR = 15 dB.</caption><thead><tr><th>Method</th><th>GRE (%)</th><th>FRE (%)</th><th>CC (%)</th><th>t<sub>e</sub> (s)</th><th>m<sub>f</sub> (MiB)</th></tr></thead><tbody><tr><td>FDBF</td><td>104</td><td>100</td><td>0</td><td>0.61</td><td>498</td></tr><tr><td>FDBF (fictitious situation)</td><td>10.89</td><td>10.89</td><td>99.47</td><td>0.65</td><td>570</td></tr><tr><td>R-FDBF</td><td>10.06</td><td>10.06</td><td>99.59</td><td>0.67</td><td>577</td></tr><tr><td>RVR</td><td>30.93</td><td>21.98</td><td>89.34</td><td>0.47</td><td>390</td></tr><tr><td>ℓ<sub>q</sub>-reg.</td><td>41.03</td><td>23.67</td><td>88.70</td><td>1.48</td><td>983</td></tr><tr><td>ℓ<sub>p,q</sub>-reg.</td><td>9.29</td><td>9.27</td><td>99.63</td><td>0.89</td><td>676</td></tr></tbody></table>

Fig. 12(a) presents the evolution of the GRE indicator as the frequency resolution varies from 0.5 Hz to 5 Hz. The results show that the GRE obtained with the three methods considered is constant and similar up to a frequency resolution of 3 Hz. Above this frequency resolution, the (forward) FDBF fails to compute a good estimate of the input force vector, while the results obtained by the R-FDBF and the $ℓ_{p,q}$-regularization remain similar. For a deeper insight, the mean SNR estimated by the method proposed in Section 3.2.3 is plotted in Fig. 12(b). This figure shows that the target SNR is correctly estimated for frequency resolution below 1 Hz. For a frequency resolution comprises between 1 Hz and 2.5 Hz, the actual noise level is underestimated, while above this range, it is overestimated above this range. This observation highlights the sensitivity of the proposed noise variance estimation procedure to frequency resolution. This sensitivity has a direct impact on the accuracy of the initial solution computed for the FDBF and R-FDBF, as emphasized in Section 3.2.2, and subsequently on the entire filtering process. Fig. 13 illustrates that the calculation
<figure><img src="image_16.png" /><figcaption>Fig. 12. Influence of the frequency resolution on (a) the GRE indicator and (b) the SNR of the measured data estimated by the proposed Bayesian approach.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Influence of the frequency resolution on the accuracy of the computed initial point for the FDBF and the R-FDBF for a frequency resolution of 5 Hz.</figcaption></figure>

<table><caption>Table 6
Performance metrics – Frequency resolution of 5 Hz.</caption><thead><tr><th>Method</th><th>GRE (%)</th><th>FRE (%)</th><th>CC (%)</th><th>t<sub>e</sub> (s)</th><th>m<sub>f</sub> (MiB)</th></tr></thead><tbody><tr><td>FDBF</td><td>104</td><td>100</td><td>0</td><td>0.19</td><td>51</td></tr><tr><td>R-FDBF</td><td>8.58</td><td>8.58</td><td>99.69</td><td>0.13</td><td>61</td></tr><tr><td>RVR</td><td>34.94</td><td>25.28</td><td>85.06</td><td>0.06</td><td>34</td></tr><tr><td>l<sub>q</sub>-reg.</td><td>31.70</td><td>10.42</td><td>97.12</td><td>0.18</td><td>127</td></tr><tr><td>l<sub>p,q</sub>-reg.</td><td>6.86</td><td>4.71</td><td>99.90</td><td>0.17</td><td>160</td></tr></tbody></table>

of the initial solution is less sensitive to the accuracy of the noise variance estimation at higher frequencies. Since the same causes have the same effects, the observed consequences on the performance of the proposed method are the same as those described in the previous section.

Finally, for the sake of the completeness, the results obtained for all the methods and a frequency resolution of 5 Hz are summarized in Table 6. In this case, the results show that the R-FDBF and the $l_{p,q}$-regularization are the most accurate methods. They are also less sensitive to the overestimation of the measurement noise level, unlike the RVR and the $l_q$-regularization, which solve the inverse problem at each frequency in a non-recursive manner. Again, the R-FDBF achieves the best trade-off between the computation time/resources and the solution accuracy.
<table><caption>**Table 7**
Performance metrics – Under-determined and non-collocated configuration.</caption><thead><tr><th>Method</th><th>GRE (%)</th><th>FRE (%)</th><th>CC (%)</th><th>t<sub>e</sub> (s)</th><th>m<sub>f</sub> (MiB)</th></tr></thead><tbody><tr><td>FDBF</td><td>106</td><td>100</td><td>0</td><td>0.31</td><td>217</td></tr><tr><td>R-FDBF</td><td>11.31</td><td>11.31</td><td>99.48</td><td>0.37</td><td>280</td></tr><tr><td>RVR</td><td>51.85</td><td>31.45</td><td>78.87</td><td>0.82</td><td>460</td></tr><tr><td>ℓ'<sub>q</sub>-reg.</td><td>3077</td><td>160</td><td>7.64</td><td>0.71</td><td>393</td></tr><tr><td>ℓ'<sub>p,q</sub>-reg.</td><td>4550</td><td>229</td><td>8.15</td><td>4.21</td><td>2305</td></tr></tbody></table>

### 4.2.5. Case 5 – Sensor configuration

The results obtained for this particular configuration are summarized in Table 7. In particular, we can see that the FDBF, the $ℓ'_q$-regularization and the $ℓ'_{p,q}$-regularization fail to provide a good estimate of the excitation field. The RVR seems to behave reasonably well, while the R-FDBF outperforms all the previous methods in terms of accuracy and computational performance. A closer analysis of the memory footprint of the RVR and the $ℓ'_q$ and $ℓ'_{p,q}$-regularization methods points out convergence difficulties of the respective resolution algorithms. In addition, as observed earlier, the difference between the forward FDBF and its reverse version lies in the accuracy of the initial solutions (they are not shown for the sake of the brevity). For a proper interpretation of the results, the excitation field estimated by each method is proposed in Fig. 14.

## 4.3. Summary

As suggested in Section 3.2 and further developed in this section, the performance of the proposed method is determined by several parameters. As for any Kalman-like algorithm, the accuracy of the estimates and the efficiency of the resolution algorithm are conditioned to:

* the sensor configuration (equal-determined or under-determined configurations),

* the noise level in the measured data,

* the availability of an accurate initial solution,

* a proper tuning of the covariance matrices **Q** and **R**,

* the way of introducing prior information of the excitation field to identify.

More specifically, it has been shown here that the proposed method is sensitive to the accuracy of the initial solution computed at the first frequency step. This sensitivity is particularly noticeable when the frequency resolution is high or the SNR is low. In these cases, the FDBF may fail to provide a correct estimate of the excitation field, while the R-FDBF is successful because it benefits from a more accurate initial solution. The need for an accurate initial solution can also be considered from an algorithmic point of view, since the solution at $k-1$ serves as an initial guess for computing the solution at $k$. Consequently, the initial solution is propagated over the frequencies. Second, the proposed method provides the best trade-off between solution accuracy and computational efficiency. This can be explained by the fact that the proposed method is based on a recursive filtering algorithm, which is able to take advantage of the frequency continuity of the excitation spectrum while solving the problem at each frequency step by taking into account the spatial sparsity of the excitation field. This property has a significant impact on the computational performance of the proposed method, since the recursive nature of the algorithm allows the iterative process used to estimate the excitation vector at a given frequency to converge faster and require less memory. This is not the case with the other approaches, which either solve the problem at each frequency independently and may suffer from convergence issues, or use the whole measured spectra to estimate the solution leading to a higher memory footprint.

# 5. Experimental validation

This section aims to validate the findings of the numerical experiment by conducting tests in a real-world environment. The primary goal is to evaluate how well the FDBF performs under real operating conditions and real data sets. By applying the filter to a real-world scenario, such as identifying a broadband excitation source in the frequency domain, we can determine the effectiveness and practicality of the proposed filter.

## 5.1. Description of the experimental set-up

The experimental setup consists of a thin aluminum plate, measuring 420 mm long, 360 mm wide, and 3 mm thick, which is mounted on a specially designed support structure that closely mimics simply supported boundary conditions (see Fig. 15(a)). The edges of the plate are glued to thin blades, which are then screwed to a rigid frame [65]. These blades are carefully selected to provide an optimum balance between flexibility and stiffness. They must be sufficiently flexible to allow rotational degrees of freedom, yet stiff enough to support the weight of the plate and limit the normal translational motion of the edges.

The plate is instrumented with 5 lightweight IEPE accelerometers (0.5 g – PCB Piezotronics 352C22) distributed as shown in Fig. 15(b). This sensor configuration is used to measure up to 2 kHz with a frequency resolution of 0.125 Hz, the Frequency Response Functions (FRFs) needed to construct the transfer function matrix $H_k$, and the acceleration spectra that form the output vector $y_k$.
<figure><img src="image_18.png" /><figcaption>Fig. 14. Estimated force distribution for the under-determined and non-collocated configuration.</figcaption></figure>

This sensor configuration has been chosen to avoid redundancy in the measurements and to ensure sufficient acceleration levels in the frequency range of interest by not placing the sensors in areas where the vibration level is negligible such as on the modal lines. In doing so, each sensor brings information about the excitation and the measured transfer functions matrix is of rank 5. Regardless of the type of the data considered (FRF or spectrum), the excitation is provided by a Maul-Theet vImpact-63 automatic impact hammer, which is an impact controller equipped with a miniature instrumented impact hammer (4.8 g – PCB Piezotronics 086E80). The hammer provides a bounce-free impact, which is able to excite a structure up to 60 kHz with a maximal amplitude of
<figure><img src="image_19.png" /><figcaption>Fig. 15. Experimental setup – (a) Simply-supported aluminum thin plate and (b) Sensor distribution on the plate.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 16. Typical force and acceleration spectra measured with the automatic impact hammer and the lightweight accelerometers.</figcaption></figure>

<table><caption>**Table 8**
Performance metrics – Collocated experimental configuration.</caption><thead><tr><th>Method</th><th>GRE (%)</th><th>FRE (%)</th><th>CC (%)</th><th>t<sub>e</sub> (s)</th><th>m<sub>f</sub> (MiB)</th></tr></thead><tbody><tr><td>FDBF</td><td>40.75</td><td>17.93</td><td>96.36</td><td>0.90</td><td>312</td></tr><tr><td>R-FDBF</td><td>15.35</td><td>6.5</td><td>99.17</td><td>0.89</td><td>319</td></tr><tr><td>RVR</td><td>54.92</td><td>10.96</td><td>73.90</td><td>0.37</td><td>165</td></tr><tr><td>ℓ<sub>q</sub>-reg.</td><td>44.16</td><td>9.06</td><td>95.31</td><td>1.14</td><td>670</td></tr><tr><td>ℓ<sub>p,q</sub>-reg.</td><td>55.42</td><td>11.0</td><td>73.81</td><td>0.65</td><td>354</td></tr></tbody></table>

220 N pk. A typical force spectrum is shown in Fig. 16(a). It could be noticed that this spectrum is relatively flat over the frequency range of interest.

The remainder of this section discusses two specific configurations. The first one corresponds to a collocated configuration where the excitation is applied at the location of one of the acceleration sensors. The second one corresponds to a non-collocated configuration where the excitation is applied outside of the sensor locations. For both configurations, the force spectrum resulting from the hammer impact is reconstructed between 400 Hz and 2 kHz. In this frequency range, the plate has approximately 20 modes (see Fig. 16(b)). Finally, the process noise covariance matrix is calibrated as described in Section 3.2.4 and the measurement noise covariance matrix is estimated from the collected data following the procedure described in Section 3.2.3.

## 5.2. Collocated configuration

In this first configuration, the impact excitation is applied at the location of sensor $S_3$ (see Fig. 15(b)). The results obtained for this configuration are presented in Figs. 17 and 18, while the corresponding performance metrics are summarized in Table 8.

The results show that the FDBF and the R-FDBF provide the best trade-off between solution accuracy and memory consumption. While the RVR is the fastest and lightest method, it is also one of the least accurate with the $ℓ_{p,q}$-regularization, which has similar
<figure><img src="image_21.png" /><figcaption>Fig. 17. Estimated force distribution for the collocated experimental configuration.</figcaption></figure>

performance. These poor results are due to large reconstruction errors at some frequencies, corresponding to low-level acceleration responses. Finally, the $ℓ_q$-regularization provides reasonable performance metrics, despite high reconstruction errors at non-excited locations. The larger memory footprint is due to the number of iterations required to achieve convergence of the resolution algorithm.

### 5.3. Non-colocated configuration

In this second configuration, the impact excitation is applied at a location outside of the sensor locations. As shown in Fig. 19, the number of sensors remains the same as in the collocated configuration, meaning that the problem is also slightly under-determined.
<figure><img src="image_22.png" /><figcaption>Fig. 18. Estimated force spectrum for the collocated experimental configuration.</figcaption></figure>

The results obtained for this configuration are shown in Figs. 20 and 21, while the corresponding performance metrics are summarized in Table 9. It should be noted that the results for the $\ell_q$-regularization are not reported here, since the resolution algorithm failed to provide a solution due to some singularity problems. Regarding the other methods, the RVR provides the most accurate solution, while the FDBF and the $\ell_{p,q}$-regularization have the worst performance. Again, the R-FDBF provides the best trade-off between accuracy and computational resources.
<figure><img src="image_23.png" /><figcaption>Fig. 19. Non-collocated configuration – (S) Sensor locations and (E) Excitation point.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 20. Estimated force distribution for the non-collocated experimental configuration.</figcaption></figure>

6. Conclusion

This paper introduced an original frequency-domain Bayesian filter that can operate in both forward and reverse modes. The basic motivation of this work was to propose a lightweight and efficient procedure for the estimation of sparse and broadband
<figure><img src="image_25.png" /><figcaption>Fig. 21. Estimated force spectrum for the non-collocated experimental configuration.</figcaption></figure>

<table><caption>Table 9
Performance metrics – Non-collocated experimental configuration.</caption><thead><tr><th>Method</th><th>GRE (%)</th><th>FRE (%)</th><th>CC (%)</th><th>t<sub>e</sub> (s)</th><th>m<sub>f</sub> (MiB)</th></tr></thead><tbody><tr><td>FDBF</td><td>242</td><td>40.82</td><td>35.64</td><td>0.89</td><td>359</td></tr><tr><td>R-FDBF</td><td>10.19</td><td>5.82</td><td>99.68</td><td>0.86</td><td>358</td></tr><tr><td>RVR</td><td>8.72</td><td>2.22</td><td>99.91</td><td>1.64</td><td>480</td></tr><tr><td>ℓ<sub>q</sub>-reg.</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>ℓ<sub>p,q</sub>-reg.</td><td>661.85</td><td>123.72</td><td>13.52</td><td>1.78</td><td>1150</td></tr></tbody></table>

sources acting on a linear mechanical structure. The proposed methodology is based on the definition of an appropriate state space representation of the dynamic behavior of the structure in the frequency domain and on the derivation of a recursive algorithm, rooted in a Bayesian framework, that promotes both the spatial sparsity of the excitation sources and the continuity of their frequency spectrum. The proposed algorithm has been validated both numerically and experimentally and compared with some regularization strategies available in the literature. The numerical and real-world experiments conducted have demonstrated the ability of the proposed strategy to accurately estimate the space-frequency characteristics of the excitation sources under different operating conditions (excitation signal, frequency resolution, measurement noise level and sensor configuration) with a low computational cost and a small memory footprint. The results presented in the paper also show that the proposed approach provides the best trade-off between computation time/resources and accuracy with respect to the regularization strategies implemented for comparison purposes. Finally, although it was not the initial motivation of this work, the proposed methodology offers the possibility of estimating the state of the system, which here corresponds to the noiseless vibration field. This observation opens up new perspectives for the development of new virtual sensing strategies.
CRediT authorship contribution statement

**M. Aucejo:** Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**O. De Smet:** Writing – original draft, Validation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

References

[1] B.J. Dobson, E. Rider, A review of the indirect calculation of excitation forces from measured structural response data, Proc. Inst. Mech. Eng. Part C: Mech. Eng. Sci. 204 (2) (1990) 69–75.

[2] S.R. Ibrahim, A. Fregolent, A. Sestieri, Structural force identification at unmeasured locations, in: Proceedings of the 14th International Modal Analysis Conference, Dearborn, MI, United States, 1996.

[3] S.S. Law, T.H.T. Chan, Moving force identification: A time domain method, J. Sound Vib. 201 (1) (1997) 1–22.

[4] D. Goutaudier, D. Gendre, V. Kehr-Candille, R. Ohayon, Single-sensor approach for impact localization and force reconstruction by using discriminating vibration modes, Mech. Syst. Signal Process. 138 (2020) 106534.

[5] J.M. Zhou, L. Dong, W. Guan, J. Yan, Impact load identification of nonlinear structures using deep Recurrent Neural Network, Mech. Syst. Signal Process. 133 (2019) 106292.

[6] C. Wang, D. Chen, J. Chen, X. Lai, T. He, Deep regression adaptation networks with model-based transfer learning for dynamic load identification in the frequency domain, Eng. Appl. Artif. Intell. 102 (2021) 104244.

[7] H. Yang, J. Jiang, G. Chen, J. Zhao, Dynamic load identification based on deep convolution neural network, Mech. Syst. Signal Process. 185 (2023) 109757.

[8] R. Zhou, B. Qiao, J. Liu, W. Cheng, X. Chen, Impact force localization and reconstruction via gated temporal convolutional network, Aerosp. Sci. Technol. 144 (2024) 108819.

[9] R. Zhou, B. Qiao, L. Jiang, W. Cheng, X. Yang, Y. Wang, X. Chen, A model-based deep learning approach to interpretable impact force localization and reconstruction, Mech. Syst. Signal Process. 224 (2025) 111977.

[10] G.R. Tondo, I. Kavrakov, G. Morgenthal, Efficient dynamic modal load reconstruction using physics-informed Gaussian processes based on frequency-sparse Fourier basis functions, Mech. Syst. Signal Process. 225 (2025) 112295.

[11] M.H.A. Janssens, J.W. Verheij, The use of an equivalent forces method for the experimental quantification of structural sound transmission in ships, J. Sound Vib. 226 (2) (1999) 305–328.

[12] A.N. Tite, D.J. Thompson, The quantification of structure-borne transmission paths by inverse meth-ods. part 2: Use of regularization techniques, J. Sound Vib. 264 (2003) 433–451.

[13] J. Ortega Almiron, F. Bianciardi, P. Corbeels, R. Ullmann, W. Desmet, Mount characterization analysis in the context of FBS for component-based TPA on a wiper system, in: Proceedings of Forum Acusticum 2020, Lyon, France, 2020.

[14] F. Ablitzer, C. Pézerat, J.-M. Génevaux, J. Begué, Identification of stiffness and damping properties of plates by using the local equation of motion, J. Sound Vib. 333 (2014) 2454–2468.

[15] O. Grosset, J.-H. Thomas, F. Ablitzer, C. Pézerat, Identification of the acoustic wall pressure of turbulent flow exciting structures, in: Proceedings of the 23rd International Congress on Sound and Vibration, Athens, Greece, 2016.

[16] Y. Liu, L. Wang, A two-step weighting regularization method for stochastic excitation identification under multi-source uncertainties based on response superposition-decomposition principle, Mech. Syst. Signal Process. 182 (2023) 109565.

[17] K. Janssens, B. Forrier, R. Cumbo, E. Risalti, B. Cornelis, T. Tamarozzi, W. Desmet, Virtual sensing application cases exploiting various degrees of model complexity, in: Proceedings of Asia-Pacific Vibration Conference 2019, 2019.

[18] C. Pézerat, J.-L. Guyader, Two inverse methods for localization of external sources exciting a beam, Acta Acust. 3 (1995) 1–10.

[19] Q. Leclère, C. Pézerat, Vibration source identification using corrected finite difference schemes, J. Sound Vib. 331 (2012) 1366–1377.

[20] J. Carpentier, J.-H. Thomas, C. Pézerat, Development of the corrected force analysis technique in polarcoordinates. Application on car window excited by a turbulent air, J. Sound Vib. 571 (2024) 118105.

[21] A. Berry, O. Robin, F. Pierron, Identification of dynamic loading on a bending plate usingthe Virtual Fields Method, J. Sound Vib. 333 (2014) 7151–7164.

[22] A. Berry, O. Robin, Identification of spatially correlated excitations on a bending plate using the Virtual Fields Method, J. Sound Vib. 375 (2016) 76–91.

[23] N. Madinier, A. Berry, Q. Leclère, K. Ege, Comparison of two force identification methods from full field vibration measurements performed with optical deflectometry, in: Proceedings of 30th International Conference on Noise and Vibration Engineering, ISMA 2022, Leuven, Belgium, 2022.

[24] N. Totaro, Selective structural source identification, J. Sound Vib. 420 (2018) 114–128.

[25] Y. Liu, W.S. Shepard Jr., An improved method for the reconstruction of a distributed force acting on a vibrating structure, J. Sound Vib. 291 (2006) 369–387.

[26] G. Chardon, L. Daudet, Narrowband source localization in an unknown reverberant environment using wavefield sparse decomposition, in: Proceedings of the 37th International Conference on Acoustics, Speech, and Signal Processing, ICASSP, Kyoto, Japan, 2012.

[27] M. Aucejo, Structural source identification using a generalized Tikhonov regularization, J. Sound Vib. 333 (5693–5707) (2014).

[28] M. Aucejo, O. De Smet, Bayesian source identification using local priors, Mech. Syst. Signal Process. 66–67 (2016) 120–136.

[29] M. Aucejo, O. De Smet, A multiplicative regularization for force reconstruction, Mech. Syst. Signal Process. 85 (2017) 730–745.

[30] E. Zhang, J. Antoni, P. Feissel, Bayesian force reconstruction with an uncertain model, J. Sound Vib. 331 (2012) 798–814.

[31] C. Faure, F. Ablitzer, J. Antoni, C. Pézerat, Empirical and fully Bayesian approaches for the identification of vibration sources from transverse displacement measurements, Mech. Syst. Signal Process. 94 (2017) 180–201.

[32] M. Aucejo, O. De Smet, On a full Bayesian inference for force reconstruction, Mech. Syst. Signal Process. 104 (2018) 36–59.

[33] M. Aucejo, O. De Smet, An optimal Bayesian regularization for force reconstructionproblems, Mech. Syst. Signal Process. 126 (2019) 98–115.
[34] C. Faure, J. Antoni, F. Ablitzer, C. Pézerat, A fully Bayesian approach based on Bernoulli-Gaussian prior for the identification of sparse vibratory sources from displacement measurements, J. Sound Vib. 524 (2022) 116726.

[35] A. Rezayat, V. Nassiri, B. De Pauw, J. Ertveldt, S. Vanlanduit, P. Guillaume, Identification of dynamic forces using group-sparsity in frequency domain, Mech. Syst. Signal Process. 70–71 (2016) 756–768.

[36] M. Aucejo, O. De Smet, A space-frequency multiplicative regularization for force reconstruction problems, Mech. Syst. Signal Process. 104 (2018) 1–18.

[37] J. Wambacq, K. Maes, A. Rezayat, P. Guillaume, G. Lombaert, Localization of dynamic forces on structures with an interior point method using group sparsity, Mech. Syst. Signal Process. 115 (2019) 593–606.

[38] M. Aucejo, O. De Smet, A novel algorithm for solving multiplicative mixed-norm regularization problems, Mech. Syst. Signal Process. 144 (2020) 106887.

[39] E. Lourens, E. Reynders, G. De Roeck, G. Degrande, G. Lombaert, An augmented Kalman filter for force identification in structural dynamics, Mech. Syst. Signal Process. 27 (2012) 446–460.

[40] S. Eftekhar Azam, E. Chatzi, C. Papadimitriou, A dual Kalman filter approach for state estimation via output-only acceleration measurements, Mech. Syst. Signal Process. 60–61 (2015) 866–886.

[41] O. Sedehi, C. Papadimitriou, D. Teymouri, L.S. Katafigyotis, Sequential Bayesian estimation of state and input in dynamical systems using output-only measurements, Mech. Syst. Signal Process. 131 (2019) 659–688.

[42] J. Ghibaudo, M. Aucejo, O. De Smet, Performances of a family of new sequential Bayesian filters for input estimation, Mech. Syst. Signal Process. 204 (2023) 110794.

[43] G. Enzner, P. Vary, Frequency-domain adaptive Kalman Filter for acoustic echo control in hands-free telephones, Signal Process. 86 (6) (2006) 1140–1156.

[44] F. Yang, G. Enzner, J. Yang, Frequency-domain adaptive Kalman filter with fast recovery of abrupt echo-path changes, IEEE Signal Process. Lett. 24 (12) (2017) 1778–1782.

[45] W. Fan, K. Chen, J. Lu, J. Tao, Effective improvement of under-modeling frequency-domain Kalman Filter, IEEE Signal Process. Lett. 26 (2) (2019) 342–346.

[46] T. Kabzinski, P. Jax, A unified perspective on time-domain and frequency-domain kalman filters for acoustic system identification, in: Proceedings of the 30th European Signal Processing Conference, EUSIPCO, Belgrade, Serbia, 2022.

[47] L.J.L. Nordström, A dynamic programming algorithm for input estimation on linear time-variant systems, Comput. Methods Appl. Mech. Eng. 195 (2006) 6407–6427.

[48] D. Bernal, Optimal discrete to continuous transfer for band limited inputs, J. Eng. Mech. 133 (12) (2007) 1370–1377.

[49] J. Ghibaudo, M. Aucejo, O. De Smet, A Sparse adaptive Bayesian filter for input estimation problems, Mech. Syst. Signal Process. 180 (2022) 109416.

[50] M.E. Tipping, Sparse bayesian learning and the relevance vector machine, J. Mach. Learn. Res. 1 (2001) 241–244.

[51] S. Samagassi, A. Khamlichi, A. Driouach, E. Jacquelin, Reconstruction of multiple impact forces by wavelet relevance vector machine approach, J. Sound Vib. 359 (2015) 56–67.

[52] W. Feng, Q. Li, Q. Lu, C. Li, B. Wang, Element-wise Bayesian regularization for fast and adaptive force reconstruction, J. Sound Vib. 490 (2020) 115713.

[53] M. AUCEJO, O. De Smet, Kalmanfreq.jl, Zenodo, http://dx.doi.org/10.5281/zenodo.15165899.

[54] P.C. Hansen, D.P. O'Leary, The use of the L-curve in the regularization of discrete ill-posed problems, SIAM J. Sci. Comput. 14 (6) (1993) 1487–1503.

[55] R. Mehra, On the identification of variances and adaptive Kalman filtering, IEEE Trans. Autom. Control 15 (2) (1970) 175–184.

[56] L. Zhang, D. Sidoti, A. Bienkowski, K.R. Pattipati, Y. Bar-Shalom, D.L. Kleinman, On the identification of noise covariances and adaptive Kalman filtering: A new look at a 50 year-old problem, IEEE Access 8 (2020) 59362–59388.

[57] M. Song, R. Astroza, H. Ebrahimian, B. Moaveni, C. Papadimitriou, Adaptive Kalman filters for nonlinear finite element model updating, Mech. Syst. Signal Process. 143 (2020) 106837.

[58] D. Teymouri, O. Sedehi, L.S. Katafygiotis, C. Papadimitriou, A Bayesian Expectation-Maximization (BEM) methodology for joint input-state estimation and virtual sensing of structures, Mech. Syst. Signal Process. 169 (2022) 108602.

[59] A. Pereira, J. Antoni, Q. Leclère, Empirical regularization of the inverse acoustic problem, Appl. Acoust. 97 (2015) 11–29.

[60] D. Garcia, Robust smoothing of gridded data in one and higher dimensions with missing values, Comput. Statist. Data Anal. 54 (2010) 1167–1178.

[61] G. Strang, The Discrete cosine transform, SIAM Rev. 41 (1) (1999) 135–147.

[62] W.C. Yueh, Eigenvalues of several tridiagonal matrices, Appl. Math. E-Notes 5 (2005) 66–74.

[63] R.P. Brent, Algorithms for Minimization Without Derivatives, Prentice-Hall, Englewood Cliffs, NJ, USA, 1973.

[64] R. Nayek, S. Chakraborty, S. Narasimhan, A Gaussian process latent force model for joint input-state estimation in linear structural systems, Mech. Syst. Signal Process. 128 (2019) 497–530.

[65] O. Robin, J.-D. Chazot, R. Boulandet, M. Michau, A. Berry, N. Atalla, A plane and thin Panel with representative simply supported boundary conditions for laboratory vibroacoustic tests, Acta Acust. United Acust. 102 (1) (2016) 170–182.