Batch Automatically find peaks
Ander
I'm probably being really stupid here, but I have over 1200 curves oscillating over a few cycles that are all collected on the same timebase. I want to identify all of the peaks in the data as a batch - both X (time of peak) and Y (amplitude of peak).
I've been able to get nice peak recognition using the Automatically Detect Peaks package, but this only allows me to identify peaks in one curve at a time. Ideally, I want to automate this and output the X and Y coordinates to separate tables with each column representing each individual curve. Does anyone have any ideas how to automate this?
Many Thanks,
Ander
It is handy that Igor is completely programmable, so you can automate the procedure yourself, possibly with some help here.
To get that help, you'd need to upload an example of your data, the code you've already got (if you have any), and a clear description of what you want.
That "clear description" is the hard part; it is easy to believe that what one writes is necessarily exactly what the reader perceives.
--Jim Prouty
Software Engineer, WaveMetrics, Inc.
February 22, 2017 at 07:28 pm - Permalink
* Write the code you need to load a data file.
* Write the code that does the work required to get the answers from that one data file.
* Put the two pieces of code in to a function with a loop over all the files that you want to process.
Is this the direction of the hints that you are after?
--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH
February 22, 2017 at 07:34 pm - Permalink
February 23, 2017 at 08:12 am - Permalink
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
February 23, 2017 at 12:44 pm - Permalink
I once wrote a wrapper function for WM's PeakFind. Probably not the nicest solution but take a look at PeakFind() in https://github.com/ukos-git/igor-common-utilities/blob/master/app/utili…
Function/Wave PeakFind(wavInput, [wvXdata, noiselevel, smoothingFactor, minPeakPercent, maxPeaks, verbose])
Wave wavInput, wvXdata
Variable noiselevel, smoothingFactor, minPeakPercent, maxPeaks
variable verbose
if (ParamIsDefault(verbose))
verbose = 0
endif
Variable pBegin, pEnd
Variable/C estimates
Variable peaksFound
String newName
if(ParamIsDefault(wvXdata))
Make/N=(DimSize(wavInput, 0))/Free wvXdata = DimOffset(wavInput, 0) + DimDelta(wavInput, 0) * p
endif
Wave wavYdata = wavInput
Make/FREE/N=(maxPeaks,7) wavOutput
// label columns of wave for readability
SetDimLabel 1, 0, wavelength, wavOutput
SetDimLabel 1, 1, height, wavOutput
SetDimLabel 1, 2, width, wavOutput
SetDimLabel 1, 3, positionY, wavOutput
SetDimLabel 1, 4, positionX, wavOutput
SetDimLabel 1, 5, widthL, wavOutput // from WM: somewhat free width left and right
SetDimLabel 1, 6, widthR, wavOutput
// input parameters for WM's AutoFindPeaksNew
pBegin = 0
pEnd = DimSize(wavYdata, 0) - 1
if(ParamIsDefault(maxPeaks))
maxPeaks = 10
endif
if(ParamIsDefault(minPeakPercent))
minPeakPercent = 5
endif
if(ParamIsDefault(noiselevel))
noiselevel = 1
endif
estimates = EstPeakNoiseAndSmfact(wavYdata, pBegin, pEnd)
noiselevel *= real(estimates)
if(ParamIsDefault(smoothingFactor))
smoothingFactor = imag(estimates)
endif
if(!(noiselevel>0))
noiselevel = 0.01
endif
peaksFound = AutoFindPeaksNew(wavYdata, pBegin, pEnd, noiseLevel, smoothingFactor, maxPeaks)
WAVE W_AutoPeakInfo // output of AutoFindPeaksNew
if(verbose > 2)
print "== peakFind =="
printf "peaks: \t %d\r", peaksFound
printf "smooth: \t %.1f\r", smoothingFactor
printf "noise: \t %.4f\r", noiseLevel
endif
// Remove too-small peaks
if(peaksFound > 0)
peaksFound = TrimAmpAutoPeakInfo(W_AutoPeakInfo,minPeakPercent/100)
endif
// process peaks
if(peaksFound > 0)
// Redimension to number of peaks
Redimension/N=(peaksFound, -1) wavOutput
// save peak positions in input wave
wavOutput[][%positionX] = W_AutoPeakInfo[p][0]
wavOutput[][%positionY] = wavYdata[wavOutput[p][%positionX]]
// The x values in W_AutoPeakInfo are still actually points, not X
AdjustAutoPeakInfoForX(W_AutoPeakInfo, wavYdata, wvXdata)
wavOutput[][%wavelength] = W_AutoPeakInfo[p][0]
// save all data from WM procedure
wavOutput[][%width] = W_AutoPeakInfo[p][1]
wavOutput[][%height] = W_AutoPeakInfo[p][2]
wavOutput[][%widthL] = W_AutoPeakInfo[p][3]
wavOutput[][%widthR] = W_AutoPeakInfo[p][4]
endif
return wavOutput
End
You could use this in a for loop over all of your waves. Just an example for fitting one wave:
variable center
Make/FREE myPeak
setscale x, -10, 10, myPeak
myPeak = gauss(x,center,1) + gnoise(0.02)
WAVE peakresults = peakfind(myPeak)
print peakresults[0][%wavelength], peakresults[0][%height]
End
Automatically Find Peaks is using the second derivative for peakfind. You should do a fit in addition after the "initial values" were found using peakfind. If you consider a PeakFit take a look at https://github.com/ukos-git/igor-common-utilities/blob/master/app/utili…. I implemented the Gauss Fit from the MultiPeak Fit package.
February 25, 2017 at 08:14 am - Permalink