/Go4Fit/f_find_peaks.c Source File

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