NcDSP.cxx

Go to the documentation of this file.

 
#include "NcDSP.h"
#include "Riostream.h"
 
ClassImp(NcDSP); // Class implementation to enable ROOT I/O
 
NcDSP::NcDSP(const char* name,const char* title) : TNamed(name,title)
{
 
 fProc=0;
 Reset();
 fSample=0;
 fKeepOutput=kFALSE;
}

NcDSP::~NcDSP()
{
 
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
}

NcDSP::NcDSP(const NcDSP& q) : TNamed(q)
{
 
 fProc=0;
 fN=q.fN;
 fNwf=q.fNwf;
 fReIn=q.fReIn;
 fImIn=q.fImIn;
 fReOut=q.fReOut;
 fImOut=q.fImOut;
 fHisto=q.fHisto;
 fWaveform=q.fWaveform;
 fNorm=q.fNorm;
 fSample=q.fSample;
 fKeepOutput=q.fKeepOutput;
}

void NcDSP::Reset()
{
 
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
 fN=0;
 fNwf=0;
 fReIn.Set(0);
 fImIn.Set(0);
 if (!fKeepOutput)
 {
  fReOut.Set(0);
  fImOut.Set(0);
  fHisto.Set(0);
 }
 fNorm="NONE";
}

void NcDSP::SetSamplingFrequency(Float_t f)
{
 
 fSample=f;
}

Float_t NcDSP::GetSamplingFrequency() const
{
 
 return fSample;
}

void NcDSP::Load(Int_t n,Double_t* re,Double_t* im,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (n<1) return;
 
 fN=n;
 if (re) fReIn.Set(n);
 if (im) fImIn.Set(n);
 
 for (Int_t i=0; i<n; i++)
 {
  if (re) fReIn.SetAt(re[i],i);
  if (im) fImIn.SetAt(im[i],i);
 }
}

void NcDSP::Load(TArray* re,TArray* im,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 Int_t nre=0;
 Int_t nim=0;
 if (re) nre=re->GetSize();
 if (im) nim=im->GetSize();
 
 Int_t n=nre;
 if (!n) n=nim;
 if (nre && nim>nre) n=nre;
 if (nim && nre>nim) n=nim;
 
 if (n<1) return;
 
 fN=n;
 if (nre) fReIn.Set(n);
 if (nim) fImIn.Set(n);
 
 Double_t valre=0;
 Double_t valim=0;
 for (Int_t i=0; i<n; i++)
 {
  if (nre)
  {
   valre=re->GetAt(i);
   fReIn.SetAt(valre,i);
  }
  if (nim)
  {
   valim=im->GetAt(i);
   fImIn.SetAt(valim,i);
  }
 }
}

void NcDSP::Load(NcSample* s,Int_t i,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (!s)
 {
  cout << " *" << ClassName() <<"::Load* No input sample specified." << endl;
  return;
 }
 
 Int_t n=s->GetN();
 Int_t store=s->GetStoreMode();
 Int_t dim=s->GetDimension();
 
 if (n<1 || !store || dim<1 || i<1 || i>dim)
 {
  cout << " *" << ClassName() <<"::Load* Inconsistent input for NcSample treatment." << endl;
  cout << " Store Mode:" << store << " Entries:" << n << " Dimension:" << dim << " i:" << i << " f:" << fSample << endl; 
  return;
 }
 
 fN=n;
 fReIn.Set(n);
 
 Double_t val=0;
 for (Int_t idx=1; idx<=n; idx++)
 {
  val=s->GetEntry(idx,i);
  fReIn.SetAt(val,idx-1);
 }
}

void NcDSP::Load(NcSample* s,TString name,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (!s)
 {
  cout << " *" << ClassName() <<"::Load* No input sample specified." << endl;
  return;
 }
 
 Int_t n=s->GetN();
 Int_t store=s->GetStoreMode();
 Int_t dim=s->GetDimension();
 Int_t i=s->GetIndex(name);
 
 
 if (n<1 || !store || dim<1 || !i)
 {
  cout << " *" << ClassName() <<"::Load* Inconsistent input for NcSample treatment." << endl;
  cout << " Store Mode:" << store << " Entries:" << n << " Dimension:" << dim << " name:" << name << " f:" << fSample << endl; 
  return;
 }
 
 Load(s,i,f);
}

void NcDSP::Load(TH1* h,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (!h)
 {
  cout << " *" << ClassName() <<"::Load* No input histogram specified." << endl;
  return;
 }
 
 Int_t nbins=h->GetNbinsX();
 Double_t nentries=h->GetEntries();
 
 if (!nbins || !nentries)
 {
  cout << " *" << ClassName() <<"::Load* Inconsistent input for histogram treatment." << endl;
  cout << " Nbins:" << nbins << " Nentries:" << nentries << " f:" << fSample << endl; 
  return;
 }
 
 fN=nbins;
 fReIn.Set(nbins);
 
 Double_t val=0;
 for (Int_t i=1; i<=nbins; i++)
 {
  val=h->GetBinContent(i);
  fReIn.SetAt(val,i-1);
 }
}

void NcDSP::Load(TGraph* gr,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (!gr)
 {
  cout << " *" << ClassName() <<"::Load* No input TGraph object specified." << endl;
  return;
 }
 
 Int_t n=gr->GetN();
 
 if (!n)
 {
  cout << " *" << ClassName() <<"::Load* Inconsistent input for TGraph treatment." << endl;
  cout << " n:" << n << " f:" << fSample << endl; 
  return;
 }
 
 fN=n;
 fReIn.Set(n);
 
 gr->Sort();
 
 Double_t x=0;
 Double_t y=0;
 for (Int_t i=0; i<n; i++)
 {
  gr->GetPoint(i,x,y);
  fReIn.SetAt(y,i);
 }
}

void NcDSP::LoadResult()
{
 
 fReIn=fReOut;
 fImIn=fImOut;
 fReOut.Set(0);
 fImOut.Set(0);
 fHisto.Set(0);
}

void NcDSP::SetWaveform(Int_t n,Double_t* h,Float_t f)
{
 
 fWaveform.Set(0);
 fNwf=0;
 
 if (f>=0) fSample=f;
 
 if (!h)
 {
  cout << " *" << ClassName() <<"::SetWaveform* No data array specified." << endl;
  return;
 }
 
 if (n<1) return;
 
 fWaveform.Set(n);
 fNwf=n;
 
 for (Int_t i=0; i<n; i++)
 {
  fWaveform[i]=h[i];
 }
}

void NcDSP::SetWaveform(TArray* h,Float_t f)
{
 
 fWaveform.Set(0);
 fNwf=0;
 
 if (f>=0) fSample=f;
 
 if (!h)
 {
  cout << " *" << ClassName() <<"::SetWaveform* No data array specified." << endl;
  return;
 }
 
 Int_t n=0;
 if (h) n=h->GetSize();
 
 if (n<1) return;
 
 fWaveform.Set(n);
 fNwf=n;
 
 Double_t val=0;
 for (Int_t i=0; i<n; i++)
 {
  val=h->GetAt(i);
  fWaveform[i]=val;
 }
}

void NcDSP::SetWaveform(NcSample* s,Int_t i,Float_t f)
{
 
 fWaveform.Set(0);
 fNwf=0;
 
 if (f>=0) fSample=f;
 
 if (!s)
 {
  cout << " *" << ClassName() <<"::SetWaveform* No data sample specified." << endl;
  return;
 }
 
 Int_t n=s->GetN();
 Int_t store=s->GetStoreMode();
 Int_t dim=s->GetDimension();
 
 if (n<1 || !store || dim<1 || i<1 || i>dim)
 {
  cout << " *" << ClassName() <<"::SetWaveform* Inconsistent input for NcSample treatment." << endl;
  cout << " Store Mode:" << store << " Entries:" << n << " Dimension:" << dim << " i:" << i << " f:" << fSample << endl; 
  return;
 }
 
 fWaveform.Set(n);
 fNwf=n;
 
 Double_t val=0;
 for (Int_t idx=1; idx<=n; idx++)
 {
  val=s->GetEntry(idx,i);
  fWaveform[idx-1]=val;
 }
}

void NcDSP::SetWaveform(NcSample* s,TString name,Float_t f)
{
 
 fWaveform.Set(0);
 fNwf=0;
 
 if (f>=0) fSample=f;
 
 if (!s)
 {
  cout << " *" << ClassName() <<"::SetWaveform* No data sample specified." << endl;
  return;
 }
 
 Int_t n=s->GetN();
 Int_t store=s->GetStoreMode();
 Int_t dim=s->GetDimension();
 Int_t i=s->GetIndex(name);
 
 
 if (n<1 || !store || dim<1 || !i)
 {
  cout << " *" << ClassName() <<"::SetWaveform* Inconsistent input for NcSample treatment." << endl;
  cout << " Store Mode:" << store << " Entries:" << n << " Dimension:" << dim << " name:" << name << " f:" << fSample << endl; 
  return;
 }
 
 SetWaveform(s,i,f);
}

void NcDSP::SetWaveform(TH1* h,Float_t f)
{
 
 fWaveform.Set(0);
 fNwf=0;
 
 if (f>=0) fSample=f;
 
 if (!h)
 {
  cout << " *" << ClassName() <<"::SetWaveform* No input histogram specified." << endl;
  return;
 }
 
 Int_t nbins=h->GetNbinsX();
 Double_t nentries=h->GetEntries();
 
 if (!nbins || !nentries)
 {
  cout << " *" << ClassName() <<"::SetWaveform* Inconsistent input for histogram treatment." << endl;
  cout << " Nbins:" << nbins << " Nentries:" << nentries << " f:" << fSample << endl; 
  return;
 }
 
 fWaveform.Set(nbins);
 fNwf=nbins;
 
 Double_t val=0;
 for (Int_t i=1; i<=nbins; i++)
 {
  val=h->GetBinContent(i);
  fWaveform[i-1]=val;;
 }
}

void NcDSP::SetWaveform(TGraph* gr,Float_t f)
{
 
 fWaveform.Set(0);
 fNwf=0;
 
 if (f>=0) fSample=f;
 
 if (!gr)
 {
  cout << " *" << ClassName() <<"::SetWaveform* No input TGraph object specified." << endl;
  return;
 }
 
 Int_t n=gr->GetN();
 
 if (!n)
 {
  cout << " *" << ClassName() <<"::SetWaveform* Inconsistent input for TGraph treatment." << endl;
  cout << " n:" << n << " f:" << fSample << endl; 
  return;
 }
 
 fWaveform.Set(n);
 fNwf=n;
 
 gr->Sort();
 
 Double_t x=0;
 Double_t y=0;
 for (Int_t i=0; i<n; i++)
 {
  gr->GetPoint(i,x,y);
  fWaveform[i]=y;
 }
}

Int_t NcDSP::GetN(Int_t mode) const
{
 
 Int_t n=fN;
 
 if (mode==1) n=fNwf;
 
 return n;
}

TArrayD NcDSP::GetData(TString sel) const
{
 
 if (sel.Contains("RE") && sel.Contains("in")) return fReIn;
 if (sel.Contains("IM") && sel.Contains("in")) return fImIn;
 if (sel.Contains("RE") && sel.Contains("out")) return fReOut;
 if (sel.Contains("IM") && sel.Contains("out")) return fImOut;
 if (sel.Contains("Hist")) return fHisto;
 if (sel.Contains("Wave")) return fWaveform;
 
 TArrayD data(fN);
 Double_t pi=acos(-1.);
 Double_t re=0;
 Double_t im=0;
 Double_t amp=0;
 Double_t phi=0;
 for (Int_t i=0; i<fN; i++)
 {
  if (sel.Contains("in"))
  {
   re=fReIn.At(i);
   im=fImIn.At(i);
  }
  if (sel.Contains("out"))
  {
   re=fReOut.At(i);
   im=fImOut.At(i);
  }
  amp=sqrt(re*re+im*im);
  phi=0;
  if (im || re) phi=atan2(im,re);
 
  if (sel.Contains("AMP")) data.SetAt(amp,i);
  if (sel.Contains("PHIR")) data.SetAt(phi,i);
  if (sel.Contains("PHID")) data.SetAt(phi*180./pi,i);
 }
 
 return data;
}

void NcDSP::Fourier(TString mode,TH1* hist,TString sel)
{
 
 fReOut.Set(0);
 fImOut.Set(0);
 fHisto.Set(0);
 
 if (fN<1) return;
 
 // Construct the Fast Fourier Transform (FFT) processor
 TString opt=mode;
 // Comply with the TVirtualFFT conventions
 if (mode=="C2C") opt="C2CFORWARD";
 if (mode=="C2CI") opt="C2CBACKWARD";
 // Set the optimalisation initialisation to "estimate" and create a new processor
 opt+=" ES K";
 
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
 
 fProc=TVirtualFFT::FFT(1,&fN,opt);
 
 // Enter the input data
 Int_t nReIn=fReIn.GetSize();
 Int_t nImIn=fImIn.GetSize();
 if (mode=="R2C")
 {
  Double_t* data=fReIn.GetArray();
  fProc->SetPoints(data);
 }
 else
 {
  for (Int_t i=0; i<fN; i++)
  {
   if (nReIn && nImIn) fProc->SetPoint(i,fReIn[i],fImIn[i]);
   if (nReIn && !nImIn) fProc->SetPoint(i,fReIn[i],0);
   if (!nReIn && nImIn) fProc->SetPoint(i,0,fImIn[i]);
  }
 }
 
 // Perform the Fast Fourier Transform
 fProc->Transform();
 
 Double_t rN=fN;
 
 // Copy the resulting data in the output arrays
 fReOut.Set(fN);
 fImOut.Set(fN);
 Double_t re=0;
 Double_t im=0;
 for (Int_t i=0; i<fN; i++)
 {
  fProc->GetPointComplex(i,re,im);
  re/=sqrt(rN);
  im/=sqrt(rN);
  fReOut[i]=re;
  fImOut[i]=im;
 }
 
 Int_t imode=1;
 if (mode=="C2R" || mode=="C2CI") imode=-1;
 TString name;
 name.Form("DFT (%-s)",mode.Data());
 if (imode<0) name.Form("Inverse DFT (%-s)",mode.Data());
 
 // Check for only domain specification specification
 TString opts[8]={"RE","IM","AMP","dB","PHIR","PHID","CPHI","SPHI"};
 Bool_t found=0;
 for (Int_t i=0; i<8; i++)
 {
  if (sel.Contains(opts[i]))
  {
   found=1;
   break;
  }  
 }
 
 if (!found)
 {
  if (sel.Contains("n") || sel.Contains("t")) // Time domain
  {
   if (mode=="C2C" || mode=="C2CI") sel+=" RE";
  }
  else // Frequency domain
  {
   sel+=" RE";
  }
 }
 
 // Fill the requested histogram
 HistogramTrafoResult(name,imode,hist,sel);
}

void NcDSP::Hartley(Int_t mode,TH1* hist,TString sel)
{
 
 fReOut.Set(0);
 fImOut.Set(0);
 fHisto.Set(0);
 
 if (!mode || fN<1) return;
 
 // Construct the Fast Fourier Transform (FFT) processor
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
 
 fProc=TVirtualFFT::FFT(1,&fN,"DHT ES K");
 
 // Enter the input data
 Double_t* data=fReIn.GetArray();
 fProc->SetPoints(data);
 
 // Perform the Discrete Hartley Transform
 fProc->Transform();
 
 Double_t rN=fN;
 
 // Copy the resulting data in the output arrays
 fReOut.Set(fN);
 fImOut.Set(0);
 Double_t re=0;
 for (Int_t i=0; i<fN; i++)
 {
  re=fProc->GetPointReal(i,kFALSE);
  re/=sqrt(rN);
  fReOut.SetAt(re,i);
 }
 
 // Fill the requested histogram
 TString name="Discrete Hartley Transform (DHT)";
 if (mode<0) name="Inverse Discrete Hartley Transform (IDHT)";
 
 HistogramTrafoResult(name,mode,hist,sel);
}

void NcDSP::Cosine(Int_t type,TH1* hist,TString sel)
{
 
 fReOut.Set(0);
 fImOut.Set(0);
 fHisto.Set(0);
 
 if (abs(type)<1 || abs(type)>4 || fN<1 || (abs(type)==1 && fN<2)) return;
 
 // Convert negative type specifications to the corresponding "normal" ones
 Int_t type2=type;
 if (type==-1 || type==-4) type2=abs(type);
 if (type==-2) type2=3;
 if (type==-3) type2=2;
 
 // Comply with the TVirtualFFT conventions
 Int_t kind=type2-1;
 
 // Construct the Fast Fourier Transform (FFT) processor
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
 
 fProc=TVirtualFFT::SineCosine(1,&fN,&kind,"ES");
 
 // Enter the input data
 Double_t* data=fReIn.GetArray();
 fProc->SetPoints(data);
 
 // Perform the Discrete Cosine Transform
 fProc->Transform();
 
 Double_t rN=fN;
 
 // Copy the resulting data in the output arrays
 fReOut.Set(fN);
 fImOut.Set(0);
 Double_t re=0;
 for (Int_t i=0; i<fN; i++)
 {
  re=fProc->GetPointReal(i,kFALSE);
  if (type2==1)
  {
   re/=sqrt(2.*(rN-1.));
  }
  else
  {
   re/=sqrt(2.*rN);
  }
  fReOut[i]=re;
 }
 
 TString name="Discrete Cosine Transform (DCT-";
 if (type<0) name="Inverse Discrete Cosine Transform (IDCT-";
 if (abs(type)==1) name+="I)";
 if (abs(type)==2) name+="II)";
 if (abs(type)==3) name+="III)";
 if (abs(type)==4) name+="IV)";
 
 // Fill the requested histogram
 HistogramTrafoResult(name,type,hist,sel);
}

void NcDSP::Sine(Int_t type,TH1* hist,TString sel)
{
 
 fReOut.Set(0);
 fImOut.Set(0);
 fHisto.Set(0);
 
 if (abs(type)<1 || abs(type)>4 || fN<1 || (abs(type)==1 && fN<2)) return;
 
 // Convert negative type specifications to the corresponding "normal" ones
 Int_t type2=type;
 if (type==-1 || type==-4) type2=abs(type);
 if (type==-2) type2=3;
 if (type==-3) type2=2;
 
 // Comply with the TVirtualFFT conventions
 Int_t kind=type2+3;
 
 // Construct the Fast Fourier Transform (FFT) processor
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
 
 fProc=TVirtualFFT::SineCosine(1,&fN,&kind,"ES K");
 
 // Enter the input data
 Double_t* data=fReIn.GetArray();
 fProc->SetPoints(data);
 
 // Perform the Discrete Cosine Transform
 fProc->Transform();
 
 Double_t rN=fN;
 
 // Copy the resulting data in the output arrays
 fReOut.Set(fN);
 fImOut.Set(0);
 Double_t re=0;
 for (Int_t i=0; i<fN; i++)
 {
  re=fProc->GetPointReal(i,kFALSE);
  if (type2==1)
  {
   re/=sqrt(2.*(rN+1.));
  }
  else
  {
   re/=sqrt(2.*rN);
  }
  fReOut[i]=re;
 }
 
 TString name="Discrete Sine Transform (DST-";
 if (type<0) name="Inverse Discrete Sine Transform (IDST-";
 if (abs(type)==1) name+="I)";
 if (abs(type)==2) name+="II)";
 if (abs(type)==3) name+="III)";
 if (abs(type)==4) name+="IV)";
 
 // Fill the requested histogram
 HistogramTrafoResult(name,type,hist,sel);
}

void NcDSP::Hilbert(Int_t mode,TH1* hist,TString sel)
{
 
 fReOut.Set(0);
 fImOut.Set(0);
 fHisto.Set(0);
 
 if (!mode || fN<1) return;
 
 // Save the original input array
 TArrayD temp=fReIn;
 
 // Perform the DFT on X[] to obtain Q[]
 Fourier("R2C");
 
 // Perform the Hilbert transformation in the frequency domain
 Double_t a=0;
 Double_t b=0;
 Double_t re=0;
 Double_t im=0;
 for (Int_t k=0; k<fN; k++)
 {
  a=fReOut[k];
  b=fImOut[k];
  re=0;
  im=0;
  if (k && k<=fN/2)
  {
   re=b;
   im=-a;
  }
  if (k>fN/2)
  {
   re=-b;
   im=a;
  }
  fReOut[k]=re;
  fImOut[k]=im;
 }
 
 // Perform the inverse DFT on H[k] to obtain HT[n]
 TArrayD ReZ=fReIn;  // The real part of the analytic signal
 LoadResult();
 Fourier("C2R");
 TArrayD ImZ=fReOut; // The imaginary part of the analytic signal
 
 // Provide the analytic signal Z[n]=X[n]+i*HT[n] as intermediate output result
 if (mode>0) // Forward transformation
 {
  fReOut=ReZ;
  fImOut=ImZ;
 }
 else // Backward transformation
 {
  fImOut=ReZ;
  for (Int_t n=0; n<fN; n++)
  {
   fReOut[n]=-ImZ[n];
  } 
 }
 
 // Restore the original input array
 fReIn=temp;   
 
 // Check for only "n" or "t" specification
 TString opts[8]={"X","HT","AMP","dB","PHIR","PHID","CPHI","OMEGA"};
 Bool_t found=0;
 for (Int_t i=0; i<8; i++)
 {
  if (sel.Contains(opts[i]))
  {
   found=1;
   break;
  }  
 }
 
 if (!found)
 {
  if (mode>0) // Forward transformation
  {
   sel+=" HT";
  }
  else // Inverse transformation
  {
   sel+=" X";
  }
 }
 
 // Convert the "X"and "HT" selections to the corresponding "RE" and "IM'.
 sel.ReplaceAll("X","RE");
 sel.ReplaceAll("HT","IM");
 
 // Fill the requested histogram
 TString name="Discrete Hilbert Transform (DHIT)";
 if (mode<0) name="Inverse Discrete Hilbert Transform (IDHIT)";
 
 HistogramTrafoResult(name,mode,hist,sel);
 
 // Store the correct data array as real output
 if (mode>0) fReOut=ImZ; // Forward transformation
}

void NcDSP::HistogramTrafoResult(TString name,Int_t mode,TH1* hist,TString sel)
{
 
 fHisto.Set(0);
 
 if (!mode || fN<1 || !hist) return;
 
 Int_t n=1+fN/2;
 Float_t maxfrac=0.5;
 if (sel.Contains("n") || sel.Contains("t") || sel.Contains("2"))
 {
  n=fN;
  maxfrac=1;
 }
 
 if ((sel.Contains("Hz") || sel.Contains("t") || sel.Contains("OMEGA")) && fSample<=0) return;
 
 if (sel.Contains("OMEGA") && !(sel.Contains("n")) && !(sel.Contains("t"))) return; 
 
 hist->Reset();
 
 // Initialize the histogram title
 TString title=name;
 title+=" ";
 
 TString title2="";
 
 // Define and fill the requested result histogram
 if (sel.Contains("k"))
 {
  hist->SetBins(n,0,n-1);
  title+="index frequency domain";
  if (mode<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Index k";
  if (sel.Contains("RE")) title+=";Re(Q[k])";
  if (sel.Contains("IM")) title+=";Im(Q[k])";
  if (sel.Contains("AMP")) title+=";Amplitude |Q[k]|";
  if (sel.Contains("dB")) title+=";Amplitude |Q[k]| in dB";
  if (sel.Contains("PHIR")) title+=";Phase #varphi[k] (rad)";
  if (sel.Contains("PHID")) title+=";Phase #varphi[k] (deg)";
  if (sel.Contains("CPHI")) title+=";cos(#varphi[k])";
  if (sel.Contains("SPHI")) title+=";sin(#varphi[k])";
  if (!(title.Contains("[k]")))
  {
   sel+=" RE";
   title+=";Value Q[k]";
  }
 }
 if (sel.Contains("f"))
 {
  hist->SetBins(n,0,maxfrac);
  title+="fractional frequency domain";
  if (mode<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Fraction f of sampling rate";
  if (sel.Contains("RE")) title+=";Re(Q[f])";
  if (sel.Contains("IM")) title+=";Im(Q[f])";
  if (sel.Contains("AMP")) title+=";Amplitude |Q[f]|";
  if (sel.Contains("dB")) title+=";Amplitude |Q[f]| in dB";
  if (sel.Contains("PHIR")) title+=";Phase #varphi[f] (rad)";
  if (sel.Contains("PHID")) title+=";Phase #varphi[f] (deg)";
  if (sel.Contains("CPHI")) title+=";cos(#varphi[f])";
  if (sel.Contains("SPHI")) title+=";sin(#varphi[f])";
  if (!(title.Contains("[f]")))
  {
   sel+=" RE";
   title+=";Value Q[f]";
  }
 }
 if (sel.Contains("Hz"))
 {
  hist->SetBins(n,0,maxfrac*fSample);
  title+="actual frequency domain";
  if (mode<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Frequency #nu (Hz)";
  if (sel.Contains("RE")) title+=";Re(Q[#nu])";
  if (sel.Contains("IM")) title+=";Im(Q[#nu])";
  if (sel.Contains("AMP")) title+=";Amplitude |Q[#nu]|";
  if (sel.Contains("dB")) title+=";Amplitude |Q[#nu]| in dB";
  if (sel.Contains("PHIR")) title+=";Phase #varphi[#nu] (rad)";
  if (sel.Contains("PHID")) title+=";Phase #varphi[#nu] (deg)";
  if (sel.Contains("CPHI")) title+=";cos(#varphi[#nu])";
  if (sel.Contains("SPHI")) title+=";sin(#varphi[#nu])";
  if (!(title.Contains("[#nu]")))
  {
   sel+=" RE";
   title+=";Value Q[#nu]";
  }
 }
 
 if (sel.Contains("n"))
 {
  hist->SetBins(fN,0,fN);
  title+="sampling time domain";
  if (mode>0 && !(name.Contains("Hilbert"))) title+=" input";
  if (fSample>0) title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Sample number n";
  if (name.Contains("Hilbert"))
  {
   if (sel.Contains("RE")) title+=";Value X[n]";
   if (sel.Contains("IM")) title+=";Value HT[n]";
   if (sel.Contains("AMP")) title+=";Amplitude |Z[n]|";
   if (sel.Contains("dB")) title+=";Amplitude |Z[n]| in dB";
  }
  else
  {
   if (sel.Contains("RE")) title+=";Re(X[n])";
   if (sel.Contains("IM")) title+=";Im(X[n])";
   if (sel.Contains("AMP")) title+=";Amplitude |X[n]|";
   if (sel.Contains("dB")) title+=";Amplitude |X[n]| in dB";
  }
  if (sel.Contains("PHIR")) title+=";Phase #varphi[n] (rad)";
  if (sel.Contains("PHID")) title+=";Phase #varphi[n] (deg)";
  if (sel.Contains("CPHI")) title+=";cos(#varphi[n])";
  if (sel.Contains("SPHI")) title+=";sin(#varphi[n])";
  if (sel.Contains("OMEGA")) title+=";#omega[n] (rad/s)";
  if (!(title.Contains("[n]")))
  {
   sel+=" RE";
   title+=";Value X[n]";
  }
 }
 if (sel.Contains("t"))
 {
  hist->SetBins(fN,0,double(fN)/fSample);
  title+="actual time domain";
  if (mode>0 && !(name.Contains("Hilbert"))) title+=" input";
  if (fSample>0) title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Time t (seconds)";
  if (name.Contains("Hilbert"))
  {
   if (sel.Contains("RE")) title+=";Value X[t]";
   if (sel.Contains("IM")) title+=";Value HT[t]";
   if (sel.Contains("AMP")) title+=";Amplitude |Z[t]|";
   if (sel.Contains("dB")) title+=";Amplitude |Z[t]| in dB";
  }
  else
  {
   if (sel.Contains("RE")) title+=";Re(X[t])";
   if (sel.Contains("IM")) title+=";Im(X[t])";
   if (sel.Contains("AMP")) title+=";Amplitude |X[t]|";
   if (sel.Contains("dB")) title+=";Amplitude |X[t]| in dB";
  }
  if (sel.Contains("PHIR")) title+=";Phase #varphi[t] (rad)";
  if (sel.Contains("PHID")) title+=";Phase #varphi[t] (deg)";
  if (sel.Contains("CPHI")) title+=";cos(#varphi[t])";
  if (sel.Contains("SPHI")) title+=";sin(#varphi[t])";
  if (sel.Contains("OMEGA")) title+=";#omega[t] (rad/s)";
  if (!(title.Contains("[t]")))
  {
   sel+=" RE";
   title+=";Value X[t]";
  }
 }
 
 hist->SetTitle(title);
 
 Int_t nReIn=fReIn.GetSize();
 Int_t nImIn=fImIn.GetSize();
 Int_t nReOut=fReOut.GetSize();
 Int_t nImOut=fImOut.GetSize();
 Double_t pi=acos(-1.);
 Double_t amp=0;
 Double_t phi=0;
 Double_t cphi=0;
 Double_t sphi=0;
 Double_t omega=0;
 Double_t re=0;
 Double_t im=0;
 Double_t re2=0;
 Double_t im2=0;
 Double_t phi2=0;
 fHisto.Set(n);
 
 // Determine stepsize in fractional sampling frequency for the Cosine and Sine transformation
 Double_t fstep=1./(2.*double(fN));
 if (name.Contains("DCT-I") && fN>1) fstep=1./(2.*double(fN-1));
 if (name.Contains("DST-I")) fstep=1./(2.*double(fN+1));
 
 Double_t x=-1;
 for (Int_t i=0; i<n; i++)
 {
  if (name.Contains("DCT-"))
  {
   x=i;
   if (mode==3 || mode==-2 || abs(mode)==4) x+=0.5;
  }
  if (name.Contains("DST-"))
  {
   x=i+1;
   if (mode==3 || mode==-2 || abs(mode)==4) x-=0.5;
  }
  x*=fstep;
  if (sel.Contains("Hz")) x*=fSample;
 
  if (sel.Contains("n") || sel.Contains("t")) // Time domain data requested
  {
   if (mode>0) // Forward transformation
   {
    if (name.Contains("Hilbert"))
    {
     if (nReOut)
     {
      re=fReOut[i];
      if (i<n-1) re2=fReOut[i+1];
     }
     if (nImOut)
     {
      im=fImOut[i];
      if (i<n-1) im2=fImOut[i+1];
     }
    }
    else
    {
     if (nReIn) re=fReIn[i];
     if (nImIn) im=fImIn[i];
    }
   }
   if (mode<0) // Inverse transformation
   {
    if (nReOut) re=fReOut[i];
    if (nImOut) im=fImOut[i];
    if (name.Contains("Hilbert") && i<n-1)
    {
     if (nReOut) re2=fReOut[i+1];
     if (nImOut) im2=fImOut[i+1];
    }
   }
  }
  else // Frequency domain data requested
  {
   if (mode<0) // Inverse transformation
   {
    if (nReIn) re=fReIn[i];
    if (nImIn) im=fImIn[i];
   }
   else // Forward transformation
   {
    if (nReOut) re=fReOut[i];
    if (nImOut) im=fImOut[i];
   }
  }
  if (name.Contains("DCT-") || name.Contains("DST-"))
  {
   if (sel.Contains("f") || sel.Contains("Hz"))
   {
    hist->Fill(x,re);
   }
   else
   {
    hist->SetBinContent(i+1,re);
   }
  }
  else
  {
   amp=sqrt(re*re+im*im);
   phi=0;
   if (im || re) phi=atan2(im,re);
   cphi=cos(phi);
   sphi=sin(phi);
   phi2=0;
   if (fSample>0 && i<n-1 && (im2 || re2))
   {
    phi2=atan2(im2,re2);
    if ((phi2*phi)>=0) omega=fabs(phi2-phi)*fSample;
   }
   if (sel.Contains("RE")) hist->SetBinContent(i+1,re);
   if (sel.Contains("IM")) hist->SetBinContent(i+1,im);
   if (sel.Contains("AMP")) hist->SetBinContent(i+1,amp);
   if (amp<=0) amp=hist->GetMinimum();
   if (sel.Contains("dB"))
   {
    // Check for rounding errors
    if (amp<=0)
    {
     amp=hist->GetMinimum();
     hist->SetBinContent(i+1,amp);
    }
    else
    {
     hist->SetBinContent(i+1,20.*log10(amp));
    }
   }
   if (sel.Contains("PHIR")) hist->SetBinContent(i+1,phi);
   if (sel.Contains("PHID")) hist->SetBinContent(i+1,phi*180./pi);
   if (sel.Contains("CPHI")) hist->SetBinContent(i+1,cphi);
   if (sel.Contains("SPHI")) hist->SetBinContent(i+1,sphi);
   if (sel.Contains("OMEGA")) hist->SetBinContent(i+1,omega);
  }
 }
 
 // Fill the fHisto array with the histogram contents
 for (Int_t i=0; i<n; i++)
 {
  fHisto[i]=hist->GetBinContent(i+1);
 }
}

TArrayD NcDSP::Convolve(TH1* hist,Int_t* i1,Int_t* i2,Int_t shift)
{
 
 TArrayD y(0);
 if (hist) hist->Reset();
 
 TArrayD x=fReIn;
 TArrayD h=fWaveform;
 Int_t nx=x.GetSize();
 Int_t nh=h.GetSize();
 
 if (nh<1 || nx<1)
 {
  printf(" *%-s::Convolve* Input or Waveform data are missing. Ninput=%-i Nwaveform=%-i \n",ClassName(),nx,nh);
  return y;
 }
 
 if (shift<0 || shift>2)
 {
  printf(" *%-s::Convolve* Invalid input argument shift=%-i \n",ClassName(),shift);
  return y;
 }
 
 Int_t ny=nx+nh-1;
 y.Set(ny);
 
 // The normalization data for correlation studies
 TArrayD x1(ny);
 TArrayD x2(ny);
 TArrayD h1(ny);
 TArrayD h2(ny);
 TArrayI ns(ny);
 
 // Convolution from the input signal viewpoint
 for (Int_t ix=0; ix<nx; ix++)
 {
  for (Int_t ih=0; ih<nh; ih++)
  {
   // Perform the unnormalized convolution
   y[ix+ih]=y[ix+ih]+x[ix]*h[ih];
 
   if (fNorm=="NONE" || fNorm=="GNCC") continue;
 
   // Gather the normalization data for the "NCC" and "ZNCC" cross-correlation studies
   x1[ix+ih]+=x[ix];
   x2[ix+ih]+=pow(x[ix],2);
   h1[ix+ih]+=h[ih];
   h2[ix+ih]+=pow(h[ih],2);
   ns[ix+ih]+=1;
  }
 }
 
 // Normalize values for cross-correlation studies if requested
 if (fNorm!="NONE")
 {
  Double_t val=0;
  if (fNorm=="GNCC") // Global vector length normalization
  {
   Double_t xnorm=0;
   Double_t hnorm=0;
   for (Int_t i=0; i<nx; i++)
   {
    xnorm+=pow(x[i],2);
   }
   for (Int_t i=0; i<nh; i++)
   {
    hnorm+=pow(h[i],2);
   }
   xnorm=sqrt(xnorm);
   hnorm=sqrt(hnorm);
   for (Int_t i=0; i<ny; i++)
   {
    val=xnorm*hnorm;
    if (val>0) y[i]=y[i]/val;
   }
  }
  else if (fNorm=="NCC") // NCC normalization
  {
   for (Int_t i=0; i<ny; i++)
   {
    val=sqrt(x2[i]*h2[i]);
    if (val>0) y[i]=y[i]/val;
   }
  }
  else // ZNCC normalization
  {
   y.Reset();
   Double_t x1mean=0;
   Double_t h1mean=0;
   Double_t x2mean=0;
   Double_t h2mean=0;
   Double_t rn=0;
   Double_t sx=0;
   Double_t sh=0;
   Double_t sxsh=0;
   Int_t iy=0;
   for (Int_t ix=0; ix<nx; ix++)
   {
    for (Int_t ih=0; ih<nh; ih++)
    { 
     iy=ix+ih;
 
     if (!ns[iy]) continue;
 
     rn=ns[iy];
     x1mean=x1[iy]/rn;
     x2mean=x2[iy]/rn;
     h1mean=h1[iy]/rn;
     h2mean=h2[iy]/rn;
     sx=sqrt(x2mean-pow(x1mean,2));
     sh=sqrt(h2mean-pow(h1mean,2));
     sxsh=sx*sh;
     val=0;
     if (sxsh>0) val=(x[ix]-x1mean)*(h[ih]-h1mean)/(rn*sxsh);
     y[iy]+=val;
    }
   }
  }
 }
 
 Int_t ilow=nh-1;
 Int_t iup=ny-nh;
 
 Int_t ioffset=0;
 if (shift==2)
 {
  ioffset=nh/2;
  TArrayD temp(nx);
  for (Int_t i=0; i<nx; i++)
  {
   temp[i]=y[i+ioffset];
  }
  y=temp;
  ny=nx;
  ilow=ilow-ioffset;
  iup=iup-ioffset;
 }
 
 if (i1) *i1=ilow;
 if (i2) *i2=iup;
 
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   title.Form("NcDSP Convolution result (%-g samples/sec);Time in seconds;Value",fSample);
   if (shift==0 || shift==2)
   {
    hist->SetBins(ny,0,double(ny)/fSample);
   }
   else
   {
    hist->SetBins(ny,double(-nh+1)/fSample,double(nx)/fSample);
   }
  }
  else
  {
   title.Form("NcDSP Convolution result;Sample number;Value");
   if (shift==0 || shift==2)
   {
    hist->SetBins(ny,0,ny);
   }
   else
   {
    hist->SetBins(ny,-nh+1,nx);
   }
  }
 
  hist->SetTitle(title);
  hist->SetMarkerStyle(20);
  for (Int_t ibin=1; ibin<=ny; ibin++)
  {
   hist->SetBinContent(ibin,y[ibin-1]);
  }
 
  Double_t ymin=hist->GetMinimum();
  Double_t ymax=hist->GetMaximum();
 
  Double_t xlow=0;
  Double_t xup=0;
  if (i1)  xlow=hist->GetBinLowEdge(ilow+1);
  if (i2)
  {
   xup=hist->GetBinLowEdge(iup+1);
   xup+=hist->GetBinWidth(1);
  }
 
  TLine* vline1=0;
  TLine* vline2=0;
 
  if (i1)
  {
   vline1=new TLine(xlow,ymin,xlow,ymax);
   vline1->SetLineStyle(2); // Dashed line
   vline1->SetLineWidth(2);
   vline1->SetLineColor(4); // Blue color
  }
  if (i2)
  {
   vline2=new TLine(xup,ymin,xup,ymax);
   vline2->SetLineStyle(2); // Dashed line
   vline2->SetLineWidth(2);
   vline2->SetLineColor(4); // Blue color
  }
 
   TList* hlist=hist->GetListOfFunctions();
   if (vline1) hlist->Add(vline1);
   if (vline2) hlist->Add(vline2);
 }
 
 return y;
}

TArrayD NcDSP::Correlate(TH1* hist,Int_t* i1,Int_t* i2,Double_t* peak,TString norm)
{
 
 TArrayD y(0);
 if (hist) hist->Reset();
 Int_t nx=fReIn.GetSize();
 Int_t nh=fWaveform.GetSize();
 
 if (nh<1 || nx<1)
 {
  printf(" *%-s::Correlate* Input or Waveform data are missing. Ninput=%-i Nwaveform=%-i \n",ClassName(),nx,nh);
  return y;
 }
 
 norm.ToUpper(); // Convert norm to uppercase characters
 if (norm!="NONE" && norm!="GNCC" && norm!="NCC" && norm!="ZNCC")
 {
  printf(" *%-s::Correlate* Unsupported normalization mode : norm=%-s \n",ClassName(),norm.Data());
  return y;
 }
 
 // The temporary "flipped" waveform
 TArrayD store=fWaveform;
 TArrayD temp(nh);
 for (Int_t i=1; i<=nh; i++)
 {
  temp[i-1]=fWaveform[nh-i];
 }
 
 fWaveform=temp;
 fNorm=norm; // Normalize the correlation values if requested
 y=Convolve(hist,i1,i2,1);
 fNorm="NONE"; // Reset the normalization indicator for convolution studies
 
 // Get the index of the maximum in the y[] array
 Int_t imax=0;
 Double_t ymax=0;
 Int_t ny=y.GetSize();
 if (ny) ymax=y[0];
 for (Int_t i=1; i<ny; i++)
 {
  if (y[i]>ymax)
  {
   ymax=y[i];
   imax=i;
  }
 }
 
 // Convert to the index matching the x[] convention
 Int_t idx=imax-nh+1;
 Double_t xpeak=idx;
 if (fSample>0) xpeak=xpeak/fSample;
 
 if (peak) *peak=xpeak;
 
 // Put the correct histogram title
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   if (norm=="NCC")
   {
    title.Form("%-s Normalized Cross-Correlation (NCC): max. at lag=%-g sec. (%-g samples/sec)",ClassName(),xpeak,fSample);
   }
   else if (norm=="ZNCC")
   {
    title.Form("%-s Zero-Normalized Cross-Correlation (ZNCC): max. at lag=%-g sec. (%-g samples/sec)",ClassName(),xpeak,fSample);
   }
   else if (norm=="GNCC")
   {
    title.Form("%-s Globally Normalized Cross-Correlation (GNCC): max. at lag=%-g sec. (%-g samples/sec)",ClassName(),xpeak,fSample);
   }
   else
   {
    title.Form("%-s Unnormalized Cross-Correlation: max. at lag=%-g sec. (%-g samples/sec)",ClassName(),xpeak,fSample);
   }
   title+=";Time lag in seconds;Cross-Correlation value";
  }
  else
  {
   if (norm=="NCC")
   {
    title.Form("%-s Normalized Cross-Correlation (NCC): max. at lag=%-i samplings",ClassName(),idx);
   }
   else if (norm=="ZNCC")
   {
    title.Form("%-s Zero-Normalized Cross-Correlation (ZNCC): max. at lag=%-i samplings",ClassName(),idx);
   }
   else if (norm=="GNCC")
   {
    title.Form("%-s Gobally Normalized Cross-Correlation (GNCC): max. at lag=%-i samplings",ClassName(),idx);
   }
   else
   {
    title.Form("%-s Unnormalized Cross-Correlation: max. at lag=%-i samplings",ClassName(),idx);
   }
   title+=";Lag in samplings;Cross-Correlation value";
  }
  hist->SetTitle(title);
 }
 
 // Restore the original waveform
 fWaveform=store;
 
 return y;
}

TArrayD NcDSP::Digitize(Int_t nbits,Double_t vcal,Int_t mode,TH1* hist,Double_t* stp,Double_t* scale) const
{
 
 TArrayD arrdig(0);
 if (hist) hist->Reset();
 
 if (mode<0 || mode>3)
 {
  cout << " *" << ClassName() << "::Digitize* Inconsistent input mode=" << mode << endl;
  return arrdig;
 }
 
 if (!nbits || nbits>60 || nbits<-60 || fabs(vcal)<=0)
 {
  cout << " *" << ClassName() << "::Digitize* Inconsistent input nbits=" << nbits << " vcal=" << vcal << endl;
  return arrdig;
 }
 
 if ((mode==1 || mode==3) && nbits==1)
 {
  cout << " *" << ClassName() << "::Digitize* Inconsistent input nbits=" << nbits << " mode=" << mode << endl;
  return arrdig;
 }
 
 Bool_t logmode=kFALSE;
 if (nbits<0) logmode=kTRUE;
 
 nbits=abs(nbits);
 
 Long_t nlevels=pow(2,nbits);
 Double_t range=fabs(vcal);
 if (mode==1 || mode==3) range*=2.;
 Double_t step=fabs(vcal);
 if (mode==0) step=range/double(nlevels-1);
 if (mode==1) step=range/double(nlevels-2);
 
 if (step<=0) return arrdig;
 
 Long_t nstepsmin=0;
 Long_t nstepsmax=nlevels-1;
 if ((mode==0 || mode==2) && vcal<0)
 {
  nstepsmin=-nlevels+1;
  nstepsmax=0;
 }
 if (mode==1 || mode==3)
 {
  nstepsmin=-nlevels/2;
  nstepsmax=nlevels/2-1;
 }
 
 Double_t digvalmin=double(nstepsmin)*step;
 Double_t digvalmax=double(nstepsmax)*step;
 
 cout << " *" << ClassName() << "::Digitize* Parameters of the " << nbits << "-bits";
 if (logmode) cout << " Log10";
 cout << " ADC digitization." << endl;
 TString s="linear";
 if (logmode) s="Log10";
 if (mode==0 || mode==2)
 {
  cout << " Digitized " << s << " full scale range : [" << digvalmin << "," << digvalmax << "]  Step size : " << step << endl;
 }
 if (mode==1 || mode==3)
 {
  cout << " Digitized " << s << " full scale range : [" << (digvalmin+step) << "," << digvalmax << "]  Step size : " << step;
  cout << "  Actual range : [" << digvalmin << "," << digvalmax << "]" << endl;
 }
 
 if (logmode)
 {
  Double_t linvalmin=pow(10,digvalmin);
  Double_t linvalmax=pow(10,digvalmax);
  if (mode==0 || mode==2)
  {
   cout << " Digitized linear full scale range : [" << linvalmin << "," << linvalmax << "]" << endl;
  }
  if (mode==1 || mode==3)
  {
   cout << " Digitized linear full scale range : [" << (linvalmin*pow(10,step)) << "," << linvalmax << "]";
   cout << "  Actual range : [" << linvalmin << "," << linvalmax << "]" << endl;
  }
 }
 
 if (stp) *stp=step;
 if (scale)
 {
  *scale=fabs(vcal);
  if (mode==2)
  {
   if (vcal<0) *scale=digvalmin;
   if (vcal>0) *scale=digvalmax;
  }
  if (mode==3) *scale=digvalmax;
 }
 
 if (fNwf<1)
 {
  printf(" === No waveform data present: Only listing of ADC specs.===\n");
  return arrdig;
 }
 
 arrdig.Set(fNwf);
 
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   hist->SetBins(fNwf,0,double(fNwf)/fSample);
   title.Form("NcDSP Digitize result (%-g samples/sec);Time in seconds;Value",fSample);
  }
  else
  {
   title.Form("NcDSP Digitize result;Sample number;Value");
   hist->SetBins(fNwf,0,fNwf);
  }
  hist->SetMarkerStyle(20);
  hist->SetTitle(title);
 }
 
 Double_t val=0;
 Long_t nsteps=0;
 Double_t digval=0;
 for (Int_t j=0; j<fNwf; j++)
 {
  val=fWaveform[j];
  
  if (logmode)
  {
   if (val>0)
   {
    val=log10(val);
   }
   else
   {
    cout << endl;
    cout << " *" << ClassName() << "::Digitize* Error: Non-positive input value encountered for Log10 ADC." << endl;
    arrdig.Set(0);
    return arrdig;
   }
  }
  nsteps=val/step;
 
  if (nsteps<nstepsmin) nsteps=nstepsmin;
  if (nsteps>nstepsmax) nsteps=nstepsmax;
 
  digval=double(nsteps)*step;
 
  if (logmode) digval=pow(10,digval);
 
  arrdig[j]=digval;
 
  if (hist) hist->SetBinContent(j+1,digval);
 }
 
 return arrdig;
}

TArrayL64 NcDSP::ADC(Int_t nbits,Double_t range,Double_t Vbias,TArray* Vsig,TH1* hist,Int_t B,Int_t C) const
{
 
 TArrayL64 arradc(0);
 if (hist) hist->Reset();
 
 if (!Vsig) Vsig=(TArray*)&fWaveform;
 Int_t nVsig=Vsig->GetSize();
 
 if (nbits<=0 || nbits>60 || !range || fabs(Vbias)>fabs(range) || B<0 || (B && C<1))
 {
  printf("\n *%-s::ADC* Inconsistent input nbits=%-i range=%-g Vbias=%-g B=%-i C=%-i \n",ClassName(),nbits,range,Vbias,B,C);
  return arradc;
 }
 
 NcMath math;
 
 Long64_t N=pow(2,nbits); // The number of quantization levels
 Long64_t adcmin=0;
 Long64_t adcmax=N-1;
 Double_t Vref=0;
 Double_t Vfs=0;
 Double_t LSB=0;
 
 // Floating point version of some parameters
 Double_t rN=N;
 Double_t rB=B;
 if (B==1) rB=exp(1);
 Double_t rC=C;
 Double_t radcmax=adcmax;
 
 if (range<0) // |range| represents Vref
 {
  Vref=fabs(range);
  LSB=Vref/rN;
  Vfs=Vref-LSB;
 }
 else // range represents Vfs
 {
  Vfs=range;
  LSB=Vfs/radcmax;
  Vref=Vfs+LSB;
 }
 
 Long64_t ped=0;
 Double_t frac=0;
 if (B) // Log ADC
 {
  if (range<0) Vfs=pow(rB,-rC)*pow(rB,radcmax*rC/rN)*Vref;
  if (range>0) Vref=pow(rB,rC)*pow(rB,-radcmax*rC/rN)*Vfs;
  LSB=Vref*pow(rB,-rC)*(pow(rB,rC/rN)-1.);
  frac=Vbias/Vref;
  ped=0;
  if (frac>0) ped=Long64_t(rN*(math.Log(rB,frac)+rC)/rC);
 }
 
 if (LSB<=0 || Vfs<=0)
 {
  printf("\n *%-s::ADC* Inconsistent parameters : LSB=%-g Vfs=%-g \n",ClassName(),LSB,Vfs);
  return arradc;
 }
 
 if (!B) ped=Vbias/LSB; // Pedestal value for a linear ADC
 
 Double_t DR=20.*log10(Vfs/LSB);
 
 TString mode="linear";
 if (B==1) mode="Log_e";
 if (B>1)
 {
  mode="Log_";
  mode+=B;
 }
 
 if (!nVsig)
 {
  printf(" *%-s::ADC* No input data have been provided --> Only the value of adc(Vbias) is returned. \n",ClassName());
  if (!B)
  {
   printf(" Parameters for a %-i-bits %-s ADC with adc=[%-lli,%-lli] : \n",nbits,mode.Data(),adcmin,adcmax);
   printf(" Vref=%-g Vfs=%-g LSB=%-g DR=%-g (dB) Vbias=%-g adc(Vbias)=%-lli \n",Vref,Vfs,LSB,DR,Vbias,ped);
  }
  else
  {
   printf(" Parameters for a %-i-bits %-s ADC with adc=[%-lli,%-lli] and a code efficiency factor of %-i: \n",nbits,mode.Data(),adcmin,adcmax,C);
   printf(" Vref=%-g Vfs=%-g LSB=%-g DR=%-g (dB) Vbias=%-g adc(Vbias)=%-lli \n",Vref,Vfs,LSB,DR,Vbias,ped);
  }
  arradc.Set(1);
  arradc[0]=ped;
  return arradc;
 }
 
 arradc.Set(nVsig);
 
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   hist->SetBins(nVsig,0,double(nVsig)/fSample);
   title.Form("%-s %-i-bits %-s ADC with Vfs=%-g LSB=%-g (%-g samples/sec);Time in seconds;ADC counts",ClassName(),nbits,mode.Data(),Vfs,LSB,fSample);
  }
  else
  {
   title.Form("%-s %-i-bits %-s ADC with Vfs=%-g LSB=%-g;Sample number;ADC counts",ClassName(),nbits,mode.Data(),Vfs,LSB);
   hist->SetBins(nVsig,0,nVsig);
  }
  hist->SetMarkerStyle(20);
  hist->SetTitle(title);
 }
 
 Double_t val=0;
 Double_t radcval=0;
 Long64_t adcval=0;
 for (Int_t j=0; j<nVsig; j++)
 {
  val=Vbias+Vsig->GetAt(j);
  
  if (B) // Log ADC
  {
   radcval=0;
   frac=val/Vref;
   if (frac>0) radcval=(rN/rC)*(math.Log(rB,frac)+rC);
  }
  else // Linear ADC
  {
   radcval=val/LSB;
  }
 
  adcval=Long64_t(radcval);
 
  // Check for saturation
  if (adcval<adcmin) adcval=adcmin;
  if (adcval>adcmax) adcval=adcmax;
 
  arradc[j]=adcval;
 
  if (hist) hist->SetBinContent(j+1,adcval);
 }
 
 return arradc;
}

TArrayD NcDSP::DAC(Int_t nbits,Double_t range,Double_t Vbias,TArray* adcs,TArray* peds,TH1* hist,Int_t B,Int_t C) const
{
 
 TArrayD arrdac(0);
 if (hist) hist->Reset();
 
 if (!adcs) adcs=(TArray*)&fWaveform;
 Int_t nadcs=adcs->GetSize();
 
 Int_t npeds=0;
 if (peds) npeds=peds->GetSize();
 
 if (nbits<=0 || nbits>60 || !range || fabs(Vbias)>fabs(range) || (peds && npeds<nadcs) || B<0 || (B && C<1))
 {
  printf("\n *%-s::DAC* Inconsistent input nbits=%-i range=%-g Vbias=%-g nadcs=%-i npeds=%-i B=%-i C=%-i \n",ClassName(),nbits,range,Vbias,nadcs,npeds,B,C);
  return arrdac;
 }
 
 NcMath math;
 
 nbits=abs(nbits);
 
 Long64_t N=pow(2,nbits); // The number of quantization levels
 Long64_t adcmin=0;
 Long64_t adcmax=N-1;
 Double_t Vref=0;
 Double_t Vfs=0;
 Double_t LSB=0;
 
 // Floating point version of some parameters
 Double_t rN=N;
 Double_t rB=B;
 if (B==1) rB=exp(1);
 Double_t rC=C;
 Double_t radcmax=adcmax;
 
 if (range<0) // |range| represents Vref
 {
  Vref=fabs(range);
  LSB=Vref/rN;
  Vfs=Vref-LSB;
 }
 else // range represents Vfs
 {
  Vfs=range;
  LSB=Vfs/radcmax;
  Vref=Vfs+LSB;
 }
 
 Long64_t ped=0;
 Double_t frac=0;
 if (B) // Log DAC
 {
  if (range<0) Vfs=pow(rB,-rC)*pow(rB,radcmax*rC/rN)*Vref;
  if (range>0) Vref=pow(rB,rC)*pow(rB,-radcmax*rC/rN)*Vfs;
  LSB=Vref*pow(rB,-rC)*(pow(rB,rC/rN)-1.);
  frac=Vbias/Vref;
  ped=0;
  if (frac>0) ped=Long64_t(rN*(math.Log(rB,frac)+rC)/rC);
 }
 
 if (LSB<=0 || Vfs<=0)
 {
  printf("\n *%-s::DAC* Inconsistent parameters : LSB=%-g Vfs=%-g \n",ClassName(),LSB,Vfs);
  return arrdac;
 }
 
 if (!B) ped=Vbias/LSB; // Pedestal value for a linear DAC
 
 Double_t DR=20.*log10(Vfs/LSB);
 
 TString mode="linear";
 if (B==1) mode="Log_e";
 if (B>1)
 {
  mode="Log_";
  mode+=B;
 }
 
 if (!nadcs)
 {
  printf(" *%-s::DAC* No input data have been provided --> Only the value of adc(Vbias) is returned. \n",ClassName());
  if (!B)
  {
   printf(" Parameters for a %-i-bits %-s DAC with adc=[%-lli,%-lli] : \n",nbits,mode.Data(),adcmin,adcmax);
   printf(" Vref=%-g Vfs=%-g LSB=%-g DR=%-g (dB) Vbias=%-g adc(Vbias)=%-lli \n",Vref,Vfs,LSB,DR,Vbias,ped);
  }
  else
  {
   printf(" Parameters for a %-i-bits %-s DAC with adc=[%-lli,%-lli] and a code efficiency factor of %-i: \n",nbits,mode.Data(),adcmin,adcmax,C);
   printf(" Vref=%-g Vfs=%-g LSB=%-g DR=%-g (dB) Vbias=%-g adc(Vbias)=%-lli \n",Vref,Vfs,LSB,DR,Vbias,ped);
  }
  arrdac.Set(1);
  arrdac[0]=ped;
  return arrdac;
 }
 
 arrdac.Set(nadcs);
 
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   hist->SetBins(nadcs,0,double(nadcs)/fSample);
   title.Form("%-s %-i-bits %-s DAC with Vfs=%-g LSB=%-g (%-g samples/sec);Time in seconds;Amplitude",ClassName(),nbits,mode.Data(),Vfs,LSB,fSample);
  }
  else
  {
   title.Form("%-s %-i-bits %-s DAC with Vfs=%-g LSB=%-g;Sample number;Amplitude",ClassName(),nbits,mode.Data(),Vfs,LSB);
   hist->SetBins(nadcs,0,nadcs);
  }
  hist->SetMarkerStyle(20);
  hist->SetTitle(title);
 }
 
 Long64_t adc=0;
 Double_t radc=0;
 Double_t dacval=0;
 for (Int_t j=0; j<nadcs; j++)
 {
  adc=adcs->GetAt(j);
  radc=adc;
 
  if (B) // Log DAC
  {
   dacval=Vref*pow(rB,-rC)*pow(rB,radc*rC/rN);
  }
  else // Linear DAC
  {
   if (peds)
   {
    ped=peds->GetAt(j);
    adc=adc-ped;
    radc=adc;
   }
   dacval=radc*LSB;
  }
 
  // Correct for bias voltage
  if (B || (!B && !peds))
  {
   dacval=dacval-Vbias;
  }
 
  arrdac[j]=dacval;
 
  if (hist) hist->SetBinContent(j+1,dacval);
 }
 
 return arrdac;
}

TArrayD NcDSP::Transmit(Int_t nbits,Double_t range,Double_t Vbias,TArray* Vsig,TArray* peds,TH1* hist,Int_t B,Int_t C) const
{
 
 TArrayL64 adcarr;
 TArrayD dacarr;
 if (hist) hist->Reset();
 
 // Provide the ADC and DAC specs in case no input is provided
 if (!Vsig && fNwf<1)
 {
  printf(" *%-s::Transmit* Specifications for the ADC-DAC transmission chain. \n",ClassName());
  adcarr=ADC(nbits,range,Vbias,0,0,B,C);
  dacarr.Set(1);
  dacarr[0]=adcarr[0];
  return dacarr;
 }
 
 // Perform the digitization via the ADC processeor
 adcarr=ADC(nbits,range,Vbias,Vsig,hist,B,C);
 
 // Convert the digital data into analog signals via the DAC processor
 dacarr=DAC(nbits,range,Vbias,&adcarr,peds,hist,B,C);
 
 if (hist)
 {
  TString title=hist->GetTitle();
  title.ReplaceAll("DAC","ADC-DAC (Transmit)");
  hist->SetTitle(title);
 }
 
 return dacarr;
}

TArrayD NcDSP::SampleAndHold(TF1 f,Double_t step,Double_t vmin,Double_t vmax,TH1* hist,Int_t loc) const
{
 
 TArrayD data(0);
 if (hist) hist->Reset();
 
 if (step<=0 || vmax<=vmin)
 {
  cout << " *" << ClassName() << "::SampleAndHold* Error : Invalid input step=" << step << " vmin=" << vmin << " vmax=" << vmax << endl;
  return data;
 }
 
 // The number of samples
 Int_t n=(vmax-vmin)/step;
 data.Set(n);
 
 if (hist)
 {
  hist->SetBins(n,vmin,vmax);
  TString sloc="start";
  if (!loc) sloc="center";
  if (loc>0) sloc="end";
  TString title;
  title.Form("NcDSP SampleAndHold for Function %-s in steps of %-.3g;X value;F(x) at the %-s of each step",(f.GetExpFormula("p")).Data(),step,sloc.Data());
  hist->SetTitle(title);
  hist->SetMarkerStyle(20);
 }
 
 // Enter the sampled data into the output array
 Double_t xlow=vmin;
 Double_t xval=0;
 Double_t yval=0;
 Int_t i=0;
 while (xlow<=vmax)
 {
  if (i>=n) break;
 
  if (loc<0) xval=xlow;
  if (!loc) xval=xlow+step/2.;
  if (loc>0) xval=xlow+step;
  if (xval>vmax) xval=vmax;
 
  yval=f.Eval(xval);
  data[i]=yval;
 
  if (hist) hist->SetBinContent(i+1,data[i]);
 
  xlow+=step;
  i++;
 }
 
 return data;
}

TArrayD NcDSP::SampleAndHold(Int_t n,TH1* hist,Int_t loc,Int_t jmin,Int_t jmax) const
{
 
 TArrayD data(0);
 if (hist) hist->Reset();
 
 if (fNwf<1)
 {
  cout << " *" << ClassName() << "::SampleAndHold* Error : No waveform present." << endl;
  return data;
 }
 
 if (jmax<=jmin)
 {
  jmin=0;
  jmax=fNwf-1;
 }
 
 if (n<=0 || jmin<0 || jmin>=fNwf || jmax<0 || jmax>=fNwf)
 {
  cout << " *" << ClassName() << "::SampleAndHold* Invalid input n=" << n << " jmin=" << jmin << " jmax=" << jmax << endl;
  return data;
 }
 
 // Fill the data array
 Int_t ndata=fNwf/n;
 if (fNwf%n) ndata++;
 data.Set(ndata);
 for (Int_t i=0; i<ndata; i++)
 {
  data[i]=0;
 }
 
 Double_t val=0;
 Int_t j=0;
 Int_t k=0;
 for (Int_t i=0; i<ndata; i++)
 {
  if (j>jmax) break;
 
  if (loc<0) k=j;
  if (!loc) k=j+n/2;
  if (loc>0) k=j+n;
 
  if (k<jmin) continue;
 
  if (k>jmax) k=jmax;
 
  val=fWaveform[k];
  data[i]=val;
  j+=n;
 }
 
 if (hist)
 {
  TString sloc="start";
  if (!loc) sloc="center";
  if (loc>0) sloc="end";
  TString title;
  if (fSample>0)
  {
   Double_t fnew=fSample/double(n);
   title.Form("NcDSP SampleAndHold over %-i original samplings (New: %-g samples/sec);Time in seconds;Value at the %-s of each new sample",n,fnew,sloc.Data());
   hist->SetBins(ndata,double(jmin)/fSample,double(jmax)/fSample);
  }
  else
  {
   title.Form("NcDSP SampleAndHold over %-i original samplings;New sample number;Value at the %-s of each new sample",n,sloc.Data());
   hist->SetBins(ndata,jmin,jmax);
  }
  hist->SetTitle(title);
  hist->SetMarkerStyle(20);
  for (Int_t i=1; i<=ndata; i++)
  {
   hist->SetBinContent(i,data[i-1]);
  }
 }
 
 return data;
}

TArrayD NcDSP::SampleAndSum(TF1 f,Double_t step,Double_t vmin,Double_t vmax,TH1* hist) const
{
 
 TArrayD data(0);
 if (hist) hist->Reset();
 
 if (step<=0 || vmax<=vmin)
 {
  cout << " *" << ClassName() << "::SampleAndSum* Error : Invalid input step=" << step << " vmin=" << vmin << " vmax=" << vmax << endl;
  return data;
 }
 
 // The number of samples
 Int_t n=(vmax-vmin)/step;
 data.Set(n);
 for (Int_t i=0; i<n; i++)
 {
  data[i]=0;
 }
 
 if (hist)
 {
  hist->SetBins(n,vmin,vmax);
  TString title;
  title.Form("NcDSP SampleAndSum for Function %-s with sampling steps of %-.3g;X value;Integral over each step",(f.GetExpFormula("p")).Data(),step);
  hist->SetTitle(title);
  hist->SetMarkerStyle(20);
 }
 
 // Enter the sampled data into the output array
 Double_t xlow=vmin;
 Double_t xup=0;
 Double_t yval=0;
 Int_t i=0;
 while (xlow<=vmax)
 {
  if (i>=n) break;
 
  xup=xlow+step;
  if (xup>vmax) xup=vmax;
 
  yval=f.Integral(xlow,xup);
  data[i]=yval;
 
  if (hist) hist->SetBinContent(i+1,data[i]);
 
  xlow+=step;
  i++;
 }
 
 return data;
}

TArrayD NcDSP::SampleAndSum(Int_t n,TH1* hist,Int_t jmin,Int_t jmax) const
{
 
 TArrayD data(0);
 if (hist) hist->Reset();
 
 if (fNwf<1)
 {
  cout << " *" << ClassName() << "::SampleAndSum* Error : No waveform present." << endl;
  return data;
 }
 
 if (jmax<=jmin)
 {
  jmin=0;
  jmax=fNwf-1;
 }
 
 if (n<=0 || jmin<0 || jmin>=fNwf || jmax<0 || jmax>=fNwf)
 {
  cout << " *" << ClassName() << "::SampleAndSum* Invalid input n=" << n << " jmin=" << jmin << " jmax=" << jmax << endl;
  return data;
 }
 
 // Fill the data array
 Int_t ndata=fNwf/n;
 if (fNwf%n) ndata++;
 data.Set(ndata);
 for (Int_t i=0; i<ndata; i++)
 {
  data[i]=0;
 }
 
 Double_t sum=0;
 Int_t j=0;
 for (Int_t i=0; i<ndata; i++)
 {
  sum=0;
  for (Int_t k=0; k<n; k++)
  {
   if (j<jmin) continue;
   if (j>jmax) break;
   sum+=fWaveform[j];
   j++;
  }
  data[i]=sum;
 }
 
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   Double_t fnew=fSample/double(n);
   title.Form("NcDSP SampleAndSum over %-i original samplings;Time in seconds;Summed value per %-.3g seconds",n,1./fnew);
   hist->SetBins(ndata,double(jmin)/fSample,double(jmax)/fSample);
 }
  else
  {
   title.Form("NcDSP SampleAndSum over %-i original samplings;New sample number;Summed value in each new sample",n);
   hist->SetBins(ndata,jmin,jmax);
  }
  hist->SetTitle(title);
  hist->SetMarkerStyle(20);
  for (Int_t i=1; i<=ndata; i++)
  {
   hist->SetBinContent(i,data[i-1]);
  }
 }
 
 return data;
}

TArrayD NcDSP::FilterMovingAverage(Int_t n,TString mode,TH1* hist,Int_t* i1,Int_t* i2,TH1* hisf,Bool_t dB)
{
 
 Int_t ny=0;
 TArrayD y(ny);
 if (hist) hist->Reset();
 if (hisf) hisf->Reset();
 Int_t nx=fReIn.GetSize();
 
 if (nx<1)
 {
  printf(" *%-s::FilterMovingAverage* No loaded input data present. \n",ClassName());
  return y;
 }
 
 if (n<1 || n>nx || (mode!="rec" && mode!="conv"))
 {
  printf(" *%-s::FilterMovingAverage* Inconsistent input n=%-i for x[%-i] and mode=%-s \n",ClassName(),n,nx,mode.Data());
  return y;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 if (mode=="conv") // Convolution mode
 {
  // The filter kernel
  TArrayD h=GetMovingAverageKernel(n);
  SetWaveform(&h);
 
  // Perform the convolution
  y=Convolve(hist,i1,i2);
 }
 else // Recursive mode
 {
  TArrayD x=fReIn;
  ny=nx-n+1;
  if (i1 || i2) ny=nx;
 
  if (i1) *i1=0;
  if (i2) *i2=nx-n;
 
  y.Set(ny);
  y.Reset();
 
  // Calculate the first y-value summation
  for (Int_t i=0; i<n; i++)
  {
   y[0]+=x[i];
  }
 
  // The recursive summation
  Double_t add=0;
  for (Int_t k=1; k<ny; k++)
  {
   add=0;
   if ((k+n-1)<nx) add=x[k+n-1];
   y[k]=y[k-1]-x[k-1]+add;
  }
 
  Double_t rn=n;
  // Calculate the average values
  for (Int_t i=0; i<ny; i++)
  {
   y[i]=y[i]/rn;
  }
 }
 
 TString title;
 if (hist)
 {
  if (mode=="rec")
  {
   if (fSample>0)
   {
    title.Form("%-s;Time in seconds;Value",ClassName());
    hist->SetBins(ny,0,double(ny)/fSample);
   }
   else
   {
    title.Form("%-s;Sample number;Value",ClassName());
    hist->SetBins(ny,0,ny);
   }
 
   // Set histogram axis labels and provisonal title
   hist->SetTitle(title);
 
   // Fill the histogram
   for (Int_t ibin=1; ibin<=ny; ibin++)
   {
    hist->SetBinContent(ibin,y[ibin-1]);
   }
 
   // Indicate the values of i1 and i2 (if requested) by vertical blue dashed lines
   if (i1 || i2)
   {
    Double_t ymin=hist->GetMinimum();
    Double_t ymax=hist->GetMaximum();
 
    Double_t xlow=0;
    Double_t xup=0;
    if (i1)  xlow=hist->GetBinLowEdge(*i1+1);
    if (i2)
    {
     xup=hist->GetBinLowEdge(*i2+1);
     xup+=hist->GetBinWidth(1);
    }
 
    TLine* vline1=0;
    TLine* vline2=0;
 
    if (i1)
    {
     vline1=new TLine(xlow,ymin,xlow,ymax);
     vline1->SetLineStyle(2); // Dashed line
     vline1->SetLineWidth(2);
     vline1->SetLineColor(4); // Blue color
    }
    if (i2)
    {
     vline2=new TLine(xup,ymin,xup,ymax);
     vline2->SetLineStyle(2); // Dashed line
     vline2->SetLineWidth(2);
     vline2->SetLineColor(4); // Blue color
    }
 
     TList* hlist=hist->GetListOfFunctions();
     if (vline1) hlist->Add(vline1);
     if (vline2) hlist->Add(vline2);
   }
  }
 
  // Set the appropriate histogram title
  title.Form("%-s Moving Average (%-s mode) Filter: Time domain result averaged over %-i samples",ClassName(),mode.Data(),n);
  hist->SetTitle(title);
 }
 
 // Fill the filtered frequency domain histogram
 if (hisf)
 {
  // Obtain the frequency domain histogram
  HistogramFilterFFT(&y,hisf,dB,kFALSE);
 
  // Set the appropriate histogram title
  title.Form("%-s Moving Average (%-s mode) Filter: Frequency domain result (%-i sample averaging in time domain)",ClassName(),mode.Data(),n);
  hisf->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return y;
}

TArrayD NcDSP::FilterLowPass(Double_t fcut,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Int_t* i1,Int_t* i2,Bool_t adaptn)
{
 
 TArrayD y(0);
 if (hisf) hisf->Reset();
 if (hist) hist->Reset();
 Int_t nx=fReIn.GetSize();
 
 if (nx<1)
 {
  printf(" *%-s::FilterLowPass* No loaded input data present. \n",ClassName());
  return y;
 }
 
 if (n<1 || fcut<=0 || fcut>0.5)
 {
  printf(" *%-s::FilterLowPass* Invalid input fcut=%-g n=%-i \n",ClassName(),fcut,n);
  return y;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 // The filter kernel
 TArrayD h=GetLowPassKernel(fcut,n,0,0,0,adaptn);
 
 // Perform the convolution
 SetWaveform(&h);
 y=Convolve(hist,i1,i2,2);
 
 // Set title for the filtered time domain histogram
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   Float_t nucut=fcut*fSample;
   title.Form("%-s Low Pass Filter: Time domain result with #nu-cut=%-.3g Hz (%-i-point kernel)",ClassName(),nucut,n);
  }
  else
  {
   title.Form("%-s Low Pass Filter: Time domain result with #nu-cut/#nu-sample=%-.3g (%-i-point kernel)",ClassName(),fcut,n);
  }
  hist->SetTitle(title);
 }
 
 // Fill the filtered frequency domain histogram
 if (hisf)
 {
  // Obtain the frequency domain histogram
  HistogramFilterFFT(&y,hisf,dB,kFALSE);
 
  // Set the appropriate histogram title
  TString title;
  if (fSample>0)
  {
   Float_t nucut=fcut*fSample;
   title.Form("%-s Low Pass Filter: Frequency domain result with #nu-cut=%-.3g Hz (%-i-point time domain kernel)",ClassName(),nucut,n);
  }
  else
  {
   title.Form("%-s Low Pass Filter: Frequency domain result with #nu-cut/#nu-sample=%-.3g (%-i-point time domain kernel)",ClassName(),fcut,n);
  }
  hisf->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return y;
}

TArrayD NcDSP::FilterHighPass(Double_t fcut,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Int_t* i1,Int_t* i2,Bool_t adaptn)
{
 
 TArrayD y(0);
 if (hisf) hisf->Reset();
 if (hist) hist->Reset();
 Int_t nx=fReIn.GetSize();
 
 if (nx<1)
 {
  printf(" *%-s::FilterHighPass* No loaded input data present. \n",ClassName());
  return y;
 }
 
 if (n<1 || fcut<=0 || fcut>0.5)
 {
  printf(" *%-s::FilterHighPass* Invalid input fcut=%-g n=%-i \n",ClassName(),fcut,n);
  return y;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 // The filter kernel
 TArrayD h=GetHighPassKernel(fcut,n,0,0,0,adaptn);
 
 // Perform the convolution
 SetWaveform(&h);
 y=Convolve(hist,i1,i2,2);
 
 // Set title for the filtered time domain histogram
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   Float_t nucut=fcut*fSample;
   title.Form("NcDSP High Pass Filter: Time domain result with #nu-cut=%-.3g Hz (%-i-point kernel)",nucut,n);
  }
  else
  {
   title.Form("NcDSP High Pass Filter: Time domain result with #nu-cut/#nu-sample=%-.3g (%-i-point kernel)",fcut,n);
  }
  hist->SetTitle(title);
 }
 
 // Fill the filtered frequency domain histogram
 if (hisf)
 {
  // Obtain the frequency domain histogram
  HistogramFilterFFT(&y,hisf,dB,kFALSE);
 
  // Set the appropriate histogram title
  TString title;
  if (fSample>0)
  {
   Float_t nucut=fcut*fSample;
   title.Form("NcDSP High Pass Filter: Frequency domain result with #nu-cut=%-.3g Hz (%-i-point time domain kernel)",nucut,n);
  }
  else
  {
   title.Form("NcDSP High Pass Filter: Frequency domain result with #nu-cut/#nu-sample=%-.3g (%-i-point time domain kernel)",fcut,n);
  }
  hisf->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return y;
}

TArrayD NcDSP::FilterBandPass(Double_t f1,Double_t f2,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Int_t* i1,Int_t* i2,Bool_t adaptn)
{
 
 TArrayD y(0);
 if (hisf) hisf->Reset();
 if (hist) hist->Reset();
 Int_t nx=fReIn.GetSize();
 
 if (nx<1)
 {
  printf(" *%-s::FilterBandPass* No loaded input data present. \n",ClassName());
  return y;
 }
 
 if (n<1 || f1<=0 || f1>0.5 || f2<=0 || f2>0.5 || f2<=f1)
 {
  printf(" *%-s::FilterBandPass* Invalid input f1=%-g f2=%-g n=%-i \n",ClassName(),f1,f2,n);
  return y;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 // The filter kernel
 TArrayD h=GetBandPassKernel(f1,f2,n,0,0,0,adaptn);
 
 // Perform the convolution
 SetWaveform(&h);
 y=Convolve(hist,i1,i2,2);
 
 // Set title for the filtered time domain histogram
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   Float_t nu1=f1*fSample;
   Float_t nu2=f2*fSample;
   title.Form("NcDSP Band Pass Filter: Time domain result for [%-.3g,%-.3g] Hz (%-i-point kernel)",nu1,nu2,n);
  }
  else
  {
   title.Form("NcDSP Band Pass Filter: Time domain result for #nu/#nu-sample=[%-.3g,%-.3g] (%-i-point kernel)",f1,f2,n);
  }
  hist->SetTitle(title);
 }
 
 // Fill the filtered frequency domain histogram
 if (hisf)
 {
  // Obtain the frequency domain histogram
  HistogramFilterFFT(&y,hisf,dB,kFALSE);
 
  // Set the appropriate histogram title
  TString title;
  if (fSample>0)
  {
   Float_t nu1=f1*fSample;
   Float_t nu2=f2*fSample;
   title.Form("NcDSP Band Pass Filter: Frequency domain result for [%-.3g,%-.3g] Hz (%-i-point time domain kernel)",nu1,nu2,n);
  }
  else
  {
   title.Form("NcDSP Band Pass Filter: Frequency domain result for #nu/#nu-sample=[%-.3g,%-.3g] (%-i-point time domain kernel)",f1,f2,n);
  }
  hisf->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return y;
}

TArrayD NcDSP::FilterBandReject(Double_t f1,Double_t f2,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Int_t* i1,Int_t* i2,Bool_t adaptn)
{
 
 TArrayD y(0);
 if (hisf) hisf->Reset();
 if (hist) hist->Reset();
 Int_t nx=fReIn.GetSize();
 
 if (nx<1)
 {
  printf(" *%-s::FilterBandReject* No loaded input data present. \n",ClassName());
  return y;
 }
 
 if (n<1 || f1<=0 || f1>0.5 || f2<=0 || f2>0.5 || f2<=f1)
 {
  printf(" *%-s::FilterBandReject* Invalid input f1=%-g f2=%-g n=%-i \n",ClassName(),f1,f2,n);
  return y;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 // The filter kernel
 TArrayD h=GetBandRejectKernel(f1,f2,n,0,0,0,adaptn);
 
 // Perform the convolution
 SetWaveform(&h);
 y=Convolve(hist,i1,i2,2);
 
 // Set title for the filtered time domain histogram
 if (hist)
 {
  TString title;
  if (fSample>0)
  {
   Float_t nu1=f1*fSample;
   Float_t nu2=f2*fSample;
   title.Form("NcDSP Band Reject Filter: Time domain result for [%-.3g,%-.3g] Hz (%-i-point kernel)",nu1,nu2,n);
  }
  else
  {
   title.Form("NcDSP Band Reject Filter: Time domain result for #nu/#nu-sample=[%-.3g,%-.3g] (%-i-point kernel)",f1,f2,n);
  }
  hist->SetTitle(title);
 }
 
 // Fill the filtered frequency domain histogram
 if (hisf)
 {
  // Obtain the frequency domain histogram
  HistogramFilterFFT(&y,hisf,dB,kFALSE);
 
  // Set the appropriate histogram title
  TString title;
  if (fSample>0)
  {
   Float_t nu1=f1*fSample;
   Float_t nu2=f2*fSample;
   title.Form("NcDSP Band Reject Filter: Frequency domain result for [%-.3g,%-.3g] Hz (%-i-point time domain kernel)",nu1,nu2,n);
  }
  else
  {
   title.Form("NcDSP Band Reject Filter: Frequency domain result for #nu/#nu-sample=[%-.3g,%-.3g] (%-i-point time domain kernel)",f1,f2,n);
  }
  hisf->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return y;
}

TArrayD NcDSP::FilterMultiBand(TArray& freqs,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Int_t* i1,Int_t* i2,Bool_t adaptn)
{
 
 TArrayD y(0);
 if (hisf) hisf->Reset();
 if (hist) hist->Reset();
 Int_t nx=fReIn.GetSize();
 
 if (nx<1)
 {
  printf(" *%-s::FilterMultiBand* No loaded input data present. \n",ClassName());
  return y;
 }
 
 Int_t arrsize=freqs.GetSize();
 Int_t oddsize=arrsize%2;
 Int_t nbands=arrsize/2;
 if (nbands<1 || n<1 || oddsize)
 {
  printf(" *%-s::FilterMultiBand* Invalid input array size=%-i nbands=%-i n=%-i \n",ClassName(),arrsize,nbands,n);
  return y;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 // The filter kernel
 TArrayD h=GetMultiBandKernel(freqs,n,0,0,0,adaptn);
 
 // Convolve the composite kernel with the original time domain data
 Load(&xsave);
 SetWaveform(&h);
 y=Convolve(hist,i1,i2,2);
 
 // Determine the number of actually specified bands
 Double_t flow=0;
 Double_t fup=0;
 Int_t index=0;
 Int_t neff=0; // The number of off actually specified bands
 // Loop over the specified frequency bands
 for (Int_t jband=1; jband<=nbands; jband++)
 {
  index=2*(jband-1);
  flow=freqs.GetAt(index);
  fup=freqs.GetAt(index+1);
 
  // Skip empty entries in the "freqs" array
  if (!flow || !fup) continue;
 
  neff++;
 }
 
 TString title;
 // Set title for the filtered time domain histogram
 if (hist)
 {
  title.Form("%-s MultiBand Filter: Time domain result for %-i bands (%-i-point kernel each)",ClassName(),neff,n);
  hist->SetTitle(title);
 }
 
 // Fill the filtered frequency domain histogram
 if (hisf)
 {
  // Obtain the frequency domain histogram
  HistogramFilterFFT(&y,hisf,dB,kFALSE);
 
  // Set the appropriate histogram title
  title.Form("%-s MultiBand Filter: Frequency domain result for %-i bands (%-i-point time domain kernel each)",ClassName(),neff,n);
  hisf->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return y;
}

TArrayD NcDSP::GetMovingAverageKernel(Int_t n,TH1* hisf,Bool_t dB,TH1* hist)
{
 
 TArrayD h(0);
 if (hisf) hisf->Reset();
 if (hist) hist->Reset();
 
 if (n<1)
 {
  printf(" *%-s::GetMovingAverageKernel* Invalid input n=%-i \n",ClassName(),n);
  return h;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 h.Set(n);
 for (Int_t i=0; i<n; i++)
 {
  h[i]=1./double(n);
 }
 
 HistogramFilterFFT(&h,hisf,dB,kTRUE,hist);
 
 // Set the appropriate histogram titles
 TString title;
 if (hisf)
 {
  title.Form("NcDSP Moving Average Filter: Frequency domain kernel (%-i-point time domain kernel)",n);
  hisf->SetTitle(title);
 }
 
 if (hist)
 {
  title.Form("NcDSP Moving Average Filter: Time domain kernel (%-i-point zero padded)",n);
  hist->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return h;
}

TArrayD NcDSP::GetLowPassKernel(Double_t fcut,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Bool_t adaptn)
{
 
 TArrayD h(0);
 if (hisf) hisf->Reset();
 
 if (n<1 || fcut<=0 || fcut>0.5)
 {
  printf(" *%-s::GetLowPassKernel* Invalid input fcut=%-g n=%-i \n",ClassName(),fcut,n);
  return h;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 Double_t twopi=2.*acos(-1.);
 h.Set(n);
 Int_t m=n-1;
 Double_t rm=m;
 Double_t ri=0;
 Double_t sum=0;
 for (Int_t i=0; i<=m; i++)
 {
  if (i==m/2)
  {
   h[i]=twopi*fcut;
  }
  else
  {
   ri=i;
   h[i]=sin(twopi*fcut*double(i-m/2))*(0.42-0.5*cos(twopi*ri/rm)+0.08*cos(2.*twopi*ri/rm))/double(i-m/2);
  }
  sum+=h[i];
 }
 
 // Normalize the filter kernel to 1
 for (Int_t i=0; i<n; i++)
 {
  h[i]=h[i]/sum;
 }
 
 HistogramFilterFFT(&h,hisf,dB,kTRUE,hist);
 
 // Set the appropriate histogram titles
 TString title;
 if (hisf)
 {
  if (fSample>0)
  {
   Float_t nucut=fcut*fSample;
   title.Form("NcDSP Low Pass Filter: Frequency domain kernel with #nu-cut=%-.3g Hz (%-i-point time domain kernel)",nucut,n);
  }
  else
  {
   title.Form("NcDSP Low Pass Filter: Frequency domain kernel with #nu-cut/#nu-sample=%-.3g (%-i-point time domain kenel)",fcut,n);
  }
  hisf->SetTitle(title);
 }
 
 if (hist)
 {
  if (fSample>0)
  {
   Float_t nucut=fcut*fSample;
   title.Form("NcDSP Low Pass Filter: Time domain kernel (%-i-point zero padded) with #nu-cut=%-.3g Hz",n,nucut);
  }
  else
  {
   title.Form("NcDSP Low Pass Filter: Time domain kernel (%-i-point zero padded) with #nu-cut/#nu-sample=%-.3g",n,fcut);
  }
  hist->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return h;
}

TArrayD NcDSP::GetHighPassKernel(Double_t fcut,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Bool_t adaptn)
{
 
 TArrayD h(0);
 if (hisf) hisf->Reset();
 
 if (n<1 || fcut<=0 || fcut>0.5)
 {
  printf(" *%-s::GetHighPassKernel* Invalid input fcut=%-g n=%-i \n",ClassName(),fcut,n);
  return h;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 // The corresponding Low Pass filter kernel
 h=GetLowPassKernel(fcut,n,0,0,0,adaptn);
 
 // Spectrally invert the filter kernel to obtain the High Pass kernel
 for (Int_t i=0; i<n; i++)
 {
  h[i]=-h[i];
 }
 h[n/2]=h[n/2]+1.;
 
 HistogramFilterFFT(&h,hisf,dB,kTRUE,hist);
 
 // Set the appropriate histogram title
 TString title;
 if (hisf)
 {
  if (fSample>0)
  {
   Float_t nucut=fcut*fSample;
   title.Form("NcDSP High Pass Filter: Frequency domain kernel with #nu-cut=%-.3g Hz (%-i-point time domain kernel)",nucut,n);
  }
  else
  {
   title.Form("NcDSP High Pass Filter: Frequency domain kernel with #nu-cut/#nu-sample=%-.3g (%-i-point time domain kernel)",fcut,n);
  }
  hisf->SetTitle(title);
 }
 
 if (hist)
 {
  if (fSample>0)
  {
   Float_t nucut=fcut*fSample;
   title.Form("NcDSP High Pass Filter: Time domain kernel (%-i-point zero padded) with #nu-cut=%-.3g Hz",n,nucut);
  }
  else
  {
   title.Form("NcDSP High Pass Filter: Time domain kernel (%-i-point zero padded) with #nu-cut/#nu-sample=%-.3g",n,fcut);
  }
  hist->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return h;
}

TArrayD NcDSP::GetBandPassKernel(Double_t f1,Double_t f2,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Bool_t adaptn)
{
 
 TArrayD h(0);
 if (hisf) hisf->Reset();
 
 if (n<1 || f1<=0 || f1>0.5 || f2<=0 || f2>0.5 || f2<=f1)
 {
  printf(" *%-s::GetBandPassKernel* Invalid input f1=%-g f2=%-g n=%-i \n",ClassName(),f1,f2,n);
  return h;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 // Get the corresponding Band Reject kernel
 h=GetBandRejectKernel(f1,f2,n,0,0,0,adaptn); 
 
 // Spectrally invert the Band Reject filter kernel to obtain the Band Pass kernel
 for (Int_t i=0; i<n; i++)
 {
  h[i]=-h[i];
 }
 h[n/2]=h[n/2]+1.;
 
 HistogramFilterFFT(&h,hisf,dB,kTRUE,hist);
 
 // Set the appropriate histogram title
 TString title;
 if (hisf)
 {
  if (fSample>0)
  {
   Float_t nu1=f1*fSample;
   Float_t nu2=f2*fSample;
   title.Form("NcDSP Band Pass Filter: Frequency domain kernel for [%-.3g,%-.3g] Hz (%-i-point time domain kernel)",nu1,nu2,n);
  }
  else
  {
   title.Form("NcDSP Band Pass Filter: Frequency domain kernel for #nu/#nu-sample=[%-.3g,%-.3g] (%-i-point time domain kernel)",f1,f2,n);
  }
  hisf->SetTitle(title);
 }
 
 if (hist)
 {
  if (fSample>0)
  {
   Float_t nu1=f1*fSample;
   Float_t nu2=f2*fSample;
   title.Form("NcDSP Band Pass Filter: Time domain kernel (%-i-point zero padded) for [%-.3g,%-.3g] Hz",n,nu1,nu2);
  }
  else
  {
   title.Form("NcDSP Band Pass Filter: Time domain kernel (%-i-point zero padded) for #nu/#nu-sample=[%-.3g,%-.3g]",n,f1,f2);
  }
  hist->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return h;
}

TArrayD NcDSP::GetBandRejectKernel(Double_t f1,Double_t f2,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Bool_t adaptn)
{
 
 TArrayD h(0);
 if (hisf) hisf->Reset();
 
 if (n<1 || f1<=0 || f1>0.5 || f2<=0 || f2>0.5 || f2<=f1)
 {
  printf(" *%-s::GetBandRejectKernel* Invalid input f1=%-g f2=%-g n=%-i \n",ClassName(),f1,f2,n);
  return h;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 // The Low Pass filter kernel for f1
 TArrayD hlow=GetLowPassKernel(f1,n,0,0,0,adaptn);
 
 // The High Pass filter kernel for f2
 TArrayD hhigh=GetHighPassKernel(f2,n,0,0,0,adaptn);
 
 // Add the Low Pass and High Pass kernels to obtain a Band Reject kernel
 h.Set(n);
 for (Int_t i=0; i<n; i++)
 {
  h[i]=hlow[i]+hhigh[i];
 }
 
 HistogramFilterFFT(&h,hisf,dB,kTRUE,hist);
 
 // Set the appropriate histogram title
 TString title;
 if (hisf)
 {
  if (fSample>0)
  {
   Float_t nu1=f1*fSample;
   Float_t nu2=f2*fSample;
   title.Form("NcDSP Band Reject Filter: Frequency domain kernel for [%-.3g,%-.3g] Hz (%-i-point time domain kernel)",nu1,nu2,n);
  }
  else
  {
   title.Form("NcDSP Band Reject Filter: Frequency domain kernel for #nu/#nu-sample=[%-.3g,%-.3g] (%-i-point time domain kernel)",f1,f2,n);
  }
  hisf->SetTitle(title);
 }
 
 if (hist)
 {
  if (fSample>0)
  {
   Float_t nu1=f1*fSample;
   Float_t nu2=f2*fSample;
   title.Form("NcDSP Band Reject Filter: Time domain kernel (%-i-point zero padded) for [%-.3g,%-.3g] Hz",n,nu1,nu2);
  }
  else
  {
   title.Form("NcDSP Band Reject Filter: Time domain kernel (%-i-point zero padded) for #nu/#nu-sample=[%-.3g,%-.3g]",n,f1,f2);
  }
  hist->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return h;
}

TArrayD NcDSP::GetMultiBandKernel(TArray& freqs,Int_t n,TH1* hisf,Bool_t dB,TH1* hist,Bool_t adaptn)
{
 
 TArrayD h(0); // The convolution of the various filter kernels
 if (hisf) hisf->Reset();
 
 Int_t arrsize=freqs.GetSize();
 Int_t oddsize=arrsize%2;
 Int_t nbands=arrsize/2;
 if (nbands<1 || n<1 || oddsize)
 {
  printf(" *%-s::GetMultiBandKernel* Invalid input array size=%-i nbands=%-i n=%-i \n",ClassName(),arrsize,nbands,n);
  return h;
 }
 
 // Save the current stored data
 TArrayD xsave=fReIn;
 TArrayD wfsave=fWaveform;
 
 // Adapt "n" to an odd value if needed
 Int_t odd=n%2;
 if (!odd && adaptn) n++;
 
 Double_t flow=0;
 Double_t fup=0;
 Int_t index=0;
 TArrayD hj; // The filter kernel for a single band
 Bool_t first=kTRUE;
 Int_t neff=0; // The number of off actually specified bands
 // Loop over the specified frequency bands
 for (Int_t jband=1; jband<=nbands; jband++)
 {
  index=2*(jband-1);
  flow=freqs.GetAt(index);
  fup=freqs.GetAt(index+1);
 
  // Skip empty entries in the "freqs" array
  if (!flow || !fup) continue;
 
  neff++;
 
  if (flow>0 && fup>0) hj=GetBandPassKernel(flow,fup,n,0,0,0,adaptn);
 
  if (flow<0 && fup<0) hj=GetBandRejectKernel(fabs(flow),fabs(fup),n,0,0,0,adaptn);
 
  if (flow<0 && fup>0) hj=GetLowPassKernel(fup,n,0,0,0,adaptn);
 
  if (flow>0 && fup<0) hj=GetHighPassKernel(flow,n,0,0,0,adaptn);
 
  SetWaveform(&hj);
 
  if (first)
  {
   h=hj;
   first=kFALSE;
  }
  else
  {
    Load(&h);
    h=Convolve();
  }
 }
 
 HistogramFilterFFT(&h,hisf,dB,kTRUE,hist);
 
 // Set the appropriate histogram titles
 TString title;
 if (hisf)
 {
  title.Form("%-s MultiBand Filter: Frequency domain kernel for %-i bands (%-i-point time domain kernel each)",ClassName(),neff,n);
  hisf->SetTitle(title);
 }
 
 if (hist)
 {
  title.Form("%-s MultiBand Filter: Time domain kernel (zero padded) for %-i bands (%-i-point kernel each)",ClassName(),neff,n);
  hist->SetTitle(title);
 }
 
 // Restore the original input data
 Load(&xsave);
 SetWaveform(&wfsave);
 
 // Let user invokations of Load() reset the output data again
 fKeepOutput=kFALSE;
 
 return h;
}

void NcDSP::HistogramFilterFFT(TArray* h,TH1* hisf,Bool_t dB,Bool_t kernel,TH1* hist)
{
 
 if (hisf) hisf->Reset("M");
 if (hist) hist->Reset("M");
 
 if (!h) return;
 
 Int_t nh=h->GetSize();
 
 // Create a temporary array to contain a zero padded copy
 // of the input time domain array "h".
 // This allows to get a better view in the time domain filter kernel
 // and a better resolution in the DFT produced frequency spectrum.
 Int_t npad=0;   // The total number of padded zeros
 Int_t nfront=0; // The number of padded zeros at the front of the array
 Int_t narr=0;   // The length of the zero padded array
 TArrayD arr(narr);
 
 // The time domain kernel histogram
 if (hist && kernel)
 {
  // The temporary array will be made twice the length of "h"
  // and symmetrically zero padded to display the kernel in the center of the histogram
  narr=2*nh;
  nfront=nh/2;
  arr.Set(narr);
  arr.Reset();
  for (Int_t i=0; i<nh; i++)
  {
   arr[nfront+i]=h->GetAt(i);
  }
 
  if (fSample>0)
  {
   hist->SetTitle("NcDSP HistogramFilterFFT time domain kernel (zero padded);Time in seconds;Value");
   hist->SetBins(narr,0,double(narr)/fSample);
  }
  else
  {
   hist->SetTitle("NcDSP HistogramFilterFFT time domain kernel (zero padded);Sample number;Value");
   hist->SetBins(narr,0,narr);
  }
  for (Int_t i=1; i<=narr; i++)
  {
   hist->SetBinContent(i,arr[i-1]);
  }
  hist->SetLineWidth(2);
  hist->SetLineColor(kBlue);
 }
 
 if (!hisf) return;
 
 // The frequency domain kernel or filter result histogram.
 // Create a temporary zero-padded time domain array with a length of narr=2^k for the FFT,
 // such that the frequency spectrum contains (narr/2)+1 samples.
 // The minimum array length will be 1024 samples.
 Int_t k=log(nh)/log(2);
 k++;
 if (k<10) k=10; // Minimal 1024 samples
 narr=pow(2,k);
 npad=narr-nh;
 nfront=npad/2;
 
 arr.Set(narr);
 arr.Reset();
 for (Int_t i=0; i<nh; i++)
 {
  arr[nfront+i]=h->GetAt(i);
 }
 
 // Load the temporary zero-padded time domain (kernel) data for Fourier transformation
 Load(&arr);
 
 TString sel="dB";
 if (!dB) sel="AMP";
 if (fSample>0)
 {
  sel+=" Hz";
 }
 else
 {
  sel+=" f";
 }
 
 // Perform the Fourier transform
 Fourier("R2C",hisf,sel);
 
 // Perform the Fourier transform with the original unpadded time domain data
 // in order to store Fourier coefficients that allow retrieving the original
 // unpadded time domain filter result via an inverse Fourier transform from
 // the frequency domain filter result that was produced here.
 
 // Load the original time domain input data
 Load(h);
 
 // Perform the Fourier transform
 Fourier("R2C");
 
 // Prevent following internal Load() invokations to reset the output
 fKeepOutput=kTRUE;
 
 if (!kernel) return;
 
 // Normalize the maximum amplitude of "hisf" to 1 (or 0 dB)
 Double_t max=hisf->GetMaximum();
 if (!dB)
 {
  if (max) hisf->Scale(1./fabs(max));
 }
 else
 {
  for (Int_t i=1; i<=hisf->GetNbinsX(); i++)
  {
   hisf->AddBinContent(i,-max);
  }
 }
 hisf->SetLineWidth(2);
 hisf->SetLineColor(kBlue);
}

TGraph NcDSP::Periodogram(TString tu,Double_t Tmin,Double_t Tmax,Int_t n,TArray& t,TArray& y,TArray* dy,TF1* Z,NcDevice* results) const
{
 
  if (results) results->Reset();
 
  TGraph gr;
 
  if (tu!="ps" && tu!="ns" && tu!="sec" && tu!="hour" && tu!="day")
  {
   printf("*%-s::Periodogram* Unsupported time unit : %-s \n",ClassName(),tu.Data());
   return gr;
  } 
 
  if (Tmin<=0 || (Tmax>Tmin && Tmax<=0) || n<1)
  {
   printf("*%-s::Periodogram* Inconsistent input : Tmin=%-g Tmax=%-g n=%-i \n",ClassName(),Tmin,Tmax,n);
   return gr;
  }
 
  // Check the compatibility of the input array sizes
  Int_t nt=t.GetSize();
  Int_t ny=y.GetSize();
  Int_t ndy=ny;
  if (dy) ndy=dy->GetSize();
 
  if ((nt+ny+ndy)!=3*nt)
  {
   if (dy)
   {
    printf("*%-s::Periodogram* Incompatible input array sizes : nt=%-i ny=%-i ndy=%-i \n",ClassName(),nt,ny,ndy);
   }
   else
   {
    printf("*%-s::Periodogram* Incompatible input array sizes : nt=%-i ny=%-i \n",ClassName(),nt,ny);
   }
   return gr;
  }
 
  // Determine the frequency bounds and step size
  Float_t fmax=1./Tmin;
  Float_t fmin=0;
  if (Tmax>Tmin) fmin=1./Tmax;
  Float_t fstep=(fmax-fmin)/float(n);
 
  // Create compact internal arrays for efficient processing 
  Int_t na=0;
  Double_t ti,yi,dyi2,wi,zi;
  TArrayD tarr(nt);
  TArrayD yarr(nt);
  TArrayD warr(nt);
  TArrayD zarr(nt);
  Double_t Y=0;
  Double_t Y2=0;
  for (Int_t i=0; i<nt; i++)
  {
   ti=t.GetAt(i);
   yi=y.GetAt(i);
   wi=1;
   zi=1;
   if (dy)
   {
    dyi2=pow(dy->GetAt(i),2);
    if (!dyi2) continue;
    wi=1./dyi2;
   }
   if (Z) zi=Z->Eval(ti);
   tarr[na]=ti;
   yarr[na]=yi;
   warr[na]=wi;
   zarr[na]=zi;
   Y+=wi*yi;
   Y2+=wi*pow(yi,2);
   na++;
  }
 
  tarr.Set(na);
  yarr.Set(na);
  warr.Set(na);
  Double_t W=warr.GetSum();
 
  if (!na || W<=0)
  {
   printf("*%-s::Periodogram* No valid data available. \n",ClassName());
   return gr;
  }
 
  Y=Y/W;
  Y2=Y2/W;
 
  // Loop over the frequencies
  Double_t twopi=2.*acos(-1.);
  Double_t omega=0;
  Double_t zcosi=0;
  Double_t zsini=0;
  Double_t C=0;
  Double_t S=0;
  Double_t CS=0;
  Double_t C2=0;
  Double_t S2=0;
  Double_t YC=0;
  Double_t YS=0;
  Double_t D=0;
  Double_t P=0;
  Double_t a=0;
  Double_t b=0;
  Double_t c=0;
  Int_t ip=0;
  NcSignal sx;
  sx.AddNamedSlot("f");
  sx.AddNamedSlot("a");
  sx.AddNamedSlot("b");
  sx.AddNamedSlot("c");
  sx.AddNamedSlot("P");
  Float_t f=fmin;
  fmax+=0.1*fstep; // To prevent accuracy issues
  while (f<=fmax)
  {
   omega=twopi*f;
   C=0;
   S=0;
   CS=0;
   C2=0;
   S2=0;
   YC=0;
   YS=0;
   // Loop over the data points
   for (Int_t i=0; i<na; i++)
   {
    ti=tarr[i];
    yi=yarr[i];
    wi=warr[i]/W;
    zi=zarr[i];
    zcosi=zi*cos(omega*ti);
    zsini=zi*sin(omega*ti);
    C+=wi*zcosi;
    S+=wi*zsini;
    CS+=wi*zcosi*zsini;
    C2+=wi*pow(zcosi,2);
    S2+=wi*pow(zsini,2);
    YC+=wi*yi*zcosi;
    YS+=wi*yi*zsini;
   }
 
   CS-=C*S;
   Y2-=pow(Y,2);
   C2-=pow(C,2);
   S2-=pow(S,2);
   YC-=Y*C;
   YS-=Y*S;
 
   D=C2*S2-pow(CS,2);
 
   P=0;
   if (Y2*D) P=(S2*pow(YC,2)+C2*pow(YS,2)-2.*CS*YC*YS)/(Y2*D);
 
   a=0;
   b=0;
   if (D)
   {
    a=(YC*S2-YS*CS)/D;
    b=(YS*C2-YC*CS)/D;
   }
   c=Y-a*C-b*S;
 
   // Record this (f,a,b,c,P) dataset as a hit in the (optional) NcDevice "results"
   if (results)
   {
    sx.SetSignal(f,"f");
    sx.SetSignal(a,"a");
    sx.SetSignal(b,"b");
    sx.SetSignal(c,"c");
    sx.SetSignal(P,"P");
    results->AddHit(sx);
   }
 
   // Record this entry in the periodogram graph
   gr.SetPoint(ip,f,P);
 
   f+=fstep;
   ip++;
  }
 
 TString unit="Hz";
 if (tu=="ps") unit="THz";
 if (tu=="ns") unit="GHz";
 if (tu=="hour") unit="1/hour";
 if (tu=="day") unit="1/day";
 TString title;
 title.Form("%-s Periodogram;Frequency in %-s;Normalized Power",ClassName(),unit.Data());
 gr.SetNameTitle("Periodogram",title);
 gr.SetMarkerStyle(8);
 gr.SetMarkerSize(0.5);
 gr.SetMarkerColor(kBlue);
 gr.SetLineWidth(2);
 
 return gr;
}

TGraph NcDSP::Periodogram(TString tu,Double_t Tmin,Double_t Tmax,Int_t n,NcSample s,Int_t it,Int_t iy,Int_t idy,TF1* Z,NcDevice* results) const
{
 
 TGraph gr;
 
 if (!s.GetStoreMode())
 {
  printf("*%-s::Periodogram* Error: NcSample storage mode is not activated. \n",ClassName());
  return gr;
 }
 
 Int_t ndim=s.GetDimension();
 Int_t ns=s.GetN();
 
 if (ns<1 || it<1 || it>ndim || iy<1 || iy>ndim || idy<0 || idy>ndim)
 {
  printf("*%-s::Periodogram* Inconsistent NcSample input: it=%-i iy=%-i idy=%-i #entries=%-i \n",ClassName(),it,iy,idy,ns);
  return gr;
 }
 
 TArrayD t(ns);
 TArrayD y(ns);
 TArrayD dy(ns);
 
 for (Int_t i=1; i<=ns; i++)
 {
  t[i-1]=s.GetEntry(i,it);
  y[i-1]=s.GetEntry(i,iy);
  if (idy) dy[i-1]=s.GetEntry(i,idy);
 }
 
 if (idy) // Uncertainties will be used as weights
 {
  gr=Periodogram(tu,Tmin,Tmax,n,t,y,&dy,Z,results);
 }
 else // No weights will be used
 {
  gr=Periodogram(tu,Tmin,Tmax,n,t,y,0,Z,results);
 }
 
 return gr;
}

TGraph NcDSP::Periodogram(TString tu,Double_t Tmin,Double_t Tmax,Int_t n,NcSample s,TString namet,TString namey,TString namedy,TF1* Z,NcDevice* results) const
{
 
 TGraph gr;
 
 // Get the indices corresponding to the specified variable names
 Int_t it=s.GetIndex(namet);
 Int_t iy=s.GetIndex(namey);
 Int_t idy=s.GetIndex(namedy);
 
 // Check for unsupported namedy
 if (!idy && namedy!="-") idy=-1;
 
 Int_t ns=s.GetN();
 
 if (ns<1 || it<1 || iy<1 || idy<0)
 {
  printf("*%-s::Periodogram* Inconsistent NcSample input: namet=%-s namey=%-s namedy=%-s #entries=%-i \n",ClassName(),namet.Data(),namey.Data(),namedy.Data(),ns);
  return gr;
 }
 
 gr=Periodogram(tu,Tmin,Tmax,n,s,it,iy,idy,Z,results);
 
 return gr;
}

TGraph NcDSP::Periodogram(TString tu,Double_t Tmin,Double_t Tmax,Int_t n,TH1& h,Bool_t err,TF1* Z,NcDevice* results) const
{
 
 TGraph gr;
 
 Int_t nbins=h.GetNbinsX();
 Int_t nen=h.GetEntries();
 
 if (nbins<1 || nen<1)
 {
  printf("*%-s::Periodogram* Inconsistent histogram input: nbins=%-i #entries=%-i \n",ClassName(),nbins,nen);
  return gr;
 }
 
 TArrayD t(nbins);
 TArrayD y(nbins);
 TArrayD dy(nbins);
 
 for (Int_t i=1; i<=nbins; i++)
 {
  t[i-1]=h.GetBinCenter(i);
  y[i-1]=h.GetBinContent(i);
  if (err) dy[i-1]=h.GetBinError(i);
 }
 
 if (err) // Uncertainties will be used as weights
 {
  gr=Periodogram(tu,Tmin,Tmax,n,t,y,&dy,Z,results);
 }
 else // No weights will be used
 {
  gr=Periodogram(tu,Tmin,Tmax,n,t,y,0,Z,results);
 }
 
 return gr;
}

TGraph NcDSP::Periodogram(TString tu,Double_t Tmin,Double_t Tmax,Int_t n,TGraph& g,Bool_t err,TF1* Z,NcDevice* results) const
{
 
 TGraph gr;
 
 if (err && !(g.InheritsFrom("TGraphErrors")))
 {
  printf("*%-s::Periodogram* Error: For weighted fitting the input graph has to be a TGraphErrors object. \n",ClassName());
  return gr;
 }
 
 Int_t np=g.GetN();
 
 if (np<1)
 {
  printf("*%-s::Periodogram* Error: Input graph is empty. \n",ClassName());
  return gr;
 }
 
 TArrayD t(np);
 TArrayD y(np);
 TArrayD dy(np);
 
 Double_t gx=0;
 Double_t gy=0;
 for (Int_t i=0; i<np; i++)
 {
  g.GetPoint(i,gx,gy);
  t[i]=gx;
  y[i]=gy;
  if (err) dy[i]=g.GetErrorY(i);
 }
 
 if (err) // Uncertainties will be used as weights
 {
  gr=Periodogram(tu,Tmin,Tmax,n,t,y,&dy,Z,results);
 }
 else // No weights will be used
 {
  gr=Periodogram(tu,Tmin,Tmax,n,t,y,0,Z,results);
 }
 
 return gr;
}

TObject* NcDSP::Clone(const char* name) const
{
 
 NcDSP* q=new NcDSP(*this);
 if (name)
 {
  if (strlen(name)) q->SetName(name);
 }
 return q;
}