An Iterative Wiener Filter Based on a Fourth-Order Tensor Decomposition

Abstract

This work focuses on linear system identification problems in the framework of the Wiener filter. Specifically, it addresses the challenging identification of systems characterized by impulse responses of long length, which poses significant difficulties due to the existence of large parameter space. The proposed solution targets a dimensionality reduction of the problem by involving the decomposition of a fourth-order tensor, using low-rank approximations in conjunction with the nearest Kronecker product. In addition, the rank of the tensor is controlled and limited to a known value without involving any approximation technique. The final estimate is obtained based on a combination of four (shorter) optimal filters, which are alternatively iterated. As a result, the designed iterative Wiener filter outperforms the traditional counterpart, being more robust to the accuracy of the statistics’ estimates and/or noisy conditions. In addition, simulations performed in the context of acoustic echo cancellation indicate that the proposed iterative Wiener filter that exploits this fourth-order tensor decomposition achieves better performance as compared to some previously developed solutions based on lower decomposition levels. This study could further lead to the development of computationally efficient tensor-based adaptive filtering algorithms.

1. Introduction

In the context of system identification [1], the overall difficulty increases when dealing with the identification of long-length impulse responses, which raises significant challenges in terms of the computational complexity and accuracy of the solution. A well-known example is related to the acoustic echo cancellation scenarios [2], where the acoustic echo paths are modeled as finite impulse response filters that can reach thousands of coefficients. In this context, the desired performance criteria are difficult to achieve when dealing with such a large parameter space. To improve the overall performance of the identification methods when facing long-length impulse responses, it is natural to exploit the intrinsic properties of the systems to be identified. These may refer to sparseness characteristics [3,4], symmetry features [5,6], low-rank properties [7,8], among others.

Several recent approaches have focused on decomposition-based techniques that involve low-rank approximations and the nearest Kronecker product (NKP), e.g., see [9,10] (with the references therein). These methods exploit the singular value decomposition (SVD) of the corresponding matrix of the reshaped impulse response (obtained by inverse vectorization) and its low-rank feature. Therefore, the original system identification problem (that relies on a single long length filter) is reformulated using a combination of two (much) shorter filters. This could lead to important performance gains, so that the decomposition-based technique has been successfully applied in many practical applications, e.g., [11,12,13,14,15] (including their references). Nevertheless, the NKP-based technique used in these works relies on the second-order decomposition, exploiting the SVD of the (low-rank) corresponding matrix associated with the impulse response [9]. Extending this technique to a higher-order decomposition faces the difficulty of dealing with a higher-rank tensor, where its rank results from a sum of rank-1 tensors. In the context of the second-order decomposition, we rely on the fact that the rank of a matrix is limited by the number of its columns. On the other hand, for higher-order decompositions, the rank of a higher-order tensor could exceed its largest dimension; therefore, approximating its rank would be a sensitive issue [16,17,18,19,20,21].

Very recently, an alternative method to extend the NKP technique to a third-order decomposition level has been presented in [22], targeting a higher dimensionality reduction (i.e., a combination of three filters of shorter lengths). In this previous work, the rank of a third-order tensor is controlled in terms of a matrix rank. The resulting iterative Wiener filter developed in [22] has proved to outperform the traditional Wiener filter and a previously designed version that relies on the second-order NKP-based decomposition, which was proposed in [9]. Nevertheless, the extension of the third-order decomposition approach to a higher-order level cannot be performed in a straightforward manner.

Developing efficient and robust signal processing methods is an important target when dealing with applications where large amounts of data would be required but are not always available or easy to handle. In this context, tensor decomposition and modeling techniques have also been exploited in the framework of big data [23]. However, there is the so-called “curse of dimensionality” [23], which limits the order of the tensors that can be handled. This can be partially alleviated by using decompositions to represent the tensor instead of using the tensor itself. Nevertheless, most decomposition-based algorithms require full tensors, which is infeasible when dealing with many data sets. On the other hand, there are many important fields nowadays that rely on the concept of big data, such as artificial intelligence, ambient assisted living, cognitive neuroscience, machine learning, quantum information theory, and scientific computing, among others, e.g., see [24,25,26] and the references therein. Consequently, there is a high interest in developing methods that “extract” value from big data and operate with a reduced size of the data set.

The dimensionality reduction of a system identification problem with a large parameter space was initially addressed in the framework of linearly separable systems, which rely on the modeling of rank-1 tensors [27,28,29,30]. The main goal was to transform a high-dimension system identification problem into “smaller” problems (i.e., shorter filters) that can be further combined to match the original purpose. This approach is based on a Kronecker product decomposition of the global impulse response, and it is related to the processing of multilinear forms. Nevertheless, in such a framework, the global impulse response of the spatiotemporal system is considered to be perfectly separable. As a consequence, the resulting algorithms are applicable only for a particular (quasi-periodic) form of the impulse response. Nevertheless, the multilinear forms have found applications in several frameworks, such as channel equalization [31], nonlinear acoustic echo cancellation [32], multiple-input multiple-output wireless communication systems [33,34], and separable linearly constrained beamformers for large-scale antenna systems [35]. Overall, the main advantage of this decomposition-based technique that involves the modeling of rank-1 tensors is related to the dimensionality reduction (i.e., from the product of the components’ lengths to the sum of them). On the other hand, a significant limitation of this method is that it can be successfully applied only for linearly separable systems, which is not the case in most real-world applications. Therefore, the NKP-based decomposition techniques in conjunction with low-rank approximations represent better fits for the identification of realistic impulse responses of long length.

In this paper, motivated toward this direction, we develop an iterative Wiener filter using fourth-order tensor decomposition. As compared to the traditional version of the Wiener filter, which results as the solution of a single (high-dimensional) set of Wiener–Hopf equations, the advantage of this iterative technique is basically twofold. First, the final estimate is obtained by combining the solutions provided by four shorter sets of Wiener–Hopf equations, which result in “extracting” the components of the fourth-order tensor. These solutions are alternatively iterated, similar to a block coordinate descent method [36]. Such a technique fits very well when the cost functions have a partially decomposable structure, with respect to the optimization variables [37]. Second, an important feature of the proposed decomposition-based approach is that the tensor’s rank is controlled and limited to a known value without requiring any approximation technique. As a result, this decomposition method is suitable for the identification of long-length impulse responses since it operates with much smaller structures of data, as compared to the traditional Wiener solution. Moreover, it owns improved robustness to the accuracy of the statistics’ estimates and the signal-to-noise ratio (SNR). Thus, it could provide reliable solutions when only a limited amount (or incomplete set) of data samples is available and/or in low SNR conditions. In such cases, it could also outperform the previously designed iterative version that relies on third-order decomposition. The main novelty of this work with respect to the previous counterparts developed in [9,22] consists in the extension of the NKP-based decomposition to a higher-order level, which is a challenging task due to the specific difficulties associated with higher-order/higher-rank tensors, as outlined before.

Following this introduction, Section 2 presents the proposed fourth-order decomposition technique of the system impulse response. Next, the signal model for system identification and the overall problem formulation are described in Section 3. Then, Section 4 is dedicated to the development of the iterative Wiener filter based on this approach. Supporting simulation results are obtained in an acoustic echo cancellation scenario, as will be shown in Section 5. In the end, Section 6 outlines the main conclusions and presents several perspectives for future research.

The notation involved in this paper follows a common convention for data structures, using italic letters for scalars [e.g., L, $x\left(t\right)$, etc.], lowercase bold letters for vectors [e.g., $\mathbf{h}$, $\mathbf{x}\left(t\right)$, etc.], uppercase bold letters for matrices (e.g., $\mathbf{H}$, ${\mathbf{R}}_{\mathbf{x}}$, etc.), and uppercase calligraphic bold letters for tensors (e.g., $\mathcal{H}$). In addition, all the vectors are considered as column vectors.

2. Fourth-Order Decomposition of the System Impulse Response

Let us consider the framework of a linear system identification problem operating in a single-input single-output (SISO) scenario. The basic goal is to obtain a reliable estimate of the unknown impulse response, having available the input sequence and a reference signal, which is the output of the unknown system [1,5]. In this context, let

$$\begin{array}{cc}\hfill \mathbf{h}& ={\left[\begin{array}{cccc}{h}_0& h_1& \cdots & h_{L-1}\end{array}\right]}^{\top }\hfill \end{array}$$

(1)

be an impulse response with L real-valued coefficients; here, the superscript ^⊤ notation stands for the transpose operator. It is assumed that $L=L_1{L}_2$, with $L_1\ge L_2$, and $L_1$ and $L_2$ are of the same order. Then, as shown in [7,9], this system impulse response can be represented in the following form:

$$\begin{array}{cc}\hfill \mathbf{h}& =\sum _{i=1}^{L_2}{\mathbf{h}}_2^i\otimes {\mathbf{h}}_1^i,\hfill \end{array}$$

(2)

where ${\mathbf{h}}_1^i$ and ${\mathbf{h}}_2^i$ are two impulse responses of shorter lengths (as compared to $\mathbf{h}$), having $L_1$ and $L_2$ coefficients, respectively; here, ⊗ denotes the Kronecker product [38].

Next, let us consider that $L_1=L_{11}{L}_{12}$, where $L_{11}\ge L_{12}$; also, let $L_2=L_{21}{L}_{22}$, with $L_{21}\ge L_{22}$. We can apply the same decomposition in (2) to ${\mathbf{h}}_1^i$ and ${\mathbf{h}}_2^i$, i.e.,

$$\begin{array}{cc}\hfill {\mathbf{h}}_1^i& =\sum _{j=1}^{L_{12}}{\mathbf{h}}_{12}^{ij}\otimes {\mathbf{h}}_{11}^{ij},\hfill \end{array}$$

(3)

where ${\mathbf{h}}_{11}^{ij}$ and ${\mathbf{h}}_{12}^{ij}$ are two impulse responses of shorter lengths (as compared to ${\mathbf{h}}_1^i$), with $L_{11}$ and $L_{12}$ coefficients, respectively, while

$$\begin{array}{cc}\hfill {\mathbf{h}}_2^i& =\sum _{k=1}^{L_{22}}{\mathbf{h}}_{22}^{ik}\otimes {\mathbf{h}}_{21}^{ik},\hfill \end{array}$$

(4)

where ${\mathbf{h}}_{21}^{ik}$ and ${\mathbf{h}}_{22}^{ik}$ are two impulse responses of shorter lengths (as compared to ${\mathbf{h}}_2^i$), with $L_{21}$ and $L_{22}$ coefficients, respectively. Clearly, we must have $L=L_1{L}_2=L_{11}{L}_{12}{L}_{21}{L}_{22}$, with $L_1\ge L_2$, $L_{11}\ge L_{12}$, and $L_{21}\ge L_{22}$.

Substituting (3) and (4) into (2), we get

$$\begin{array}{cc}\hfill \mathbf{h}& =\sum _{i=1}^{L_2}\left(\sum _{k=1}^{L_{22}}{\mathbf{h}}_{22}^{ik}\otimes {\mathbf{h}}_{21}^{ik}\right)\otimes \left(\sum _{j=1}^{L_{12}}{\mathbf{h}}_{12}^{ij}\otimes {\mathbf{h}}_{11}^{ij}\right)\hfill \\ \hfill & =\sum _{i=1}^{L_2}\sum _{j=1}^{L_{12}}\sum _{k=1}^{L_{22}}{\mathbf{h}}_{22}^{ik}\otimes {\mathbf{h}}_{21}^{ik}\otimes {\mathbf{h}}_{12}^{ij}\otimes {\mathbf{h}}_{11}^{ij}.\hfill \end{array}$$

(5)

Now, let us assume that ${\mathbf{h}}_1^i$ and ${\mathbf{h}}_2^i$ are low rank [7], i.e.,

$$\begin{array}{cc}\hfill {\mathbf{h}}_1^i& =\sum _{p=1}^P{\mathbf{h}}_{12}^{ip}\otimes {\mathbf{h}}_{11}^{ip},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\mathbf{h}}_2^i& =\sum _{q=1}^Q{\mathbf{h}}_{22}^{iq}\otimes {\mathbf{h}}_{21}^{iq},\hfill \end{array}$$

(7)

where $P\le L_{12}$ and $Q\le L_{22}$. Then, (5) can be expressed as

$$\begin{array}{cc}\hfill \mathbf{h}& =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q{\mathbf{h}}_{22}^{lq}\otimes {\mathbf{h}}_{21}^{lq}\otimes {\mathbf{h}}_{12}^{lp}\otimes {\mathbf{h}}_{11}^{lp}.\hfill \end{array}$$

(8)

Following (2) and (8), equivalent representations in matrix and fourth-order tensor forms, respectively, result in

$$\begin{array}{cc}\hfill \mathbf{H}& =\sum _{l=1}^{L_2}{\mathbf{h}}_1^l{\left({\mathbf{h}}_2^l\right)}^{\top }\hfill \\ \hfill & =\sum _{l=1}^{L_2}{\mathbf{h}}_1^l\circ {\mathbf{h}}_2^l\hfill \end{array}$$

(9)

and

$$\begin{array}{cc}\hfill \mathcal{H}& =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q{\mathbf{h}}_{11}^{lp}\circ {\mathbf{h}}_{12}^{lp}\circ {\mathbf{h}}_{21}^{lq}\circ {\mathbf{h}}_{22}^{lq}\hfill \\ & =\sum _{p=1}^P\sum _{q=1}^Q\left(\sum _{l=1}^{L_2}{\mathbf{h}}_{11}^{lp}\circ {\mathbf{h}}_{12}^{lp}\circ {\mathbf{h}}_{21}^{lq}\circ {\mathbf{h}}_{22}^{lq}\right)\hfill \\ & =\sum _{p=1}^P\sum _{q=1}^Q{\mathcal{H}}^{:pq},\hfill \end{array}$$

(10)

where ∘ is the outer product and ${\mathcal{H}}^{:pq}={\sum }_{l=1}^{L_2}{\mathbf{h}}_{11}^{lp}\circ {\mathbf{h}}_{12}^{lp}\circ {\mathbf{h}}_{21}^{lq}\circ {\mathbf{h}}_{22}^{lq}\in {\mathbb{R}}^{L_{11}\times L_{12}\times L_{21}\times L_{22}}$. It can be noticed that ${\mathcal{H}}^{:pq}$ corresponds to the canonical polyadic (CP) decomposition, i.e., the factorization of a fourth-order tensor into a sum of fourth-order tensors of rank-1 [39]. Accordingly, $L_2$ is the minimum number that generates ${\mathcal{H}}^{:pq}$. As a result, the rank of ${\mathcal{H}}^{:pq}$ is also $L_2$, while $\mathcal{H}$ results as a double-sum of $PQ$ fourth-order tensors (each one of rank $L_2$).

In order to exploit (8) in practical system identification problems, the next task is to write this expression in four alternative ways for the purpose of extracting each block (of short impulse responses), thus leading to a more efficient identification. First, we focus on ${\mathbf{h}}_{22}^{lq}, l=1,2,\dots ,L_2, q=1,2,\dots ,Q$. We have

$$\begin{array}{cc}\hfill \mathbf{h}& =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q\left({\mathbf{I}}_{L_{22}}\otimes {\mathbf{h}}_{21}^{lq}\otimes {\mathbf{h}}_{12}^{lp}\otimes {\mathbf{h}}_{11}^{lp}\right){\mathbf{h}}_{22}^{lq}\hfill \\ & =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q{\mathbf{H}}_{-22}^{lpq}{\mathbf{h}}_{22}^{lq}\hfill \\ & =\sum _{q=1}^Q\sum _{p=1}^P{\mathbf{H}}_{-22}^{:pq}{\mathbf{h}}_{22}^{:q}\hfill \\ & =\sum _{q=1}^Q{\mathbf{H}}_{-22}^{:+:q}{\mathbf{h}}_{22}^{:q}={\underline{\mathbf{H}}}_{-22}{\underline{\mathbf{h}}}_{22},\hfill \end{array}$$

(11)

where ${\mathbf{I}}_{L_{22}}$ is the $L_{22}\times L_{22}$ identity matrix, ${\mathbf{H}}_{-22}^{lpq}={\mathbf{I}}_{L_{22}}\otimes {\mathbf{h}}_{21}^{lq}\otimes {\mathbf{h}}_{12}^{lp}\otimes {\mathbf{h}}_{11}^{lp}$ is an $L\times L_{22}$ matrix,

$$\begin{array}{cc}\hfill {\mathbf{H}}_{-22}^{:pq}& =\left[\begin{array}{cccc}{\mathbf{H}}_{-22}^{1pq}& {\mathbf{H}}_{-22}^{2pq}& \cdots & {\mathbf{H}}_{-22}^{L_{2}pq}\end{array}\right]\hfill \end{array}$$

(12)

is an $L\times L_2{L}_{22}$ matrix,

$$\begin{array}{cc}\hfill {\mathbf{h}}_{22}^{:q}& ={\left[\begin{array}{cccc}{\left({\mathbf{h}}_{22}^{1q}\right)}^{\top }& {\left({\mathbf{h}}_{22}^{2q}\right)}^{\top }& \cdots & {\left({\mathbf{h}}_{22}^{L_{2}q}\right)}^{\top }\end{array}\right]}^{\top }\hfill \end{array}$$

(13)

is a vector of length $L_2{L}_{22}$, ${\mathbf{H}}_{-22}^{:+:q}={\sum }_{p=1}^P{\mathbf{H}}_{-22}^{:pq}$ is an $L\times L_2{L}_{22}$ matrix,

$$\begin{array}{cc}\hfill {\underline{\mathbf{H}}}_{-22}& =\left[\begin{array}{cccc}{\mathbf{H}}_{-22}^{:+:1}& {\mathbf{H}}_{-22}^{:+:2}& \cdots & {\mathbf{H}}_{-22}^{:+:Q}\end{array}\right]\hfill \end{array}$$

(14)

is an $L\times Q{L}_2{L}_{22}$ matrix, and

$$\begin{array}{cc}\hfill {\underline{\mathbf{h}}}_{22}& ={\left[\begin{array}{cccc}{\left({\mathbf{h}}_{22}^{:1}\right)}^{\top }& {\left({\mathbf{h}}_{22}^{:2}\right)}^{\top }& \cdots & {\left({\mathbf{h}}_{22}^{:Q}\right)}^{\top }\end{array}\right]}^{\top }\hfill \end{array}$$

(15)

is a column vector with $Q{L}_2{L}_{22}$ coefficients.

Second, we express (8) by considering ${\mathbf{h}}_{21}^{lq}, l=1,2,\dots ,L_2, q=1,2,\dots ,Q$. We have

$$\begin{array}{cc}\hfill \mathbf{h}& =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q\left({\mathbf{h}}_{22}^{lq}\otimes {\mathbf{I}}_{L_{21}}\otimes {\mathbf{h}}_{12}^{lp}\otimes {\mathbf{h}}_{11}^{lp}\right){\mathbf{h}}_{21}^{lq}\hfill \\ \hfill & =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q{\mathbf{H}}_{-21}^{lpq}{\mathbf{h}}_{21}^{lq}\hfill \\ \hfill & =\sum _{q=1}^Q\sum _{p=1}^P{\mathbf{H}}_{-21}^{:pq}{\mathbf{h}}_{21}^{:q}\hfill \\ \hfill & =\sum _{q=1}^Q{\mathbf{H}}_{-21}^{:+:q}{\mathbf{h}}_{21}^{:q}={\underline{\mathbf{H}}}_{-21}{\underline{\mathbf{h}}}_{21},\hfill \end{array}$$

(16)

where ${\mathbf{I}}_{L_{21}}$ is the $L_{21}\times L_{21}$ identity matrix, ${\mathbf{H}}_{-21}^{lpq}={\mathbf{h}}_{22}^{lq}\otimes {\mathbf{I}}_{L_{21}}\otimes {\mathbf{h}}_{12}^{lp}\otimes {\mathbf{h}}_{11}^{lp}$ is an $L\times L_{21}$ matrix,

$$\begin{array}{cc}\hfill {\mathbf{H}}_{-21}^{:pq}& =\left[\begin{array}{cccc}{\mathbf{H}}_{-21}^{1pq}& {\mathbf{H}}_{-21}^{2pq}& \cdots & {\mathbf{H}}_{-21}^{L_{2}pq}\end{array}\right]\hfill \end{array}$$

(17)

is an $L\times L_2{L}_{21}$ matrix,

$$\begin{array}{cc}\hfill {\mathbf{h}}_{21}^{:q}& ={\left[\begin{array}{cccc}{\left({\mathbf{h}}_{21}^{1q}\right)}^{\top }& {\left({\mathbf{h}}_{21}^{2q}\right)}^{\top }& \cdots & {\left({\mathbf{h}}_{21}^{L_{2}q}\right)}^{\top }\end{array}\right]}^{\top }\hfill \end{array}$$

(18)

is a column vector with $L_2{L}_{21}$ coefficients, ${\mathbf{H}}_{-21}^{:+:q}={\sum }_{p=1}^P{\mathbf{H}}_{-21}^{:pq}$ is an $L\times L_2{L}_{21}$ matrix,

$$\begin{array}{cc}\hfill {\underline{\mathbf{H}}}_{-21}& =\left[\begin{array}{cccc}{\mathbf{H}}_{-21}^{:+:1}& {\mathbf{H}}_{-21}^{:+:2}& \cdots & {\mathbf{H}}_{-21}^{:+:Q}\end{array}\right]\hfill \end{array}$$

(19)

is an $L\times Q{L}_2{L}_{21}$ matrix, and

$$\begin{array}{cc}\hfill {\underline{\mathbf{h}}}_{21}& ={\left[\begin{array}{cccc}{\left({\mathbf{h}}_{21}^{:1}\right)}^{\top }& {\left({\mathbf{h}}_{21}^{:2}\right)}^{\top }& \cdots & {\left({\mathbf{h}}_{21}^{:Q}\right)}^{\top }\end{array}\right]}^{\top }\hfill \end{array}$$

(20)

is a vector of length $Q{L}_2{L}_{21}$.

Third, (8) can be developed by extracting ${\mathbf{h}}_{12}^{lp}, l=1,2,\dots ,L_2, p=1,2,\dots ,P$. We have

$$\begin{array}{cc}\hfill \mathbf{h}& =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q\left({\mathbf{h}}_{22}^{lq}\otimes {\mathbf{h}}_{21}^{lq}\otimes {\mathbf{I}}_{L_{12}}\otimes {\mathbf{h}}_{11}^{lp}\right){\mathbf{h}}_{12}^{lp}\hfill \\ \hfill & =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q{\mathbf{H}}_{-12}^{lpq}{\mathbf{h}}_{12}^{lp}\hfill \\ \hfill & =\sum _{p=1}^P\sum _{q=1}^Q{\mathbf{H}}_{-12}^{:pq}{\mathbf{h}}_{12}^{:p}\hfill \\ \hfill & =\sum _{p=1}^P{\mathbf{H}}_{-12}^{:p:+}{\mathbf{h}}_{12}^{:p}={\underline{\mathbf{H}}}_{-12}{\underline{\mathbf{h}}}_{12},\hfill \end{array}$$

(21)

where ${\mathbf{I}}_{L_{12}}$ is the $L_{12}\times L_{12}$ identity matrix, ${\mathbf{H}}_{-12}^{lpq}={\mathbf{h}}_{22}^{lq}\otimes {\mathbf{h}}_{21}^{lq}\otimes {\mathbf{I}}_{L_{12}}\otimes {\mathbf{h}}_{11}^{lp}$ is an $L\times L_{12}$ matrix,

$$\begin{array}{cc}\hfill {\mathbf{H}}_{-12}^{:pq}& =\left[\begin{array}{cccc}{\mathbf{H}}_{-12}^{1pq}& {\mathbf{H}}_{-12}^{2pq}& \cdots & {\mathbf{H}}_{-12}^{L_{2}pq}\end{array}\right]\hfill \end{array}$$

(22)

is an $L\times L_2{L}_{12}$ matrix,

$$\begin{array}{cc}\hfill {\mathbf{h}}_{12}^{:p}& ={\left[\begin{array}{cccc}{\left({\mathbf{h}}_{12}^{1p}\right)}^{\top }& {\left({\mathbf{h}}_{12}^{2p}\right)}^{\top }& \cdots & {\left({\mathbf{h}}_{12}^{L_{2}p}\right)}^{\top }\end{array}\right]}^{\top }\hfill \end{array}$$

(23)

is a column vector with $L_2{L}_{12}$ coefficients, ${\mathbf{H}}_{-12}^{:p:+}={\sum }_{q=1}^Q{\mathbf{H}}_{-12}^{:pq}$ is an $L\times L_2{L}_{12}$ matrix,

$$\begin{array}{cc}\hfill {\underline{\mathbf{H}}}_{-12}& =\left[\begin{array}{cccc}{\mathbf{H}}_{-12}^{:1:+}& {\mathbf{H}}_{-12}^{:2:+}& \cdots & {\mathbf{H}}_{-12}^{:P:+}\end{array}\right]\hfill \end{array}$$

(24)

is an $L\times P{L}_2{L}_{12}$ matrix, and

$$\begin{array}{cc}\hfill {\underline{\mathbf{h}}}_{12}& ={\left[\begin{array}{cccc}{\left({\mathbf{h}}_{12}^{:1}\right)}^{\top }& {\left({\mathbf{h}}_{12}^{:2}\right)}^{\top }& \cdots & {\left({\mathbf{h}}_{12}^{:P}\right)}^{\top }\end{array}\right]}^{\top }\hfill \end{array}$$

(25)

is a vector of length $P{L}_2{L}_{12}$.

Finally, we may express (8) by considering ${\mathbf{h}}_{11}^{lp}, l=1,2,\dots ,L_2, p=1,2,\dots ,P$. We have

$$\begin{array}{cc}\hfill \mathbf{h}& =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q\left({\mathbf{h}}_{22}^{lq}\otimes {\mathbf{h}}_{21}^{lq}\otimes {\mathbf{h}}_{12}^{lp}\otimes {\mathbf{I}}_{L_{11}}\right){\mathbf{h}}_{11}^{lp}\hfill \\ \hfill & =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q{\mathbf{H}}_{-11}^{lpq}{\mathbf{h}}_{11}^{lp}\hfill \\ \hfill & =\sum _{p=1}^P\sum _{q=1}^Q{\mathbf{H}}_{-11}^{:pq}{\mathbf{h}}_{11}^{:p}\hfill \\ \hfill & =\sum _{p=1}^P{\mathbf{H}}_{-11}^{:p:+}{\mathbf{h}}_{11}^{:p}={\underline{\mathbf{H}}}_{-11}{\underline{\mathbf{h}}}_{11},\hfill \end{array}$$

(26)

where ${\mathbf{I}}_{L_{11}}$ is the $L_{11}\times L_{11}$ identity matrix, ${\mathbf{H}}_{-11}^{lpq}={\mathbf{h}}_{22}^{lq}\otimes {\mathbf{h}}_{21}^{lq}\otimes {\mathbf{h}}_{12}^{lp}\otimes {\mathbf{I}}_{L_{11}}$ is an $L\times L_{11}$ matrix,

$$\begin{array}{cc}\hfill {\mathbf{H}}_{-11}^{:pq}& =\left[\begin{array}{cccc}{\mathbf{H}}_{-11}^{1pq}& {\mathbf{H}}_{-11}^{2pq}& \cdots & {\mathbf{H}}_{-11}^{L_{2}pq}\end{array}\right]\hfill \end{array}$$

(27)

is an $L\times L_2{L}_{11}$ matrix,

$$\begin{array}{cc}\hfill {\mathbf{h}}_{11}^{:p}& ={\left[\begin{array}{cccc}{\left({\mathbf{h}}_{11}^{1p}\right)}^{\top }& {\left({\mathbf{h}}_{11}^{2p}\right)}^{\top }& \cdots & {\left({\mathbf{h}}_{11}^{L_{2}p}\right)}^{\top }\end{array}\right]}^{\top }\hfill \end{array}$$

(28)

is a column vector with $L_2{L}_{11}$ coefficients, ${\mathbf{H}}_{-11}^{:p:+}={\sum }_{q=1}^Q{\mathbf{H}}_{-11}^{:pq}$ is an $L\times L_2{L}_{11}$ matrix,

$$\begin{array}{cc}\hfill {\underline{\mathbf{H}}}_{-11}& =\left[\begin{array}{cccc}{\mathbf{H}}_{-11}^{:1:+}& {\mathbf{H}}_{-11}^{:2:+}& \cdots & {\mathbf{H}}_{-11}^{:P:+}\end{array}\right]\hfill \end{array}$$

(29)

is an $L\times P{L}_2{L}_{11}$ matrix, and

$$\begin{array}{cc}\hfill {\underline{\mathbf{h}}}_{11}& ={\left[\begin{array}{cccc}{\left({\mathbf{h}}_{11}^{:1}\right)}^{\top }& {\left({\mathbf{h}}_{11}^{:2}\right)}^{\top }& \cdots & {\left({\mathbf{h}}_{11}^{:P}\right)}^{\top }\end{array}\right]}^{\top }\hfill \end{array}$$

(30)

is a vector of length $P{L}_2{L}_{11}$.

3. Signal Model for System Identification and Problem Formulation

In order to exploit the previous approach in the framework of system identification problems, let us consider a linear SISO system with the impulse response $\mathbf{h}$ of length L. At the discrete-time index t, the output of this system, denoted by $d\left(t\right)$, is given by the convolution between the samples of the input signal, $x\left(t\right)$, and the coefficients of the system impulse response from (1), usually corrupted by additive noise, i.e.,

$$\begin{array}{cc}\hfill d\left(t\right)& =\sum _{\tau =0}^{L-1}{h}_{\tau }x(t-\tau )+w\left(t\right)\hfill \\ \hfill & ={\mathbf{h}}^{\top }\mathbf{x}\left(t\right)+w\left(t\right),\hfill \end{array}$$

(31)

where

$$\begin{array}{cc}\hfill \mathbf{x}\left(t\right)& ={\left[\begin{array}{cccc}x\left(t\right)& x(t-1)& \cdots & x(t-L+1)\end{array}\right]}^{\top }\hfill \end{array}$$

(32)

is a vector that groups the most recent L time samples of the zero-mean input signal, and $w\left(t\right)$ is the zero-mean additive noise. In this framework, $x\left(t\right)$ and $w\left(t\right)$ are real-valued signals, which are uncorrelated.

From (31), we can evaluate the SNR of the SISO system. Using the notation $y\left(t\right)={\mathbf{h}}^{\top }\mathbf{x}\left(t\right)$, the SNR results in

$$\begin{array}{cc}\hfill \mathrm{SNR}& =\frac{{\sigma }_y^2}{{\sigma }_w^2}\hfill \\ \hfill & =\frac{E\left[{\mathbf{h}}^{\top }\mathbf{x}\left(t\right){\mathbf{x}}^{\top }\left(t\right)\mathbf{h}\right]}{{\sigma }_w^2}\hfill \\ \hfill & =\frac{{\mathbf{h}}^{\top }{\mathbf{R}}_{\mathbf{x}}\mathbf{h}}{{\sigma }_w^2},\hfill \end{array}$$

(33)

where ${\sigma }_y^2$ and ${\sigma }_w^2$ are the variances of $y\left(t\right)$ and $w\left(t\right)$, respectively, which are evaluated as ${\sigma }_y^2=E\left[y^2\left(t\right)\right]$ and ${\sigma }_w^2=E\left[w^2\left(t\right)\right]$, with $E[·]$ denoting mathematical expectation. In addition, ${\mathbf{R}}_{\mathbf{x}}=E\left[\mathbf{x}\left(t\right){\mathbf{x}}^{\top }\left(t\right)\right]$ represents the covariance matrix of the input signal.

The traditional Wiener filter results by solving the Wiener–Hopf equations [5], which lead to the solution

$$\begin{array}{cc}\hfill {\widehat{\mathbf{h}}}_{\mathrm{W}}& ={\mathbf{R}}_{\mathbf{x}}^{-1}{\mathbf{r}}_{\mathbf{x}d},\hfill \end{array}$$

(34)

where ${\mathbf{r}}_{\mathbf{x}d}=E\left[\mathbf{x}\left(t\right)d\left(t\right)\right]$ represents the cross-correlation vector between the input signal $\mathbf{x}\left(t\right)$ and the reference (or desired) sequence $d\left(t\right)$. In practice, this optimal filter depends on several factors, such as the accuracy of the statistics’ estimates and the SNR. In other words, there are two main issues that limit its performance: (i) a low amount of available data [i.e., samples of $x\left(t\right)$ and $d\left(t\right)$], which influence the estimates of ${\mathbf{R}}_{\mathbf{x}}$ and ${\mathbf{r}}_{\mathbf{x}d}$, and (ii) a low SNR, which also affects the estimation error. Under these circumstances, the normalized Euclidean distance between the true impulse response, $\mathbf{h}$, and the conventional Wiener filter, ${\widehat{\mathbf{h}}}_{\mathrm{W}}$, could be very high. Nevertheless, we can improve the estimation of $\mathbf{h}$ by exploiting its structure and involving the decomposition technique, as shown in the following.

Based on the decompositions presented in Section 2, the reference signal from (31) can be rewritten as

$$\begin{array}{cc}\hfill d\left(t\right)& ={\underline{\mathbf{h}}}_{22}^{\top }{\underline{\mathbf{H}}}_{-22}^{\top }\mathbf{x}\left(t\right)+w\left(t\right)\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill \hspace{1em} & ={\underline{\mathbf{h}}}_{21}^{\top }{\underline{\mathbf{H}}}_{-21}^{\top }\mathbf{x}\left(t\right)+w\left(t\right)\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \hspace{1em} & ={\underline{\mathbf{h}}}_{12}^{\top }{\underline{\mathbf{H}}}_{-12}^{\top }\mathbf{x}\left(t\right)+w\left(t\right)\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill \hspace{1em} & ={\underline{\mathbf{h}}}_{11}^{\top }{\underline{\mathbf{H}}}_{-11}^{\top }\mathbf{x}\left(t\right)+w\left(t\right).\hfill \end{array}$$

(38)

In the following, we show how to use the previous expressions in order to derive an iterative Wiener filter, which leads to better performance than the traditional Wiener filter in harsh conditions.

4. Proposed Iterative Wiener Filter for Linear System Identification

Let $\mathbf{g}$ be an estimate of $\mathbf{h}$, i.e., a linear filter of length L, with real-valued coefficients. Furthermore, we can construct ${\underline{\mathbf{g}}}_{22}$, ${\underline{\mathbf{g}}}_{21}$, ${\underline{\mathbf{g}}}_{12}$, ${\underline{\mathbf{g}}}_{11}$, and ${\underline{\mathbf{G}}}_{-22}$, ${\underline{\mathbf{G}}}_{-21}$, ${\underline{\mathbf{G}}}_{-12}$, ${\underline{\mathbf{G}}}_{-11}$, in the same way we constructed ${\underline{\mathbf{h}}}_{22}$, ${\underline{\mathbf{h}}}_{21}$, ${\underline{\mathbf{h}}}_{12}$, ${\underline{\mathbf{h}}}_{11}$, and ${\underline{\mathbf{H}}}_{-22}$, ${\underline{\mathbf{H}}}_{-21}$, ${\underline{\mathbf{H}}}_{-12}$, ${\underline{\mathbf{H}}}_{-11}$. In this context, the error signal results as the difference between the desired signal, $d\left(t\right)$, and $y_{\mathrm{E}}\left(t\right)={\mathbf{g}}^{\top }\mathbf{x}\left(t\right)$, which is the estimated signal. Consequently,

$$\begin{array}{cc}\hfill e\left(t\right)& =d\left(t\right)-y_{\mathrm{E}}\left(t\right),\hfill \end{array}$$

(39)

which can be further rewritten in four alternative manners:

$$\begin{array}{cc}\hfill e\left(t\right)& =d\left(t\right)-{\underline{\mathbf{g}}}_{22}^{\top }{\underline{\mathbf{G}}}_{-22}^{\top }\mathbf{x}\left(t\right)\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill \hspace{1em} & =d\left(t\right)-{\underline{\mathbf{g}}}_{21}^{\top }{\underline{\mathbf{G}}}_{-21}^{\top }\mathbf{x}\left(t\right)\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill \hspace{1em} & =d\left(t\right)-{\underline{\mathbf{g}}}_{12}^{\top }{\underline{\mathbf{G}}}_{-12}^{\top }\mathbf{x}\left(t\right)\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill \hspace{1em} & =d\left(t\right)-{\underline{\mathbf{g}}}_{11}^{\top }{\underline{\mathbf{G}}}_{-11}^{\top }\mathbf{x}\left(t\right).\hfill \end{array}$$

(43)

From (40) and considering that ${\underline{\mathbf{g}}}_{21}$, ${\underline{\mathbf{g}}}_{12}$, and ${\underline{\mathbf{g}}}_{11}$ are fixed, the mean-squared error (MSE) criterion can be expressed as

$$\begin{array}{cc}\hfill & \mathcal{J}\left({\underline{\mathbf{g}}}_{22}|{\underline{\mathbf{g}}}_{21},{\underline{\mathbf{g}}}_{12},{\underline{\mathbf{g}}}_{11}\right)=E\left[{\left|d\left(t\right)-{\underline{\mathbf{g}}}_{22}^{\top }{\underline{\mathbf{G}}}_{-22}^{\top }\mathbf{x}\left(t\right)\right|}^2\right]\hfill \\ \hfill & ={\sigma }_d^2-2{\underline{\mathbf{g}}}_{22}^{\top }{\underline{\mathbf{G}}}_{-22}^{\top }{\mathbf{r}}_{\mathbf{x}d}+{\underline{\mathbf{g}}}_{22}^{\top }{\underline{\mathbf{G}}}_{-22}^{\top }{\mathbf{R}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-22}{\underline{\mathbf{g}}}_{22},\hfill \end{array}$$

(44)

where ${\sigma }_d^2$ represents the variance of the desired signal, i.e., ${\sigma }_d^2=E\left[d^2\left(t\right)\right]$. The minimization of (44) leads to the Wiener filter for ${\underline{\mathbf{g}}}_{22}$, i.e.,

$$\begin{array}{cc}\hfill {\underline{\mathbf{g}}}_{22,\mathrm{W}}& ={\left({\underline{\mathbf{G}}}_{-22}^{\top }{\mathbf{R}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-22}\right)}^{-1}{\underline{\mathbf{G}}}_{-22}^{\top }{\mathbf{r}}_{\mathbf{x}d}.\hfill \end{array}$$

(45)

In a similar manner, following (41) and considering that ${\underline{\mathbf{g}}}_{22}$, ${\underline{\mathbf{g}}}_{12}$, and ${\underline{\mathbf{g}}}_{11}$ are fixed, the MSE criterion can be written as

$$\begin{array}{cc}\hfill & \mathcal{J}\left({\underline{\mathbf{g}}}_{21}|{\underline{\mathbf{g}}}_{22},{\underline{\mathbf{g}}}_{12},{\underline{\mathbf{g}}}_{11}\right)=E\left[{\left|d\left(t\right)-{\underline{\mathbf{g}}}_{21}^{\top }{\underline{\mathbf{G}}}_{-21}^{\top }\mathbf{x}\left(t\right)\right|}^2\right]\hfill \\ \hfill & ={\sigma }_d^2-2{\underline{\mathbf{g}}}_{21}^{\top }{\underline{\mathbf{G}}}_{-21}^{\top }{\mathbf{r}}_{\mathbf{x}d}+{\underline{\mathbf{g}}}_{21}^{\top }{\underline{\mathbf{G}}}_{-21}^{\top }{\mathbf{R}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-21}{\underline{\mathbf{g}}}_{21},\hfill \end{array}$$

(46)

which leads to the Wiener filter for ${\underline{\mathbf{g}}}_{21}$, i.e.,

$$\begin{array}{cc}\hfill {\underline{\mathbf{g}}}_{21,\mathrm{W}}& ={\left({\underline{\mathbf{G}}}_{-21}^{\top }{\mathbf{R}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-21}\right)}^{-1}{\underline{\mathbf{G}}}_{-21}^{\top }{\mathbf{r}}_{\mathbf{x}d}.\hfill \end{array}$$

(47)

Next, based on (42) and considering that ${\underline{\mathbf{g}}}_{22}$, ${\underline{\mathbf{g}}}_{21}$, and ${\underline{\mathbf{g}}}_{11}$ are fixed, the MSE criterion becomes

$$\begin{array}{cc}\hfill & \mathcal{J}\left({\underline{\mathbf{g}}}_{12}|{\underline{\mathbf{g}}}_{22},{\underline{\mathbf{g}}}_{21},{\underline{\mathbf{g}}}_{11}\right)=E\left[{\left|d\left(t\right)-{\underline{\mathbf{g}}}_{12}^{\top }{\underline{\mathbf{G}}}_{-12}^{\top }\mathbf{x}\left(t\right)\right|}^2\right]\hfill \\ \hfill & ={\sigma }_d^2-2{\underline{\mathbf{g}}}_{12}^{\top }{\underline{\mathbf{G}}}_{-12}^{\top }{\mathbf{r}}_{\mathbf{x}d}+{\underline{\mathbf{g}}}_{12}^{\top }{\underline{\mathbf{G}}}_{-12}^{\top }{\mathbf{R}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-12}{\underline{\mathbf{g}}}_{12}.\hfill \end{array}$$

(48)

The minimization of (48) results in the Wiener filter for ${\underline{\mathbf{g}}}_{12}$, which is

$${\underline{\mathbf{g}}}_{12,\mathrm{W}}={\left({\underline{\mathbf{G}}}_{-12}^{\top }{\mathbf{R}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-12}\right)}^{-1}{\underline{\mathbf{G}}}_{-12}^{\top }{\mathbf{r}}_{\mathbf{x}d}.$$

(49)

Finally, from (43) and considering that ${\underline{\mathbf{g}}}_{22}$, ${\underline{\mathbf{g}}}_{21}$, and ${\underline{\mathbf{g}}}_{12}$ are fixed, the MSE criterion is

$$\begin{array}{cc}\hfill & \mathcal{J}\left({\underline{\mathbf{g}}}_{11}|{\underline{\mathbf{g}}}_{22},{\underline{\mathbf{g}}}_{21},{\underline{\mathbf{g}}}_{12}\right)=E\left[{\left|d\left(t\right)-{\underline{\mathbf{g}}}_{11}^{\top }{\underline{\mathbf{G}}}_{-11}^{\top }\mathbf{x}\left(t\right)\right|}^2\right]\hfill \\ \hfill & ={\sigma }_d^2-2{\underline{\mathbf{g}}}_{11}^{\top }{\underline{\mathbf{G}}}_{-11}^{\top }{\mathbf{r}}_{\mathbf{x}d}+{\underline{\mathbf{g}}}_{11}^{\top }{\underline{\mathbf{G}}}_{-11}^{\top }{\mathbf{R}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-11}{\underline{\mathbf{g}}}_{11}.\hfill \end{array}$$

(50)

The minimization of the cost function from (50) results in the Wiener filter for ${\underline{\mathbf{g}}}_{11}$, so that

$$\begin{array}{cc}\hfill {\underline{\mathbf{g}}}_{11,\mathrm{W}}& ={\left({\underline{\mathbf{G}}}_{-11}^{\top }{\mathbf{R}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-11}\right)}^{-1}{\underline{\mathbf{G}}}_{-11}^{\top }{\mathbf{r}}_{\mathbf{x}d}.\hfill \end{array}$$

(51)

The solutions ${\underline{\mathbf{g}}}_{22,\mathrm{W}}$, ${\underline{\mathbf{g}}}_{21,\mathrm{W}}$, ${\underline{\mathbf{g}}}_{12,\mathrm{W}}$, and ${\underline{\mathbf{g}}}_{11,\mathrm{W}}$ can be alternatively iterated (as shown in [9,30]) so that the four optimal filters can be combined to finally converge to the Wiener solution after a certain number of iterations. In this way, denoting by ${\mathbf{g}}_{\mathrm{W}}$ the global optimal filter, this estimate of $\mathbf{h}$ result in

$$\begin{array}{cc}\hfill {\mathbf{g}}_{\mathrm{W}}& =\sum _{l=1}^{L_2}\sum _{p=1}^P\sum _{q=1}^Q{\mathbf{g}}_{22,\mathrm{W}}^{lq}\otimes {\mathbf{g}}_{21,\mathrm{W}}^{lq}\otimes {\mathbf{g}}_{12,\mathrm{W}}^{lp}\otimes {\mathbf{g}}_{11,\mathrm{W}}^{lp},\hfill \end{array}$$

(52)

where ${\mathbf{g}}_{22,\mathrm{W}}^{lq}$, ${\mathbf{g}}_{21,\mathrm{W}}^{lq}$, ${\mathbf{g}}_{12,\mathrm{W}}^{lp}$, and ${\mathbf{g}}_{11,\mathrm{W}}^{lp}$ are extracted from ${\underline{\mathbf{g}}}_{22,\mathrm{W}}$, ${\underline{\mathbf{g}}}_{21,\mathrm{W}}$, ${\underline{\mathbf{g}}}_{12,\mathrm{W}}$, and ${\underline{\mathbf{g}}}_{11,\mathrm{W}}$, respectively. The resulting algorithm is summarized in

Table 1. The Main Steps of the Proposed Iterative Wiener Filter.

$$\begin{array}{c}\hspace{1em}\hspace{1em} \mathit{Statistics}: {\widehat{\mathbf{R}}}_{\mathbf{x}}, {\widehat{\mathbf{r}}}_{\mathbf{x}d} [\mathrm{estimated} \mathrm{based} \mathrm{on} (53) \mathrm{and} (54), \mathrm{respectively}], \delta > 0 (\mathrm{regularization} \mathrm{constant})\hfill \\ \hspace{1em}\hspace{1em} \mathit{Setup}: L=L_1{L}_2=L_{11}{L}_{12}{L}_{21}{L}_{22}, L_1\ge L_2, L_{11}\ge L_{12}, L_{21}\ge L_{22}, P\le L_{12}, Q\le L_{22}\hfill \\ \hspace{1em}\hspace{1em} \mathit{Initialization}: {\mathbf{g}}_{11,\mathrm{W}}^{lp}={\left[\begin{array}{cc} 1& 0_{L_{11}-1}^{\top } \end{array}\right]}^{\top }, {\mathbf{g}}_{12,\mathrm{W}}^{lp}={\left[\begin{array}{cc} 1& 0_{L_{12}-1}^{\top } \end{array}\right]}^{\top }, l=1,2,\dots ,L_2, p=1,2,\dots ,P\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathbf{g}}_{21,\mathrm{W}}^{lq}={\left[\begin{array}{cc} 1& 0_{L_{21}-1}^{\top }\end{array}\right]}^{\top }, l=1,2,\dots ,L_2, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \mathit{Algorithm}\mathit{:}\hfill \\ \hspace{1em}\hspace{1em} \underline{\mathbf{Step} \mathbf{1}\mathbf{:}} {\mathbf{G}}_{-22}^{lpq}={\mathbf{I}}_{L_{22}}\otimes {\mathbf{g}}_{21,\mathrm{W}}^{lq}\otimes {\mathbf{g}}_{12,\mathrm{W}}^{lp}\otimes {\mathbf{g}}_{11,\mathrm{W}}^{lp}, l=1,2,\dots ,L_2, p=1,2,\dots ,P, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hookrightarrow {\mathbf{G}}_{-22}^{:pq}=\left[\begin{array}{cccc}{\mathbf{G}}_{-22}^{1pq}& {\mathbf{G}}_{-22}^{2pq}& \cdots & {\mathbf{G}}_{-22}^{L_{2}pq}\end{array}\right], p=1,2,\dots ,P, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hookrightarrow {\mathbf{G}}_{-22}^{:+:q}={\displaystyle \sum _{p=1}^P}{G}_{-22}^{:pq}, q=1,2,\dots ,Q ⟶ {\underline{\mathbf{G}}}_{-22}=\left[\begin{array}{cccc}{\mathbf{G}}_{-22}^{:+:1}& {\mathbf{G}}_{-22}^{:+:2}& \cdots & {\mathbf{G}}_{-22}^{:+:Q}\end{array}\right]\hfill \\ \hspace{1em}\hspace{1em} {\underline{\mathbf{g}}}_{22,\mathrm{W}}={\left({\underline{\mathbf{G}}}_{-22}^{\top }{\widehat{\mathbf{R}}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-22}+\delta {\mathbf{I}}_{Q{L}_2{L}_{22}}\right)}^{-1}{\underline{\mathbf{G}}}_{-22}^{\top }{\widehat{\mathbf{r}}}_{\mathbf{x}d}={\left[\begin{array}{cccc}{\left({\mathbf{g}}_{22,\mathrm{W}}^{:1}\right)}^{\top }& {\left({\mathbf{g}}_{22,\mathrm{W}}^{:2}\right)}^{\top }& \cdots & {\left({\mathbf{g}}_{22,\mathrm{W}}^{:Q}\right)}^{\top }\end{array}\right]}^{\top }\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hookrightarrow {\mathbf{g}}_{22,\mathrm{W}}^{:q}={\left[\begin{array}{cccc}{\left({\mathbf{g}}_{22,\mathrm{W}}^{1q}\right)}^{\top }& {\left({\mathbf{g}}_{22,\mathrm{W}}^{2q}\right)}^{\top }& \cdots & {\left({\mathbf{g}}_{22,\mathrm{W}}^{L_{2}q}\right)}^{\top }\end{array}\right]}^{\top }, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \underline{\mathbf{Step} \mathbf{2}\mathbf{:}} {\mathbf{G}}_{-21}^{lpq}={\mathbf{g}}_{22,\mathrm{W}}^{lq}\otimes {\mathbf{I}}_{L_{21}}\otimes {\mathbf{g}}_{12,\mathrm{W}}^{lp}\otimes {\mathbf{g}}_{11,\mathrm{W}}^{lp}, l=1,2,\dots ,L_2, p=1,2,\dots ,P, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hookrightarrow {\mathbf{G}}_{-21}^{:pq}=\left[\begin{array}{cccc}{\mathbf{G}}_{-21}^{1pq}& {\mathbf{G}}_{-21}^{2pq}& \cdots & {\mathbf{G}}_{-21}^{L_{2}pq}\end{array}\right], p=1,2,\dots ,P, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hookrightarrow {\mathbf{G}}_{-21}^{:+:q}={\displaystyle \sum _{p=1}^P}{\mathbf{G}}_{-21}^{:pq}, q=1,2,\dots ,Q ⟶ {\underline{\mathbf{G}}}_{-21}=\left[\begin{array}{cccc}{\mathbf{G}}_{-21}^{:+:1}& {\mathbf{G}}_{-21}^{:+:2}& \cdots & {\mathbf{G}}_{-21}^{:+:Q}\end{array}\right]\hfill \\ \hspace{1em}\hspace{1em} {\underline{\mathbf{g}}}_{21,\mathrm{W}}={\left({\underline{\mathbf{G}}}_{-21}^{\top }{\widehat{\mathbf{R}}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-21}+\delta {\mathbf{I}}_{Q{L}_2{L}_{21}}\right)}^{-1}{\underline{\mathbf{G}}}_{-21}^{\top }{\widehat{\mathbf{r}}}_{\mathbf{x}d}={\left[\begin{array}{cccc}{\left({\mathbf{g}}_{21,\mathrm{W}}^{:1}\right)}^{\top }& {\left({\mathbf{g}}_{21,\mathrm{W}}^{:2}\right)}^{\top }& \cdots & {\left({\mathbf{g}}_{21,\mathrm{W}}^{:Q}\right)}^{\top }\end{array}\right]}^{\top }\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hookrightarrow {\mathbf{g}}_{21,\mathrm{W}}^{:q}={\left[\begin{array}{cccc}{\left({\mathbf{g}}_{21,\mathrm{W}}^{1q}\right)}^{\top }& {\left({\mathbf{g}}_{21,\mathrm{W}}^{2q}\right)}^{\top }& \cdots & {\left({\mathbf{g}}_{21,\mathrm{W}}^{L_{2}q}\right)}^{\top }\end{array}\right]}^{\top }, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \underline{\mathbf{Step} \mathbf{3}\mathbf{:}} {\mathbf{G}}_{-12}^{lpq}={\mathbf{g}}_{22,\mathrm{W}}^{lq}\otimes {\mathbf{g}}_{21,\mathrm{W}}^{lq}\otimes {\mathbf{I}}_{L_{12}}\otimes {\mathbf{g}}_{11,\mathrm{W}}^{lp}, l=1,2,\dots ,L_2, p=1,2,\dots ,P, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hookrightarrow {\mathbf{G}}_{-12}^{:pq}=\left[\begin{array}{cccc}{\mathbf{G}}_{-12}^{1pq}& {\mathbf{G}}_{-12}^{2pq}& \cdots & {\mathbf{G}}_{-12}^{L_{2}pq}\end{array}\right], p=1,2,\dots ,P, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hookrightarrow {\mathbf{G}}_{-12}^{:p:+}={\displaystyle \sum _{q=1}^Q}{\mathbf{G}}_{-12}^{:pq}, p=1,2,\dots ,P ⟶ {\underline{\mathbf{G}}}_{-12}=\left[\begin{array}{cccc}{\mathbf{G}}_{-12}^{:1:+}& {\mathbf{G}}_{-12}^{:2:+}& \cdots & {\mathbf{G}}_{-12}^{:P:+}\end{array}\right]\hfill \\ \hspace{1em}\hspace{1em} {\underline{\mathbf{g}}}_{12,\mathrm{W}}={\left({\underline{\mathbf{G}}}_{-12}^{\top }{\widehat{\mathbf{R}}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-12}+\delta {\mathbf{I}}_{P{L}_2{L}_{12}}\right)}^{-1}{\underline{\mathbf{G}}}_{-12}^{\top }{\widehat{\mathbf{r}}}_{\mathbf{x}d}={\left[\begin{array}{cccc}{\left({\mathbf{g}}_{12,\mathrm{W}}^{:1}\right)}^{\top }& {\left({\mathbf{g}}_{12,\mathrm{W}}^{:2}\right)}^{\top }& \cdots & {\left({\mathbf{g}}_{12,\mathrm{W}}^{:P}\right)}^{\top }\end{array}\right]}^{\top }\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hookrightarrow {\mathbf{g}}_{12,\mathrm{W}}^{:p}={\left[\begin{array}{cccc}{\left({\mathbf{g}}_{12,\mathrm{W}}^{1p}\right)}^{\top }& {\left({\mathbf{g}}_{12,\mathrm{W}}^{2p}\right)}^{\top }& \cdots & {\left({\mathbf{g}}_{12,\mathrm{W}}^{L_{2}p}\right)}^{\top }\end{array}\right]}^{\top },p=1,2,\dots ,P\hfill \\ \hspace{1em}\hspace{1em} \underline{\mathbf{Step} \mathbf{4}\mathbf{:}} {\mathbf{G}}_{-11}^{lpq}={\mathbf{g}}_{22,\mathrm{W}}^{lq}\otimes {\mathbf{g}}_{21,\mathrm{W}}^{lq}\otimes {\mathbf{g}}_{12,\mathrm{W}}^{lp}\otimes {\mathbf{I}}_{L_{11}}, l=1,2,\dots ,L_2, p=1,2,\dots ,P, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hookrightarrow {\mathbf{G}}_{-11}^{:pq}=\left[\begin{array}{cccc}{\mathbf{G}}_{-11}^{1pq}& {\mathbf{G}}_{-11}^{2pq}& \cdots & {\mathbf{G}}_{-11}^{L_{2}pq}\end{array}\right], p=1,2,\dots ,P, q=1,2,\dots ,Q\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hookrightarrow {\mathbf{G}}_{-11}^{:p:+}={\displaystyle \sum _{q=1}^Q}{\mathbf{G}}_{-11}^{:pq}, p=1,2,\dots ,P ⟶ {\underline{\mathbf{G}}}_{-11}=\left[\begin{array}{cccc}{\mathbf{G}}_{-11}^{:1:+}& {\mathbf{G}}_{-11}^{:2:+}& \cdots & {\mathbf{G}}_{-11}^{:P:+}\end{array}\right]\hfill \\ \hspace{1em}\hspace{1em} {\underline{\mathbf{g}}}_{11,\mathrm{W}}={\left({\underline{\mathbf{G}}}_{-11}^{\top }{\widehat{\mathbf{R}}}_{\mathbf{x}}{\underline{\mathbf{G}}}_{-11}+\delta {\mathbf{I}}_{P{L}_2{L}_{11}}\right)}^{-1}{\underline{\mathbf{G}}}_{-11}^{\top }{\widehat{\mathbf{r}}}_{\mathbf{x}d}={\left[\begin{array}{cccc}{\left({\mathbf{g}}_{11,\mathrm{W}}^{:1}\right)}^{\top }& {\left({\mathbf{g}}_{11,\mathrm{W}}^{:2}\right)}^{\top }& \cdots & {\left({\mathbf{g}}_{11,\mathrm{W}}^{:P}\right)}^{\top }\end{array}\right]}^{\top }\hfill \\ \hspace{1em}\hspace{1em} \hspace{1em}\hookrightarrow {\mathbf{g}}_{11,\mathrm{W}}^{:p}={\left[\begin{array}{cccc}{\left({\mathbf{g}}_{11,\mathrm{W}}^{1p}\right)}^{\top }& {\left({\mathbf{g}}_{11,\mathrm{W}}^{2p}\right)}^{\top }& \cdots & {\left({\mathbf{g}}_{11,\mathrm{W}}^{L_{2}p}\right)}^{\top }\end{array}\right]}^{\top }, p=1,2,\dots ,P\hfill \\ \hspace{1em}\hspace{1em} \underline{\mathbf{Optimal} \mathbf{filter}\mathbf{:}} {\mathbf{g}}_{\mathrm{W}}={\displaystyle \sum _{l=1}^{L_2}}{\displaystyle \sum _{p=1}^P}{\displaystyle \sum _{q=1}^Q}{\mathbf{g}}_{22,\mathrm{W}}^{lq}\otimes {\mathbf{g}}_{21,\mathrm{W}}^{lq}\otimes {\mathbf{g}}_{12,\mathrm{W}}^{lp}\otimes {\mathbf{g}}_{11,\mathrm{W}}^{lp}\hfill \end{array}$$

, in a step-by-step manner. In practice, a regularization constant ($\delta >0$) should be added to the diagonal of the matrices to be inverted in (45), (47), (49), and (51), in order to prevent any potential numerical issue related to the inversion of badly-scaled or ill-conditioned matrices. Additional insights about different regularization techniques that could be used for the Wiener filter can be found in [40]. However, the analysis of the influence of the regularization parameter on the overall performance of the proposed iterative Wiener filter is not relevant for the purpose of this paper.

5. Simulation Results

The experimental framework considered in this section is acoustic echo cancellation [2]. In this context, the main goal is to identify the acoustic impulse response that results due to the coupling between the loudspeaker and microphone (e.g., when using a hands-free communication device). This system is modeled using a long-length finite impulse response filter with thousands of coefficients. The excessive length in time is due to the slow speed of sound in the air, together with multiple reflections caused by the environment. In our simulations, we use a measured acoustic impulse response $\mathbf{h}$ with $L=4096$ coefficients, while the sampling rate is 16 kHz. The fourth-order tensor decomposition of this impulse response is performed using $L_1=L_2=64$ and $L_{11}=L_{12}=L_{21}=L_{22}=8$.

The acoustic impulse response is depicted in Figure 1a. In addition, the singular values (normalized to the largest one) corresponding to the SVD of the matrix $\mathbf{H}$ from (9) are shown in Figure 1b, while their logarithmic scale representation is provided in Figure 1c. These plots illustrate the low-rank property of the acoustic impulse response since most of the singular values are close to zero. The low-rank property can also be noticed in the case of the component impulse responses from (2), indicated in Figure 1d–g. In these subfigures, the singular values corresponding to the SVD of the matrix associated to ${\mathbf{h}}_2^i$ (of size $L_{21}\times L_{22}$) are shown, for different values of i (1, $L_2/4$, $L_2/2$, and $L_2$).

The input signal $x\left(t\right)$ is a first-order autoregressive process. Such a sequence is the result of filtering a white Gaussian noise through a system with the pole at 0.9. The output of the echo path, $y\left(t\right)$, is corrupted by white Gaussian noise, $w\left(t\right)$, with different SNRs, thus resulting in the desired signal $d\left(t\right)$ from (31). The estimates of the required statistics are obtained by averaging across N data samples of $x\left(t\right)$ and $d\left(t\right)$. Consequently, they result in

$$\begin{array}{c}\hfill {\widehat{\mathbf{R}}}_{\mathbf{x}}=\sum _{t=1}^N\mathbf{x}\left(t\right){\mathbf{x}}^{\top }\left(t\right),\end{array}$$

(53)

$$\begin{array}{c}\hfill {\widehat{\mathbf{r}}}_{\mathbf{x}d}=\sum _{t=1}^N\mathbf{x}\left(t\right)d\left(t\right),\end{array}$$

(54)

which are used instead of ${\mathbf{R}}_{\mathbf{x}}$ and ${\mathbf{r}}_{\mathbf{x}d}$, respectively. The accuracy of these estimates depends on the value of N. In the following, we consider that $N=KL$, with $K>0$. A large value of K leads to a reliable set of statistics estimates. However, in practice, there are various scenarios in which we deal with a low quantity of data, which affects the quality of the estimates from (53) and (54). In this case, the accuracy of the Wiener filter is also affected. The statistics’ estimates are also influenced by the SNR, so it is challenging to obtain a reliable solution for the Wiener filter in noisy environments. Nevertheless, since the proposed iterative Wiener filter involves smaller structures of data, it is expected to be more robust against such adverse conditions.

The performance measure used in all the experiments is the normalized misalignment, which basically relies on an Euclidian distance. It is evaluated (in dB) as $20{\log }_{10}\left[∥\mathbf{h}-\widehat{\mathbf{h}}∥/∥\mathbf{h}∥\right]$, where $∥·∥$ is the Euclidean norm and $\widehat{\mathbf{h}}$ generally denotes an estimate of $\mathbf{h}$.

The first set of results illustrate the comparison between the conventional version of the Wiener filter (WF) and the designed iterative Wiener filter that exploits a fourth-order tensor decomposition, which is referred as WF-FOT. Different amounts of data samples are used (i.e., $N=KL$), by varying the value of K. In addition, different SNR conditions are considered. The conventional WF is obtained based on (34), i.e., ${\widehat{\mathbf{h}}}_{\mathrm{W}}$, as the solution vector (of length $L=4096$) that results by directly solving the Wiener–Hopf equations [5]. On the other hand, the estimate provided by the WF-FOT, i.e., ${\mathbf{g}}_{\mathrm{W}}$, results based on (52), using a combination of shorter sets of coefficients.

In Figure 2, the settings are $K=5$ and $\mathrm{SNR}=20$ dB, which represent good conditions for obtaining a reliable set of estimates from (53) and (54). Thus, the conventional WF provides a reasonable attenuation of the misalignment. Nevertheless, the WF-FOT is able to reach lower misalignment levels, for different settings of its decomposition parameters, P and Q, which are much smaller as compared to their maximum limits, $L_{12}$ and $L_{22}$, respectively. Moreover, the solution obtained for $P=Q=1$ provides not only a good attenuation of the misalignment, but also the fewest sets of coefficients to be combined based on (52).

The performance gain of the WF-FOT is more apparent when a smaller amount of data can be used for the statistics’ estimates from (53) and (54). Such a scenario is considered in Figure 3, where $K=2$. As we can notice, the difference between the proposed WF-FOT and the conventional WF is more significant in this case.

The case $K=1$ (or $N=L$) represents a difficult scenario for the conventional WF, due to the limited amount of data available for the estimation of the statistics. This can be noticed in Figure 4, where the conventional WF cannot provide a reliable estimate, while the proposed WF-FOT is still able to achieve some attenuation of the misalignment.

The SNR is also an important factor that influences the accuracy of the Wiener solution. This aspect is supported in Figure 5, where a sufficient amount of data is available (i.e., $N=5L$), but $\mathrm{SNR}=10$ dB. The noisy conditions considered in this simulation lead to higher misalignment levels for all the algorithms. Nevertheless, the proposed WF-FOT still outperforms its conventional counterpart, especially for the case when $P=Q=1$, which is also advantageous in terms of the parameter space (i.e., the lowest number of coefficients to be estimated).

The second set of experiments focuses on the comparison between the proposed WF-FOT and two previously developed solutions that also exploit the impulse response decomposition. First, we consider the iterative Wiener filter using the nearest Kronecker product (namely WF-NKP), which was presented in [9], in the framework of a second-order decomposition level. In other words, if an impulse response $\mathbf{h}$ of length $L=L_1^{\prime }{L}_2^{\prime }$ (with $L_2^{\prime }\le L_1^{\prime }$) is reshaped as a matrix $\mathbf{H}$ of size $L_1^{\prime }\times L_2^{\prime }$ and the rank of this corresponding matrix is $P^{\prime }<L_2^{\prime }$, then the system to be identified owns the low-rank property. Hence, using the NKP, the impulse response is expressed as $\mathbf{h}={\sum }_{p^{\prime }=1}^{P^{\prime }}{\mathbf{h}}_{2,p^{\prime }}\otimes {\mathbf{h}}_{1,p^{\prime }}$, where ${\mathbf{h}}_{1,p^{\prime }}$ and ${\mathbf{h}}_{2,p^{\prime }}$, with $p^{\prime }=1,2,\dots ,P^{\prime }$, are (shorter) impulse responses of lengths $L_1^{\prime }$ and $L_2^{\prime }$, respectively. Basically, the WF-NKP operates with two filters of lengths $P^{\prime }{L}_1^{\prime }$ and $P^{\prime }{L}_2^{\prime }$, which group the sets of coefficients ${\mathbf{h}}_{1,p^{\prime }}$ and ${\mathbf{h}}_{2,p^{\prime }}$, with $p^{\prime }=1,2,\dots ,P^{\prime }$. In our experimental scenario (with $L=L_1^{\prime }{L}_2^{\prime }=4096$), the decomposition setup of the WF-NKP is $L_1^{\prime }=L_2^{\prime }=64$ and $P^{\prime }<L_2^{\prime }$. For example, according to Figure 1b,c, we can choose $P^{\prime }=20$. Second, the iterative Wiener filter that relies on a third-order tensor decomposition (referred as WF-TOT) from [22] is considered in comparisons. In this case, the length of impulse response is factorized as $L=L_{11}^{″}{L}_{12}^{″}{L}_2^{″}=4096$, with $L_{11}^{″}{L}_{12}^{″}\gg L_2^{″}$, so that the decomposition setup of the WF-TOT is $L_{11}^{″}=L_{12}^{″}=32$, $L_2^{″}=4$, and $P^{″}<L_{12}^{″}$. For the identification of acoustic impulse responses, $P^{″}$ can be chosen in the range $L_{12}^{″}/3$ to $L_{12}^{″}/2$. In the current setup, we choose $P^{″}=10,12,$ or 14, as indicated in the captions of the following figures.

In Figure 6, the number of available data samples is $N=5L$ and $\mathrm{SNR}=20$ dB. It can be noticed that the three comparing algorithms achieve similar performance in terms of their accuracy. However, the WF-FOT reaches the lowest misalignment level among its counterparts, with a slightly slower convergence rate.

Next, in Figure 7, the number of available data samples is reduced to $N=2L$, while keeping the same SNR conditions. As we can notice, the gain of the WF-FOT becomes more apparent in this scenario, especially as compared to WF-NKP.

A similar behavior can be observed in Figure 8, where $N=5L$ and $\mathrm{SNR}=10$ dB. The lower SNR level influences the performances of all the algorithms. Nevertheless, the WF-FOT still outperforms both its counterparts. In the current setup, the WF-NKP [9] involves two filters of lengths $P^{\prime }{L}_1^{\prime }$ and $P^{\prime }{L}_1^{\prime }$, so that a total of $2\times 20\times 64=2560$ coefficients. In addition, the WF-TOT [22] uses three filters of lengths $P^{″}{L}_2^{″}{L}_{11}^{″}$, $P^{″}{L}_2^{″}{L}_{11}^{″}$, and $L_2^{″2}$, which result in $2\times 4\times 10\times 32+4^2=2576$ coefficients. The smallest parameter space is obtained by the WF-FOT, which requires four filters of lengths $P{L}_2{L}_{11}$, $P{L}_2{L}_{12}$, $Q{L}_2{L}_{21}$, and $Q{L}_2{L}_{22}$, i.e., a total of $4\times 8\times 64=2048$ coefficients. Finally, we should outline that all the decomposition-based iterative Wiener filters involve a much smaller parameter space as compared to the conventional WF, which needs to estimate $L=4096$ coefficients, corresponding to the full length of the impulse response.

6. Conclusions and Perspectives

In this paper, we have designed an iterative Wiener filter based on a fourth-order tensor decomposition associated with the impulse response. Due to the dimensionality reduction, the proposed WF-FOT is suitable for system identification problems with a large parameter space since it combines the estimates provided by four shorter filters. The algorithm has been tested in an acoustic echo cancellation scenario, targeting the identification of a long-length acoustic impulse response with $L=4096$ coefficients. Based on the obtained results, the gain is basically twofold. First, since the WF-FOT operates with smaller data structures, it shows improved robustness to the accuracy of the statistics’ estimates and in very noisy environments, compared to the conventional Wiener filter. As indicated by the experiments, it can provide a reliable solution when only a limited amount of data (e.g., N close to L) is available to estimate the statistics. Second, due to the higher-order decomposition, the WF-FOT outperforms the previously designed iterative Wiener filters that exploit the second-order and third-order decomposition levels, i.e., WF-NKP [9] and WF-TOT [22], respectively. In addition, for the reported settings, it also owns the smallest parameter space (i.e., $L/2$ coefficients) among the decomposition-based algorithms involved in comparisons. Another important advantage is related to the rank of the fourth-order tensor, which is not obtained by any approximation technique, but it is controlled and limited to a known value.

Finally, we should note that there are some inherent limitations of the Wiener filter, even in this iterative form. First, it relies on the statistics’ estimates, which are challenging to evaluate, especially in real-time applications. Second, it is based on the stationary assumption, which is biased in time-varying environments. Third, it involves matrix inversion operations to solve the associated Wiener–Hopf equations, which could raise numerical and complexity issues. Consequently, in perspective, it would be beneficial to explore two main directions for future research. First, we aim to design adaptive filtering algorithms following the fourth-order tensor decomposition approach. Such adaptive algorithms would address two of the previously mentioned limitations concerning the statistics’ estimates (which would not be required anymore) and the operation in nonstationary conditions. Second, to address the matrix inversion issues, the proposed WF-FOT can be reformulated using efficient iterative methods for solving the four sets of Wiener–Hopf equations, such as the conjugate gradient method and different coordinate descent techniques [41].