Batch Automatically find peaks

Share icon

Hi!

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
Everyone wants to do this differently, in my experience, so no one solution has presented itself.

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
I might suggest these steps ...

* 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
I haven't used the Automatically Detect Peaks package, but the Multipeak Fitting 2 package under the Analysis-Packages menu may be able to do what you need. First check if it can do a good job without user intervention on a single data set just by trying "auto-locate peaks now" followed by "OK" and "do fit". If that result looks reasonable, take a look in the package's help where there's a batch fitting procedure near the end called MultiPeakBatchFit(datafoldername,ywavelist,xwavelist). It'll probably require some modification for your application.

DisplayHelpTopic("Multi-Peak Fitting for Igor Programmers")

February 23, 2017 at 08:12 am - Permalink
Hi.
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…

#include <Peak AutoFind>

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:

Function testfit(center)
	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