FFT - again

ROOT
1

Hi,

I have another question regarding FFT. I am trying to use it as a tool for deconvolution. So if I transform the system response, s -> S, and the ouput after convolution, h -> H, I can recover the transform of the input distribution I with:

I = H / S.

So, for example, if h and s were the histograms in time space:
TH1F hTrans = new TH1F(“h in frequecy space”, “h in frequecy space”, iter, 0, iter);
hTrans = (TH1F)h->FFT(hTrans, “MAG R2C M”);
hTrans->Draw();
hTrans->SetStats(kFALSE);

Double_t H = new Double_t[2((iter)/2+1)];

TVirtualFFT *current = TVirtualFFT::GetCurrentTransform();
current->GetPoints(H);

TH1F STrans = new TH1F(“h in frequecy space”, “h in frequecy space”, iter, 0, iter);
STrans = (TH1F)s->FFT(STrans, “MAG R2C M”);
STrans->Draw();
STrans->SetStats(kFALSE);

Double_t S = new Double_t[2((iter)/2+1)];

TVirtualFFT *current = TVirtualFFT::GetCurrentTransform();
current->GetPoints(S);

But when I transform to get H and S, the imaginary parts are also included so I can not simply do :

for (int i=0; i< NoBins; i++)
I[i] = H[i] / S[i]

Infact this breaks down because some values of S[n] = 0. Do you have any ideas how this can be acheived?

Many thanks.

Ben

2

Hi Ben,

If you have complex numbers, you should divide them as complex numbers. Your input distribution, that you are trying to compute, is also complex in frequency domain.
Re(I)=(Re(H)*Re(S)+Im(H)*Im(S))/(Re(H)*Re(H)+Re(S)*Re(S));
Im(I)= (Im(H)*Re(S)-Re(H)*Im(S))/(Re(H)*Re(H)+Re(S)*Re(S));

In the output you get from TVirtualFFT, real and imaginary parts are written one after the other, so it’s (re0, im0, re1, im1,…)

If some of your S points are zero, it means that the input signal at these points can’t be reconstructed from the information you have. You could try to use some filtering, but I’m not a DSP specialist to recommend something in particular. Maybe this page http://www.dspguide.com/ch9/3.htm can give you more info (it also has good sections on using DFT).

Cheers,
Anna