Multi Peak Fit 2.0 - Adding user defined functions

Share icon

Mon, 05/19/2014 - 11:33 pm

Hello everyone,

I just "moved in" from O...never mind to Igor and I must say - I love Igor!
However, as any beginner I seek the advise and help of code experts.

I'm trying to add a sum frequency generation (SFG) curve fitting based on two (as a simple case) Lorentzian peaks that also have a coherence effect on one another.
The equation as follows is:

yw = ((abs(A1/((xw-w1)+G1*sqrtJ) + A2/(x-w2)+G2*sqrtJ)))^2

Where: A1 and A2 are the amplitudes of the two peaks.
w1 and w2 are the peak centers (usually in cm^-1)
G1 and G2 are the phases, respectfully.
sqrtJ is the imaginary number = i.

I load a csv file, than I use:
SetScale/P x,2800,5, '14051607'
// In order to set the x scale to 2800 - 3100 reciprocal cm.

I wrote the following code based on the help file of Multi-Peak fitting 2
#pragma rtGlobals=3     // Use modern global access method and strict wave access.

Function/S DoubleSFG_PeakFuncInfo(InfoDesired)
    Variable InfoDesired
   
    String info=""
   
    Switch (InfoDEsired)
        case PeakFuncInfo_ParamNames:
            info = "w;w1;w2;A1;A2;G1;G2" // w = wavenumber or "x", w1 = peak center, A1 = peak's amplitude,  
            break;                  //G1 = the peak's phase. 2 refres to the second peak.
        case PeakFuncInfo_PeakFName:
            info = "DoubleSFG_Peak"
            break;
        case PeakFuncInfo_GaussConvFName:
            info = "SimpleLorentzianGuess"
            break;
        case PeakFuncInfo_ParameterFunc:
            info = "DoubleSFG_PeakParams"
            break;
        case PeakFuncInfo_DerivedParamNames:
            info = "Location;Height;Area;FWHM;Wavenumber;Amplitude;Phase"
            break;
        default:
            break;
    endSwitch

    return info
end

Function SimpleLorentzianGuess(w) //Initial guess?
    Wave w
   
    Redimension/N=3 w
    // width of Lorenzian needs to be modified from the estimated Gaussian width
    w[2] = w[1]*w[2]*sqrt(pi)
    return 0
end



Function DoubleSFG_Peak(w, yw, xw, A1, A2, G1, G2, w1, w2) // w = wavenumber or "x", w1 = peak center, A1 = peak's amplitude,  
                                                                   //G1 = the peak's phase.
    Wave w
    Wave yw, xw, A1, A2, G1, G2, w1, w2
   
    Variable/C sqrtJ = cmplx(0,1)
    yw = ((abs(A1/((xw-w1)+G1*sqrtJ) + A2/(x-w2)+G2*sqrtJ)))^2
end

Function DoubleSFG_PeakParams(wW, wA, wP, outWave) //wW=Wavenuber, wA=Amplitude, wP=Phase initial guesses.
    Wave wW, wA, wP, outWave
   
   
    //Location
    outWave[0][0] = wW[0]
    outWave[0][1] = NaN
   
    Variable amp = wA[2]
   
    //Amplitude
    outWave[1][0] = amp
    outWave[1][1] = NaN
   
    //Phase
    outWave[2][0] = wP[2]
    outWave[2][1] = NaN
end


I am well aware that my code is full of holes yet I don't understand the following error when I hit the FIT button:
"Multi-Fit Failed:
Number of dimensions must match data wave".

I have two questions:
1) To what dimensions should I pay attention to - parameters, raw data, fit?
2) How to fix it?

Thanks,
Y.H

 14051607.csv DoubleSFG_Peak.ipf
As you have already guessed, there are some holes in your code. ;) But most of the things are actually simple misinterpretations of the inner workings of MPF. You put things into the Function definitions which are not supposed to be there. But let's look at it step by step:

0) You don't need w as a parameter, i.e.,
    case PeakFuncInfo_ParamNames:
        info = "w1;w2;A1;A2;G1;G2" // w1 = peak center, A1 = peak's amplitude,  

Also, you could use some more meaningful names (e.g., "Center 1,Center 2, Amplitude 1, Amplitude 2, Phase 1, Phase 2"), as only the order is important here. But that's to everybody's taste.

1) Your fitting function: The parameters don't get fed into the function as variables. That's what 'w' is for, the parameter wave.
At this stage the parameter wave should be properly resized (i.e. should have the size of 6 entries) and should be filled with initial guesses. Look how to use 'w' for your calculation (assuming above order of the parameters):
Function DoubleSFG_Peak(w, yw, xw)
    Wave w, yw, xw
    Variable/C sqrtJ = cmplx(0,1)
    yw = ((abs(w[2]/((xw-w[0])+w[4]*sqrtJ) + w[3]/(x-w[1])+w[5]*sqrtJ)))^2
end


2) Initial guesses: You need to change the Redimension command to match your number of parameters (6). Then fill the rows 0-5 with what is meaningful for your fit. Note that the wave is already pre-filled with parameters from the automatic or "by-hand" guess from the side of MPF assuming a possibly asymmetric Gaussian. Parameters w[0],w[1],w[2], are corresponding to position, width, and height respectively, which should give you something to work out your A1, w1, ... etc. You can of course also use constant values here, if you know some numbers which work for the majority of cases (e.g., for an asymmetric peak the asymmetry is difficult to guess and one could always start small like w[...] = 0.1 or something)

3) Output: This function converts the fit results to some meaningful output parameters, and takes exactly three waves. The output wave is two-dimensional to hold the value and some associated error value (which goes to the second column). The output wave size needs to match the number of parameters defined in PeakFuncInfo_DerivedParamNames. I guess you also want to think a bit more about what you actually want to have as output and how to calculate this. If you just want to have your six parameters directly, it's easy. But it seems you are shooting for some 'overall' amplitude etc. Do you know how to calculate this from your given parameters?
Function DoubleSFG_PeakParams(cw, sw, outwave)
    Wave cw, sw, outwave // cw=parameters from fit, sw = standard deviations for the parameters
    Redimension/N=(...,2) outWave   //adjust your output wave

    outWave[0][0] = cw[0]   // just outputs parameter 1 here
    outWave[0][1] = sqrt(sw[0][0])   // just the standard error for parameter 1
    // you can do any calculation here to get your desired "output", but you have to think about how to do this using the parameters in 'cw'
    // ...
 end


Hope that helps a bit. If not you have to wait until the real experts take over.

May 20, 2014 at 12:22 am - Permalink

Thank you very much Chozo!
I now realize that I didn't understand the parameter wave (w[]) role at all.
I've already begun to work on it and hope to post a solution (or another question) soon.

Thanks again,

May 20, 2014 at 10:02 am - Permalink

You may have trouble fitting a double-peak function into MPF2, which is really oriented toward fitting a single peak per peak. But if you are willing to put up with some manual labor, it may work.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 20, 2014 at 11:05 am - Permalink

It's possible that you should write it as a stand-alone user-defined fitting function. How many of these double peaks do you typically have?

Wait- I just looked at the CSV file you attached to the original post. I don't see double peaks. Can you describe exactly what's in that file, and what you expect to get out of this?

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 20, 2014 at 04:42 pm - Permalink

Hello everyone,

I apologize for my rather late response.
I have succeeded in writing some code for two Lorentzian peaks with a coherent interference between them. So for peak-A and peak-B we get the general equation of:

|A^2 + B^2 + 2A*B|^2.

If the two peaks centers are clearly separated, i.e., no overlapping can occur (the distance between the two peak centers is greater than their FWHM value) then there is no problem to treat them as simple Lorentzians. In general, as the molecule at hand has more carbon hydrogen groups (C-H2, C-H3) the number of coherent interference increases - there are triplets, quadruplets and so forth.

However, due to the physical nature of SFG, any overlapping peaks have also a coherent interference contribution (that is 2A*B, generally speaking) that arrises from the theoretical model of summing up two frequencies in order to get a third..

John, I'm looking for a better data representation for a ""double peak" sample, this will take time since it is not my data. Nevertheless, if you use the multi peak package on the data in 14051607.csv, the peak assigned as "1" has this mild interference.

I've added the code that I've written so far.
It seems to create an ""exploding peak" (higher amplitude than possible) each time I use it (for now only on the added .csv file) and I receive the following error message:
"The fitting function returned INF for at least one X value."

Looking forward to your insights :)
Y.H
DoubleSFG_Peak_0.ipf DoubleSFG_Peak_0.ipf DoubleSFG_Peak_0.ipf

May 23, 2014 at 07:31 am - Permalink

Which of the many data sets in that file are you fitting?

Could you post an Igor experiment file with a graph annotated with the features that you believe are the doubled peaks?

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 23, 2014 at 02:10 pm - Permalink

Hi John,

Thank you for your patience.
I have attached a .pxp file as you have requested with one data set (SFG double peak spectrum.pxp).
Please note that the "peaks" are actually local minima ("deeps") due to the non-resonant background of the substrate.

In accordance with the above remark, should I write a procedure that will "flip" the Y-axis data so a minimum will now be a maximum when I use the multi-peak fit? Can it detect "deeps"?

Thank you,
Y.H
SFG double peak spectrum1.pxp

May 23, 2014 at 03:35 pm - Permalink

Thanks- that's helpful.

Bear with us for a bit longer here...

So the "peaks" are actually "deeps" (what I would call a "negative peak")? Are all the features you want to fit negative peaks? And a few (one in the data set you just posted) are doubled due to the processes that you explained earlier?

There's a broad "deep" from about 2850 to about 2925. Is that a doubled peak?

The DoubleSFG_Peak function has separate parameters for center, amp, and width for each peak of the doubled peak. It effectively just fits a sum of two independent peaks. How about just using the built-in peak shapes? I have attached a picture of a fit that I did to your data set, using just negative Lorenzian peaks. I clipped off the right-most part because it looks like a peak (or two) that lacks the right tail, making it difficult to fit. Lacking the tail on one side will cause ambiguity in the width that can sink a curve fit.

One thing is that the automatic peak finder finds positive peaks by default. You can ask it to find negative peaks instead.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
MultipeakFit_Set2.png

May 23, 2014 at 05:26 pm - Permalink

Hi John,

In general, a positive peak or a negative one depend on the sample that I'm probing. Can you please tell me how you assigned the negative peaks? Can I incorporate it into the Multi-fit algorithm?
You are right that a simple set of Lorentzians describes this spectrum, however, while mathematically it sometimes works - physically it's not correct. Later on I expect to encounter triple peaks and even quadruple, moreover, most likely that I will be the first to try to assign any peak profile so I can't make any further simplifications.

Despite all of the above, it's possible to use the Lorentzian fit WITH a correcting factor (the 2A*B, i.e., the coherent interference part) AND a good background description - but one thing at a time.

Cheers,
Y.H

May 23, 2014 at 09:43 pm - Permalink

I just looked again at your DoubleSFG_Peak() function. Last week I failed to notice (Friday afternoon problem?) that the whole expression is squared. That changes things just a bit...

Let me try again.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 27, 2014 at 11:18 am - Permalink

OK- I just read (a small amount) about SFG and realized that in your expression abs() needs to be cabs(). I could tell you whole nitty-gritty story about why we have to have a different function for plain absolute value and for complex absolute value...

Even so, I think you need to debug your expression somehow. With the change to cabs() your function now gives one negative and one positive peak. At least they are finite at the center!

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 27, 2014 at 11:39 am - Permalink

Hi John,

Thanks for your comment and the extra effort in reading about SFG.
I've changed the abs() to cabs() - it looks way better!
When you write that the expression needs further debugging do you mean that adding another coherent interference, i.e., adding another peak is needed until the best fit is achieved?
Like: yw = ((Cabs(A1/((xw-w1)+G1*sqrtJ) + A2/(x-w2)+G2*sqrtJ) +A3/((x-w3)+G3*sqrtJ) ,etc....))^2

Cheers,
Y.H

May 28, 2014 at 12:12 pm - Permalink

No, I meant that from your previous descriptions I expected both peaks to be either positive or negative. But I didn't read the SFG stuff carefully, and it's not my field...

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 28, 2014 at 02:18 pm - Permalink

Hello everyone,

Indeed as John had suggested my function was missing something...a bracket :)
The code is near completion but before I publish it I'd like to resolve a final issue.
Once I achieve a "perfect" fit (it's perfect because I generated it according to the fit function) and press "Peaks results" an error message appears:
"...Index out of range for wave "Peak 0 Coefs". - the full message is attached.

I think that I have generated a coefficient wave with six parameters consisting of three variables (Center, Amplitude, FWHM) for each peak.
I seem to be wrong, hence the error message - How can I fix this?
Moreover, for connivance I'd like that each set of variables will be in a row of it's own instead of one very long row, much like for the case of fitting two different peaks, for example:
  • Center1 , Amplitude1, FWHM1

    • Center2 , Amplitude2, FWHM2

    Do you have any idea how should I do it?

    I have attached a .pxp file in which both data and fit are present. I hope this will help.

    Thank you,
    Y.H
Function Execution Error.png Y.H SFG profile test.pxp

May 30, 2014 at 10:58 am - Permalink

Is the error happening in your code? The best way to figure out what's wrong is to right-click in any procedure window and select Enable Debugger. Then right-click again and select Debug on Error. I also recommend the NVAR, SVAR, Wave Checking.

When an error occurs, the symbolic debugger will pop up, with the little yellow arrow (most likely) pointing at the line *after* the one causing the error.
Quote:
Function DoubleSFG_Peak(w, yw, xw, A1, A2, G1, G2, w1, w2) // w = waven

That's your fit function, right? You can't have all those extra input parameters beyond xw. The format must be that of an all-at-once fitting function.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 30, 2014 at 02:06 pm - Permalink

Hi John,

I took your advice. The error is not in my code. It is somewhere in the Multi-Fit code.
I haven't yet got the skill to modify it.
The good news are that if I select "Peak Results" in the format of "Report in Notebook"
I can easily copy-paste the relevant data.
However, it still bugs me - why the Fit Function Parameters are arranged according to their affiliation, e.g., Center2 = 2960.6 [1/cm]
while the first lines of the report, for example, PeakCenter2 equals gibberish? (kindly refer to Notebook).

This is a minor nuisance but I believe that I can learn from this type of error in my future code writing adventures.
Here is my code so far:
#pragma rtGlobals=3     // Use modern global access method and strict wave access.

Function/S DoubleSFG_PeakFuncInfo(InfoDesired)
    Variable InfoDesired
   
    String info=""
   
    Switch (InfoDEsired)
        case PeakFuncInfo_ParamNames:
            info = "Center1;Amplitude1;Width1;Center2;Amplitude2;Width2;NR-A;NR-B;NR-C;NR-Phase" // Center1 = peak center, Amp1 = peak's amplitude,  
            break;                          //Width1 = the peak's FWHM which is also refered to as "Gama". 2's refer to the second peak. Name the parameters as you like
        case PeakFuncInfo_PeakFName:        // BUT you must maintain their order of apperance in the parameter w[] wave!
            info = "DoubleSFG_Peak"
            break;
        case PeakFuncInfo_GaussConvFName:
            info = "SimpleLorentzianGuess"
            break;
        case PeakFuncInfo_ParameterFunc:
            info = "DoubleSFG_PeakParams"
            break;
        case PeakFuncInfo_DerivedParamNames:
            info = "PeakCenter1;Amplitude1;Area1;PeakWidth1;PeakCenter2;Amplitude2;Area2;PeakWidth2;NR-A;NR-B;NR-C;NR-Phase" // PeakWidth = FWHM = gama...
            break;
        default:
            break;
    endSwitch

    return info
end

Function SimpleLorentzianGuess(w) // Let's assume that a simple Lorentzian is a good guess to start with...
    Wave w
    //Note that the wave is already pre-filled with parameters from the automatic or "by-hand" guess from the side of MPF
    //assuming a possibly asymmetric Gaussian.
    //Parameters w[0],w[1],w[2], are corresponding to position, width, and height respectively,
   
    Redimension/N=10w
    // w[0] = Position1==>Center1               w[3] = Center2
    // w[1] = Width1==>Width1               w[4] = Amplitude2
    // w[2] = Height1==>Amplitude1          w[5] = Width2
    //*****************Non Resonat Background*******************************
    // sqrt(x*A^2 +Bx+ C)*exp(J*Phi)
    // w[6] = a, order of magnitude of SFG intensity
    // w[7] = b, three orders of magnitude smaller than SFG Inst.
    // w[8] = c, six orders of magnitude smaller than SFG Inst.
    // w[9] = Phi, the pahse due to substarte (delay or electronic effect) in [Degrees].
   
    // 1st peak: width of Lorenzian needs to be modified from the estimated Gaussian width
    w[2] = w[1]*w[2]*sqrt(pi)
   
    // 2nd peak: width of Lorenzian needs to be modified from the estimated Gaussian width
    w[3] = w[4]*w[5]*sqrt(pi)
    return 0
   
   
   
end



Function DoubleSFG_Peak(w, yw, xw)  // The w wave is a 10 object wave that has the fit parameters. Use w[0] to w[6], i.e.;
                                     // w[0] = Center1
                                     // w[1] = Amp1
                                       // w[2] = Width1
                                     // w[3] = Center2
                                     // w[4] = Amp2
                                     // w[5] = Width2
                                    // w[6] = a, order of magnitude of SFG intensity
                                    // w[7] = b, three orders of magnitude smaller than SFG Int.
                                    // w[8] = c, six orders of magnitude smaller than SFG Int.
                                    // w[9] = Phi, the phase in degrees!
    Wave w, yw, xw
   
    Variable/C sqrtJ = cmplx(0,1)
   
    yw =( cabs( sqrt(w[6]^2*xw+w[7]*xw+w[8])*exp(sqrtJ*(w[9]*pi/180)) + w[1]/((xw-w[0])+w[2]*sqrtJ) + w[4]/((xw-w[3])+w[5]*sqrtJ) ) )^2//Cahnged abs() to CABS()
end

Function DoubleSFG_PeakParams(cw, sw, outWave) //cw is the wave parameter from the fit,
    Wave cw, sw, outWave                         //sw is the standard deviations for the parameters - both are 12 wave parameters as defined in "PeakFuncInfo_DerivedParamNames"
   
   
    //Location (central position) 1st peak - treated as a simple Lorenztian
    outWave[0][0] = cw[0]           // This is actually the peak center of a Lorentian.
    outWave[0][1] = sqrt(sw[0][0])  // The error is treated accordingly.
   
    //Amplitude 1st peak - treated as a simple Lorenztian
    outWave[1][0] = cw[1]
    outWave[2][1] =  NaN // sqrt( (outWave[1][0]^2)*((sw[2][2]/cw[2]^2) + (sw[1][1]/cw[1]^2) - 2*sw[1][2]/(cw[1]*cw[2])) )
   
    //Area 1st peak = Intensity or oscillator strength
    outWave[2][0] =  cw[1]/sqrt(cw[2])  // Amplitude/sqrt(Gama)
    outWave[2][1] = NaN
   
    //FWHM 1st peak - this is "Gama", i.e., the width of the peak...
    outWave[3][0] = cw[2]
    outWave[3][1] = NaN
   
    //***********************************Second Peak*************************************//
   
    //Location 2nd peak
    outWave[4][0] = cw[4]
    outWave[4][1] = sqrt(sw[4][4])
   
    //Amplitude 2nd peak
    outWave[5][0] = cw[5] // 2*cw[5]/(pi*cw[5])
    outWave[5][1] = NaN //sqrt( (outWave[5][0]^2)*((sw[5][5]/cw[5]^2) + (sw[4][4]/cw[4]^2) - 2*sw[4][5]/(cw[4]*cw[5])) )
   
    //Area 2nd peak
    outWave[6][0] = cw[5]/sqrt(cw[6])
    outWave[6][1] = NaN
   
    //FWHM 2nd peak
    outWave[7][0] = cw[6]
    outWave[7][1] = NaN // sqrt(sw[4][4])
   
    //**********************************Non Resonant Background*************************************//
    //sqrt(a^2*xw+b*xw+c)*exp(sqrtJ*Phase*pi/180)
    //The a coefficient
    outWave[8][0] = cw[6]
    outWave[8][1] = NaN
   
    // The b coefficient
    outWave[9][0] = cw[7]
    outWave[9][1] = Nan
   
    // The c coefficient
    outWave[10][0] = cw[8]
    outWave[10][1] = Nan
   
    // The non-resonant pahse of the substarte
    outWave[11][0] = cw[9]
    outWave[11][1] = Nan
end

I've also attached a .pxp file (of actual data fitting!)
Cheers,
Y.H
Test my SFG profile.pxp

May 31, 2014 at 11:28 pm - Permalink

I looked at your experiment a bit. There are several things to improve in the code, but I guess this is sort of an in-between version so I just give a quick comment. First, I don't get an error when pressing 'peak results' (win7, Igor 6.34A). Did this problem resolved itself in the meantime? Regarding your output, I see that you assignment of the results in your DoubleSFG_PeakParams has some misplaced numbers which mixes up some output parameters. Here's a tip how to work with the output easier. Go to the folder Packages => MultiPeakFit2 => MPF_SetFolder_2 and display the wave 'Peak 0 Coefs' in a table. There you can see what parameter is at which position, and with that information you can fix your assignment of 'cw' positions.

June 1, 2014 at 06:02 am - Permalink

Hi Chozo,

I've done as you have suggested and indeed it worked out perfectly.
I'm using a MacBook Pro (OS X 10.9.3) Igor 6.34A.
Whenever I execute my code I get several different error messages (like string index exceeded, etc.), however, if I "force run" it using the arrow in the debugger I still get my results.
Is this due to OS issues or is it time to refine my code?

Cheers and thanks for your help,
Y.H

June 1, 2014 at 09:12 pm - Permalink

The errors could be both. Your function is rather fiddly and complex (e.g., it has a lot of parameters), so it could very well be that you are just triggering some buggy edge cases in MPF. Since the newer Igor versions have a rather strict checking for data boundaries, probably most errors are due to this which got just overlooked in earlier versions (which you are emulating now by pressing the green arrow continuously).
That you get an overall "good" looking result by this method doesn't necessary mean that everything is all right now. The fit seems not too stable considering the many parameters. Also you have several orders of magnitude between some of the parameters (like peak position and NR stuff), which could lead to rounding problems.
You can try to optimize/complete your code further and see how it goes. As for the MPF problems, I guess that's better looked over by John.

June 2, 2014 at 12:47 am - Permalink

If you think it's real bug in MPF2, let's move this to a support incident. Send an e-mail to support@wavemetrics.com. Attach a copy of your experiment, making sure to include all your code, and instructions for me to reproduce your problem.

It's possible that the error is induced by a bug in your code. I could look for that, too. It's much more convenient for me to deal with such things as a tech support incident.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

June 2, 2014 at 12:09 pm - Permalink