some more testing towards fourier

src/tests/fourier/Makefile.fourierPhaseTest

@@ -0,0 +1,94 @@
+#---------------------------------------------------
+# Makefile.fourierPhaseTest
+#
+# Author: Andreas Suter
+# e-mail: andreas.suter@psi.ch
+#
+# $Id$
+#---------------------------------------------------
+
+#---------------------------------------------------
+# get compilation and library flags from root-config
+
+ROOTCFLAGS    = $(shell $(ROOTSYS)/bin/root-config --cflags)
+ROOTLIBS      = $(shell $(ROOTSYS)/bin/root-config --libs)
+ROOTGLIBS     = $(shell $(ROOTSYS)/bin/root-config --glibs)
+
+#---------------------------------------------------
+# depending on the architecture, choose the compiler,
+# linker, and the flags to use
+#
+
+OSTYPE = linux
+
+ifeq ($(OSTYPE),linux)
+OS = LINUX
+endif
+ifeq ($(OSTYPE),linux-gnu)
+OS = LINUX
+endif
+ifeq ($(OSTYPE),darwin)
+OS = DARWIN
+endif
+
+# -- Linux
+ifeq ($(OS),LINUX)
+CXX           = g++
+CXXFLAGS      = -g -Wall -fPIC
+PMUSRPATH     = ./include
+MNPATH        = $(ROOTSYS)/include
+INCLUDES      = -I $(PMUSRPATH) -I $(MNPATH)
+LD            = g++
+LDFLAGS       = -g
+endif
+
+# -- Darwin
+ifeq ($(OS),DARWIN)
+CXX           = g++
+CXXFLAGS      = -g -Wall -fPIC
+INCLUDES      = -I./
+LD            = g++
+LDFLAGS       = -g
+endif
+
+# the output from the root-config script:
+CXXFLAGS      += $(ROOTCFLAGS)
+LDFLAGS       +=
+
+# the ROOT libraries (G = graphic)
+LIBS          = $(ROOTLIBS) -lfftw3
+GLIBS         = $(ROOTGLIBS) -lfftw3
+
+INSTALLPATH = $(HOME)/analysis/bin
+
+EXEC = fourierPhaseTest
+
+# some definitions: headers, sources, objects,...
+OBJS =
+OBJS += $(EXEC).o
+
+# make the executable:
+#
+all:  $(EXEC)
+
+$(EXEC): $(OBJS)
+	@echo "---> Building $(EXEC) ..."
+	/bin/rm -f $(SHLIB)
+	$(LD) $(OBJS) -o $(EXEC) $(GLIBS)
+	@echo "done"
+
+# clean up: remove all object file (and core files)
+# semicolon needed to tell make there is no source
+# for this target!
+#
+clean:; 	@rm -f $(OBJS)
+	@echo "---> removing $(OBJS)"
+
+#
+$(OBJS): %.o: %.cpp
+	$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
+
+install: all
+	@echo "copy $(EXEC) to $(INSTALLPATH)"
+	cp -p $(EXEC) $(INSTALLPATH)
+

src/tests/fourier/PMusrFourier.cpp

@@ -107,6 +107,10 @@ PMusrFourier::PMusrFourier(int dataType, PDoubleVector &data, double timeResolut
         rebinCounter = 0;
       }
     }
+  } else { // no rebinning, just copy data
+    for (unsigned int i=0; i<fData.size(); i++) {
+      fDataRebinned.push_back(fData[i]);
+    }
   }
 
 cout << endl << "dB = " << 1.0/(2.0 * F_GAMMA_BAR_MUON * (fEndTime-fStartTime)) << " (G), Bmax = " << 1.0/(2.0 * F_GAMMA_BAR_MUON * fTimeResolution) << " (G)" << endl;
@@ -192,54 +196,50 @@ void PMusrFourier::Transform(unsigned int apodizationTag, unsigned int filterTag
     PrepareFFTwInputData(apodizationTag, filterTag);
   }
 
-// for test only
+  fftw_execute(fFFTwPlan);
-// keep data
-fftw_complex *data; 
-data = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins);
-for (unsigned int i=0; i<fNoOfBins; i++) {
-  data[i][0] = fIn[i][0];
-  data[i][1] = 0.0;
-}
 
+// for test only
 // loop over the phase
 double sum;
-TH1F sumHist("sumHist", "sumHist", 361, -180.5, 180.5);
+double corr_phase;
+double min, min_phase;;
+TH1F sumHist("sumHist", "sumHist", 181, -90.5, 90.5);
 double dB   = 1.0/(2.0 * F_GAMMA_BAR_MUON * (fEndTime-fStartTime));
 double Bmax = 1.0/(2.0 * F_GAMMA_BAR_MUON * fTimeResolution);
 TH1F re("re", "re", fNoOfBins/2+1, -dB/2.0, Bmax+dB/2.0);
 TH1F im("im", "im", fNoOfBins/2+1, -dB/2.0, Bmax+dB/2.0);
-for (int p=-180; p<180; p++) {
+for (int p=-90; p<90; p++) {
-  for (unsigned int i=0; i<fNoOfBins; i++) {
+  // calculate sum of the rotated imaginary part of fOut
-    // recalculate fIn including the phase
-    fIn[i][0] = data[i][0]*cos(p/180.0*PI);
-    fIn[i][1] = data[i][0]*sin(p/180.0*PI);
-  }
-  fftw_execute(fFFTwPlan);
-
-  if (p == 7) {
-    for (unsigned int j=0; j<fNoOfBins/2; j++) {
-      re.SetBinContent(j+1, fOut[j][0]);
-      im.SetBinContent(j+1, fOut[j][1]);
-    }
-  }
-
-  // calculate sum of the imaginary part of fOut
   sum = 0.0;
+  corr_phase = (double)p/180.0*PI;
   for (unsigned int i=0; i<fNoOfBins/2; i++) {
-    sum += fOut[i][1];
+    sum += fOut[i][0]*sin(corr_phase) + fOut[i][1]*cos(corr_phase);
   }
-  if (p == 7) {
+cout << endl << "-> sum = " << fabs(sum) << ", min = " << min << ", min_phase = " << min_phase/PI*180.0;
-    cout << endl << ">> sum = " << sum << endl;
+  if (p==-90) {
+    min = fabs(sum);
+    min_phase = corr_phase;
+cout << endl << "!> min = " << min << ", min_phase = " << min_phase/PI*180.0;
   }
-  sumHist.SetBinContent(p+181, fabs(sum));
+  if (fabs(sum) < min) {
+    min = fabs(sum);
+    min_phase = corr_phase;
+cout << endl << "-> min = " << min << ", min_phase = " << min_phase/PI*180.0;
   }
+  sumHist.SetBinContent(p+91, fabs(sum));
+}
+cout << endl << ">> min = " << min << ", min_phase = " << min_phase/PI*180.0;
+
+for (unsigned int i=0; i<fNoOfBins/2; i++) {
+  re.SetBinContent(i+1, fOut[i][0]*cos(min_phase) - fOut[i][1]*sin(min_phase));
+  im.SetBinContent(i+1, fOut[i][0]*sin(min_phase) + fOut[i][1]*cos(min_phase));
+}
+
 TFile f("test_out.root", "RECREATE");
 re.Write();
 im.Write();
 sumHist.Write();
 f.Close();
-
-fftw_free(data);
 }
 
 //--------------------------------------------------------------------------

src/tests/fourier/fourierPhaseTest.cpp

@@ -0,0 +1,153 @@
+/***************************************************************************
+
+  fourierPhaseTest.cpp
+
+  Author: Andreas Suter
+  e-mail: andreas.suter@psi.ch
+
+  $Id$
+
+***************************************************************************/
+
+#include <iostream>
+using namespace std;
+
+#include "fftw3.h"
+
+#include <TFile.h>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TGraph.h>
+#include <TColor.h>
+
+int main(int argc, char *argv[])
+{
+  if (argc != 2) {
+    cout << endl << "usage: fourierPhaseTest <root-output-file-name>";
+    cout << endl;
+    return -1;
+  }
+
+  UInt_t   noOfBins = (UInt_t)pow(2, 11);
+  Double_t freq = TMath::TwoPi() * 10.0 / (Double_t)noOfBins;
+  Double_t phase = 23.7 * TMath::Pi() / 180.0;
+  cout << endl << ">> noOfBins = " << noOfBins;
+  cout << endl << ">> freq     = " << freq;
+  cout << endl << ">> phase    = " << phase;
+
+  fftw_plan fftwPlan;
+  fftw_complex *in;
+  fftw_complex *out;
+
+  // allocate necessary memory
+  in  = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*noOfBins);
+  out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*noOfBins);
+
+  // generate plan
+  fftwPlan = fftw_plan_dft_1d(noOfBins, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
+  cout << endl << ">> fftwPlan = " << fftwPlan;
+  if (!fftwPlan) {
+    cout << endl << "**ERROR** Problems to establish a FFTW plan ...";
+    cout << endl;
+    return -1;
+  }
+
+  // fill data
+  cout << endl << ">> fill data ...";
+  for (UInt_t i=0; i<noOfBins; i++) {
+    in[i][0] = cos(freq*i + phase)*exp(-(Double_t)i/(Double_t)noOfBins/0.1);
+    in[i][1] = 0.0;
+  }
+
+  TCanvas *canv_in = new TCanvas("canv_in", "phase test ..");
+
+  TGraph *gin = new TGraph(noOfBins);
+  for (UInt_t i=0; i<noOfBins; i++) {
+    gin->SetPoint(i, i, in[i][0]);
+  }
+  gin->Draw("ap");
+
+  fftw_execute(fftwPlan);
+
+  // find peak
+  Double_t dval, peak = 0.0;
+  UInt_t peak_pos = 0;
+  for (UInt_t i=0; i<noOfBins/2; i++) {
+    dval = out[i][0]*out[i][0]+out[i][1]*out[i][1];
+    if (peak < dval) {
+      peak = dval;
+      peak_pos = i;
+    }
+  }
+  cout << endl << ">> peak at peak_pos = " << peak_pos;
+  if (2*peak_pos > noOfBins/2)
+    peak_pos = noOfBins/4;
+
+  TCanvas *canv_out = new TCanvas("canv_out", "phase test ..");
+
+  TGraph *goutRe = new TGraph(noOfBins/2);
+  TGraph *goutIm = new TGraph(noOfBins/2);
+
+  // estimate phase
+  Double_t sum;
+  Double_t corr_phase;
+  Double_t min, min_phase;
+  for (Int_t k=-90; k<90; k++) {
+    sum = 0.0;
+    corr_phase = (Double_t)k/180.0*TMath::Pi();
+    for (UInt_t i=0; i<2*peak_pos; i++) {
+      sum += out[i][0]*sin(corr_phase) + out[i][1]*cos(corr_phase);
+    }
+    if (k == -90) {
+      min = fabs(sum);
+      min_phase = 0.0;
+    }
+    if (fabs(sum) < fabs(min)) {
+      min = fabs(sum);
+      min_phase = (Double_t)k;
+    }
+    cout << endl << ">> phase = " << k << ", sum = " << fabs(sum);
+  }
+
+  cout << endl << ">> min = " << min << ", min phase = " << min_phase;
+  cout << endl;
+
+  min_phase = -23.7*TMath::Pi()/180.0;
+  for (UInt_t i=0; i<noOfBins/2; i++) {
+    goutRe->SetPoint(i, i, out[i][0]*cos(min_phase) - out[i][1]*sin(min_phase));
+    goutIm->SetPoint(i, i, out[i][0]*sin(min_phase) + out[i][1]*cos(min_phase));
+  }
+  goutRe->SetMarkerStyle(20); // bullet
+  goutRe->SetMarkerColor(TColor::GetColor(0,0,0)); // black
+  goutRe->SetLineColor(TColor::GetColor(0,0,0));   // black
+  goutRe->SetFillColor(TColor::GetColor(255,255,255)); // white
+  goutRe->SetTitle("Re");
+  goutIm->SetMarkerStyle(29); // filled star
+  goutIm->SetMarkerColor(TColor::GetColor(255,0,0)); // red
+  goutIm->SetLineColor(TColor::GetColor(255,0,0));   // red
+  goutIm->SetFillColor(TColor::GetColor(255,255,255)); // white
+  goutIm->SetTitle("Im");
+  goutRe->Draw("apl");
+  goutIm->Draw("plsame");
+
+  // write data to file
+  TFile f(argv[1], "RECREATE");
+  canv_in->Write();
+  canv_out->Write();
+  f.Close();
+
+  // clean up
+  delete canv_in;
+  delete gin;
+  delete canv_out;
+  delete goutRe;
+  delete goutIm;
+
+  fftw_destroy_plan(fftwPlan);
+  fftw_free(in);
+  fftw_free(out);
+
+  cout << endl;
+
+  return 0;
+}