File Exchange

## Fast Volterra Filtering

version 1.7.0.0 (2.76 KB) by José Goulart

### José Goulart (view profile)

Fast algorithm for computing the response of a discrete Volterra model given a input sequence.

Updated 12 Oct 2011

Implementation of the algorithm proposed in the paper:

Morhac, M., 1991 - A Fast Algorithm of Nonlinear Volterra Filtering

Basically, the output of each n-homogeneous subsystem is decomposed into a summation of convolutions. Each of these convolutions is performed in the frequency domain, which improves the efficiency of the computation.

The Volterra model must be provided to the algorithm as a structure containing the kernels in the triangular form, i.e., such that

h(t1,...,tn) ~= 0 only if t1<=t2<=...<=tn

A example is given below:

% Building the model
kernels = {k1 k2 k3 ... kP};

% Defining the memory extension of
% each kernel (in samples)
memories = [N1 N2 N3 ... NP];

% Given an input u, the output
% is computed with
y = fastVMcell(u, kernels, memories);

% y is a P x M vector where the p-th
% row is the output of the
% p-th homogeneous nonlinear system,
% i.e., the multidimensional convolution
% of the input with the kernel kp.

% To get the overall output of the
% Volterra filter, use
yP = sum(y,1);

Note: The code provided *does not* compute Volterra kernels for a model. It only computes the output of a Volterra filter, given its kernels and the input.

### Cite As

José Goulart (2020). Fast Volterra Filtering (https://www.mathworks.com/matlabcentral/fileexchange/32248-fast-volterra-filtering), MATLAB Central File Exchange. Retrieved .

Hira

### Hira (view profile)

Can anyone guide on how to evaluate volterra kernels of an ADC please ?

Gauth

### Gauth (view profile)

Very useful routine. Just what I needed,
saved me a lot of work.

José Goulart

### José Goulart (view profile)

Maher,

This is because the kernels in the paper are assumed to be symmetric. So, this constant basically multiplies each term by the number of all possible equivalent permutations of its kernels indexes. Since I am supposing that the kernels provided to the function are in the triangular form, this constant is not necessary. In fact, to convert a symmetric kernel to a triangular kernel, one simply multiplies the terms with indexes such that k1<=k2<=...<=kp by their respective S constants.

Maher

### Maher (view profile)

Goulart,

I have read the paper which you implemented and also I have followed your script. I am wonder about the constant ( S in the paper) which multiple by the kernel of Volterra series. I didn't notice this constant in the script.

José Goulart

### José Goulart (view profile)

Goryn,

Each row p of y corresponds to the output of the p-th order homogeneous system output.

That is, the p-th row is the multidimensional convolution of the input with the kernel kp.

About the zeros in the result, this is just as expected, because your kernels are all null. To get a nonzero result, you must provide the correct kernels of your Volterra filter.

Goryn

### Goryn (view profile)

Ok. Now I use:

u = data; %1000x1 matrix

N1 = 100;
k1 = zeros(1,N1);
N2 = 100;
k2 = zeros(N2,N2);
N3 = 100;
k3 = zeros(N3,N3,N3);

kernels = {k1 k2 k3};
memories = [N1 N2 N3];

y = fastVMcell(u, kernels, memories);

Now it works! But the result 'y' is 3x1099 matrix with zeros. No filtered data in output, just three rows of zeros.

José Goulart

### José Goulart (view profile)

Goryn,

In fact, the previous version was still wrong. I am really sorry, I have been using a new version lately and haven't still updated the files because I was having no time to do it.

In this new version, which I have just uploaded (and will typically take one day to be available for download) the model is specified to the function as a cell containing tensors that correspond to the kernels.

This update adds flexibility, making it possible to use a model where the kernels have different memory lengths. Also, i have prepared the function to use the parallel processing toolbox if its available. Finally, it is now simpler to use.

So, about your question, let us suppose you have three kernels: k1 being a 1xN1 vector, k2 a N2xN2 matrix and k3 a N3xN3xN3 tensor. Also, let us assume your input is in a vector called u. Then you should do the following:

% Building the model
kernels = {k1 k2 k3};
memories = [N1 N2 N3];

% Computing the output
y = fastVMcell(u, kernels, memories);

Goryn

### Goryn (view profile)

Now it gives:

??? Error using ==> plus
Matrix dimensions must agree.

Error in ==> fastVM at 61
y = y + frconv(hs,us);

Jose,
could you please provide an example code how to filter 1000x1 matrix of nonstationary chatic datasample from scratch. Thank you.

José Goulart

### José Goulart (view profile)

Goryn, there's nothing wrong with your code.

I'm sorry, this function (mysub2ind) was missing in the previous submission. I have just updated the file, in a few minutes it should be available for download. Thank you.

Goryn

### Goryn (view profile)

Ok. That's the code that I use:

% data - 1000x1 matrix that I want to filter

u = data;

k1 = zeros(100);
k2 = zeros(100,100);
kernels.k(1).k = k1;
kernels.k(2).k = k2;
kernels.orders = [1 2];

y = fastVM(u,kernels)

ERROR:

??? Undefined function or method 'mysub2ind' for input arguments of type 'double'.

Error in ==> fastVM at 58
hinds = mysub2ind(siz,hinds+1);

Whats wrong with the code?

José Goulart

### José Goulart (view profile)

Each kp, p=1,...,P, must be a p-dimensional tensor with size M in each dimension.

For example, to build a kernel with a memory of 100 samples, then:

k1 = zeros(100);
% compute elements of k1
k2 = zeros(100,100);
% compute elements of k2
% and so on.

Goryn

### Goryn (view profile)

Thank you, Jose. Is there any way to do something like this:
I just have a data (chaotic, nonstationary) that I whant to filter. I just copy an example code where U is my data and Y is my filtered result. All other params like k1,k2,...,kP, M, P and all others are predefined in example code. To see the result of filter.

José Goulart

### José Goulart (view profile)

Sorry, Goryn.

I included an example of use in the description. Please tell me if there is any additional detail that needs explanation.

José Goulart

### José Goulart (view profile)

That is, i included an example of use in the description of the file.

Goryn

### Goryn (view profile)

How to use it? Could anyone give an example of code using this filter?