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

Last change on this file since 3074 was 3074, checked in by gaug, 21 years ago
*** empty log message ***
File size: 20.8 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, TArrayF &data)
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(data[j-1],data[i-1]);
114 Swap(data[j],data[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*data[j-1] - wi*data[j];
149 tempi = wr*data[j] + wi*data[j-1];
150 data[j-1] = data[i-1] - tempr;
151 data[j] = data[i] - tempi;
152 data[i-1] += tempr;
153 data[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, TArrayD &data)
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(data[j-1],data[i-1]);
186 Swap(data[j],data[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*data[j-1] - wi*data[j];
221 tempi = wr*data[j] + wi*data[j-1];
222 data[j-1] = data[i-1] - tempr;
223 data[j] = data[i] - tempi;
224 data[i-1] += tempr;
225 data[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,fDataF);
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,fDataF);
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,fDataD);
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,fDataD);
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// Power Spectrum Density calculation for TArrayF
652//
653TArrayF* MFFT::PowerSpectrumDensity(const TArrayF *array)
654{
655
656 fDim = array->GetSize();
657 CheckDim(fDim);
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] = array->At(i);
665
666 RealFTF(1);
667
668 const Int_t dim2 = fDim*fDim;
669 const Int_t dim05 = fDim/2;
670 Float_t c02;
671 Float_t ck2;
672 Float_t cn2;
673
674 TArrayF *newarray = new TArrayF(dim05);
675
676 //
677 // Fill the new histogram:
678 //
679 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
680 //
681 c02 = (fDataF[0]*fDataF[0]);
682 // newarray->AddAt(c02/dim2,0);
683 //
684 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
685 //
686 for (Int_t k=1;k<dim05-1;k++)
687 {
688 const Int_t k2 = k+k;
689 ck2 = (fDataF[k2]*fDataF[k2] + fDataF[k2+1]*fDataF[k2+1]);
690 newarray->AddAt(ck2/dim2,k);
691 }
692 //
693 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
694 //
695 cn2 = (fDataF[1]*fDataF[1]);
696 // newarray->AddAt(cn2,dim05-1);
697
698 return newarray;
699}
700
701
702TArrayD* MFFT::PowerSpectrumDensity(const TArrayD *array)
703{
704
705 fDim = array->GetSize();
706 CheckDim(fDim);
707
708 fDataD.Set(fDim);
709 //
710 // Copy the hist into an array
711 //
712 for (Int_t i=0;i<fDim;i++)
713 fDataD[i] = array->At(i);
714
715 RealFTD(1);
716
717 const Int_t dim2 = fDim*fDim;
718 const Int_t dim05 = fDim/2;
719 Float_t c02;
720 Float_t ck2;
721 Float_t cn2;
722
723 TArrayD *newarray = new TArrayD(dim05);
724
725 //
726 // Fill the new histogram:
727 //
728 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
729 //
730 c02 = (fDataD[0]*fDataD[0]);
731 // newarray->AddAt(c02/dim2,0);
732 //
733 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
734 //
735 for (Int_t k=1;k<dim05-1;k++)
736 {
737 const Int_t k2 = k+k;
738 ck2 = (fDataD[k2]*fDataD[k2] + fDataD[k2+1]*fDataD[k2+1]);
739 newarray->AddAt(ck2/dim2,k);
740 }
741 //
742 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
743 //
744 cn2 = (fDataD[1]*fDataD[1]);
745 // newarray->AddAt(cn2,dim05-1);
746
747 return newarray;
748}
749
750
751//
752// Power Spectrum Density calculation for TH1
753//
754TH1F* MFFT::PowerSpectrumDensity(const TH1 *hist)
755{
756
757 TH1F *newhist = (TH1F*)CheckHist(hist,0);
758
759 fDataF.Set(fDim);
760 //
761 // Copy the hist into an array
762 //
763 for (Int_t i=0;i<fDim;i++)
764 fDataF[i] = hist->GetBinContent(i);
765
766 RealFTF(1);
767
768 Int_t dim2 = fDim*fDim;
769 Float_t c02;
770 Float_t ck2;
771 Float_t cn2;
772 //
773 // Fill the new histogram:
774 //
775 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
776 //
777 c02 = (fDataF[0]*fDataF[0]);
778 newhist->Fill(0.,c02/dim2);
779 //
780 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
781 //
782 for (Int_t k=2;k<fDim;k+=2)
783 {
784 ck2 = (fDataF[k]*fDataF[k] + fDataF[k+1]*fDataF[k+1]);
785 newhist->Fill(k/2.,ck2/dim2);
786 }
787 //
788 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
789 //
790 cn2 = (fDataF[1]*fDataF[1]);
791 newhist->Fill(fDim/2.-1.,cn2/dim2);
792
793 return newhist;
794}
795
796
797//
798// Power Spectrum Density calculation for TH1I
799//
800TH1F* MFFT::PowerSpectrumDensity(const TH1F *hist)
801{
802 return PowerSpectrumDensity((TH1*)hist);
803}
804
805//
806// Power Spectrum Density calculation for TH1I
807//
808TH1F* MFFT::PowerSpectrumDensity(const TH1I *hist)
809{
810 return PowerSpectrumDensity((TH1*)hist);
811}
812
813
814void MFFT::CheckDim(Int_t a)
815{
816
817 // If even number, return 0
818 if (a==2) return;
819
820 // If odd number, return the closest power of 2
821 if (a & 1)
822 {
823 Int_t b = 1;
824 while (b < fDim/2+1)
825 b <<= 1;
826
827 fDim = b;
828 // gLog << warn << "Dimension of Data is not a multiple of 2, will take only first "
829 // << fDim << " entries! " << endl;
830 return;
831 }
832
833 CheckDim(a>>1);
834}
835
836TH1* MFFT::CheckHist(const TH1 *hist, const Int_t flag)
837{
838
839 // number of entries
840 fDim = hist->GetNbinsX();
841 CheckDim(fDim);
842
843 // Step width
844 Double_t delta = hist->GetBinWidth(1);
845
846 // Nyquist frequency
847 Axis_t fcrit = 1./(2.*delta);
848 Axis_t low = -0.5;
849 Axis_t up = fcrit;
850
851 switch (flag)
852 {
853 case 0:
854 return new TH1F(Form("%s%s",hist->GetName()," PSD"),
855 Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
856 fDim/2,low,up);
857 break;
858 case 1:
859 return new TH1D(Form("%s%s",hist->GetName()," PSD"),
860 Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
861 fDim/2,low,up);
862 break;
863 default:
864 return new TH1F(Form("%s%s",hist->GetName()," PSD"),
865 Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
866 fDim/2,low,up);
867 break;
868 }
869}
870
871//
872// Real function spectrum with data windowing
873//
874TArrayF* MFFT::RealFunctionSpectrum(const TArrayF *data)
875{
876
877 fDim = data->GetSize();
878 CheckDim(fDim);
879
880 fDataF.Set(fDim);
881 //
882 // Copy the hist into an array
883 //
884 for (Int_t i=0;i<fDim;i++)
885 fDataF[i] = data->At(i);
886
887 fWindowF.Set(fDim);
888
889 Int_t dim2 = fDim/2;
890
891 TArrayF *power = new TArrayF(dim2);
892
893 //
894 // Start program spctrm from NR's
895 //
896 Float_t w, facp, facm, sumw=0.;
897
898 facm = dim2;
899 facp = 1./dim2;
900
901 for (Int_t j=0;j<dim2;j++)
902 {
903 Int_t j2 = j+j;
904 w = ApplyWindow(j,facm,facp);
905 sumw += w*w;
906 fWindowF[j2] = fDataF[j]*w;
907 fWindowF[j2+1] = fDataF[dim2+j]*w;
908 }
909
910 TransformF(1,fWindowF);
911
912 power->AddAt(fWindowF[0]*fWindowF[0] + fWindowF[1]*fWindowF[1],0);
913
914 // power->AddAt(fWindowF[0]*fWindowF[0],0);
915 // power->AddAt(fWindowF[1]*fWindowF[1],dim2-1);
916
917
918 for (Int_t j=1;j<dim2;j++)
919 // for (Int_t j=1;j<dim2;j++)
920 {
921 Int_t j2 = j+j;
922 Float_t buf = fWindowF[j2+1] *fWindowF[j2+1]
923 + fWindowF[j2 ] *fWindowF[j2 ]
924 + fWindowF[fDim-j2+1]*fWindowF[fDim-j2+1]
925 + fWindowF[fDim-j2 ]*fWindowF[fDim-j2 ] ;
926 power->AddAt(buf/sumw/(fDim+fDim),j);
927 }
928
929 return power;
930
931}
932
Note: See TracBrowser for help on using the repository browser.