_ROOT Version:6.14/04
Hello everyone,
I try to use the Deconvolution of the TSpectrum class and I get results I cannot understand.
This deconvolution process is based on the gold deconvolution algorithm. You can see explanations in the following website :
Gold deconvolution process.
Deconvolution(double *source,double *response, int size, int numberIterations, int numberRepetitions)
What I understood is that we set the source in argument is the y, the response is the h and after the process, the y is transformed to be x the solution of the system.
I tried to use this to solve a convolution between a gaussian signal (the h) and a random x (that I chose). After deconvolution, I cannot get the x at all. Here is my very simple code :
void Deconv_test() {
/*** Fill the x and h ***/
// x
int N=20;
double *x = new double[N];
x[0]=1.0; x[1]=33.0; x[2]=5.0; x[3]=4.0; x[4]=10.1; x[5]=22.3; x[6]=6.0; x[7]=8.7; x[8]=9.3; x[9]=1.1;
x[10]=0.2; x[11]=0.3; x[12]=1.05; x[13]=5.7; x[14]=1.2; x[15]=0.2; x[16]=3.6; x[17]=7.0; x[18]=1.0; x[19]=0.1;
// elec
int Nr=14.0;
double sigma = 1.5;
double *elec=new double[Nr];
for(int i = 0 ; i < Nr ; i++) {
elec[i] = exp(-0.5 * (i - Nr/2) * (i - Nr/2) / sigma / sigma);
elec[i] = elec[i] / sqrt(2 * M_PI) / sigma;}
double sum=0.0;
for (int i=0 ; i < Nr ; i++) sum += elec[i];
for (int i = 0 ; i < Nr ; i++) elec[i] = elec[i] * 100. / sum;
// h (to be at the same format as the Deconvolution method) OR response !
double *h=new double[N];
for(int i=0;i<N;i++) h[i]=0.0;
for(int i=0;i<N;i++) h[i]=elec[i];
/*** Convolution ***/
double** H = new double*[2*N];
for (int i=0;i<2*N;i++) H[i]=new double[N];
for (int i = 0; i < 2*N; i++) {
for (int j = 0; j < N; j++) H[i][j] = 0.0; }
for (int j = 0; j < N; j++) {
for (int i = j; i < j+N; i++) H[i][j] = h[i-j]; }
double *source=new double[2*N]; // or y
for(int k=0;k<2*N;k++) {
source[k]=0;
for(int j=0;j<N;j++)
source[k]+=H[k][j]*x[j]; }
/*** Deconvolution ***/
TSpectrum *s = new TSpectrum();
s->Deconvolution(source,h,N,1000,1,1);
for(int i=0;i<N;i++) cout << x[i] <<" "<< source[i]<<endl;
delete[] source;
delete[] h;
delete[] x;
delete[] elec;
delete[] H;
}
I printed values of the source after the process and the values of x, they are completely different. There should be something I did not understand in the use of this method.
Thank you in advance for your answer,
Kali95