source: trunk/Mars/mtools/MFFT.cc@ 19594

Last change on this file since 19594 was 4958, checked in by gaug, 20 years ago
*** empty log message ***
File size: 26.9 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// Fourier-Transformation:
67// =======================
68
69// (taken from http://www.parasitaere-kapazitaeten.net/Pd/ft.htm)
70//
71// The Fourier-Transformation is a mathematical function that breaks
72// down a signal (like sound) into its frequency-spectrum as a set of
73// sinusoidal components, converting it from the Time Domain to the
74// Frequency Domain.
75//
76// In the Time Domain the signal x[ ] consists of N samples, labeled
77// from 0 to N-1. In the Frequency Domain the RFFT produces two signals
78// (signalvectors), treated as complex numbers representing the Real Part:
79// Re X[ ] and the Imaginary Part: Im X[ ]. They are seen as the Cosine-
80// und Sine-Components of the base frequencies. Each of these two signals
81// contains one more sample than the half of the original signal: N/2 + 1
82// samples. (this results from the fact, that the sine-components of the
83// first frequency (0) and the last (nyquist, N/2) are always 0). With the
84// complex Fourier-Transformation N complexe values are transformed to N
85// new complex values. For both it applies to: the Frequency Domain
86// contains exactly the same information as the Time-Domain.
87//
88// A Real FFT over 64 samples produces values for 33 cosine- and 33
89// sine-wave-amplitudes with the frequencies 0, 1, 2, 3, ..., 30, 31, 32.
90// The first value (frequency 0) is the DC (direct current), the other
91// values have to be seen in practice as factors of a
92// fundamental-frequency which can be calculated by dividing samplerate by
93// windowsize. The highest frequency is the nyquist-frequency
94// (samplerate/2).
95//
96//////////////////////////////////////////////////////////////////////////////
97
98#include "MFFT.h"
99
100#include <TMath.h>
101
102#include "MLog.h"
103#include "MLogManip.h"
104
105#include "MArrayD.h"
106#include "MArrayF.h"
107#include "MArrayI.h"
108
109ClassImp(MFFT);
110
111using namespace std;
112
113// ---------------------------------------------------------------------------
114//
115// Default Constructor
116// Initializes random number generator and default variables
117//
118MFFT::MFFT() : fDim(0)
119{
120}
121
122// --------------------------------------------------------------------------
123//
124// Destructor.
125//
126MFFT::~MFFT()
127{
128}
129
130void MFFT::TransformF(const Int_t isign, TArrayF &data)
131{
132
133 UInt_t n,mmax,m,j,istep,i;
134 Float_t wtemp,wr,wpr,wpi,wi,theta;
135 Float_t tempr,tempi;
136
137 Int_t nn = fDim/2;
138 n = nn << 1;
139
140 //
141 // The bit-reversal section of the routine:
142 // Exchange the two complex numbers
143 //
144 j=1;
145 for (i=1;i<n;i+=2) {
146 if (j > i) {
147 Swap(data[j-1],data[i-1]);
148 Swap(data[j],data[i]);
149 }
150 m=nn;
151 while (m >= 2 && j > m) {
152 j -= m;
153 m >>= 1;
154 }
155 j += m;
156 }
157 //
158 // Here begins the Danielson-Lanczos section of the routine
159 //
160 mmax=2;
161 while (n > mmax) { // Outer loop executed log_2(nn) times
162
163 istep = mmax << 1;
164 //
165 // Initialize the trigonometric recurrence:
166 //
167 theta = isign*(6.28318530717959/mmax);
168
169 wtemp = TMath::Sin(0.5*theta);
170 wpr = -2.0*wtemp*wtemp;
171 wpi = TMath::Sin(theta);
172
173 wr=1.0;
174 wi=0.0;
175
176 for (m=1; m<mmax; m+=2) {
177 for (i=m; i<=n; i+=istep) {
178 //
179 // The Danielson-Lanczos formula:
180 //
181 j = i+mmax;
182 tempr = wr*data[j-1] - wi*data[j];
183 tempi = wr*data[j] + wi*data[j-1];
184 data[j-1] = data[i-1] - tempr;
185 data[j] = data[i] - tempi;
186 data[i-1] += tempr;
187 data[i] += tempi;
188 }
189
190 //
191 // Trigonometric recurrence
192 //
193 wr = (wtemp=wr)*wpr - wi*wpi+wr;
194 wi = wi*wpr + wtemp*wpi+wi;
195
196 }
197 mmax=istep;
198 }
199}
200
201
202void MFFT::TransformD(const Int_t isign, TArrayD &data)
203{
204
205 UInt_t n,mmax,m,j,istep,i;
206 Double_t wtemp,wr,wpr,wpi,wi,theta;
207 Double_t tempr,tempi;
208
209 Int_t nn = fDim/2;
210 n = nn << 1;
211
212 //
213 // The bit-reversal section of the routine:
214 // Exchange the two complex numbers
215 //
216 j=1;
217 for (i=1;i<n;i+=2) {
218 if (j > i) {
219 Swap(data[j-1],data[i-1]);
220 Swap(data[j],data[i]);
221 }
222 m=nn;
223 while (m >= 2 && j > m) {
224 j -= m;
225 m >>= 1;
226 }
227 j += m;
228 }
229 //
230 // Here begins the Danielson-Lanczos section of the routine
231 //
232 mmax=2;
233 while (n > mmax) { // Outer loop executed log_2(nn) times
234
235 istep = mmax << 1;
236 //
237 // Initialize the trigonometric recurrence:
238 //
239 theta = isign*(6.28318530717959/mmax);
240
241 wtemp = TMath::Sin(0.5*theta);
242 wpr = -2.0*wtemp*wtemp;
243 wpi = TMath::Sin(theta);
244
245 wr=1.0;
246 wi=0.0;
247
248 for (m=1; m<mmax; m+=2) {
249 for (i=m; i<=n; i+=istep) {
250 //
251 // The Danielson-Lanczos formula:
252 //
253 j = i+mmax;
254 tempr = wr*data[j-1] - wi*data[j];
255 tempi = wr*data[j] + wi*data[j-1];
256 data[j-1] = data[i-1] - tempr;
257 data[j] = data[i] - tempi;
258 data[i-1] += tempr;
259 data[i] += tempi;
260 }
261
262 //
263 // Trigonometric recurrence
264 //
265 wr = (wtemp=wr)*wpr - wi*wpi+wr;
266 wi = wi*wpr + wtemp*wpi+wi;
267
268 }
269 mmax=istep;
270 }
271}
272
273//
274// Calculates the Fourier transform of a set of n real-valued data points.
275// Replaces this data (which is stored in array data[1..n]) by the positive
276// frequency half of its complex Fourier transform. The real-valued first
277// and last components of the complex transform are returned as elements
278// data[1] and data[2], respectively. n must be a power of 2. This routine
279// also calculates the inverse transform of a complex data array if it is
280// the transform of real data. (Result in this case mus be multiplied by
281// 2/n.). From NUMERICAL RECIPES IN C++.
282//
283void MFFT::RealFTF(const Int_t isign)
284{
285
286 Int_t i,i1,i2,i3,i4;
287 Float_t c1=0.5,c2,h1r,h1i,h2r,h2i;
288 Float_t wr,wi,wpr,wpi,wtemp,theta;
289
290 //
291 // Initialize the recurrence
292 //
293 theta = TMath::Pi() / (Double_t)(fDim>>1);
294
295 if(isign==1) // forward transform
296 {
297 c2 = -0.5;
298 TransformF(1,fDataF);
299 }
300 else // set up backward transform
301 {
302 c2 = 0.5;
303 theta = -theta;
304 }
305
306 wtemp = TMath::Sin(0.5*theta);
307 wpr = -2.0*wtemp*wtemp;
308 wpi = TMath::Sin(theta);
309
310 wr = 1.0 + wpr;
311 wi = wpi;
312
313 for(i=1;i<(fDim>>2);i++) // case i=0 done separately below
314 {
315
316 i2 = 1 + (i1 = i+i);
317 i4 = 1 + (i3 = fDim-i1);
318
319 //
320 // The two separate transforms are separated out of the data
321 //
322 h1r = c1*(fDataF[i1]+fDataF[i3]);
323 h1i = c1*(fDataF[i2]-fDataF[i4]);
324 h2r = -c2*(fDataF[i2]+fDataF[i4]);
325 h2i = c2*(fDataF[i1]-fDataF[i3]);
326
327 //
328 // Here, they are recombined to from the true transform
329 // of the orginal real data
330 //
331 fDataF[i1] = h1r + wr*h2r - wi*h2i;
332 fDataF[i2] = h1i + wr*h2i + wi*h2r;
333 fDataF[i3] = h1r - wr*h2r + wi*h2i;
334 fDataF[i4] = -h1i + wr*h2i + wi*h2r;
335
336 //
337 // The recurrence
338 //
339 wr = (wtemp=wr)*wpr - wi*wpi + wr;
340 wi = wi*wpr + wtemp*wpi + wi;
341 }
342
343 //
344 // Squeeze the first and last data together to get them all
345 // within the original array
346 //
347 if(isign==1)
348 {
349 fDataF[0] = (h1r=fDataF[0]) + fDataF[1];
350 fDataF[1] = h1r - fDataF[1];
351 }
352 else
353 {
354
355 fDataF[0] = c1*((h1r=fDataF[0]) + fDataF[1]);
356 fDataF[1] = c1*(h1r-fDataF[1]);
357
358 //
359 // The inverse transform for the case isign = -1
360 //
361 TransformF(-1,fDataF);
362
363 //
364 // normalize correctly (not done in original NR's)
365 //
366 for(i=1;i<=fDim;i++)
367 fDataF[i] *= (2./fDim);
368 }
369}
370void MFFT::RealFTD(const Int_t isign)
371{
372
373 Int_t i,i1,i2,i3,i4;
374 Float_t c1=0.5,c2,h1r,h1i,h2r,h2i;
375 Double_t wr,wi,wpr,wpi,wtemp,theta;
376
377 //
378 // Initialize the recurrence
379 //
380 theta=3.141592653589793/(Double_t) (fDim>>1);
381
382 if(isign==1) // forward transform
383 {
384 c2 = -0.5;
385 TransformD(1,fDataD);
386 }
387 else // set up backward transform
388 {
389 c2 = 0.5;
390 theta = -theta;
391 }
392
393 wtemp = TMath::Sin(0.5*theta);
394 wpr = -2.0*wtemp*wtemp;
395 wpi = TMath::Sin(theta);
396
397 wr = 1.0 + wpr;
398 wi = wpi;
399
400 for(i=1;i<(fDim>>2);i++) // case i=0 done separately below
401 {
402
403 i2 = 1 + (i1 = i+i);
404 i4 = 1 + (i3 = fDim-i1);
405
406 //
407 // The two separate transforms are separated out of the data
408 //
409 h1r = c1*(fDataD[i1]+fDataD[i3]);
410 h1i = c1*(fDataD[i2]-fDataD[i4]);
411 h2r = -c2*(fDataD[i2]+fDataD[i4]);
412 h2i = c2*(fDataD[i1]-fDataD[i3]);
413
414 //
415 // Here, they are recombined to from the true transform
416 // of the orginal real data
417 //
418 fDataD[i1] = h1r + wr*h2r - wi*h2i;
419 fDataD[i2] = h1i + wr*h2i + wi*h2r;
420 fDataD[i3] = h1r - wr*h2r + wi*h2i;
421 fDataD[i4] = -h1i + wr*h2i + wi*h2r;
422
423 //
424 // The recurrence
425 //
426 wr = (wtemp=wr)*wpr - wi*wpi + wr;
427 wi = wi*wpr + wtemp*wpi + wi;
428 }
429
430 //
431 // Squeeze the first and last data together to get them all
432 // within the original array
433 //
434 if(isign==1)
435 {
436 fDataD[0] = (h1r=fDataD[0]) + fDataD[1];
437 fDataD[1] = h1r - fDataD[1];
438 }
439 else
440 {
441
442 fDataD[0] = c1*((h1r=fDataD[0]) + fDataD[1]);
443 fDataD[1] = c1*(h1r-fDataD[1]);
444
445 //
446 // The inverse transform for the case isign = -1
447 //
448 TransformD(-1,fDataD);
449
450 //
451 // normalize correctly (not done in original NR's)
452 //
453 for(i=1;i<=fDim;i++)
454 fDataD[i] *= (2./fDim);
455 }
456}
457
458
459//
460// Fast Fourier Transform for float arrays
461//
462Float_t* MFFT::RealFunctionFFT(const Int_t n, const Float_t *data)
463{
464
465 fDim = n;
466 CheckDim(n);
467
468 fDataF.Set(fDim);
469 //
470 // Clone the array
471 //
472 for (Int_t i=0;i<fDim;i++)
473 fDataF[i] = data[i];
474
475 RealFTF(1);
476
477 return fDataF.GetArray();
478
479}
480
481//
482// Fast Inverse Fourier Transform for float arrays
483//
484Float_t* MFFT::RealFunctionIFFT(const Int_t n, const Float_t *data)
485{
486
487 fDim = n;
488 CheckDim(fDim);
489
490 fDataF.Set(fDim);
491 //
492 // Clone the array
493 //
494 for (Int_t i=0;i<fDim;i++)
495 fDataF[i] = data[i];
496
497 RealFTF(-1);
498
499 return fDataF.GetArray();
500
501}
502
503//
504// Fast Fourier Transform for double arrays
505//
506Double_t* MFFT::RealFunctionFFT(const Int_t n, const Double_t *data)
507{
508
509 fDim = n;
510 CheckDim(n);
511
512 fDataD.Set(fDim);
513 //
514 // Clone the array
515 //
516 for (Int_t i=0;i<fDim;i++)
517 fDataD[i] = data[i];
518
519 RealFTD(1);
520
521 return fDataD.GetArray();
522
523}
524
525//
526// Fast Inverse Fourier Transform for double arrays
527//
528Double_t* MFFT::RealFunctionIFFT(const Int_t n, const Double_t *data)
529{
530
531 fDim = n;
532 CheckDim(fDim);
533
534 fDataD.Set(fDim);
535 //
536 // Clone the array
537 //
538 for (Int_t i=0;i<fDim;i++)
539 fDataD[i] = data[i];
540
541 RealFTD(-1);
542
543 return fDataD.GetArray();
544
545}
546
547//
548// Fast Fourier Transform for TArrayF's
549//
550TArrayF* MFFT::RealFunctionFFT(const TArrayF *data)
551{
552
553 fDim = data->GetSize();
554 CheckDim(fDim);
555
556 fDataF.Set(fDim);
557 //
558 // Clone the array
559 //
560 for (Int_t i=0;i<fDim;i++)
561 fDataF[i] = data->At(i);
562
563 RealFTF(1);
564
565 return new TArrayF(fDim,fDataF.GetArray());
566
567}
568
569//
570// Inverse Fast Fourier Transform for TArrayF's
571//
572TArrayF* MFFT::RealFunctionIFFT(const TArrayF *data)
573{
574
575 fDim = data->GetSize();
576 CheckDim(fDim);
577
578 fDataF.Set(fDim);
579 //
580 // Clone the array
581 //
582 for (Int_t i=0;i<fDim;i++)
583 fDataF[i] = data->At(i);
584
585 RealFTF(-1);
586
587 return new TArrayF(fDim,fDataF.GetArray());
588}
589
590
591//
592// Fast Fourier Transform for TArrayD's
593//
594TArrayD* MFFT::RealFunctionFFT(const TArrayD *data)
595{
596
597 fDim = data->GetSize();
598 CheckDim(fDim);
599
600 fDataD.Set(fDim);
601 //
602 // Clone the array
603 //
604 for (Int_t i=0;i<fDim;i++)
605 fDataD[i] = data->At(i);
606
607 RealFTD(1);
608
609 return new TArrayD(fDim,fDataD.GetArray());
610
611}
612
613//
614// Inverse Fast Fourier Transform for TArrayD's
615//
616TArrayD* MFFT::RealFunctionIFFT(const TArrayD *data)
617{
618
619 fDim = data->GetSize();
620 CheckDim(fDim);
621
622 fDataD.Set(fDim);
623 //
624 // Clone the array
625 //
626 for (Int_t i=0;i<fDim;i++)
627 fDataD[i] = data->At(i);
628
629 RealFTD(-1);
630
631 return new TArrayD(fDim,fDataD.GetArray());
632}
633
634//----------------------------------------------------------
635//
636// Power Spectrum Density Calculation
637//
638TH1D* MFFT::PowerSpectrumDensity(const TH1D *hist)
639{
640
641 TH1D *newhist = (TH1D*)CheckHist(hist,1);
642
643 fDataD.Set(fDim);
644 //
645 // Copy the hist into an array
646 //
647 for (Int_t i=0;i<fDim;i++)
648 fDataD[i] = hist->GetBinContent(i);
649
650 RealFTD(1);
651
652 Int_t dim2 = fDim*fDim;
653 Double_t c02;
654 Double_t ck2;
655 Double_t cn2;
656 //
657 // Fill the new histogram:
658 //
659 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
660 // (stored in fData{0])
661 //
662 c02 = fDataD[0]*fDataD[0];
663 newhist->Fill(c02/dim2);
664 //
665 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)| + |C(N-k)|*|C(N-k)|)
666 //
667 for (Int_t k=2;k<fDim-2;k+=2)
668 {
669
670 Int_t ki = k+1;
671 ck2 = (fDataD[k]*fDataD[k] + fDataD[ki]*fDataD[ki]);
672 newhist->Fill(ck2/dim2);
673 }
674 //
675 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
676 // (stored in fData[1])
677 //
678 cn2 = (fDataD[1]*fDataD[1]);
679 newhist->Fill(cn2/dim2);
680
681 return newhist;
682}
683
684// -------------------------------------------------
685//
686// Power Spectrum Density calculation for TArrayF
687//
688TArrayF* MFFT::PowerSpectrumDensity(const TArrayF *array)
689{
690
691 fDim = array->GetSize();
692 CheckDim(fDim);
693
694 fDataF.Set(fDim);
695 //
696 // Copy the hist into an array
697 //
698 for (Int_t i=0;i<fDim;i++)
699 fDataF[i] = array->At(i);
700
701 RealFTF(1);
702
703 const Int_t dim2 = fDim*fDim;
704 const Int_t dim05 = fDim/2;
705 Float_t c02;
706 Float_t ck2;
707 Float_t cn2;
708
709 TArrayF *newarray = new TArrayF(dim05);
710
711 //
712 // Fill the new histogram:
713 //
714 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
715 //
716 c02 = (fDataF[0]*fDataF[0]);
717 newarray->AddAt(c02/dim2,0);
718 //
719 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
720 //
721 for (Int_t k=1;k<dim05-1;k++)
722 {
723 const Int_t k2 = k+k;
724 ck2 = (fDataF[k2]*fDataF[k2] + fDataF[k2+1]*fDataF[k2+1]);
725 newarray->AddAt(ck2/dim2,k);
726 }
727 //
728 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
729 //
730 cn2 = (fDataF[1]*fDataF[1]);
731 newarray->AddAt(cn2,dim05-1);
732
733 return newarray;
734}
735
736// -------------------------------------------------
737//
738// Power Spectrum Density calculation for TArrayI
739//
740TArrayF* MFFT::PowerSpectrumDensity(const TArrayI *array)
741{
742
743 fDim = array->GetSize();
744 CheckDim(fDim);
745
746 fDataF.Set(fDim);
747 //
748 // Copy the hist into an array
749 //
750 for (Int_t i=0;i<fDim;i++)
751 fDataF[i] = (Float_t)array->At(i);
752
753 RealFTF(1);
754
755 const Int_t dim2 = fDim*fDim;
756 const Int_t dim05 = fDim/2;
757 Float_t c02;
758 Float_t ck2;
759 Float_t cn2;
760
761 TArrayF *newarray = new TArrayF(dim05);
762
763 //
764 // Fill the new histogram:
765 //
766 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
767 //
768 c02 = (fDataF[0]*fDataF[0]);
769 newarray->AddAt(c02/dim2,0);
770 //
771 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
772 //
773 for (Int_t k=1;k<dim05-1;k++)
774 {
775 const Int_t k2 = k+k;
776 ck2 = (fDataF[k2]*fDataF[k2] + fDataF[k2+1]*fDataF[k2+1]);
777 newarray->AddAt(ck2/dim2,k);
778 }
779 //
780 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
781 //
782 cn2 = (fDataF[1]*fDataF[1]);
783 newarray->AddAt(cn2,dim05-1);
784
785 return newarray;
786}
787
788
789// -------------------------------------------------
790//
791// Power Spectrum Density calculation for TArrayD
792//
793TArrayD* MFFT::PowerSpectrumDensity(const TArrayD *array)
794{
795
796 fDim = array->GetSize();
797 CheckDim(fDim);
798
799 fDataD.Set(fDim);
800 //
801 // Copy the hist into an array
802 //
803 for (Int_t i=0;i<fDim;i++)
804 fDataD[i] = array->At(i);
805
806 RealFTD(1);
807
808 const Int_t dim2 = fDim*fDim;
809 const Int_t dim05 = fDim/2;
810 Float_t c02;
811 Float_t ck2;
812 Float_t cn2;
813
814 TArrayD *newarray = new TArrayD(dim05);
815
816 //
817 // Fill the new histogram:
818 //
819 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
820 //
821 c02 = (fDataD[0]*fDataD[0]);
822 newarray->AddAt(c02/dim2,0);
823 //
824 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
825 //
826 for (Int_t k=1;k<dim05-1;k++)
827 {
828 const Int_t k2 = k+k;
829 ck2 = (fDataD[k2]*fDataD[k2] + fDataD[k2+1]*fDataD[k2+1]);
830 newarray->AddAt(ck2/dim2,k);
831 }
832 //
833 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
834 //
835 cn2 = (fDataD[1]*fDataD[1]);
836 newarray->AddAt(cn2,dim05-1);
837
838 return newarray;
839}
840
841// -------------------------------------------------
842//
843// Power Spectrum Density calculation for MArrayF
844// The difference to the TArrayF versions is that
845// the resulting array has two entries less, namely
846// the first and last one are skipped!
847//
848MArrayF* MFFT::PowerSpectrumDensity(const MArrayF *array)
849{
850
851 fDim = array->GetSize();
852 CheckDim(fDim);
853
854 fDataF.Set(fDim);
855 //
856 // Copy the hist into an array
857 //
858 for (Int_t i=0;i<fDim;i++)
859 fDataF[i] = array->At(i);
860
861 RealFTF(1);
862
863 const Int_t dim2 = fDim*fDim;
864 const Int_t dim05 = fDim/2;
865 Float_t ck2;
866
867 MArrayF *newarray = new MArrayF(dim05-2);
868
869 //
870 // Fill the new histogram:
871 //
872 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
873 //
874 // c02 = (fDataF[0]*fDataF[0]);
875 // newarray->AddAt(c02/dim2,0);
876 //
877 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
878 //
879 for (Int_t k=1;k<dim05-1;k++)
880 {
881 const Int_t k2 = k+k;
882 ck2 = (fDataF[k2]*fDataF[k2] + fDataF[k2+1]*fDataF[k2+1]);
883 newarray->AddAt(ck2/dim2,k-1);
884 }
885 //
886 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
887 //
888 // cn2 = (fDataF[1]*fDataF[1]);
889 // newarray->AddAt(cn2,dim05-1);
890
891 return newarray;
892}
893
894//-----------------------------------------------------
895//
896// Power Spectrum Density calculation for MArrayI
897// The difference to the TArrayI versions is that
898// the resulting array has two entries less, namely
899// the first and last one are skipped!
900//
901MArrayF* MFFT::PowerSpectrumDensity(const MArrayI *array)
902{
903
904 fDim = array->GetSize();
905 CheckDim(fDim);
906
907 fDataF.Set(fDim);
908 //
909 // Copy the hist into an array
910 //
911 for (Int_t i=0;i<fDim;i++)
912 fDataF[i] = (Float_t)array->At(i);
913
914 RealFTF(1);
915
916 const Int_t dim2 = fDim*fDim;
917 const Int_t dim05 = fDim/2;
918 Float_t ck2;
919
920 MArrayF *newarray = new MArrayF(dim05-2);
921
922 //
923 // Fill the new histogram:
924 //
925 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
926 //
927 // c02 = (fDataF[0]*fDataF[0]);
928 // newarray->AddAt(c02/dim2,0);
929 //
930 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
931 //
932 for (Int_t k=1;k<dim05-1;k++)
933 {
934 const Int_t k2 = k+k;
935 ck2 = (fDataF[k2]*fDataF[k2] + fDataF[k2+1]*fDataF[k2+1]);
936 newarray->AddAt(ck2/dim2,k-1);
937 }
938 //
939 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
940 //
941 // cn2 = (fDataF[1]*fDataF[1]);
942 // newarray->AddAt(cn2,dim05-1);
943
944 return newarray;
945}
946
947// -------------------------------------------------
948//
949// Power Spectrum Density calculation for MArrayD
950// The difference to the TArrayI versions is that
951// the resulting array has two entries less, namely
952// the first and last one are skipped!
953//
954MArrayD* MFFT::PowerSpectrumDensity(const MArrayD *array)
955{
956
957 fDim = array->GetSize();
958 CheckDim(fDim);
959
960 fDataD.Set(fDim);
961 //
962 // Copy the hist into an array
963 //
964 for (Int_t i=0;i<fDim;i++)
965 fDataD[i] = array->At(i);
966
967 RealFTD(1);
968
969 const Int_t dim2 = fDim*fDim;
970 const Int_t dim05 = fDim/2;
971 Float_t ck2;
972
973 MArrayD *newarray = new MArrayD(dim05-2);
974
975 //
976 // Fill the new histogram:
977 //
978 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
979 //
980 // c02 = (fDataD[0]*fDataD[0]);
981 // newarray->AddAt(c02/dim2,0);
982 //
983 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
984 //
985 for (Int_t k=1;k<dim05-1;k++)
986 {
987 const Int_t k2 = k+k;
988 ck2 = (fDataD[k2]*fDataD[k2] + fDataD[k2+1]*fDataD[k2+1]);
989 newarray->AddAt(ck2/dim2,k-1);
990 }
991 //
992 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
993 //
994 // cn2 = (fDataD[1]*fDataD[1]);
995 // newarray->AddAt(cn2,dim05-1);
996
997 return newarray;
998}
999
1000// -----------------------------------------------
1001//
1002// Power Spectrum Density calculation for TH1
1003//
1004TH1F* MFFT::PowerSpectrumDensity(const TH1 *hist)
1005{
1006
1007 TH1F *newhist = (TH1F*)CheckHist(hist,0);
1008
1009 fDataF.Set(fDim);
1010 //
1011 // Copy the hist into an array
1012 //
1013 for (Int_t i=0;i<fDim;i++)
1014 fDataF[i] = hist->GetBinContent(i);
1015
1016 RealFTF(1);
1017
1018 Int_t dim2 = fDim*fDim;
1019 Float_t c02;
1020 Float_t ck2;
1021 Float_t cn2;
1022 //
1023 // Fill the new histogram:
1024 //
1025 // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
1026 //
1027 c02 = (fDataF[0]*fDataF[0]);
1028 newhist->Fill(0.,c02/dim2);
1029 //
1030 // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
1031 //
1032 for (Int_t k=2;k<fDim;k+=2)
1033 {
1034 ck2 = (fDataF[k]*fDataF[k] + fDataF[k+1]*fDataF[k+1]);
1035 newhist->Fill(k/2.,ck2/dim2);
1036 }
1037 //
1038 // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
1039 //
1040 cn2 = (fDataF[1]*fDataF[1]);
1041 newhist->Fill(fDim/2.-1.,cn2/dim2);
1042
1043 return newhist;
1044}
1045
1046
1047//
1048// Power Spectrum Density calculation for TH1I
1049//
1050TH1F* MFFT::PowerSpectrumDensity(const TH1F *hist)
1051{
1052 return PowerSpectrumDensity((TH1*)hist);
1053}
1054
1055//
1056// Power Spectrum Density calculation for TH1I
1057//
1058TH1F* MFFT::PowerSpectrumDensity(const TH1I *hist)
1059{
1060 return PowerSpectrumDensity((TH1*)hist);
1061}
1062
1063
1064void MFFT::CheckDim(Int_t a)
1065{
1066
1067 // If even number, return 0
1068 if (a==2) return;
1069
1070 // If odd number, return the closest power of 2
1071 if (a & 1)
1072 {
1073 Int_t b = 1;
1074 while (b < fDim/2+1)
1075 b <<= 1;
1076
1077 fDim = b;
1078 // gLog << warn << "Dimension of Data is not a multiple of 2, will take only first "
1079 // << fDim << " entries! " << endl;
1080 return;
1081 }
1082
1083 CheckDim(a>>1);
1084}
1085
1086TH1* MFFT::CheckHist(const TH1 *hist, const Int_t flag)
1087{
1088
1089 // number of entries
1090 fDim = hist->GetNbinsX();
1091 CheckDim(fDim);
1092
1093 // Step width
1094 Double_t delta = hist->GetBinWidth(1);
1095
1096 // Nyquist frequency
1097 Axis_t fcrit = 1./(2.*delta);
1098 Axis_t low = -0.5;
1099 Axis_t up = fcrit;
1100
1101 switch (flag)
1102 {
1103 case 0:
1104 return new TH1F(Form("%s%s",hist->GetName()," PSD"),
1105 Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
1106 fDim/2,low,up);
1107 break;
1108 case 1:
1109 return new TH1D(Form("%s%s",hist->GetName()," PSD"),
1110 Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
1111 fDim/2,low,up);
1112 break;
1113 default:
1114 return new TH1F(Form("%s%s",hist->GetName()," PSD"),
1115 Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
1116 fDim/2,low,up);
1117 break;
1118 }
1119}
1120
1121//
1122// Real function spectrum with data windowing
1123//
1124TArrayF* MFFT::RealFunctionSpectrum(const TArrayF *data)
1125{
1126
1127 fDim = data->GetSize();
1128 CheckDim(fDim);
1129
1130 fDataF.Set(fDim);
1131 //
1132 // Copy the hist into an array
1133 //
1134 for (Int_t i=0;i<fDim;i++)
1135 fDataF[i] = data->At(i);
1136
1137 fWindowF.Set(fDim);
1138
1139 Int_t dim2 = fDim/2;
1140
1141 TArrayF *power = new TArrayF(dim2);
1142
1143 //
1144 // Start program spctrm from NR's
1145 //
1146 Float_t w, facp, facm, sumw=0.;
1147
1148 facm = dim2;
1149 facp = 1./dim2;
1150
1151 for (Int_t j=0;j<dim2;j++)
1152 {
1153 Int_t j2 = j+j;
1154 w = ApplyWindow(j,facm,facp);
1155 sumw += w*w;
1156 fWindowF[j2] = fDataF[j]*w;
1157 fWindowF[j2+1] = fDataF[dim2+j]*w;
1158 }
1159
1160 TransformF(1,fWindowF);
1161
1162 power->AddAt(fWindowF[0]*fWindowF[0] + fWindowF[1]*fWindowF[1],0);
1163
1164 // power->AddAt(fWindowF[0]*fWindowF[0],0);
1165 // power->AddAt(fWindowF[1]*fWindowF[1],dim2-1);
1166
1167
1168 for (Int_t j=1;j<dim2;j++)
1169 // for (Int_t j=1;j<dim2;j++)
1170 {
1171 Int_t j2 = j+j;
1172 Float_t buf = fWindowF[j2+1] *fWindowF[j2+1]
1173 + fWindowF[j2 ] *fWindowF[j2 ]
1174 + fWindowF[fDim-j2+1]*fWindowF[fDim-j2+1]
1175 + fWindowF[fDim-j2 ]*fWindowF[fDim-j2 ] ;
1176 power->AddAt(buf/sumw/(fDim+fDim),j);
1177 }
1178
1179 return power;
1180
1181}
Note: See TracBrowser for help on using the repository browser.