Continuous time

A continuous time-invariant system with x(t) as input and y(t) as output can be expanded in the Volterra series as

${\displaystyle y(t)=h_{0}+\sum _{n=1}^{N}\int _{a}^{b}\cdots \int _{a}^{b}h_{n}(\tau _{1},\dots ,\tau _{n})\prod _{j=1}^{n}x(t-\tau _{j})\,d\tau _{j}.}$

Here the constant term ${\displaystyle h_{0}}$ on the right side is usually taken to be zero by suitable choice of output level ${\displaystyle y}$. The function ${\displaystyle h_{n}(\tau _{1},\dots ,\tau _{n})}$ is called the n-th-order Volterra kernel. It can be regarded as a higher-order impulse response of the system. For the representation to be unique, the kernels must be symmetrical in the n variables ${\displaystyle \tau }$. If it is not symmetrical, it can be replaced by a symmetrized kernel, which is the average over the n! permutations of these n variables ${\displaystyle \tau }$.

If N is finite, the series is said to be truncated. If a, b, and N are finite, the series is called doubly finite.

Sometimes the n-th-order term is divided by n!, a convention which is convenient when taking the output of one Volterra system as the input of another ("cascading").

The causality condition: Since in any physically realizable system the output can only depend on previous values of the input, the kernels ${\displaystyle h_{n}(t_{1},t_{2},\ldots ,t_{n})}$ will be zero if any of the variables ${\displaystyle t_{1},t_{2},\ldots ,t_{n}}$ are negative. The integrals may then be written over the half range from zero to infinity. So if the operator is causal, ${\displaystyle a\geq 0}$.

Fréchet's approximation theorem: The use of the Volterra series to represent a time-invariant functional relation is often justified by appealing to a theorem due to Fréchet. This theorem states that a time-invariant functional relation (satisfying certain very general conditions) can be approximated uniformly and to an arbitrary degree of precision by a sufficiently high finite-order Volterra series. Among other conditions, the set of admissible input functions ${\displaystyle x(t)}$ for which the approximation will hold is required to be compact. It is usually taken to be an equicontinuous, uniformly bounded set of functions, which is compact by the Arzelà–Ascoli theorem. In many physical situations, this assumption about the input set is a reasonable one. The theorem, however, gives no indication as to how many terms are needed for a good approximation, which is an essential question in applications.

Discrete time

The discrete-time case is similar to the continuous-time case, except that the integrals are replaced by summations:

$${\displaystyle y(n)=h_{0}+\sum _{p=1}^{P}\sum _{\tau _{1}=a}^{b}\cdots \sum _{\tau _{p}=a}^{b}h_{p}(\tau _{1},\dots ,\tau _{p})\prod _{j=1}^{p}x(n-\tau _{j}),}$$

where ${\displaystyle P\in \{1,2,\dots \}\cup \{\infty \}.}$ Each function ${\displaystyle h_{p}}$ is called a discrete-time Volterra kernels. If P is finite, the series operator is said to be truncated. If a, b and P are finite, the series operator is called a doubly finite Volterra series. If ${\displaystyle a\geq 0}$, the operator is said to be causal.

We can always consider, without loss of the generality, the kernel ${\displaystyle h_{p}(\tau _{1},\dots ,\tau _{p})}$ as symmetrical. In fact, for the commutativity of the multiplication it is always possible to symmetrize it by forming a new kernel taken as the average of the kernels for all permutations of the variables ${\displaystyle \tau _{1},\dots ,\tau _{p}}$.

For a causal system with symmetrical kernels we can rewrite the n-th term approximately in triangular form

${\displaystyle \sum _{\tau _{1}=0}^{M}\sum _{\tau _{2}=\tau _{1}}^{M}\cdots \sum _{\tau _{p}=\tau _{p-1}}^{M}h_{p}(\tau _{1},\dots ,\tau _{p})\prod _{j=1}^{p}x(n-\tau _{j}).}$

Crosscorrelation method

This method, developed by Lee and Schetzen, orthogonalizes with respect to the actual mathematical description of the signal, i.e. the projection onto the new basis functionals is based on the knowledge of the moments of the random signal.

We can write the Volterra series in terms of homogeneous operators, as

${\displaystyle y(n)=h_{0}+\sum _{p=1}^{P}H_{p}x(n),}$

where

${\displaystyle H_{p}x(n)=\sum _{\tau _{1}=a}^{b}\cdots \sum _{\tau _{p}=a}^{b}h_{p}(\tau _{1},\dots ,\tau _{p})\prod _{j=1}^{p}x(n-\tau _{j}).}$

To allow identification orthogonalization, Volterra series must be rearranged in terms of orthogonal non-homogeneous G operators (Wiener series):

${\displaystyle y(n)=\sum _{p}H_{p}x(n)\equiv \sum _{p}G_{p}x(n).}$

The G operators can be defined by the following:

${\displaystyle E\{H_{i}x(n)G_{j}x(n)\}=0;\quad i<j,}$

${\displaystyle E\{G_{i}x(n)G_{j}x(n)\}=0;\quad i\neq j,}$

whenever ${\displaystyle H_{i}x(n)}$ is arbitrary homogeneous Volterra, x(n) is some stationary white noise (SWN) with zero mean and variance A.

Recalling that every Volterra functional is orthogonal to all Wiener functional of greater order, and considering the following Volterra functional:

${\displaystyle H_{\overline {p}}^{*}x(n)=\prod _{j=1}^{\overline {p}}x(n-\tau _{j}),}$

we can write

${\displaystyle E\left\{y(n)H_{\overline {p}}^{*}x(n)\right\}=E\left\{\sum _{p=0}^{\infty }G_{p}x(n)H_{\overline {p}}^{*}x(n)\right\}.}$

If x is SWN, ${\displaystyle \tau _{1}\neq \tau _{2}\neq \ldots \neq \tau _{P}}$ and by letting ${\displaystyle A=\sigma _{x}^{2}}$, we have

${\displaystyle E\left\{y(n)\prod _{j=1}^{\overline {p}}x(n-\tau _{j})\right\}=E\left\{G_{\overline {p}}x(n)\prod _{j=1}^{\overline {p}}x(n-\tau _{j})\right\}={\overline {p}}!A^{\overline {p}}k_{\overline {p}}(\tau _{1},\dots ,\tau _{\overline {p}}).}$

So if we exclude the diagonal elements, ${\displaystyle {\tau _{i}\neq \tau _{j},\ \forall i,j}}$, it is

${\displaystyle k_{p}(\tau _{1},\dots ,\tau _{p})={\frac {E\left\{y(n)x(n-\tau _{1})\cdots x(n-\tau _{p})\right\}}{p!A^{p}}}.}$

If we want to consider the diagonal elements, the solution proposed by Lee and Schetzen is

${\displaystyle k_{p}(\tau _{1},\dots ,\tau _{p})={\frac {E\left\{\left(y(n)-\sum \limits _{m=0}^{p-1}G_{m}x(n)\right)x(n-\tau _{1})\cdots x(n-\tau _{p})\right\}}{p!A^{p}}}.}$

The main drawback of this technique is that the estimation errors, made on all elements of lower-order kernels, will affect each diagonal element of order p by means of the summation ${\displaystyle \sum \limits _{m=0}^{p-1}G_{m}x(n)}$, conceived as the solution for the estimation of the diagonal elements themselves. Efficient formulas to avoid this drawback and references for diagonal kernel element estimation exist^[6][7]

Once the Wiener kernels were identified, Volterra kernels can be obtained by using Wiener-to-Volterra formulas, in the following reported for a fifth-order Volterra series:

${\displaystyle h_{5}=k_{5},}$

${\displaystyle h_{4}=k_{4},}$

${\displaystyle h_{3}=k_{3}-10A\sum _{\tau _{4}}k_{5}(\tau _{1},\tau _{2},\tau _{3},\tau _{4},\tau _{4}),}$

${\displaystyle h_{2}=k_{2}-6A\sum _{\tau _{3}}k_{4}(\tau _{1},\tau _{2},\tau _{3},\tau _{3}),}$

${\displaystyle h_{1}=k_{1}-3A\sum _{\tau _{2}}k_{3}(\tau _{1},\tau _{2},\tau _{2})+15A^{2}\sum _{\tau 2}\sum _{\tau _{3}}k_{5}(\tau _{1},\tau _{2},\tau _{2},\tau _{3},\tau _{3}),}$

${\displaystyle h_{0}=k_{0}-A\sum _{\tau _{1}}k_{2}(\tau _{1},\tau _{1})+3A^{2}\sum _{\tau _{1}}\sum _{\tau _{2}}k_{4}(\tau _{1},\tau _{1},\tau _{2},\tau _{2}).}$

Multiple-variance method

In the traditional orthogonal algorithm, using inputs with high ${\displaystyle \sigma _{x}}$ has the advantage of stimulating high-order nonlinearity, so as to achieve more accurate high-order kernel identification. As a drawback, the use of high ${\displaystyle \sigma _{x}}$ values causes high identification error in lower-order kernels,^[8] mainly due to nonideality of the input and truncation errors.

On the contrary, the use of lower ${\displaystyle \sigma _{x}}$ in the identification process can lead to a better estimation of lower-order kernel, but can be insufficient to stimulate high-order nonlinearity.

This phenomenon, which can be called locality of truncated Volterra series, can be revealed by calculating the output error of a series as a function of different variances of input. This test can be repeated with series identified with different input variances, obtaining different curves, each with a minimum in correspondence of the variance used in the identification.

To overcome this limitation, a low ${\displaystyle \sigma _{x}}$ value should be used for the lower-order kernel and gradually increased for higher-order kernels. This is not a theoretical problem in Wiener kernel identification, since the Wiener functional are orthogonal to each other, but an appropriate normalization is needed in Wiener-to-Volterra conversion formulas for taking into account the use of different variances. Furthermore, new Wiener to Volterra conversion formulas are needed.

The traditional Wiener kernel identification should be changed as follows:^[8]

${\displaystyle k_{0}^{(0)}=E\{y^{(0)}(n)\},}$

${\displaystyle k_{1}^{(1)}(\tau _{1})={\frac {1}{A_{1}}}E\left\{y^{(1)}(n)x^{(1)}(n-\tau _{1})\right\},}$

${\displaystyle k_{2}^{(2)}(\tau _{1},\tau _{2})={\frac {1}{2!A_{2}^{2}}}\left\{E\left\{y^{(2)}(n)\prod _{i=1}^{2}x^{(2)}(n-\tau _{i})\right\}-A_{2}k_{0}^{(2)}\delta _{\tau _{1}\tau _{2}}\right\},}$

${\displaystyle k_{3}^{(3)}(\tau _{1},\tau _{2},\tau _{3})={\frac {1}{3!A_{3}^{3}}}\left\{E\left\{y^{(3)}(n)\prod _{i=1}^{3}x^{(3)}(n-\tau _{i})\right\}-A_{3}^{2}\left[k_{1}^{(3)}(\tau _{1})\delta _{\tau _{2}\tau _{3}}+k_{1}^{(3)}(\tau _{2})\delta _{\tau _{1}\tau _{3}}+k_{1}^{(3)}(\tau _{3})\delta _{\tau _{1}\tau _{2}}\right]\right\}.}$

In the above formulas the impulse functions are introduced for the identification of diagonal kernel points. If the Wiener kernels are extracted with the new formulas, the following Wiener-to-Volterra formulas (explicited up the fifth order) are needed:

${\displaystyle h_{5}=k_{5}^{(5)},}$

${\displaystyle h_{4}=k_{4}^{(4)},}$

${\displaystyle h_{3}=k_{3}^{(3)}-10A_{3}\sum _{\tau _{4}}k_{5}^{(5)}(\tau _{1},\tau _{2},\tau _{3},\tau _{4},\tau _{4}),}$

${\displaystyle h_{2}=k_{2}^{(2)}-6A_{2}\sum _{\tau _{3}}k_{4}^{(4)}(\tau _{1},\tau _{2},\tau _{3},\tau _{3}),}$

${\displaystyle h_{1}=k_{1}^{(1)}-3A_{1}\sum _{\tau _{2}}k_{3}^{(3)}(\tau _{1},\tau _{2},\tau _{2})+15A_{1}^{2}\sum _{\tau 2}\sum _{\tau _{3}}k_{5}^{(5)}(\tau _{1},\tau _{2},\tau _{2},\tau _{3},\tau _{3}),}$

${\displaystyle h_{0}=k_{0}^{(0)}-A_{0}\sum _{\tau _{1}}k_{2}^{(2)}(\tau _{1},\tau _{1})+3A_{0}^{2}\sum _{\tau _{1}}\sum _{\tau _{2}}k_{4}^{(4)}(\tau _{1},\tau _{1},\tau _{2},\tau _{2}).}$

As can be seen, the drawback with respect to the previous formula^[7] is that for the identification of the n-th-order kernel, all lower kernels must be identified again with the higher variance. However, an outstanding improvement in the output MSE will be obtained if the Wiener and Volterra kernels are obtained with the new formulas.^[8]

Feedforward network

This method was developed by Wray and Green (1994) and utilizes the fact that a simple 2-fully connected layer neural network (i.e., a multilayer perceptron) is computationally equivalent to the Volterra series and therefore contains the kernels hidden in its architecture. After such a network has been trained to successfully predict the output based on the current state and memory of the system, the kernels can then be computed from the weights and biases of that network.

The general notation for the n-th-order Volterra kernel is given by

${\displaystyle h_{n}(\tau _{1},\dots ,\tau _{n})=\sum _{i=1}^{M}(c_{i}a_{ni}\omega _{\tau _{1}i}\dots \omega _{\tau _{n}i}),}$

where ${\displaystyle n}$ is the order, ${\displaystyle c_{i}}$ the weights to the linear output node, ${\displaystyle a_{ji}}$ the coefficients of the polynomial expansion of the output function of the hidden nodes, and ${\displaystyle \omega _{ji}}$ are the weights from the input layer to the non-linear hidden layer. It is important to note that this method allows kernel extraction up until the number of input delays in the architecture of the network. Furthermore, it is vital to carefully construct the size of the network input layer so that it represents the effective memory of the system.

Exact orthogonal algorithm

This method and its more efficient version (fast orthogonal algorithm) were invented by Korenberg.^[9] In this method the orthogonalization is performed empirically over the actual input. It has been shown to perform more precisely than the crosscorrelation method. Another advantage is that arbitrary inputs can be used for the orthogonalization and that fewer data points suffice to reach a desired level of accuracy. Also, estimation can be performed incrementally until some criterion is fulfilled.

Linear regression

Linear regression is a standard tool from linear analysis. Hence, one of its main advantages is the widespread existence of standard tools for solving linear regressions efficiently. It has some educational value, since it highlights the basic property of Volterra series: linear combination of non-linear basis-functionals. For estimation, the order of the original should be known, since the Volterra basis functionals are not orthogonal, and thus estimation cannot be performed incrementally.

Kernel method

This method was invented by Franz and Schölkopf^[10] and is based on statistical learning theory. Consequently, this approach is also based on minimizing the empirical error (often called empirical risk minimization). Franz and Schölkopf proposed that the kernel method could essentially replace the Volterra series representation, although noting that the latter is more intuitive. ^[11]

Differential sampling

This method was developed by van Hemmen and coworkers^[12] and utilizes Dirac delta functions to sample the Volterra coefficients.