Go4-Fitpackagev4.04-2: f_find_peaks.c Source File

00001 // $Id: f_find_peaks.c 478 2009-10-29 12:26:09Z linev $
00002 //-----------------------------------------------------------------------
00003 //       The GSI Online Offline Object Oriented (Go4) Project
00004 //         Experiment Data Processing at EE department, GSI
00005 //-----------------------------------------------------------------------
00006 // Copyright (C) 2000- GSI Helmholtzzentrum für Schwerionenforschung GmbH
00007 //                     Planckstr. 1, 64291 Darmstadt, Germany
00008 // Contact:            http://go4.gsi.de
00009 //-----------------------------------------------------------------------
00010 // This software can be used under the license agreements as stated
00011 // in Go4License.txt file which is part of the distribution.
00012 //-----------------------------------------------------------------------
00013 
00014 #include <stdlib.h>
00015 
00016 #define TYPE__SHORT  0
00017 #define TYPE__INT    1
00018 #define TYPE__FLOAT  2
00019 #define TYPE__DOUBLE 3
00020 
00021 #define UP 1
00022 #define DOWN -1
00023 #define HORI 0
00024 #define PRINT 0
00025 
00026 int go4fit_position(int l_len,double *pa_data,double *pr_pos,double *pr_sig, double *pr_sum);
00027 
00028 void go4fit_find_peaks(
00029         void   *pfData,     // pointer to data
00030         int     lType,      // data type: 0=w 1=l 2=f 3=d
00031         int     lFirstChan, // first channel
00032         int     lPoints,    // number of points
00033         int     lAver,      // channels to sum up
00034         double  dDeltaFact, // noise factor
00035         double  dDeltaMin,  // noise minimum
00036         double *dNoise,     // noise array
00037         int     lPeaksMax,  // Maximum number of peaks
00038         int    *plPeaks,    // returns number of peaks found
00039         int    *plMinima,   // pointer to list of minima
00040         double *pdMinima,   // pointer to list of minima values
00041         int    *plMean,   // pointer to list of maxima
00042         int    *plWidth,    // pointer to list of width
00043         double *pdIntegral, // pointer to list of integrals
00044         int    *plLeft,     // pointer to array of left RMS index
00045         int    *plMaximum,     // pointer to array of maximum index
00046         int    *plRight,     // pointer to array of right RMS index
00047         float  *pfLowNoise,   // pointer to data array of low noise
00048         float  *pfHighNoise)  // pointer to data array of high noise
00049 {
00050   short *pw;
00051   int lMaxIndex,lMinIndex,lIndexMax,lDirection,lIndex,lLastUp,lLastDown,ii,i,lMax,lLeft,lRight,
00052     *plMinimaL,*plMeanL,*plMinimaR,*plMeanR,lPrint,*pl;
00053   float *pf;
00054   double dValue,dLow,dHigh,dDelta,dAver,*pdData,*pdNoise,dMax,dMin,dPos,dSig,dSum,*pd;
00055 
00056   *plPeaks=0;
00057   pdNoise=NULL;
00058   pdData = (double *)malloc(lPoints*8);
00059   pdNoise = (double *)malloc(lPoints*8);
00060   memset(pdData,0,lPoints*8);
00061   memset(pdNoise,0,lPoints*8);
00062   lPrint=PRINT;
00063   dMax=0.;
00064   dMin=1e20;
00065   pw = ((short  *)pfData) + lFirstChan;
00066   pl = ((int    *)pfData) + lFirstChan;
00067   pf = ((float  *)pfData) + lFirstChan;
00068   pd = ((double *)pfData) + lFirstChan;
00069   if(lAver == 0)lAver=1;
00070   lPoints=lPoints/lAver;
00071   dAver = (double) lAver;
00072   // copy the data field into local double field
00073   switch(lType)
00074     {
00075     case TYPE__SHORT:  for(ii=0;ii<lPoints;ii++)
00076       {  *(pdData+ii) = (double)*(pw+lAver*ii);
00077       for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pw+lAver*ii+i);
00078       *(pdData+ii) = *(pdData+ii)/dAver;
00079       if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
00080       if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
00081     case TYPE__INT:    for(ii=0;ii<lPoints;ii++)
00082       {  *(pdData+ii) = (double)*(pl+lAver*ii);
00083       for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pl+lAver*ii+i);
00084       *(pdData+ii) = *(pdData+ii)/dAver;
00085       if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
00086       if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
00087     case TYPE__FLOAT:  for(ii=0;ii<lPoints;ii++)
00088       {  *(pdData+ii) = (double)*(pf+lAver*ii);
00089       for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pf+lAver*ii+i);
00090       *(pdData+ii) = *(pdData+ii)/dAver;
00091       if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
00092       if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
00093     case TYPE__DOUBLE: for(ii=0;ii<lPoints;ii++)
00094       {  *(pdData+ii) = (double)*(pd+lAver*ii);
00095       for(i=1;i<lAver;i++)*(pdData+ii) += (double)*(pd+lAver*ii+i);
00096       *(pdData+ii) = *(pdData+ii)/dAver;
00097       if(*(pdData+ii) > dMax) dMax=*(pdData+ii);
00098       if(*(pdData+ii) < dMin) dMin=*(pdData+ii);} break;
00099     default: printf("Data type %d not supported!\n",lType);
00100       printf("Type 0: short, 1: int, 2: float, 3: double\n");
00101       free(pdData);
00102       free(pdNoise);
00103       return;
00104     }
00105   if(dNoise == NULL)
00106     {// calculate from square root
00107       for(ii=0;ii<lPoints;ii++)
00108    {
00109      if(*(pdData+ii) == 0.0) *(pdNoise+ii) = dDeltaFact;
00110      else                    *(pdNoise+ii) = sqrt(*(pdData+ii)) * dDeltaFact;
00111      if(dDeltaMin > *(pdNoise+ii)) *(pdNoise+ii)=dDeltaMin;
00112    }
00113     }else {// sum up bins from noise input array
00114    for(ii=0;ii<lPoints;ii++)
00115      {  *(pdNoise+ii) = *(dNoise+lFirstChan+lAver*ii);
00116         for(i=1;i<lAver;i++)*(pdNoise+ii) += *(dNoise+lFirstChan+lAver*ii+i);
00117         *(pdNoise+ii) = *(pdNoise+ii)/dAver;
00118      }
00119       }
00120   lMaxIndex=0;
00121   lMinIndex=0;
00122   lIndexMax=lPoints-1;
00123   plMinimaL = (int *)malloc(lPeaksMax*16);
00124   plMeanL = (int *)(plMinimaL+lPeaksMax);
00125   plMinimaR = (int *)(plMinimaL+2*lPeaksMax);
00126   plMeanR = (int *)(plMinimaL+3*lPeaksMax);
00127   memset(plMinimaL,0,lPeaksMax*16);
00128 
00129   dDelta= *pdNoise;
00130   dHigh = *pdData + dDelta;
00131   dLow  = *pdData - dDelta;
00132   if(pfLowNoise)*pfLowNoise=dLow;
00133   if(pfHighNoise)*pfHighNoise=dHigh;
00134   //-----------------------------------------------------
00135   // loop: get points of minima and maxima
00136   lDirection=HORI;
00137   lIndex=0;
00138   lLastDown=0;
00139   lLastUp=0;
00140   while ((lIndex <= lIndexMax) & (lMaxIndex < lPeaksMax-2))
00141     {
00142       dValue = *(pdData+lIndex);
00143       if(dValue > dHigh)
00144    {// direction goes upstairs
00145      if (lDirection != UP)
00146        { // direction was downstairs: minimum found
00147          if(lPrint)printf("%5d - %5d - %5d d %8.1f\n",
00148                 lMinIndex,
00149                 lLastDown,
00150                 lIndex,
00151                 dDelta);
00152          *(plMinimaR+lMinIndex)=lIndex;
00153          *(plMinimaL+lMinIndex)=lLastDown;
00154          lMinIndex++;
00155        }
00156      lLastUp=lIndex;
00157      lDirection=UP;// direction goes upstairs
00158    }
00159       if(dValue < dLow)
00160    {// direction goes downstairs
00161      if (lDirection == UP)
00162        { // direction was upstairs: maximum found
00163          if(lPrint)printf("%5d + %5d - %5d d %8.1f\n",
00164                 lMaxIndex,
00165                 lLastUp,
00166                 lIndex,
00167                 dDelta);
00168          *(plMeanR+lMaxIndex)=lIndex;
00169          *(plMeanL+lMaxIndex)=lLastUp;
00170          lMaxIndex++;
00171        }
00172      lLastDown=lIndex;
00173      lDirection=DOWN;// direction goes downstairs
00174    }
00175       // calculate new limits if noise band has been left
00176       if((dValue < dLow)|(dValue > dHigh))
00177    {
00178           dDelta= *(pdNoise+lIndex);
00179      dHigh = dValue + dDelta;
00180      dLow  = dValue - dDelta;
00181    }
00182       lIndex++;
00183       if(pfLowNoise)*(pfLowNoise+lIndex)=dLow;
00184       if(pfHighNoise)*(pfHighNoise+lIndex)=dHigh;
00185     }
00186   if(lDirection == DOWN)
00187     {
00188       *(plMinimaL+lMinIndex)=lLastDown;
00189       *(plMinimaR+lMinIndex)=lIndexMax;
00190       if(lPrint)printf("%5d - %5d - %5d\n",
00191              lMinIndex,
00192              lLastDown,
00193              lIndexMax);
00194       lMinIndex++;
00195     }
00196   //-----------------------------------------------------
00197   if(lPrint)printf("-------------------------\n");
00198   if(lPrint)printf("Minima %d Maxima %d  \n",lMinIndex,lMaxIndex);
00199   // for each peak there must be a minimum left and right
00200   if(lMinIndex != (lMaxIndex+1))
00201     {
00202       free(plMinimaL);
00203       free(pdData);
00204       free(pdNoise);
00205       *plPeaks=0;
00206       if(lPrint)printf("Error, wrong minima\n");
00207       return;
00208     }
00209   // calculate mean values of minima L&R, mean content, peaks
00210   // start with first minimum
00211   if(lPrint)printf("  [ Left   Right ]  \n");
00212   *plMinima = (*plMinimaL + *plMinimaR)/2;
00213   *pdMinima=0.; // average content at mean minimum channel
00214   for(ii = *plMinimaL;ii <= *plMinimaR;ii++) *pdMinima += *(pdData+ii);
00215   *pdMinima = *pdMinima/(double)(*plMinimaR - *plMinimaL + 1);
00216   if(lPrint)printf("- [%6d %6d] %6d %6.1f\n",
00217          *plMinimaL,
00218          *plMinimaR,
00219          *plMinima,
00220          *pdMinima);
00221   for(lIndex=0;lIndex < lMaxIndex;lIndex++)
00222     {
00223       // simple mean peak value, overwritten later
00224       *(plMean+lIndex)=(*(plMeanL+lIndex)+*(plMeanR+lIndex))/2;
00225       if(lPrint)printf("+ [%6d %6d] %6d\n",
00226              *(plMeanL+lIndex),
00227              *(plMeanR+lIndex),
00228              *(plMean+lIndex));
00229       // mean minima values (channel and content)
00230       *(plMinima+lIndex+1) = (*(plMinimaL+lIndex+1) + *(plMinimaR+lIndex+1))/2;
00231       *(pdMinima+lIndex+1)=0.;
00232       for(ii = *(plMinimaL+lIndex+1);ii <= *(plMinimaR+lIndex+1);ii++)
00233    *(pdMinima+lIndex+1) += *(pdData+ii);
00234       *(pdMinima+lIndex+1) = *(pdMinima+lIndex+1)/(double)(*(plMinimaR+lIndex+1) - *(plMinimaL+lIndex+1)+1);
00235       if(lPrint)printf("- [%6d %6d] %6d %6.1f\n",
00236              *(plMinimaL+lIndex+1),
00237              *(plMinimaR+lIndex+1),
00238              *(plMinima+lIndex+1),
00239              *(pdMinima+lIndex+1));
00240       // Peak positions, width and sum by momenta
00241       ii = *(plMinimaL+lIndex+1) - *(plMinimaR+lIndex) + 1;
00242       if(ii > 1)
00243    {
00244      go4fit_position(ii,pdData + *(plMinimaR+lIndex),&dPos,&dSig,&dSum);
00245 
00246      *(plMean+*plPeaks)=(int)((dPos+(double)(*(plMinimaR+lIndex)))*dAver)+lFirstChan;
00247      *(plWidth+*plPeaks) =(int)(dSig*dAver + 0.5);
00248      if(*(pdMinima+lIndex) < *(pdMinima+lIndex+1))
00249        dSum=dSum-((double)ii * *(pdMinima+lIndex)+ (double)ii * (*(pdMinima+lIndex+1)-*(pdMinima+lIndex))/2.);
00250           else
00251        dSum=dSum-((double)ii * *(pdMinima+lIndex)- (double)ii * (*(pdMinima+lIndex)-*(pdMinima+lIndex+1))/2.);
00252      *(pdIntegral+*plPeaks) = dSum*dAver;
00253           (*plPeaks)++;
00254    }
00255       else
00256    {
00257      free(plMinimaL);
00258      free(pdData);
00259       free(pdNoise);
00260       *plPeaks=0;
00261      if(lPrint)printf("Error, negative interval\n");
00262      return;
00263    }
00264     }
00265   for(i=0;i<*plPeaks;i++)
00266     {
00267       dMax=0; // search local maximum in interval
00268       dMin=1000000000; // local minimum
00269       pd=pdData+*(plMinima+i);
00270       for(ii=*(plMinima+i);ii<*(plMinima+i+1);ii++){if(dMax < *pd)dMax = *pd; if(dMin > *pd)dMin = *pd; pd++;}
00271       dMin=dMin+(dMax-dMin)/2.; // threashold for RMS
00272       pd=pdData+*(plMinima+i); // search channel of local maximum
00273       for(lMax=*(plMinima+i);lMax<*(plMinima+i+1);lMax++){if(dMax == *pd)break; pd++;} // l_max is now channel of local maximum
00274       lLeft=*(plMinima+i);
00275       lRight=*(plMinima+i+1);
00276       *(plMaximum+i)=lMax;
00277       pd=pdData+lMax;
00278       for(ii=lMax;ii<lRight;ii++){if(dMin > *pd){lRight=ii;break;}pd++;}
00279       pd=pdData+lMax;
00280       for(ii=lMax;ii>lLeft;ii--){if(dMin > *pd){lLeft=ii;break;}pd--;}
00281       *(plLeft+i)=lLeft;
00282       *(plRight+i)=lRight;
00283     }
00284   *(plMinima) = (int)((double)(*(plMinima))*dAver + dAver/2.)+lFirstChan;
00285   for(i=0;i<*plPeaks;i++)
00286     {
00287       // calculate channel numbers of original data
00288       *(plMaximum+i)  = (int)((double)(*(plMaximum+i))*dAver + dAver/2.)+lFirstChan;
00289       *(plLeft+i)     = (int)((double)(*(plLeft+i))*dAver + dAver/2.)+lFirstChan;
00290       *(plRight+i)    = (int)((double)(*(plRight+i))*dAver + dAver/2.)+lFirstChan;
00291       *(plMinima+i+1) = (int)((double)(*(plMinima+i+1))*dAver + dAver/2.)+lFirstChan;
00292     }
00293   if(lPrint)printf("peaks found %d\n",*plPeaks);
00294 
00295   free(plMinimaL);
00296   free(pdData);
00297   free(pdNoise);
00298 }
00299 //**********************************************************************************
00300 int go4fit_position(int l_len,double *pa_data,double *pr_pos,double *pr_sig, double *pr_sum)
00301 {
00302 
00303 #define LOW  2
00304 #define WIDTH  8
00305 
00306   double r_sig_f     ;
00307 //SL  int I,J,K,L     ;
00308   int J;
00309   double d_sum_prod  ;
00310 //SL  double l_max_chan  ;
00311   double *pl_data,d_max    ;
00312 
00313   /* 2.0E2*SQRT(2.0E0*LOG(2.0E0))/100. */;
00314   r_sig_f      = 55.4518;
00315   r_sig_f      = 2.;
00316   *pr_sum   = 0;
00317   d_sum_prod   = 0;
00318   *pr_sig      = 0.;
00319   *pr_pos      = 0.;
00320   d_max=0.0;
00321 
00322   /* get maximum channel content and integral */
00323   /*  pl_data = pa_data;
00324       for(J = 0; J < l_len; J++)
00325       {
00326       if(d_max < *pl_data)
00327       {
00328       d_max = *pl_data;
00329       l_max_chan = J;
00330       }
00331       pl_data++;
00332       }
00333   */
00334   /* Calculate first momentum */
00335   pl_data = pa_data;
00336   for(J = 1; J <= l_len; J++)
00337     {
00338       *pr_sum = *pr_sum +     *pl_data;
00339       d_sum_prod = d_sum_prod + J * *pl_data;
00340       pl_data++;
00341     }
00342   /* Calculate second momentum */
00343   if(*pr_sum > 0)
00344     {
00345       *pr_pos = (d_sum_prod/ *pr_sum + 0.5);
00346       pl_data = pa_data;
00347       for(J = 1; J <= l_len; J++)
00348    {
00349      *pr_sig = *pr_sig + ((double)J - *pr_pos) * ((double)J - *pr_pos) * *pl_data;
00350      pl_data++;
00351    }
00352       *pr_sig = r_sig_f * sqrt(*pr_sig/ *pr_sum);
00353     }
00354   else return(1);
00355   *pr_pos = *pr_pos -1.0;
00356   return(0);
00357 }