Go4Fit/f_find_peaks.c
Go to the documentation of this file.
00001 // $Id: f_find_peaks.c 1037 2013-11-06 13:39:24Z 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 int J; 00308 double d_sum_prod ; 00309 double *pl_data; 00310 00311 /* 2.0E2*SQRT(2.0E0*LOG(2.0E0))/100. */; 00312 r_sig_f = 55.4518; 00313 r_sig_f = 2.; 00314 *pr_sum = 0; 00315 d_sum_prod = 0; 00316 *pr_sig = 0.; 00317 *pr_pos = 0.; 00318 00319 /* Calculate first momentum */ 00320 pl_data = pa_data; 00321 for(J = 1; J <= l_len; J++) 00322 { 00323 *pr_sum = *pr_sum + *pl_data; 00324 d_sum_prod = d_sum_prod + J * *pl_data; 00325 pl_data++; 00326 } 00327 /* Calculate second momentum */ 00328 if(*pr_sum > 0) 00329 { 00330 *pr_pos = (d_sum_prod/ *pr_sum + 0.5); 00331 pl_data = pa_data; 00332 for(J = 1; J <= l_len; J++) 00333 { 00334 *pr_sig = *pr_sig + ((double)J - *pr_pos) * ((double)J - *pr_pos) * *pl_data; 00335 pl_data++; 00336 } 00337 *pr_sig = r_sig_f * sqrt(*pr_sig/ *pr_sum); 00338 } 00339 else return(1); 00340 *pr_pos = *pr_pos -1.0; 00341 return(0); 00342 }
