# FFT - again

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

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