| VF_deconvolve | VD_deconvolve | VE_deconvolve |
| VF_deconvolvewEdit | VD_deconvolvewEdit | VE_deconvolvewEdit |
| VFb_deconvolve | VDb_deconvolve | VEb_deconvolve |
| VFb_deconvolvewEdit | VDb_deconvolvewEdit | VEb_deconvolvewEdit |
|
|
| Syntax C/C++ | #include <VFstd.h>
void VF_deconvolve( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size );
void VF_deconvolvewEdit( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size; fComplex thresh );
void VFb_deconvolve( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size, fVector Buf );
void VFb_deconvolvewEdit( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size; fComplex thresh, fVector Buf ); |
| C++ VecObj | #include <OptiVec.h>
void vector<T>::deconvolve( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp );
void vector<T>::deconvolvewEdit( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp, complex<T> thresh );
void vector<T>::b_deconvolve( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp, vector<T> Buf );
void vector<T>::b_deconvolvewEdit( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp, complex<T> thresh, vector<T> Buf ); |
| Pascal/Delphi | uses VFstd;
procedure VF_deconvolve( Y, Flt, X, Rsp:fVector; size:UIntSize );
procedure VF_deconvolvewEdit( Y, Flt, X, Rsp:fVector; size:UIntSize; thresh:fComplex );
procedure VFb_deconvolve( Y, Flt, X, Rsp:fVector; size:UIntSize; Buf:fVector );
procedure VFb_deconvolvewEdit( Y, Flt, X, Rsp:fVector; size:UIntSize; thresh:fComplex; Buf:fVector ); |
|
| CUDA function C/C++ | #include <cudaVFstd.h>
int cudaVF_deconvolve( fVector d_Y, fVector d_Flt, fVector d_X, fVector d_Rsp, ui size );
void VFcu_deconvolve( fVector h_Y, fVector h_Flt, fVector h_X, fVector h_Rsp, ui size );
int cudaVF_deconvolvewEdit( fVector d_Y, fVector d_Flt, fVector d_X, fVector d_Rsp, ui size, fComplex thresh );
void VFcu_deconvolvewEdit( fVector h_Y, fVector h_Flt, fVector h_X, fVector h_Rsp, ui size, fComplex thresh );
|
| CUDA function Pascal/Delphi | uses VFstd;
function cudaVF_deconvolve( d_Y, d_Flt, d_X, d_Rsp:fVector; size:UIntSize ): IntBool;
procedure VFcu_deconvolve( h_Y, h_Flt, h_X, h_Rsp:fVector; size:UIntSize );
function cudaVF_deconvolvewEdit( d_Y, d_Flt, d_X, d_Rsp:fVector; size:UIntSize; thresh:fComplex ): IntBool;
procedure VFcu_deconvolvewEdit( h_Y, h_Flt, h_X, h_Rsp:fVector; size:UIntSize; thresh:fComplex );
|
|
| Description | X is assumed to be the result of a convolution of some "true" profile with the response function Rsp; a deconvolution is attempted and stored in Y. A filter Flt is also calculated; if more than one vector is to be deconvolved with the same Rsp, use VF_deconvolve only once and use the filter Flt thus obtained to deconvolve other vectors by calling VF_filter. The response has to be stored in Rsp in wrap-around order: the elements for zero and positive times (or whatever the independent variable is) are stored as Rsp0 to Rspsize/2 and the elements for negative times as Rspsize/2+1 to Rspsize−1.
You may wish to use VF_rotate or VF_reflect to obtain the correct order when constructing the response vector.
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.
Mathematically, Flt is the inverse of the Fourier transform of Rsp. If the Fourier transform of Rsp contains elements equal to zero, all information is lost for the respective frequency and no reconstruction is possible. The best one can do in this case is to accept this loss and to deconvolve only up to those frequencies where still something is left to be reconstructed.
You are therefore advised not to use this function blindly but rather to inspect the Fourier transform of Rsp and decide what to do on the basis of your specific application. If you wish to use this function nevertheless, you may rely on the automatic editing of the filter, built into VF_deconvolve. Thereby, Flt is set to zero (instead of infinity) at those frequences where all information has been lost.
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_deconvolve and VF_convolve, 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_deconvolve or VF_convolve. Here, you have to use the variant VF_deconvolvewEdit, 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.
This deconvolution is based on the implicit assumption that X is periodic; if this is not the case, see the description of VF_convolve about how to avoid end effects.
Internally, VF_deconvolve / VF_deconvolvewEdit allocates and frees additional workspace memory. For repeated calls, this would be inefficient. In such a case, it is recommended to use VFb_deconvolve / VFb_deconvolvewEdit 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.
|
|
| Error handling | If size is not a power of 2, VF_FFT (on which VF_deconvolve is based) complains "Size must be an integer power of 2" and the program is aborted.
If, by VF_setRspEdit, you specified Trunc.Re = Trunc.Im = 0, or if you call VF_deconvolveEdit with thresh = fcplx(0,0), SING errors may occur that are handled by setting Flt to ±HUGE_VAL at the respective frequency. During multiplication with the transform of X, this may lead to unhandled floating-point overflow errors (in case your guess of Rsp was wrong and there is some information left at the frequencies where you thought it was not). |
|
|
|