Ignore:
Timestamp:
02/08/04 20:43:07 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mtools/MFFT.cc

    r2855 r3054  
    9494}
    9595
    96 void MFFT::TransformF(const Int_t isign)
     96void MFFT::TransformF(const Int_t isign, TArrayF &data)
    9797{
    9898 
    … …  
    111111  for (i=1;i<n;i+=2) {
    112112    if (j > i) {
    113       Swap(fDataF[j-1],fDataF[i-1]);
    114       Swap(fDataF[j],fDataF[i]);
     113      Swap(data[j-1],data[i-1]);
     114      Swap(data[j],data[i]);
    115115    }
    116116    m=nn;
    … …  
    146146        //
    147147        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;
     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;
    154154      }
    155155
    … …  
    166166
    167167
    168 void MFFT::TransformD(const Int_t isign)
     168void MFFT::TransformD(const Int_t isign, TArrayD &data)
    169169{
    170170 
    … …  
    183183  for (i=1;i<n;i+=2) {
    184184    if (j > i) {
    185       Swap(fDataD[j-1],fDataD[i-1]);
    186       Swap(fDataD[j],fDataD[i]);
     185      Swap(data[j-1],data[i-1]);
     186      Swap(data[j],data[i]);
    187187    }
    188188    m=nn;
    … …  
    218218        //
    219219        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;
     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;
    226226      }
    227227
    … …  
    262262    {
    263263      c2    = -0.5;
    264       TransformF(1);
     264      TransformF(1,fDataF);
    265265    }
    266266  else         // set up backward transform
    … …  
    325325      // The inverse transform for the case isign = -1
    326326      //
    327       TransformF(-1); 
     327      TransformF(-1,fDataF); 
    328328
    329329      //
    … …  
    349349    {
    350350      c2    = -0.5;
    351       TransformD(1);
     351      TransformD(1,fDataD);
    352352    }
    353353  else         // set up backward transform
    … …  
    412412      // The inverse transform for the case isign = -1
    413413      //
    414       TransformD(-1); 
     414      TransformD(-1,fDataD); 
    415415     
    416416      //
    … …  
    648648}
    649649
    650 
    651 //
    652 // Power Spectrum Density calculation
    653 //
    654 TH1F* MFFT::PowerSpectrumDensity(const TH1F *hist)
     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=0;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=0;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)
    655755{
    656756
    … …  
    676776  //
    677777  c02 = (fDataF[0]*fDataF[0]);
    678   newhist->Fill(c02/dim2);
     778  newhist->Fill(0.,c02/dim2);
    679779  //
    680780  // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
    … …  
    682782  for (Int_t k=2;k<fDim;k+=2)
    683783    {
    684       Int_t ki  = k+1;
    685       ck2 = (fDataF[k]*fDataF[k] + fDataF[ki]*fDataF[ki]);
    686       newhist->Fill(ck2/dim2);
     784      ck2 = (fDataF[k]*fDataF[k] + fDataF[k+1]*fDataF[k+1]);
     785      newhist->Fill(k/2.,ck2/dim2);
    687786    }
    688787  //
    … …  
    690789  //
    691790  cn2 = (fDataF[1]*fDataF[1]);
    692   newhist->Fill(cn2/dim2);
     791  newhist->Fill(fDim/2.-1.,cn2/dim2);
    693792 
    694793  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);
    695811}
    696812
    … …  
    721837{
    722838 
    723   TString name = hist->GetName();
    724   name += " PSD";
    725   TString title = hist->GetTitle();
    726   title += " - Power Spectrum Density";
    727 
    728839  // number of entries
    729840  fDim = hist->GetNbinsX();
    … …  
    735846  // Nyquist frequency
    736847  Axis_t fcrit = 1./(2.*delta);
    737   Axis_t low = 0.;
     848  Axis_t low = -0.5;
    738849  Axis_t up  = fcrit;
    739850
    … …  
    741852    {
    742853    case 0:
    743       return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
     854      return new TH1F(Form("%s%s",hist->GetName()," PSD"),
     855                      Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
     856                      fDim/2,low,up);
    744857      break;
    745858    case 1:
    746       return new TH1D(name.Data(),title.Data(),fDim/2+1,low,up);
     859      return new TH1D(Form("%s%s",hist->GetName()," PSD"),
     860                      Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
     861                      fDim/2,low,up);
    747862      break;
    748863    default:
    749       return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
     864      return new TH1F(Form("%s%s",hist->GetName()," PSD"),
     865                      Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
     866                      fDim/2,low,up);
    750867      break;
    751868    }
    752869}
     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 TracChangeset for help on using the changeset viewer.