Peak Finding

Finding peaks in acquired data is a bit of an art that varies depending on the type of data.

It is helpful to remove any baseline (or "trend") and possibly to smooth the data to simplify finding peaks.

For a series of only positive-going or only negative-going peaks:

spectra of positive-going peaks

You can use the PeakAreasUnipolar.ipf procedure file (part of Technical Note #020-C Unipolar Peak Areas) to automatically locate peaks:

Init Identify Peaks dialog showing peak waves and threshold values positive peaks and markers identifying peak start, maxima, and end locations

For peaks that alternate between positive and negative values:

alternating peaks (respiration data)

you can use the Peak Areas.ipf procedure file (part of Technical Note #020-B Peak Areas) to automatically find peaks, construct a baseline and compute the areas:

If you're comfortable with programming, you can use Igor's built-in FindPeak operation repeatedly in a procedure to locate peak starts and ends:

Wave w	// peak data
Make/O/N=(maxPeaks) peakPositionsX= NaN, peakPositionsY= NaN    
Variable peaksFound=0
Variable startP=0
Variable endP= DimSize(w,0)-1
do
    FindPeak/B=(smoothing)/I/M=(threshold)/P/Q/R=[startP,endP] w
    // FindPeak outputs are V_Flag, V_PeakLoc, V_LeadingEdgeLoc,
    // V_TrailingEdgeLoc, V_PeakVal, and V_PeakWidth. 
    if( V_Flag != 0 )
        break
    endif
    
    peakPositionsX[peaksFound]=pnt2x(w,V_PeakLoc)
    peakPositionsY[peaksFound]=V_PeakVal
    peaksFound += 1
    
    startP= V_TrailingEdgeLoc+1
while( peaksFound < maxPeaks )

if( peaksFound )
    Redimension/N=(peaksFound) peakPositionsX, peakPositionsY

    DoWindow/F ShowPeaks
    if(V_Flag == 0 )
        Display/N=ShowPeaks w
        AppendToGraph/W=ShowPeaks peakPositionsY vs peakPositionsX
        ModifyGraph/W=ShowPeaks mode(peakPositionsY)=3,marker(peakPositionsY)=19, rgb(peakPositionsY)=(0,0,65535)
    endif
else
    DoAlert 0, "No peaks found using threshold= "+num2str(threshold)
    KillWaves/Z peakPositionsX, peakPositionsY
endif

The Multi-Peak Fitting package uses the Peak AutoFind.ipf procedure, which uses the FindLevel operation to search smoothed derivatives of the data. To use the procedure, type:

#include <Peak AutoFind>

in your Procedure window, or choose "Multi-peak Fitting" from the Analysis->Packages submenu. Select the Peak AutoFind procedure from the Other Windows menu to see routines, and select "Multi-Peak Fitting" from the Other Windows menu to see how the routines are used.

Here is a procedure that uses the Peak AutoFind procedure file:

#pragma rtglobals=1
#include <Peak AutoFind>

Menu "Analysis"
    "Automatically Find Peaks", /Q, MyAutomaticallyFindPeaks()
EndMacro

Function MyAutomaticallyFindPeaks()

    String wname, xdata="_calculated_"
    Variable maxPeaks=100, minPeakPercent=5
    Prompt wname, "Peak Wave", popup, WaveList("*",";","DIMS:1,TEXT:0,CMPLX:0")+"_none_;"
    Prompt xdata, "X values", popup, "_calculated_;"+WaveList("*",";","DIMS:1,TEXT:0,CMPLX:0")
    Prompt maxPeaks, "Maximum Peaks"
    Prompt minPeakPercent, "Minimum Peak Amplitude (% max)"
    DoPrompt "Automatically Find Peaks", wname, xdata, maxPeaks, minPeakPercent
    if( V_Flag != 0 )
        return 0    // user cancelled
    endif
    
    WAVE/Z w=$wname
    WAVE/Z wx=$xdata
    Variable pBegin=0, pEnd= numpnts(w)-1
    Variable/C estimates= EstPeakNoiseAndSmfact(w,pBegin, pEnd)
    Variable noiselevel=real(estimates)
    Variable smoothingFactor=imag(estimates)
    
    Variable peaksFound= AutoFindPeaks(w,pBegin,pEnd,noiseLevel,smoothingFactor,maxPeaks)
    if( peaksFound > 0 )
        WAVE W_AutoPeakInfo
        // Remove too-small peaks
        peaksFound= TrimAmpAutoPeakInfo(W_AutoPeakInfo,minPeakPercent/100)
        if( peaksFound > 0 )
            // Make waves to display in a graph
            // The x values in W_AutoPeakInfo are still actually points, not X
            Make/O/N=(peaksFound) WA_PeakCentersY = w[W_AutoPeakInfo[p][0]]
            AdjustAutoPeakInfoForX(W_AutoPeakInfo,w,wx)
            Make/O/N=(peaksFound) WA_PeakCentersX = W_AutoPeakInfo[p][0]

            // Show W_AutoPeakInfo in a table, with dimension labels
            SetDimLabel 1, 0, center, W_AutoPeakInfo
            SetDimLabel 1, 1, width, W_AutoPeakInfo
            SetDimLabel 1, 2, height, W_AutoPeakInfo
            CheckDisplayed/A W_AutoPeakInfo
            if( V_Flag == 0 )
                Edit W_AutoPeakInfo.ld
            endif

            DoWindow ShowPeaks
            if( V_Flag == 0 )
                if( WaveExists(wx) )
                    Display/N=ShowPeaks w vs wx                
                else
                    Display/N=ShowPeaks w
                endif
                AppendToGraph/W=ShowPeaks WA_PeakCentersY vs WA_PeakCentersX
                ModifyGraph/W=ShowPeaks rgb(WA_PeakCentersY)=(0,0,65535)
                ModifyGraph/W=ShowPeaks mode(WA_PeakCentersY)=3
                ModifyGraph/W=ShowPeaks marker(WA_PeakCentersY)=19
            endif
        endif
    endif
    if( peaksFound < 1 )
        DoAlert 0, "No Peaks found!"
    endif
    return peaksFound
End
parameters for the MyAutomaticallyFindPeaks example code graph showing peak centers found by MyAutomaticallyFindPeaks example code

(A variation of this code has been added to Igor Pro 5.04´s Peak AutoFind procedure.)

See Also: Level Detection