| Syntax C/C++ | #include <VFstd.h>
void VF_convolve( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size );
void VF_convolvewEdit( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size; fComplex thresh );
void VFb_convolve( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size, fVector Buf );
void VFb_convolvewEdit( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size; fComplex thresh, fVector Buf ); |
| C++ VecObj | #include <OptiVec.h>
void vector<T>::convolve( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp );
void vector<T>::convolvewEdit( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp, complex<T> thresh );
void vector<T>::b_convolve( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp, vector<T> Buf );
void vector<T>::b_convolvewEdit( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp, complex<T> thresh, vector<T> Buf ); |
| Pascal/Delphi | uses VFstd;
procedure VF_convolve( Y, Flt, X, Rsp:fVector; size:UIntSize );
procedure VF_convolvewEdit( Y, Flt, X, Rsp:fVector; size:UIntSize; thresh:fComplex );
procedure VFb_convolve( Y, Flt, X, Rsp:fVector; size:UIntSize; Buf:fVector );
procedure VFb_convolvewEdit( Y, Flt, X, Rsp:fVector; size:UIntSize; thresh:fComplex; Buf:fVector ); |
|
| 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 );
int cudaVF_convolvewEdit( fVector d_Y, fVector d_Flt, fVector d_X, fVector d_Rsp, ui size, fComplex thresh );
void VFcu_convolvewEdit( fVector h_Y, fVector h_Flt, fVector h_X, fVector h_Rsp, ui size, fComplex thresh );
|
| CUDA function Pascal/Delphi | uses 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 );
function cudaVF_convolvewEdit( d_Y, d_Flt, d_X, d_Rsp:fVector; size:UIntSize; thresh:fComplex ): IntBool;
procedure VFcu_convolvewEdit( h_Y, h_Flt, h_X, h_Rsp:fVector; size:UIntSize; thresh:fComplex );
|
|
| Description | The 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.
The result of the convolution appears scaled with the sum of all elements of Rsp. Generally, 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.
A response function which attenuates some frequences so much that any information about them is lost (e.g., high frequencies for a low-cut filter) is marked by very small values in Flt for the respective frequencies. Here, "very small" means that they are in the order of round-off error, i.e. around VF_absmax(Flt, size)*epsilon). In order to minimize round-off error, VF_convolve replaces such small values in Flt by 0. The default threshold for this editing can be retrieved by VF_setRspEdit. If you wish to set a different threshold for all calls to VF_convolve and VF_deconvolve, you may use VF_setRspEdit. However, this method is not thread-safe and should not be used to set different thresholds for different calls to VF_convolve or VF_deconvolve. Here, you have to use the variant VF_convolvewEdit, which takes the desired threshold as the argument thresh (and ignores the default threshold). As Flt consists of complex numbers, and as it is sometimes desirable to treat real and imaginary parts differently, thresh is complex as well.
Internally, VF_convolve / VF_convolvewEdit allocates and frees additional workspace memory. For repeated calls, this would be inefficient. In such a case, it is recommended to use VFb_convolve / VFb_convolvewEdit instead. The size of Buf must be ≥ 2*size. Additionally, Buf must be 128-bit (P8) or 256-bit (P9) aligned. This usually means you can only take a vector allocated by the VF_vector family as Buf. |
| 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/Delphi | 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( 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.
|
|