VF_convolveVD_convolveVE_convolve
FunctionCalculates the convolution of a vector with a response function.
Syntax C/C++#include <VFstd.h>
void VF_convolve( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size );
C++ VecObj#include <OptiVec.h>
void vector<T>::convolve( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp );
Pascal/Delphiuses VFstd;
procedure VF_convolve( Y, Flt, X, Rsp:fVector; size:UIntSize );
CUDA function C/C++#include <cudaVFstd.h>
int cudaVF_convolve( fVector d_Y, fVector d_Flt, fVector d_X, fVector d_Rsp, ui size );
void VFcu_convolve( fVector h_Y, fVector h_Flt, fVector h_X, fVector h_Rsp, ui size );
CUDA function Pascal/Delphiuses VFstd;
function cudaVF_convolve( d_Y, d_Flt, d_X, d_Rsp:fVector; size:UIntSize ): IntBool;
procedure VFcu_convolve( h_Y, h_Flt, h_X, h_Rsp:fVector; size:UIntSize );
DescriptionThe convolution of X with the response function Rsp is calculated and stored in Y. A filter Flt is also calculated. If more than one vector is to be convolved with the same Rsp, use VF_convolve only once and use VF_filter for the other vectors.
The response has to be stored in Rsp in wrap-around order: the response for zero and positive times (or whatever the independent variable is) is stored in Rsp0 to Rspsize/2 and the response for negative times (beginning with the most negative time) in Rspsize/2+1 to Rspsize-1. You may wish to use VF_rotate or VF_reflect to achieve this wrap-around order and to construct the response vector.
Notice that Rsp has to be of the same size as X.

The result of the convolution appears scaled with the sum of all elements of Rsp. Normally, therefore, Rsp should be normalized to 1.0.

X, Y, Rsp, and Flt must all be of the same size, which has to be an integer power of 2. X may be overwritten by Y, Rsp may be overwritten by Flt, but X and Flt as well as Y and Rsp have to be distinct from each other.

The treatment of round-off errors in the construction of Flt may be modified by VF_setRspEdit.

Example C/C++VF_ramp( Time, 1024, 0.0, 1.0 );
VF_Gauss( Rsp, Time, 513, 10.0, 0.0, 1.0 );
      /* Response function for zero and positive times */
VF_reflect( Rsp+1, 1023 );
      /* ... and for negative times */
VF_divC( Rsp, Rsp, 1024, VF_sum( Rsp, 1024 ) );
      /* Normalisation of Rsp */
VF_convolve( X, Rsp, X, Rsp, 1024 );
      /* Convolution; X is overwritten by the desired result
         and Rsp is overwritten by the frequency filter */
VF_filter( Y, Y, Rsp, 1024 );
      /* Next convolution: instead of another call to
         VF_convolve, Y is filtered using the frequency
         filter just obtained */
 
Example Pascal/DelphiVF_ramp( Time, 1024, 0.0, 1.0 );
VF_Gauss( Rsp, Time, 513, 10.0, 0.0, 1.0 );
      (* Response function for zero and positive times *)
VF_reflect( VF_Pelement(Rsp,1), 1023 );
      (* ... and for negative times *)
VF_divC( Rsp, Rsp, 1024, VF_sum( Rsp, 1024 ) );
      (* Normalisation of Rsp *)
VF_convolve( X, Rsp, X, Rsp, 1024 );
      (* Convolution; X is overwritten by the desired result and Rsp is overwritten by the frequency filter *)
VF_filter( Y, Y, Rsp, 1024 );
      (* Next convolution: instead of another call to VF_convolve, Y is filtered using the frequency filter just obtained   *)

Mathematically, this convolution is based on the assumption that X is periodic; it still works well if X is non-periodic but converges on both ends to the same value X0 = Xsize-1. If that is not the case, the first and the last elements of Y are spoiled by "wrap-around" from elements on the other side. Extrapolate X on both sides in order to imbed the original X in a larger vector, if wrap-around is a problem. The minimum number of elements to be added equals half the width of the response function. (In the case of an asymmetric response function, it is the broader wing that counts.) After convolving the larger vector with the response function, it will be the dummy elements just added which become spoiled by wrap-around. Those elements of the result vector which correspond to the original X will represent the desired convolution of X with Rsp.
If X is smoothly converging on both sides to different values, it is not necessary to employ the procedure just described. Rather, the difference between the end points may be regarded as a linear trend. In this case, remove the trend, convolve the resultant vector and add the trend to the result.

Example C/C++d = (X[size-1] - X[0]) / (size-1);
VF_ramp( Trend, size, 0.0, d );
VF_subV( Y, X, Trend, size );
VF_convolve( Y, Flt, Y, Rsp, size );
VF_addV( Y, Y, Trend, size );
Example for treatment of end effects in Pascal/Delphi d := (VF_element(X,size-1) - X^) / (size-1);
VF_ramp( Trend, size, 0.0, d );
VF_subV( Y, X, Trend, size );
VF_convolve( Y, Flt, Y, Rsp, size );
VF_addV( Y, Y, Trend, size );

You might notice that Flt is declared as fVector rather than cfVector, although the information stored in Flt consists of complex numbers. The reason is that these numbers are stored in the packed complex format (as described for VF_FFT) which is used only in connection with Fourier-Transform operations of real vectors.

Error handlingIf size is not a power of 2, VF_FFT (on which VF_convolve is based) complains "Size must be an integer power of 2" and the program is aborted.
Return valuenone
See alsoVF_filter,   VF_deconvolve,   VF_FFT,   VF_autocorr,   VF_xcorr,   VF_spectrum

VectorLib Table of Contents  OptiVec home