My question is regarding performing operations with fourier transforms, say multiplying them, and then taking the inverse of the result. I can take the inverse Fourier transform of a transform but I am not sure how to do it after I modify it.
//here h2 and h3 are some previously defined TH1's
TH1 *hm2 =0;
TVirtualFFT::SetTransform(0);
hm2 = h2->FFT(hm2, "MAG");
TH1 *hm3 =0;
TVirtualFFT::SetTransform(0);
hm3 = h3->FFT(hm3, "MAG");
TH1 *hCombo = 0;
hCombo = hm2;
hCombo->Multiply(hm3);
TVirtualFFT::SetTransform(0);
//This is where my issue is!
//TVirtualFFT *fftC = TVirtualFFT::SetTransform(notSure);
//Use the following method to get the full output:
Double_t *re_fullC = new Double_t[n];
Double_t *im_fullC = new Double_t[n];
fftC->GetPointsComplex(re_fullC,im_fullC);
//Now let's make a backward transform:
TVirtualFFT *fft_backC = TVirtualFFT::FFT(1, &n, "C2R M K");
fft_backC->SetPointsComplex(re_fullC,im_fullC);
fft_backC->Transform();
TH1 *hbC = 0;
//Let's look at the output
hbC = TH1::TransformHisto(fft_backC,hbC,"Re");
hbC->SetTitle("The backward transform result");
delete fft_backC;
fft_backC=0;
//////////////////////////////////////////////
I have the same doubts, because I need to perform a correlation using Fourier Transform. The correlation for two functions g and h:
F[Corr(g,h)] = [F(g)]* . F(h)
the asterisk denotes the complex conjugate. My doubt is about how to calculate the complex conjugate for the first transformation and after that, how multiply them.
A way to perform operations on Fourier transforms in ROOT is to use their individual real and imaginary parts, which can be stored in
If you take the transform of two 1D histograms and store their respective values in arrays, one can then perform any mathematical operation using complex numbers such as multiplication and division.
Here is an example of how to multiply two Fourier transforms and then take the inverse using a modified FFTC.C available in the ROOT tutorials:
//h2 and h3 are filled with whatever data you like
Int_t n=2500;
TH1 *hm2 =0;
TVirtualFFT::SetTransform(0);
hm2 = h2->FFT(hm2, "MAG");
hm2->SetTitle("Magnitude of the 1st transform");
//Look at the phase of the output
TH1 *hp2 = 0;
hp2 = h2->FFT(hp2, "PH");
hp2->SetTitle("Phase of the 1st transform");
//Look at the DC component and the Nyquist harmonic:
Double_t re2, im2;
//That's the way to get the current transform object:
TVirtualFFT *fft2 = TVirtualFFT::GetCurrentTransform();
//Use the following method to get the full output:
Double_t *re_full2 = new Double_t[n];
Double_t *im_full2 = new Double_t[n];
fft2->GetPointsComplex(re_full2,im_full2);
//Now let's make a backward transform:
TVirtualFFT *fft_back2 = TVirtualFFT::FFT(1, &n, "C2R M K");
fft_back2->SetPointsComplex(re_full2,im_full2);
fft_back2->Transform();
TH1 *hb2 = 0;
//Let's look at the output
hb2 = TH1::TransformHisto(fft_back2,hb2,"Re");
hb2->SetTitle("The backward transform result");
delete fft_back2;
fft_back2=0;
//////////////////////////////////////////////
TH1 *hm3 =0;
TVirtualFFT::SetTransform(0);
hm3 = h3->FFT(hm3, "MAG");
hm3->SetTitle("Magnitude of the 2nd transform");
//Look at the phase of the output
TH1 *hp3 = 0;
hp3 = h3->FFT(hp3, "PH");
hp3->SetTitle("Phase of the 2nd transform");
//Look at the DC component and the Nyquist harmonic:
Double_t re3, im3;
//That's the way to get the current transform object:
TVirtualFFT *fft3 = TVirtualFFT::GetCurrentTransform();
//Use the following method to get the full output:
Double_t *re_full3 = new Double_t[n];
Double_t *im_full3 = new Double_t[n];
fft3->GetPointsComplex(re_full3,im_full3);
//Now let's make a backward transform:
TVirtualFFT *fft_back3 = TVirtualFFT::FFT(1, &n, "C2R M K");
fft_back3->SetPointsComplex(re_full3,im_full3);
fft_back3->Transform();
TH1 *hb3 = 0;
//Let's look at the output
hb3 = TH1::TransformHisto(fft_back3,hb3,"Re");
hb3->SetTitle("The backward transform result");
delete fft_back3;
fft_back3=0;
//////////////////////////////////////////////////////
//need to multiply the Fourier transform of h2 and h3
double *re_conv = new double[n];
double *im_conv = new double[n];
for(int i = 0; i<n; i++)
{
re_conv[i] = re_full2[i] * re_full3[i] - im_full2[i] * im_full3[i];
im_conv[i] = re_full2[i] * im_full3[i] + re_full3[i] * im_full2[i];
}
TVirtualFFT::SetTransform(0);
//Now let's make a backward transform:
TVirtualFFT *fft_back4 = TVirtualFFT::FFT(1, &n, "C2R M K");
fft_back4->SetPointsComplex(re_conv,im_conv);
fft_back4->Transform();
TH1 *hb4 = 0;
//Let's look at the output
hb4 = TH1::TransformHisto(fft_back4,hb4,"Re");
hb4->SetTitle("Inverse of the product of two Fourier transforms");
delete fft_back4;
fft_back4=0;