TVirtualFFT.SetPointsComplex problems

Hi,

I recently updated to Root v5.14 to use the Fourier Transform package. I realize all of it is very new, but I'm having trouble doing the first things with it, despite reading the "tutorial" macro and the class documentation for TH1->FFT and TVirtualFFT.

I have modified the tutorial macro to see if I can modify a FFT and get a new function. Ultimately, I hoped to "trim" the FFT and keep only the first five or eight contribtuions or so. As you can see, I am having great difficulty. My macro is attached. I start out making a function of several sines, and histogramming it. I find the corresponding FFT and historgram it. Then I hit trouble.

1. I try to copy the FFT and then reassign my trimmed arrays with SetPointsComplex. The FFT ignores my attempt and keeps its own arrays instead, as shown in histogram "hout1".
2. I try to make a brand new FFT, using what I thought was the default that would be used in the original ("R2C"), and assign my trimmed arrays with SetPointsComplex. The FFT goes completely bananas and acts like I assigned it values >4000E246, as shown in histogram "hout2".
3. I try another new FFT using "C2CFORWARD" and get Inf/NaN complaints. There is thus no output histo "hout3".

I am not at all sure what I'm doing wrong, and think what I have done should have given me plots that looked like "out_MAG" except with zeros after the 8th element. Once I get that far, I'm not sure how to use the TVirtualFFT to create a histogram of the result, or even if I can. (Maybe I'm supposed to just built my own TF1 at that point?)

Any guidance from the code authors, or any others, would be greatly appreciated.

• John

-- 

Dr. John Krane
http://jkrane.home.comcast.net

Quantitative Financial Analyst





#include "TH1D.h"
#include "TVirtualFFT.h"
#include "TF1.h"
#include "TCanvas.h"

void test2()
{




  //prepare the canvas for drawing

  TCanvas *orig = new TCanvas("orig", "Original distribution", 600, 600);
  TCanvas *myc = new TCanvas("myc", "Fast Fourier Transform", 600, 600);
  myc->SetFillColor(45);
  TPad *c1_1 = new TPad("c1_1", "c1_1",0.01,0.51,0.49,0.99);
  TPad *c1_2 = new TPad("c1_2", "c1_2",0.51,0.51,0.99,0.99);
  TPad *c1_3 = new TPad("c1_3", "c1_3",0.01,0.01,0.49,0.49);
  TPad *c1_4 = new TPad("c1_4", "c1_4",0.51,0.01,0.99,0.49);
  c1_1->Draw();
  c1_2->Draw();
  c1_3->Draw();
  c1_4->Draw();
  c1_1->SetFillColor(30);
  c1_1->SetFrameFillColor(42);
  c1_2->SetFillColor(30);
  c1_2->SetFrameFillColor(42);
  c1_3->SetFillColor(30);
  c1_3->SetFrameFillColor(42);
  c1_4->SetFillColor(30);
  c1_4->SetFrameFillColor(42);
  
  orig->cd();


  //A function to sample

  TF1 *fsin = new TF1("fsin", "sin(x)+sin(2*x)+sin(0.5*x)+1", 0, 4*TMath::Pi());
  fsin->Draw();
  Int_t n=25;
  TH1D *hsin = new TH1D("hsin", "hsin", n+1, 0, 4*TMath::Pi());
  Double_t x;


  //Fill the histogram with function values

  for (Int_t i=0; i<=n; i++){
    x = (Double_t(i)/n)*(4*TMath::Pi());
    hsin->SetBinContent(i+1, fsin->Eval(x));
  }
  hsin->Draw("same");
  fsin->GetXaxis()->SetLabelSize(0.05);
  fsin->GetYaxis()->SetLabelSize(0.05);




  //Compute the transform and look at the magnitude of the output

  TH1 *hm =0;
  hm = hsin->FFT(hm, "MAG");
  c1_1->cd();
  hm->Draw();
  hm->SetStats(kFALSE);
  hm->GetXaxis()->SetLabelSize(0.05);
  hm->GetYaxis()->SetLabelSize(0.05);




  // J's mods.  Try (in vain) to modify a copy of the FFT

  Double_t *arr=new Double_t[26]; Double_t *arrcheck=new Double_t[26];
  Double_t *iarr=new Double_t[26]; Double_t *iarrcheck=new Double_t[26];
  TH1 *hm2 =0;
  hm2 = hsin->FFT(hm2, "MAG");
  TVirtualFFT *theFFT= TVirtualFFT::GetCurrentTransform();



  //

  // These will be my own versions of FFT

  //




  // This should be an exact copy

  TVirtualFFT *fft_own1= theFFT;    



  // These are made from scratch, but is R2C what I want?  Don't know for sure

  // but I think this is the defaut used by theFFT

  TVirtualFFT *fft_own2= TVirtualFFT::FFT(1, theFFT->GetN(), "R2C P K");
  TVirtualFFT *fft_own3= TVirtualFFT::FFT(1, theFFT->GetN(), "C2CFORWARD P K");



  // Trim the higher orders from expansion in theFFT for my versions

  Double_t val,ival;
  for (int i=0; i<*theFFT->GetN(); ++i) {
    theFFT->GetPointComplex(&i,val,ival);
    arrcheck[i]=val;    iarrcheck[i]=ival;
    if (i<8){
      arr[i]=val; iarr[i]=ival;
    } else {
      arr[i]=0.0; iarr[i]=0.0;
    }

    // screen output to see what we have
    cerr<<i<<"   :    "<<arrcheck[i]<<" -> "<<arr[i]
	   <<"   :    "<<iarrcheck[i]<<" -> "<<iarr[i]
	<<endl;
  }
  fft_own1->SetPointsComplex(arr,iarr);
  fft_own2->SetPointsComplex(arr,iarr);
  fft_own3->SetPointsComplex(arr,iarr);

  


  // Fill histo with truncated Fourier series

  TH1D *hout1= new TH1D("hout1", "hout1", n+1, 0, 4*TMath::Pi());
  TH1D *hout2= new TH1D("hout2", "hout2", n+1, 0, 4*TMath::Pi());
  TH1D *hout3= new TH1D("hout3", "hout3", n+1, 0, 4*TMath::Pi());
  hout1->TransformHisto(fft_own1,hout1,"MAG");
  hout2->TransformHisto(fft_own2,hout2,"MAG");
  hout3->TransformHisto(fft_own3,hout3,"MAG");

  c1_2->cd();  hout1->Draw();
  c1_3->cd();  hout2->Draw();
  c1_4->cd();  hout3->Draw();


}