Logo ROOT  
Reference Guide
FFT.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_fft
3/// \notebook -js
4/// This tutorial illustrates the Fast Fourier Transforms interface in ROOT.
5/// FFT transform types provided in ROOT:
6///
7/// - "C2CFORWARD" - a complex input/output discrete Fourier transform (DFT)
8/// in one or more dimensions, -1 in the exponent
9/// - "C2CBACKWARD"- a complex input/output discrete Fourier transform (DFT)
10/// in one or more dimensions, +1 in the exponent
11/// - "R2C" - a real-input/complex-output discrete Fourier transform (DFT)
12/// in one or more dimensions,
13/// - "C2R" - inverse transforms to "R2C", taking complex input
14/// (storing the non-redundant half of a logically Hermitian array)
15/// to real output
16/// - "R2HC" - a real-input DFT with output in "halfcomplex" format,
17/// i.e. real and imaginary parts for a transform of size n stored as
18/// r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1
19/// - "HC2R" - computes the reverse of FFTW_R2HC, above
20/// - "DHT" - computes a discrete Hartley transform
21///
22///
23///
24/// Sine/cosine transforms:
25///
26/// - DCT-I (REDFT00 in FFTW3 notation)
27/// - DCT-II (REDFT10 in FFTW3 notation)
28/// - DCT-III(REDFT01 in FFTW3 notation)
29/// - DCT-IV (REDFT11 in FFTW3 notation)
30/// - DST-I (RODFT00 in FFTW3 notation)
31/// - DST-II (RODFT10 in FFTW3 notation)
32/// - DST-III(RODFT01 in FFTW3 notation)
33/// - DST-IV (RODFT11 in FFTW3 notation)
34///
35/// First part of the tutorial shows how to transform the histograms
36/// Second part shows how to transform the data arrays directly
37///
38/// \macro_image
39/// \macro_output
40/// \macro_code
41///
42/// \authors Anna Kreshuk, Jens Hoffmann
43
44#include "TH1D.h"
45#include "TVirtualFFT.h"
46#include "TF1.h"
47#include "TCanvas.h"
48#include "TMath.h"
49
50void FFT()
51{
52 // Histograms
53 // =========
54 //prepare the canvas for drawing
55 TCanvas *myc = new TCanvas("myc", "Fast Fourier Transform", 800, 600);
56 myc->SetFillColor(45);
57 TPad *c1_1 = new TPad("c1_1", "c1_1",0.01,0.67,0.49,0.99);
58 TPad *c1_2 = new TPad("c1_2", "c1_2",0.51,0.67,0.99,0.99);
59 TPad *c1_3 = new TPad("c1_3", "c1_3",0.01,0.34,0.49,0.65);
60 TPad *c1_4 = new TPad("c1_4", "c1_4",0.51,0.34,0.99,0.65);
61 TPad *c1_5 = new TPad("c1_5", "c1_5",0.01,0.01,0.49,0.32);
62 TPad *c1_6 = new TPad("c1_6", "c1_6",0.51,0.01,0.99,0.32);
63 c1_1->Draw();
64 c1_2->Draw();
65 c1_3->Draw();
66 c1_4->Draw();
67 c1_5->Draw();
68 c1_6->Draw();
69 c1_1->SetFillColor(30);
70 c1_1->SetFrameFillColor(42);
71 c1_2->SetFillColor(30);
72 c1_2->SetFrameFillColor(42);
73 c1_3->SetFillColor(30);
74 c1_3->SetFrameFillColor(42);
75 c1_4->SetFillColor(30);
76 c1_4->SetFrameFillColor(42);
77 c1_5->SetFillColor(30);
78 c1_5->SetFrameFillColor(42);
79 c1_6->SetFillColor(30);
80 c1_6->SetFrameFillColor(42);
81
82 c1_1->cd();
84
85 //A function to sample
86 TF1 *fsin = new TF1("fsin", "sin(x)+sin(2*x)+sin(0.5*x)+1", 0, 4*TMath::Pi());
87 fsin->Draw();
88
89 Int_t n=25;
90 TH1D *hsin = new TH1D("hsin", "hsin", n+1, 0, 4*TMath::Pi());
91 Double_t x;
92
93 //Fill the histogram with function values
94 for (Int_t i=0; i<=n; i++){
95 x = (Double_t(i)/n)*(4*TMath::Pi());
96 hsin->SetBinContent(i+1, fsin->Eval(x));
97 }
98 hsin->Draw("same");
99 fsin->GetXaxis()->SetLabelSize(0.05);
100 fsin->GetYaxis()->SetLabelSize(0.05);
101
102 c1_2->cd();
103 //Compute the transform and look at the magnitude of the output
104 TH1 *hm =0;
106 hm = hsin->FFT(hm, "MAG");
107 hm->SetTitle("Magnitude of the 1st transform");
108 hm->Draw();
109 //NOTE: for "real" frequencies you have to divide the x-axes range with the range of your function
110 //(in this case 4*Pi); y-axes has to be rescaled by a factor of 1/SQRT(n) to be right: this is not done automatically!
111
112 hm->SetStats(kFALSE);
113 hm->GetXaxis()->SetLabelSize(0.05);
114 hm->GetYaxis()->SetLabelSize(0.05);
115 c1_3->cd();
116 //Look at the phase of the output
117 TH1 *hp = 0;
118 hp = hsin->FFT(hp, "PH");
119 hp->SetTitle("Phase of the 1st transform");
120 hp->Draw();
121 hp->SetStats(kFALSE);
122 hp->GetXaxis()->SetLabelSize(0.05);
123 hp->GetYaxis()->SetLabelSize(0.05);
124
125 //Look at the DC component and the Nyquist harmonic:
126 Double_t re, im;
127 //That's the way to get the current transform object:
129 c1_4->cd();
130 //Use the following method to get just one point of the output
131 fft->GetPointComplex(0, re, im);
132 printf("1st transform: DC component: %f\n", re);
133 fft->GetPointComplex(n/2+1, re, im);
134 printf("1st transform: Nyquist harmonic: %f\n", re);
135
136 //Use the following method to get the full output:
137 Double_t *re_full = new Double_t[n];
138 Double_t *im_full = new Double_t[n];
139 fft->GetPointsComplex(re_full,im_full);
140
141 //Now let's make a backward transform:
142 TVirtualFFT *fft_back = TVirtualFFT::FFT(1, &n, "C2R M K");
143 fft_back->SetPointsComplex(re_full,im_full);
144 fft_back->Transform();
145 TH1 *hb = 0;
146 //Let's look at the output
147 hb = TH1::TransformHisto(fft_back,hb,"Re");
148 hb->SetTitle("The backward transform result");
149 hb->Draw();
150 //NOTE: here you get at the x-axes number of bins and not real values
151 //(in this case 25 bins has to be rescaled to a range between 0 and 4*Pi;
152 //also here the y-axes has to be rescaled (factor 1/bins)
153 hb->SetStats(kFALSE);
154 hb->GetXaxis()->SetLabelSize(0.05);
155 hb->GetYaxis()->SetLabelSize(0.05);
156 delete fft_back;
157 fft_back=0;
158
159// Data array - same transform
160// ===========================
161
162 //Allocate an array big enough to hold the transform output
163 //Transform output in 1d contains, for a transform of size N,
164 //N/2+1 complex numbers, i.e. 2*(N/2+1) real numbers
165 //our transform is of size n+1, because the histogram has n+1 bins
166
167 Double_t *in = new Double_t[2*((n+1)/2+1)];
168 Double_t re_2,im_2;
169 for (Int_t i=0; i<=n; i++){
170 x = (Double_t(i)/n)*(4*TMath::Pi());
171 in[i] = fsin->Eval(x);
172 }
173
174 //Make our own TVirtualFFT object (using option "K")
175 //Third parameter (option) consists of 3 parts:
176 //- transform type:
177 // real input/complex output in our case
178 //- transform flag:
179 // the amount of time spent in planning
180 // the transform (see TVirtualFFT class description)
181 //- to create a new TVirtualFFT object (option "K") or use the global (default)
182 Int_t n_size = n+1;
183 TVirtualFFT *fft_own = TVirtualFFT::FFT(1, &n_size, "R2C ES K");
184 if (!fft_own) return;
185 fft_own->SetPoints(in);
186 fft_own->Transform();
187
188 //Copy all the output points:
189 fft_own->GetPoints(in);
190 //Draw the real part of the output
191 c1_5->cd();
192 TH1 *hr = 0;
193 hr = TH1::TransformHisto(fft_own, hr, "RE");
194 hr->SetTitle("Real part of the 3rd (array) tranfsorm");
195 hr->Draw();
196 hr->SetStats(kFALSE);
197 hr->GetXaxis()->SetLabelSize(0.05);
198 hr->GetYaxis()->SetLabelSize(0.05);
199 c1_6->cd();
200 TH1 *him = 0;
201 him = TH1::TransformHisto(fft_own, him, "IM");
202 him->SetTitle("Im. part of the 3rd (array) transform");
203 him->Draw();
204 him->SetStats(kFALSE);
205 him->GetXaxis()->SetLabelSize(0.05);
206 him->GetYaxis()->SetLabelSize(0.05);
207
208 myc->cd();
209 //Now let's make another transform of the same size
210 //The same transform object can be used, as the size and the type of the transform
211 //haven't changed
212 TF1 *fcos = new TF1("fcos", "cos(x)+cos(0.5*x)+cos(2*x)+1", 0, 4*TMath::Pi());
213 for (Int_t i=0; i<=n; i++){
214 x = (Double_t(i)/n)*(4*TMath::Pi());
215 in[i] = fcos->Eval(x);
216 }
217 fft_own->SetPoints(in);
218 fft_own->Transform();
219 fft_own->GetPointComplex(0, re_2, im_2);
220 printf("2nd transform: DC component: %f\n", re_2);
221 fft_own->GetPointComplex(n/2+1, re_2, im_2);
222 printf("2nd transform: Nyquist harmonic: %f\n", re_2);
223 delete fft_own;
224 delete [] in;
225 delete [] re_full;
226 delete [] im_full;
227}
228
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
double Double_t
Definition: RtypesCore.h:59
virtual void SetLabelSize(Float_t size=0.04)
Set size of axis labels.
Definition: TAttAxis.cxx:206
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
void SetFrameFillColor(Color_t color=1)
Definition: TAttPad.h:73
The Canvas class.
Definition: TCanvas.h:23
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
Definition: TCanvas.cxx:706
1-Dim function class
Definition: TF1.h:213
TAxis * GetYaxis() const
Get y axis of the function.
Definition: TF1.cxx:2424
void Draw(Option_t *option="") override
Draw this function with its current attributes.
Definition: TF1.cxx:1342
virtual Double_t Eval(Double_t x, Double_t y=0, Double_t z=0, Double_t t=0) const
Evaluate this function.
Definition: TF1.cxx:1454
TAxis * GetXaxis() const
Get x axis of the function.
Definition: TF1.cxx:2413
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:617
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
void SetTitle(const char *title) override
See GetStatOverflows for more information.
Definition: TH1.cxx:6708
static TH1 * TransformHisto(TVirtualFFT *fft, TH1 *h_output, Option_t *option)
For a given transform (first parameter), fills the histogram (second parameter) with the transform ou...
Definition: TH1.cxx:9194
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition: TH1.cxx:1277
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:319
virtual TH1 * FFT(TH1 *h_output, Option_t *option)
This function allows to do discrete Fourier transforms of TH1 and TH2.
Definition: TH1.cxx:3288
TAxis * GetYaxis()
Definition: TH1.h:320
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition: TH1.cxx:3070
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:9097
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8865
The most important graphics class in the ROOT system.
Definition: TPad.h:26
TVirtualPad * cd(Int_t subpadnumber=0) override
Set Current pad.
Definition: TPad.cxx:606
void Draw(Option_t *option="") override
Draw Pad in Current pad (re-parent pad if necessary).
Definition: TPad.cxx:1292
TVirtualFFT is an interface class for Fast Fourier Transforms.
Definition: TVirtualFFT.h:88
static void SetTransform(TVirtualFFT *fft)
static: set the current transfrom to parameter
virtual void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const =0
static TVirtualFFT * FFT(Int_t ndim, Int_t *n, Option_t *option)
Returns a pointer to the FFT of requested size and type.
static TVirtualFFT * GetCurrentTransform()
static: return current fgFFT
virtual void SetPoints(const Double_t *data)=0
virtual void SetPointsComplex(const Double_t *re, const Double_t *im)=0
virtual void Transform()=0
virtual void GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const =0
virtual void GetPointsComplex(Double_t *re, Double_t *im, Bool_t fromInput=kFALSE) const =0
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
constexpr Double_t Pi()
Definition: TMath.h:37