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.

9 Downloads

Updated 12 Oct 2011

View License

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.

Thanks for your effort

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?

Updates

 12 Oct 2011 1.7.0.0 Fixed the description to match the current version. 10 Oct 2011 1.6.0.0 Function updated to a version that uses a cell containing the kernels as the argument representing the model. 2 Oct 2011 1.5.0.0 Fixed error pointed out by Goryn. 19 Aug 2011 1.2.0.0 Corrected the submission. 19 Aug 2011 1.1.0.0 Included an example of use in the description.
MATLAB Release Compatibility
Created with R2009b
Compatible with any release
Platform Compatibility
Windows macOS Linux