Multipeak Fit: Fit your experimental data with ... your experimental data?!

Igor's Multipeak Fit is convenient to fit analytical peak functions to data. But there might be times when you want to include some other data (wave) into your fit. For example, to find a match between different measurements or to include reference data into the fit. This code snippet provides the framework to fit experimental data including Gaussian broadening with Multipeak Fit (MPF).

The data is prepared in numbered folders (such as 'root:myData_1:') with numbered wave names (such as 'myPeak0', 'myPeak1' ...). Since MPF can only interact via numerical values, these numbers are later used to select the folder and desired data. After loading MPF, first navigate to Analysis => Multipeak Fit => Prepare Data Peak Function. A small prompt appears where the base folder and wave names are set ("root:myData_" and "myPeak" in above example).

Then start a fit set. Now, the DataPeak fit function is available from within MPF. Here, the last two parameters control the folder and wave number of the data you have just set up. These two parameters need to be on Hold at all times, or the fit will fail. Dial in the desired value to select the data to fit (e.g., folder 1 and wave 2), then run the fit. The first three parameters control the relative x-shift, Gaussian broadening and height (scale) of the data. If you want to switch off broadening, just set and hold this parameter at zero.

You can provide data in two variants: as-is or normalized. If you provide data as-is, then the first parameter can be seen as a 'relative shift' and the third parameter is a 'relative scale'. But you might want to normalize the data before using it in the fit: This would mean to normalize the data height to a maximum of 1 and/or adjust the x-scaling to be centered at zero (so that the 'peak' in your data is placed at x=0). For normalized data, the first parameter can be seen as a 'peak position' and the third parameter is 'height'.

Some notes and caveats:

  • You might want to at least remove the baseline from your input data or this baseline will appear in (and be scaled by) the fit.
  • If the data's x-range is smaller than needed for the fit, endpoints will be repeated to extend the range.
  • The output Location, Amplitude, Area, and FWHM values are numerically determined and may be junk depending on how much your data resembles 'a peak'.
#pragma TextEncoding = "UTF-8"
#pragma rtGlobals = 3
#pragma IgorVersion = 8.00
#pragma moduleName = MPFDataFit
#pragma DefaultTab = {3,20,4}
 
static StrConstant kMPF_WorkingDir = "root:Packages:MultiPeakFit2:"
 
Menu "Analysis"
    Submenu MPFDataFit#CreateDataFitMenu()  // add menu entry after loading MPF
        "Prepare Data Peak Function", /Q, MPF_PrepareDataFit()
    End
End
 
static Function/S CreateDataFitMenu()
    return SelectString(DataFolderExists(kMPF_WorkingDir),"-","Multipeak Fit")
End
 
// +++++++++
 
Function MPF_PrepareDataFit()   // setup for data locations
    String setBaseFolder = "", setBaseName = ""
    SVAR/Z baseFolder = $(kMPF_WorkingDir+"MPF_DataFit_baseFolder") // load previous settings
    SVAR/Z baseName = $(kMPF_WorkingDir+"MPF_DataFit_baseName")
    if (SVAR_Exists(baseFolder))
        setBaseFolder = baseFolder
    endif
    if (SVAR_Exists(baseName))
        setBaseName = baseName
    endif
    
    Prompt setBaseFolder,"Set base folder name without trailing number (e.g., \"root:myFit_\"):"
    Prompt setBaseName,"Set base data name without trailing number (e.g., \"peak_\"):"
    DoPrompt "Setup MPF data fit", setBaseFolder, setBaseName
    if (V_Flag)
        return -1
    endif
    if (!strlen(setBaseFolder) || !strlen(setBaseName))
        Abort "You need to enter both a folder path and base data name."
    endif
    
    int i, j, sourceFound = 0
    if (!DataFolderExists(kMPF_WorkingDir)) // if MPF was not started yet => create basic folder structure
        DFREF saveDF = GetDataFolderDFR()
        SetDataFolder root:
        for (i = 1; i < ItemsInList(kMPF_WorkingDir,":"); i++)
            NewDataFolder/O/S $StringFromList(i,kMPF_WorkingDir,":")
        endfor
        Variable/G currentSetNumber = 0
        SetDataFolder saveDF
    endif
    setBaseFolder = RemoveEnding(setBaseFolder , ":")   // make sure no colon was entered
    String/G $(kMPF_WorkingDir+"MPF_DataFit_baseFolder") = setBaseFolder
    String/G $(kMPF_WorkingDir+"MPF_DataFit_baseName") = setBaseName
    
    for (i = 0; i < 10 && !sourceFound; i++)    // try to look for the data in the first 10 folder / waves
        for (j = 0; j < 10 && !sourceFound; j++)
            Wave/Z source = $(setBaseFolder+num2str(i)+":"+setBaseName+num2str(j))
            sourceFound = WaveExists(source)
        endfor
    endfor
    
    if (!sourceFound)   // still nothing found? => have a simple check to remind the user
        DoAlert/T="Checking folder and data" 0, "Your first folder and/or data does not seem to exist yet. Make sure to prepare the data in the format: "+setBaseFolder+"[NUM]:"+setBaseName+"[NUM]"
    endif
    return 0
End
 
//########
 
static Function/Wave getSource(Wave cw)     // get source data from saved base names
    SVAR/Z baseFolder = $(kMPF_WorkingDir+"MPF_DataFit_baseFolder")
    SVAR/Z baseName   = $(kMPF_WorkingDir+"MPF_DataFit_baseName")
    if (SVAR_Exists(baseFolder) && SVAR_Exists(baseName))
        Wave/Z source = $(baseFolder+num2str(cw[3])+":"+baseName+num2str(cw[4]))
        if (WaveExists(source) && WaveType(source,1) == 1 && WaveDims(source) == 1)
            return source
        endif
    endif
    return $""
End
 
//########
 
static Function MPF_setHoldStateForPeak(Wave cw, String holdStr)    // injects hold state for the current peak
    DFREF DFpath = GetWavesDataFolderDFR(cw)
    int peakNum
    sscanf NameOfWave(cw), "Peak %d Coefs", peakNum
    if (!DataFolderRefStatus(DFpath) || numType(peakNum) != 0 || peakNum < 0)
        return -1
    endif
    SVAR gName = DFPath:GraphName
    String pName = gName+"#MultiPeak2Panel" // the MPF panel
    if (WinType(pName) == 7)
        Wave/T HoldStrings = DFPath:HoldStrings
        Execute/Z "WMMultiPeakFit#MPF2_RefreshHoldStrings(\""+pName+"\")"   // first, grab the latest hold settings from the panel
        HoldStrings[peakNum+1] = holdStr    // inject hold string
        Execute/Z "WMMultiPeakFit#SetHoldCheckboxesFromWave(WMMultiPeakFit#GetSetNumberFromWinName(\""+pName+"\"))" // rewrite to panel
    endif
    return 0
End
 
//####### MPF peak format functions #####
 
Function/S DataPeak_PeakFuncInfo(Variable which)    // pass names to MPF
    Switch (which)
        case 0:
            return "LocShift;Broaden;Height;Folder;DataNo;"
        case 1:
            return "MPF_DataPeakFunc"
        case 2:
            return "MPF_GaussToDataPeakGuess"
        case 3:
            return "MPF_DataPeakParams"
        case 4:
            return "Location;Height;Area;FWHM;"
    EndSwitch
End
 
//+++++++
 
Function MPF_GaussToDataPeakGuess(Wave cw)  // initializes parameters, and set hold state
    // cw values are: "PosShift;GaussWidth;Height;Folder;DataNo;"
    Redimension/N=5 cw
    int i,j, sourceFound = 0
    for (i = 0; i < 10 && !sourceFound; i++)
        for (j = 0; j < 10 && !sourceFound; j++)
            cw[3] = i
            cw[4] = j
            Wave/Z source = getSource(cw)   // grabs the first valid source wave
            sourceFound = WaveExists(source)
        endfor
    endfor
    
    int normInt = 1, normPos = 1
    if (sourceFound)
        Variable wMax  = WaveMax(source)
        Variable left  = leftx(source) < rightx(source) ? leftx(source) : rightx(source)
        Variable right = leftx(source) < rightx(source) ? rightx(source) : leftx(source)
        normInt = wMax > 0.9 && wMax < 1.1 ? 0 : 1
        normPos = left < 0 && right > 0 ? 0 : 1
    endif
    
    cw[0] = normPos ? 0 : cw[0] // if data seems centered around zero, peak position can stay
    cw[1] = 0.5 // initial Gaussian broadening
    cw[2] = normInt ? 1 : cw[2] // if data is normalized, peak height can stay
 
    // this needs to be executed later or the hold state will be reset
    Execute/P/Q "MPFDataFit#MPF_setHoldStateForPeak("+GetWavesDataFolder(cw,2)+", \"00011\")"
    return 0
End
 
//+++++++
 
Function MPF_DataPeakFunc(Wave cw, Wave yw, Wave xw) : FitFunc  // the fit function
    Wave/Z source = getSource(cw)
    if (!WaveExists(source))
        yw = NaN
        if (Exists("MPF2_DisplayStatusMessage"))    // display an error in the MPF graph
            Execute/Z "MPF2_DisplayStatusMessage(\"Base folder setup wrong or data not found.\","+GetWavesDataFolder(cw,1)+")"
        endif
        return -1
    endif
    
    Duplicate/free source, yWave
    SetScale/P x, DimOffset(source,0)+cw[0], DimDelta(source,0), ywave  // prepare shifted convolution wave
    PeakConvolveGaussHelper(yWave, cw[1])   // Gauss convolution
    
    Variable minX = pnt2x(yWave,0), maxX = pnt2x(yWave,numpnts(yWave)-1)        // make sure to not run out-of-bounds (end-point data is repeated)
    yw = xw[p] >= minX ? ( xw[p] <= maxX ? cw[2]*yWave(xw[p]) : cw[2]*yWave(maxX) ) : cw[2]*yWave(minX)
    return 0
End
 
//+++++++
 
Function MPF_DataPeakParams(Wave cw, Wave sw, Wave outWave)
    Variable Amplitude = NaN, Location = NaN, FWHM = NaN, PeakArea = NaN        // set to NaN for now
    FindPeakStatsHelper(cw, Amplitude, Location, FWHM, PeakArea)    // approximately calculates all results
    // cw values are: "PosShift;GaussWidth;Height;Folder;DataNo;"
    // output values are: "Location;Height;Area;FWHM;"
    outWave[0][0] = Location        // position
    outWave[0][1] = NaN
    outWave[1][0] = Amplitude       // height
    outWave[1][1] = NaN
    outWave[2][0] = PeakArea        // area
    outWave[2][1] = NaN
    outWave[3][0] = FWHM            // FWHM
    outWave[3][1] = NaN
    return 0
End
 
//#######
 
Static Function PeakConvolveGaussHelper(Wave in, Variable width)    // convolve with a prepared Gaussian
    if (width == 0 || numtype(width) != 0)
        return 0
    endif
    
    Variable pnt, dx  = DimDelta(in,0), size = DimSize(in,0)
    width = abs(width/dx)/sqrt(2)   // width normalized by x scaling and convert to expected Gauss 'width'
    pnt = round(max(abs(10*width),11))  // create a wave 10 times the Gauss width (5 sigma on each side)
    pnt = min(pnt, 2*numpnts(in)+1)     // doesn't have to be larger than two times the full input data
    pnt = pnt+!mod(pnt, 2)  // force odd size
 
    Make/FREE/D/N=(pnt) GaussWave   
    GaussWave = Gauss(x, (pnt-1)/2, width)
    Variable A = sum(GaussWave) // make sure the Gauss area limited even for a few points
    if (A > 1)
        GaussWave /= A
    endif
    
    Make/FREE/D/N=(2*pnt+size) temp = in[0]     // enlarge the wave to get the full convolution range
    temp[pnt,size+pnt-1] = in[p-pnt]
    temp[size+pnt,size+2*pnt-1] = in[size-1]
    MatrixOP/Free temp = replaceNaNs(temp,0)    // makes sure there are no NaN values in-between
    Convolve/A GaussWave temp
    in[0,size-1] = temp[p+pnt]
End
 
//#########
 
// recalculate the (almost) full peak and extract all result values numerically
Static Function FindPeakStatsHelper(Wave cw, Variable &Amplitude, Variable &Location, Variable &FWHM, Variable &PeakArea)
    Wave/Z sw = getSource(cw)
    if (!WaveExists(sw))
        return -1
    endif
    // cw values are: "PosShift;GaussWidth;Height;Folder;DataNo;"
    Variable left = pnt2x(sw,0)+cw[0], right = pnt2x(sw,numpnts(sw)-1)+cw[0]
    Variable pntSize = DimSize(sw,0)*10 // create a big wave (fine increments)
    pntSize = pntSize > 100000 ? 100000 : pntSize   // make sure its not too big
    Make/FREE/D/N=(pntSize+1) yw, xw
    Make/FREE/D/N=2 Levels
    SetScale/I x, left, right, xw; xw = x
    //+++++ find values +++++
    MPF_DataPeakFunc(cw, yw, xw)    // recalculate the peak with the coef set
    PeakArea = AreaXY(xw, yw, left, right)      // the area is calculated for the full source range
    WaveStats/Q yw
    Amplitude = cw[2] > 0 ? V_max : V_min
    Location  = cw[2] > 0 ? xw[V_maxloc] : xw[V_minloc]
    FindLevels/Q/D=Levels yw, Amplitude/2   // find the half-maximum positions
    if (V_LevelsFound == 2)
        left  = xw[Levels[0]]
        right = xw[Levels[1]]   
        FWHM  = abs(left - right)
    endif
    if (numtype(FWHM) != 0)
        return 0
    endif
    //++++++ recalculate with higher resolution ++++++
    // choose a range just a bit bigger than the FWHM
    SetScale/I x, left - 2*FWHM/(pntsize+1), right + 2*FWHM/(pntsize+1), xw; xw = x
    MPF_DataPeakFunc(cw, yw, xw)
    WaveStats/Q yw
    Amplitude = cw[2] > 0 ? V_max : V_min
    Location  = cw[2] > 0 ? xw[V_maxloc] : xw[V_minloc]
    FindLevels/Q/D=Levels yw, Amplitude/2
    if (V_LevelsFound == 2)
        left  = xw[Levels[0]]
        right = xw[Levels[1]]   
        FWHM  = abs(left - right)
    endif
    return 0
End

 

 MPF Fit Data Peaks_0.ipf (11.82 KB)

Forum

Support

Gallery

Igor Pro 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More