I’m having an issue with the ROOT FFT. I use the TVirtualFFT ROOT class and I am using the v5r28p00n01 version of ROOT.
I have included some code below. The issue is that, for the specified function and the specified values, the imaginary part of the FFT output should be exactly zero. However, one notices that this is not the case. In the histogram of the imaginary part of the FFT result, one sees that for n = 2 and n = 4 the values are 1.11022e-16 and -1.11022e-16 respectively.
Does anyone know of a workaround for this?
[code]void ErrorPropCalc(int N = -1){
gStyle->SetOptStat(0);
gStyle->SetTitleFillColor(kWhite);
//Files
TFile *OutFile = new TFile(“ErrorPropCalcPlots.root”, “RECREATE”);
//Functions
TF1 *fsin = new TF1(“fsin”, “sin(x)”, 0, TMath::Pi());
//Histograms
TH1D *hsin_values = new TH1D(“hsin_values”, “hsin”, N, 0, TMath::Pi());
TH1 *hsinTransform = NULL;
TH1 *hsinReal = NULL;
TH1 *hsinIm = NULL;
//Canvases
TCanvas *c1 = new TCanvas(“c1”, “c1”, 1280, 800);
TCanvas *c2 = new TCanvas(“c2”, “c2”, 1280, 800);
TCanvas *c4 = new TCanvas(“c4”, “c4”, 1280, 800);
TCanvas *c5 = new TCanvas(“c5”, “c5”, 1280, 800);
//Variables
double x;
double f_R;
double f_I;
//Change to c1
c1 -> cd();
//Draw sin(x) function on canvas
fsin -> Draw();
//Fill hsin_values with pik/4 values
for (int k = 0; k < N; k++){
x = ((double(k))/(double(N)))(TMath::Pi());
hsin_values -> SetBinContent(k+1, fsin->Eval(x));
std::cout<<"x = "<<x<<std::endl;
std::cout<<"fsin->Eval(x) = "<Eval(x)<<std::endl;
}
//Draw hsin_values on same canvas as fsin
hsin_values -> Draw(“same”);
//Tidy up canvas
c1->SetFillColor(kWhite);
fsin -> SetTitle(“Sin(x) & Sin(x) with x = pi*k/4”);
fsin -> GetYaxis()-> SetTitle(“Sin(x)”);
fsin -> GetXaxis()-> SetTitle(“x”);
c1 -> Print(“Sin(x).png”);
//change to second canvas
c2 -> cd();
//Initialze FFT, perform FFT, output Magnitude
TVirtualFFT::SetTransform(0);
hsinTransform = hsin_values -> FFT(hsinTransform, “MAG”);
//Tidy up Mag histogram
c2 -> SetFillColor(kWhite);
hsinTransform -> SetTitle("|FFT[sin(x)]|");
hsinTransform->GetYaxis()->SetTitle(“FFTMag(n)”);
hsinTransform->GetXaxis()->SetTitle(“n”);
hsinTransform->SetLineWidth(2);
hsinTransform->Draw(“e”);
//change to fourth canvas
c4 -> cd();
//Initialze FFT, perform FFT, output Real part
TVirtualFFT::SetTransform(0);
hsinReal = hsin_values -> FFT(hsinReal, “RE”);
//Tidy up Real histogram
c4 -> SetFillColor(kWhite);
hsinReal -> SetTitle(“Re_FFT[sin(x)]”);
hsinReal->GetYaxis()->SetTitle(“Re_FFT(n)”);
hsinReal->GetXaxis()->SetTitle(“n”);
hsinReal->SetLineWidth(2);
hsinReal->Draw(“e”);
c4 -> Print(“RealSinX.png”);
//change to fifth canvas
c5 -> cd();
//Initialze FFT, perform FFT, output Imaginary part
TVirtualFFT::SetTransform(0);
hsinIm = hsin_values -> FFT(hsinIm, “IM P”);
//Tidy up Imaginary histogram
c5 -> SetFillColor(kWhite);
hsinIm -> SetTitle(“Im_FFT[sin(x)]”);
hsinIm->GetYaxis()->SetTitle(“Im_FFT(n)”);
hsinIm->GetXaxis()->SetTitle(“n”);
hsinIm->SetLineWidth(2);
hsinIm->Draw();
c5 -> Print(“ImSinX.png”);[/code]