source: trunk/MagicSoft/Mars/mtools/MFFT.cc@ 2815

Last change on this file since 2815 was 2808, checked in by gaug, 21 years ago
*** empty log message ***
File size: 17.1 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Markus Gaug 01/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2001-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26// //
27// Fast Fourier Transforms //
28// //
29// (Most of the code is adapted from Numerical Recipies in C++, 2nd ed., //
30// pp. 509-563) //
31// //
32// Usage: //
33// //
34// 1) Functions RealFunctionFFT: (FOURIER TRANSFORM) //
35// * Take as argument arrays of real numbers, //
36// in some cases the dimension of the array has to be given separately//
37// * Return a COMPLEX array with the following meaning: //
38// array[0]: The value of F(0) (has only real component) //
39// array[1]: The value of F(N/2) (has only real component) //
40// array[2i]: The real part of F(i) //
41// array[2i+1]: The imaginary part of F(i) //
42// * Note that F(N-i)* = F(i), therefore only the positive frequency //
43// half is stored. //
44// * The dimension MUST be an integer power of 2, //
45// otherwise, the array will be shortened!! //
46// //
47// 2) Functions RealFunctionIFFT: (INVERSER FOURIER TRANSFORM) //
48// * Take as argument a COMPLEX array //
49// of Fourier-transformed REAL numbers //
50// with the following meaning: //
51// array[0]: The value of F(0) (has only real component) //
52// array[1]: The value of F(N/2) (has only real component) //
53// array[2i]: The real part of F(i) //
54// array[2i+1]: The imaginary part of F(i) //
55// * Returns the original complex array of dimension 2N-1 //
56// //
57// 3) Functions PowerSpectrumDensity: //
58// * Return a histogram with the spectral density, i.e. //
59// P(k) = 1/(N*N) * |F(k)|*|F(k)| //
60// * The histogram is ranged between 0 and 1./(2*binwidth) //
61// * The number of bins equals N/2+1 //
62// * Note that histograms with unequal binwidth can not yet be treated! //
63// * If the PSD does NOT CONVERGE to 0 at the maximum bin, //
64// you HAVE TO sample your data finer! //
65// //
66//////////////////////////////////////////////////////////////////////////////
67
68#include "MFFT.h"
69
70#include "TMath.h"
71
72#include "MLog.h"
73#include "MLogManip.h"
74
75ClassImp(MFFT);
76
77using namespace std;
78
79// ---------------------------------------------------------------------------
80//
81// Default Constructor
82// Initializes random number generator and default variables
83//
84MFFT::MFFT() : fDim(0)
85{
86}
87
88// --------------------------------------------------------------------------
89//
90// Destructor.
91//
92MFFT::~MFFT()
93{
94}
95
96void MFFT::TransformF(const Int_t isign)
97{
98
99 UInt_t n,mmax,m,j,istep,i;
100 Float_t wtemp,wr,wpr,wpi,wi,theta;
101 Float_t tempr,tempi;
102
103 Int_t nn = fDim/2;
104 n = nn << 1;
105
106 //
107 // The bit-reversal section of the routine:
108 // Exchange the two complex numbers
109 //
110 j=1;
111 for (i=1;i<n;i+=2) {
112 if (j > i) {
113 Swap(fDataF[j-1],fDataF[i-1]);
114 Swap(fDataF[j],fDataF[i]);
115 }
116 m=nn;
117 while (m >= 2 && j > m) {
118 j -= m;
119 m >>= 1;
120 }
121 j += m;
122 }
123 //
124 // Here begins the Danielson-Lanczos section of the routine
125 //
126 mmax=2;
127 while (n > mmax) { // Outer loop executed log_2(nn) times
128
129 istep = mmax << 1;
130 //
131 // Initialize the trigonometric recurrence:
132 //
133 theta = isign*(6.28318530717959/mmax);
134
135 wtemp = TMath::Sin(0.5*theta);
136 wpr = -2.0*wtemp*wtemp;
137 wpi = TMath::Sin(theta);
138
139 wr=1.0;
140 wi=0.0;
141
142 for (m=1; m<mmax; m+=2) {
143 for (i=m; i<=n; i+=istep) {
144 //
145 // The Danielson-Lanczos formula:
146 //
147 j = i+mmax;
148 tempr = wr*fDataF[j-1] - wi*fDataF[j];
149 tempi = wr*fDataF[j] + wi*fDataF[j-1];
150 fDataF[j-1] = fDataF[i-1] - tempr;
151 fDataF[j] = fDataF[i] - tempi;
152 fDataF[i-1] += tempr;
153 fDataF[i] += tempi;
154 }
155
156 //
157 // Trigonometric recurrence
158 //
159 wr = (wtemp=wr)*wpr - wi*wpi+wr;
160 wi = wi*wpr + wtemp*wpi+wi;
161
162 }
163 mmax=istep;
164 }
165}
166
167
168void MFFT::TransformD(const Int_t isign)
169{
170
171 UInt_t n,mmax,m,j,istep,i;
172 Double_t wtemp,wr,wpr,wpi,wi,theta;
173 Double_t tempr,tempi;
174
175 Int_t nn = fDim/2;
176 n = nn << 1;
177
178 //
179 // The bit-reversal section of the routine:
180 // Exchange the two complex numbers
181 //
182 j=1;
183 for (i=1;i<n;i+=2) {
184 if (j > i) {
185 Swap(fDataD[j-1],fDataD[i-1]);
186 Swap(fDataD[j],fDataD[i]);
187 }
188 m=nn;
189 while (m >= 2 && j > m) {
190 j -= m;
191 m >>= 1;
192 }
193 j += m;
194 }
195 //
196 // Here begins the Danielson-Lanczos section of the routine
197 //
198 mmax=2;
199 while (n > mmax) { // Outer loop executed log_2(nn) times
200
201 istep = mmax << 1;
202 //
203 // Initialize the trigonometric recurrence:
204 //
205 theta = isign*(6.28318530717959/mmax);
206
207 wtemp = TMath::Sin(0.5*theta);
208 wpr = -2.0*wtemp*wtemp;
209 wpi = TMath::Sin(theta);
210
211 wr=1.0;
212 wi=0.0;
213
214 for (m=1; m<mmax; m+=2) {
215 for (i=m; i<=n; i+=istep) {
216 //
217 // The Danielson-Lanczos formula:
218 //
219 j = i+mmax;
220 tempr = wr*fDataD[j-1] - wi*fDataD[j];
221 tempi = wr*fDataD[j] + wi*fDataD[j-1];
222 fDataD[j-1] = fDataD[i-1] - tempr;
223 fDataD[j] = fDataD[i] - tempi;
224 fDataD[i-1] += tempr;
225 fDataD[i] += tempi;
226 }
227
228 //
229 // Trigonometric recurrence
230 //
231 wr = (wtemp=wr)*wpr - wi*wpi+wr;
232 wi = wi*wpr + wtemp*wpi+wi;
233
234 }
235 mmax=istep;
236 }
237}
238
239//
240// Calculates the Fourier transform of a set of n real-valued data points.
241// Replaces this data (which is stored in array data[1..n]) by the positive
242// frequency half of its complex Fourier transform. The real-valued first
243// and last components of the complex transform are returned as elements
244// data[1] and data[2], respectively. n must be a power of 2. This routine
245// also calculates the inverse transform of a complex data array if it is
246// the transform of real data. (Result in this case mus be multiplied by
247// 2/n.). From NUMERICAL RECIPES IN C++.
248//
249void MFFT::RealFTF(const Int_t isign)
250{
251
252 Int_t i,i1,i2,i3,i4;
253 Float_t c1=0.5,c2,h1r,h1i,h2r,h2i;
254 Float_t wr,wi,wpr,wpi,wtemp,theta;
255
256 //
257 // Initialize the recurrence
258 //
259 theta = TMath::Pi() / (Double_t)(fDim>>1);
260
261 if(isign==1) // forward transform
262 {
263 c2 = -0.5;
264 TransformF(1);
265 }
266 else // set up backward transform
267 {
268 c2 = 0.5;
269 theta = -theta;
270 }
271
272 wtemp = TMath::Sin(0.5*theta);
273 wpr = -2.0*wtemp*wtemp;
274 wpi = TMath::Sin(theta);
275
276 wr = 1.0 + wpr;
277 wi = wpi;
278
279 for(i=1;i<(fDim>>2);i++) // case i=0 done separately below
280 {
281
282 i2 = 1 + (i1 = i+i);
283 i4 = 1 + (i3 = fDim-i1);
284
285 //
286 // The two separate transforms are separated out of the data
287 //
288 h1r = c1*(fDataF[i1]+fDataF[i3]);
289 h1i = c1*(fDataF[i2]-fDataF[i4]);
290 h2r = -c2*(fDataF[i2]+fDataF[i4]);
291 h2i = c2*(fDataF[i1]-fDataF[i3]);
292
293 //
294 // Here, they are recombined to from the true transform
295 // of the orginal real data
296 //
297 fDataF[i1] = h1r + wr*h2r - wi*h2i;
298 fDataF[i2] = h1i + wr*h2i + wi*h2r;
299 fDataF[i3] = h1r - wr*h2r + wi*h2i;
300 fDataF[i4] = -h1i + wr*h2i + wi*h2r;
301
302 //
303 // The recurrence
304 //
305 wr = (wtemp=wr)*wpr - wi*wpi + wr;
306 wi = wi*wpr + wtemp*wpi + wi;
307 }
308
309 //
310 // Squeeze the first and last data together to get them all
311 // within the original array
312 //
313 if(isign==1)
314 {
315 fDataF[0] = (h1r=fDataF[0]) + fDataF[1];
316 fDataF[1] = h1r - fDataF[1];
317 }
318 else
319 {
320
321 fDataF[0] = c1*((h1r=fDataF[0]) + fDataF[1]);
322 fDataF[1] = c1*(h1r-fDataF[1]);
323
324 //
325 // The inverse transform for the case isign = -1
326 //
327 TransformF(-1);
328
329 //
330 // normalize correctly (not done in original NR's)
331 //
332 for(i=1;i<=fDim;i++)
333 fDataF[i] *= (2./fDim);
334 }
335}
336void MFFT::RealFTD(const Int_t isign)
337{
338
339 Int_t i,i1,i2,i3,i4;
340 Float_t c1=0.5,c2,h1r,h1i,h2r,h2i;
341 Double_t wr,wi,wpr,wpi,wtemp,theta;
342
343 //
344 // Initialize the recurrence
345 //
346 theta=3.141592653589793/(Double_t) (fDim>>1);
347
348 if(isign==1) // forward transform
349 {
350 c2 = -0.5;
351 TransformD(1);
352 }
353 else // set up backward transform
354 {
355 c2 = 0.5;
356 theta = -theta;
357 }
358
359 wtemp = TMath::Sin(0.5*theta);
360 wpr = -2.0*wtemp*wtemp;
361 wpi = TMath::Sin(theta);
362
363 wr = 1.0 + wpr;
364 wi = wpi;
365
366 for(i=1;i<(fDim>>2);i++) // case i=0 done separately below
367 {
368
369 i2 = 1 + (i1 = i+i);
370 i4 = 1 + (i3 = fDim-i1);
371
372 //
373 // The two separate transforms are separated out of the data
374 //
375 h1r = c1*(fDataD[i1]+fDataD[i3]);
376 h1i = c1*(fDataD[i2]-fDataD[i4]);
377 h2r = -c2*(fDataD[i2]+fDataD[i4]);
378 h2i = c2*(fDataD[i1]-fDataD[i3]);
379
380 //
381 // Here, they are recombined to from the true transform
382 // of the orginal real data
383 //
384 fDataD[i1] = h1r + wr*h2r - wi*h2i;
385 fDataD[i2] = h1i + wr*h2i + wi*h2r;
386 fDataD[i3] = h1r - wr*h2r + wi*h2i;
387 fDataD[i4] = -h1i + wr*h2i + wi*h2r;
388
389 //
390 // The recurrence
391 //
392 wr = (wtemp=wr)*wpr - wi*wpi + wr;
393 wi = wi*wpr + wtemp*wpi + wi;
394 }
395
396 //
397 // Squeeze the first and last data together to get them all
398 // within the original array
399 //
400 if(isign==1)
401 {
402 fDataD[0] = (h1r=fDataD[0]) + fDataD[1];
403 fDataD[1] = h1r - fDataD[1];
404 }
405 else
406 {
407
408 fDataD[0] = c1*((h1r=fDataD[0]) + fDataD[1]);
409 fDataD[1] = c1*(h1r-fDataD[1]);
410
411 //
412 // The inverse transform for the case isign = -1
413 //
414 TransformD(-1);
415
416 //
417 // normalize correctly (not done in original NR's)
418 //
419 for(i=1;i<=fDim;i++)
420 fDataD[i] *= (2./fDim);
421 }
422}
423
424
425//
426// Fast Fourier Transform for float arrays
427//
428Float_t* MFFT::RealFunctionFFT(const Int_t n, const Float_t *data)
429{
430
431 fDim = n;
432 CheckDim(n);
433
434 fDataF.Set(fDim);
435 //
436 // Clone the array
437 //
438 for (Int_t i=0;i<fDim;i++)
439 fDataF[i] = data[i];
440
441 RealFTF(1);
442
443 return fDataF.GetArray();
444
445}
446
447//
448// Fast Inverse Fourier Transform for float arrays
449//
450Float_t* MFFT::RealFunctionIFFT(const Int_t n, const Float_t *data)
451{
452
453 fDim = n;
454 CheckDim(fDim);
455
456 fDataF.Set(fDim);
457 //
458 // Clone the array
459 //
460 for (Int_t i=0;i<fDim;i++)
461 fDataF[i] = data[i];
462
463 RealFTF(-1);
464
465 return fDataF.GetArray();
466
467}
468
469//
470// Fast Fourier Transform for double arrays
471//
472Double_t* MFFT::RealFunctionFFT(const Int_t n, const Double_t *data)
473{
474
475 fDim = n;
476 CheckDim(n);
477
478 fDataD.Set(fDim);
479 //
480 // Clone the array
481 //
482 for (Int_t i=0;i<fDim;i++)
483 fDataD[i] = data[i];
484
485 RealFTD(1);
486
487 return fDataD.GetArray();
488
489}
490
491//
492// Fast Inverse Fourier Transform for double arrays
493//
494Double_t* MFFT::RealFunctionIFFT(const Int_t n, const Double_t *data)
495{
496
497 fDim = n;
498 CheckDim(fDim);
499
500 fDataD.Set(fDim);
501 //
502 // Clone the array
503 //
504 for (Int_t i=0;i<fDim;i++)
505 fDataD[i] = data[i];
506
507 RealFTD(-1);
508
509 return fDataD.GetArray();
510
511}
512
513//
514// Fast Fourier Transform for TArrayF's
515//
516TArrayF* MFFT::RealFunctionFFT(const TArrayF *data)
517{
518
519 fDim = data->GetSize();
520 CheckDim(fDim);
521
522 fDataF.Set(fDim);
523 //
524 // Clone the array
525 //
526 for (Int_t i=0;i<fDim;i++)
527 fDataF[i] = data->At(i);
528
529 RealFTF(1);
530
531 return new TArrayF(fDim,fDataF.GetArray());
532
533}
534
535//
536// Inverse Fast Fourier Transform for TArrayF's
537//
538TArrayF* MFFT::RealFunctionIFFT(const TArrayF *data)
539{
540
541 fDim = data->GetSize();
542 CheckDim(fDim);
543
544 fDataF.Set(fDim);
545 //
546 // Clone the array
547 //
548 for (Int_t i=0;i<fDim;i++)
549 fDataF[i] = data->At(i);
550
551 RealFTF(-1);
552
553 return new TArrayF(fDim,fDataF.GetArray());
554}
555
556
557//
558// Fast Fourier Transform for TArrayD's
559//
560TArrayD* MFFT::RealFunctionFFT(const TArrayD *data)
561{
562
563 fDim = data->GetSize();
564 CheckDim(fDim);
565
566 fDataD.Set(fDim);
567 //
568 // Clone the array
569 //
570 for (Int_t i=0;i<fDim;i++)
571 fDataD[i] = data->At(i);
572
573 RealFTD(1);
574
575 return new TArrayD(fDim,fDataD.GetArray());
576
577}
578
579//
580// Inverse Fast Fourier Transform for TArrayD's
581//
582TArrayD* MFFT::RealFunctionIFFT(const TArrayD *data)
583{
584
585 fDim = data->GetSize();
586 CheckDim(fDim);
587
588 fDataD.Set(fDim);
589 //
590 // Clone the array
591 //
592 for (Int_t i=0;i<fDim;i++)
593 fDataD[i] = data->At(i);
594
595 RealFTD(-1);
596
597 return new TArrayD(fDim,fDataD.GetArray());
598}
599
600
601//
602// Power Spectrum Density Calculation
603//
604TH1D* MFFT::PowerSpectrumDensity(const TH1D *hist)
605{
606
607 TH1D *newhist = (TH1D*)CheckHist(hist,1);
608
609 fDataD.Set(fDim);
610 //
611 // Copy the hist into an array
612 //
613 for (Int_t i=0;i<fDim;i++)
614 fDataD[i] = hist->GetBinContent(i);
615
616 RealFTD(1);
617
618 Int_t dim2 = fDim*fDim;
619 Double_t c02;
620 Double_t ck2;
621 Double_t cn2;
622 //
623 // Fill the new histogram:
624 //
625 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
626 // (stored in fData{0])
627 //
628 c02 = fDataD[0]*fDataD[0];
629 newhist->Fill(c02/dim2);
630 //
631 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)| + |C(N-k)|*|C(N-k)|)
632 //
633 for (Int_t k=2;k<fDim-2;k+=2)
634 {
635
636 Int_t ki = k+1;
637 ck2 = (fDataD[k]*fDataD[k] + fDataD[ki]*fDataD[ki]);
638 newhist->Fill(ck2/dim2);
639 }
640 //
641 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
642 // (stored in fData[1])
643 //
644 cn2 = (fDataD[1]*fDataD[1]);
645 newhist->Fill(cn2/dim2);
646
647 return newhist;
648}
649
650
651//
652// Power Spectrum Density calculation
653//
654TH1F* MFFT::PowerSpectrumDensity(const TH1F *hist)
655{
656
657 TH1F *newhist = (TH1F*)CheckHist(hist,0);
658
659 fDataF.Set(fDim);
660 //
661 // Copy the hist into an array
662 //
663 for (Int_t i=0;i<fDim;i++)
664 fDataF[i] = hist->GetBinContent(i);
665
666 RealFTF(1);
667
668 Int_t dim2 = fDim*fDim;
669 Float_t c02;
670 Float_t ck2;
671 Float_t cn2;
672 //
673 // Fill the new histogram:
674 //
675 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
676 //
677 c02 = (fDataF[0]*fDataF[0]);
678 newhist->Fill(c02/dim2);
679 //
680 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
681 //
682 for (Int_t k=2;k<fDim;k+=2)
683 {
684 Int_t ki = k+1;
685 ck2 = (fDataF[k]*fDataF[k] + fDataF[ki]*fDataF[ki]);
686 newhist->Fill(ck2/dim2);
687 }
688 //
689 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
690 //
691 cn2 = (fDataF[1]*fDataF[1]);
692 newhist->Fill(cn2/dim2);
693
694 return newhist;
695}
696
697
698void MFFT::CheckDim(Int_t a)
699{
700
701 // If even number, return 0
702 if (a==2) return;
703
704 // If odd number, return the closest power of 2
705 if (a & 1)
706 {
707 Int_t b = 1;
708 while (b < fDim/2+1)
709 b <<= 1;
710
711 fDim = b;
712 *fLog << warn << "Dimension of Data is not a multiple of 2, will take only first "
713 << fDim << " entries! " << endl;
714 return;
715 }
716
717 CheckDim(a>>1);
718}
719
720TH1* MFFT::CheckHist(const TH1 *hist, const Int_t flag)
721{
722
723 TString name = hist->GetName();
724 name += " PSD";
725 TString title = hist->GetTitle();
726 title += " - Power Spectrum Density";
727
728 // number of entries
729 fDim = hist->GetNbinsX();
730 CheckDim(fDim);
731
732 // Step width
733 Double_t delta = hist->GetBinWidth(1);
734
735 // Nyquist frequency
736 Axis_t fcrit = 1./(2.*delta);
737 Axis_t low = 0.;
738 Axis_t up = fcrit;
739
740 switch (flag)
741 {
742 case 0:
743 return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
744 break;
745 case 1:
746 return new TH1D(name.Data(),title.Data(),fDim/2+1,low,up);
747 break;
748 default:
749 return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
750 break;
751 }
752}
Note: See TracBrowser for help on using the repository browser.