User Defined Peak Shape in Multipeak fitting 2.0

Share icon

Hi everyone,

I was hoping to add a user defined peak shape to my multipeakfitting 2.0 toolbox, but I am struggling to understand how to do thins. I have looked into the Peakfunctions2.ipf macro that has the peak functions, but I am unsure how to go about adding my own function. If someone could walk me through doing this, or may already have this peak shape written into a macro that they are willing to share I would appreciate it.

The function I want to add is a Breit–Wigner Fano peak shape. This is an asymmetric peak that is pretty common when looking at Raman spectroscopy of amorphous carbons. The function is as follows, sorry for the horrible formatting of the equation here.

I(ω) = (A/(1+(2(ω+ω)/T)^2)) + B[1+2((ω- ω2)/qT2)]2/[1+2((ω- ω2)/T2)]2.

Any help or guidance on how to properly get this added to my fitting options would be appreciated.

Thanks
Yeah, Multipeak Fit 2 makes it possible to add your own peak shapes, not easy to do so.

I started looking at the problem and I think your equation has a typo in it. Should this: 2(ω+ω) be something more like 2(ω-ω2)? I guess ω2 is what is commonly rendered as x0.

And can you give me some info on what controls width and amplitude?

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

June 10, 2015 at 11:26 am - Permalink

Yes, it is not extremely easy, but if you understand how it works it is rather straight-forward. You basically need to write four functions:
1) a function for naming you input/output parameters, as well as telling MPF the names of the other three functions.
2) a function calculating the shape of you peak using an x & y-wave input and a parameter set
3) a function calculating useful initial guesses from the standard asymmetric gaussian (the shape you draw when adding a peak) handed from MPF, or at least plug in the ones like width, height at the right places for your peak to work
4) a function calculating useful output parameters from you fit parameters (i.e., FWHM from width etc.). If you feel lazy just hand over the fit parameters and do the rest later.
Peeking at the at the already implemented one's (PeakFunctions2.ipf) might help a bit to get a grasp how stuff works. Also, you might want to read here a bit: http://www.igorexchange.com/node/5750
I attach my definition of the Doniac-Sunjic peak-shape for reference (this one is already rather advanced, as it involves Gauss-convolution to simulate experimental broadening, so please don't feel demotivated by this).

//________________________________________________________________________________________________//
//  MultiPeak-Fit function: Doniac-Sunjic shape
//________________________________________________________________________________________________//
 
Function/S DS_PeakFuncInfo(InfoDesired)
    Variable InfoDesired
 
    String info=""
    Switch (InfoDesired)
        case PeakFuncInfo_ParamNames:
            info = "Location;Width;Height;Asymmetry;GaussW;"
            break;
        case PeakFuncInfo_PeakFName:
            info = "DSPeak"
            break;
        case PeakFuncInfo_GaussConvFName:
            info = "GausstoDSGuess"
            break;
        case PeakFuncInfo_ParameterFunc:
            info = "DSParams"
            break;
        case PeakFuncInfo_DerivedParamNames:
            info = "Location;Height;Area;FWHM;Asymmetry;GaussWidth;"
            break;
        default:
            break;
    endSwitch
    return info
End
//################################################################################################
Function DSPeak(w,yw,xw)
    Wave w
    Wave yw, xw
 
    Variable Resolution = 10                                 // the sample resolution (higher = slower)
    Variable Xdelta     = abs(xw[0] - xw[NumPnts(xw)-1]) / (Resolution * NumPnts(xw))        // calculate an high-res xdelta for scaling
    Variable yPoints    = Resolution * NumPnts(yw)           // make the wave X times bigger than the original y wave
    Make/FREE/D/N=(2*yPoints+1) yWave                        // create a wave two times as big to have space before and after the peak
    SetScale/P x xw[0]-yPoints/2*Xdelta, Xdelta, yWave           // reserve 1/2 size of space before the x wave starts (1/2 after)
 
    yWave = w[2]*cos(Pi*w[3]/2+(1-w[3])*atan((x-w[0])/w[1]))/((x-w[0])^2+w[1]^2)^((1-w[3])/2)
    
    PeakConvolveGauss(yWave, w[4])
    
    yw = yWave(xw[p])        // write the calculated values into the resultwave
End
//################################################################################################
Function GaussToDSGuess(w)
    Wave w
    
    Variable x0 = w[0]   // extract the initial guesses
    Variable width  = w[1]
    Variable height = w[2]
    Variable lwidth = w[3]
    Variable rwidth = w[4]
    
    Redimension/N=5 w        // redimension to apropriate values
    // output values are: "Location;Width;Height;Asym;GaussW;"
    w[1] = 0.8*width     // 80% into lorentzian width
    w[2] = 0.8*height    // make the peak a bit smaller
    w[3] = 5*(rwidth-lwidth)/width   // asymmetry depending on width difference
    w[4] = 0.3*width     // 30% into gaussian width
    return 0
end
//################################################################################################
Function DSParams(cw, sw, outWave)
    Wave cw, sw, outWave
 
    Variable Amplitude, Location, FWHM, PeakArea
    Peak_FindStats(cw, "DSPeak", Amplitude, Location, PeakArea, FWHM) // approximately calculates all the results from the parameters
 
    Redimension/N=(6,2) outWave
    // output values are: "Location;Height;Area;FWHM;Asym;GaussWidth"   
    outWave[0][0] = Location                 // position
    outWave[0][1] = sqrt(sw[0][0])
    outWave[1][0] = Amplitude                // height
    outWave[1][1] = NaN                  // no idea how to calculate the error
    outWave[2][0] = PeakArea                 // area (probably wrong)
    outWave[2][1] = NaN
    outWave[3][0] = FWHM                 // FWHM
    outWave[3][1] = NaN
    outWave[4][0] = cw[3]                // asymmetry value
    outWave[4][1] = sqrt(sw[3][3])
    outWave[5][0] = cw[4]                // gauss broadening width
    outWave[5][1] = sqrt(sw[3][3])
End
//###################################### Value finder helper function ###################################
Function Peak_FindStats(cw, PeakFunc, Amplitude, Location, PeakArea, FWHM) // recalculate the FULL peak and extract all result values (approximately)
    Wave cw
    String PeakFunc
    Variable &Amplitude, &Location, &PeakArea, &FWHM
    
    Funcref PCISimplePeak PeakFunction = $PeakFunc       // reference the passed peak function
    
    Variable x0 = cw[0]
    Variable width  = cw[1] + cw[4]
    Variable PosH   = cw[2] > 0              // Look if height is positive
    
    Variable Size = 3000                         // create a big wave (fine increments)
    Make/FREE/N=(Size+1) yw, xw
    xw = (p/Size-0.5) * 30 * width + x0          // create a broad x scale (to make sure the full peak fits in)
    PeakFunction(cw, yw, xw)                     // recalculate the peak with the coef set
    WaveStats/Q yw                           // find the statistics
 
    if (PosH)
        Amplitude   = V_max
        Location    = xw(V_maxloc)
    else
        Amplitude   = V_min
        Location    = xw(V_minloc)
    endif
    
    FindLevel/Q/R=(0,Size/2) yw, Amplitude/2
    Variable Left   = V_flag == 1 ? NaN : xw(V_LevelX)
    FindLevel/Q/R=(Size/2,Size) yw, Amplitude/2
    Variable Right  = V_flag == 1 ? NaN : xw(V_LevelX)  
    
    FWHM        = abs(Left -Right)
    PeakArea    = AreaXY(xw, yw)
    return 0
End

June 10, 2015 at 07:03 pm - Permalink