File Exchange

image thumbnail

Fast Volterra Filtering

version 1.7.0.0 (2.76 KB) by José Goulart
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.

Comments and Ratings (15)

Hira

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

Gauth

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

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

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

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

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.

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

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.

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

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?

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

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.

Sorry, Goryn.

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

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

Goryn

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

Updates

1.7.0.0

Fixed the description to match the current version.

1.6.0.0

Function updated to a version that uses a cell containing the kernels as the argument representing the model.

1.5.0.0

Fixed error pointed out by Goryn.

1.2.0.0

Corrected the submission.

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