[URGENT] Novice Igor pro user could really use some help. Curve fitting custom function

Share icon

Thu, 04/01/2010 - 04:38 pm

Hi there,

Very new to Igor, and I have a week until a project deadline.

What I have, is a Free Induction Decay signal from an NMR experiment. This means a sine wave, which decays exponentially.

I would like to do a fit to some oscilloscope data, so Igor interprets the x axis as time. The actual data, has two exp decaying sine waves superimposed (the other is ring down from the equipment, not important). A FFT reveals the natural frequency of the NMR signal.

I have tried to fit a function to the data, in order to extract the time constant of the exponential.

Function given to Igor:
f(time) = A + B*sin(C*time)*exp(-time/D)

Where A is y-offset, B is initial amplitude, C is frequency, & D is the time constant.

I'm having issues getting Igor to fit this function, I have suggested values but there's a problem. Substituting the FTT computed frequency in seconds^-1 for C makes the sine wave several orders of magnitude lower in frequency than the data.

When I click 'do it' the fit doesn't error, but produces a much lower amplitude wave..

Am I doing something wrong, I would really appreciate some manner of explanation of Igor's curve fitting abilities?

I hope my question isn't too novice, I haven't had time to familiarize myself with Igor yet, and I'm running on a short deadline.

Any help, or references to sources on this subject would be much appreciated.

P.S. additionally, is it possible to remove the second FFT component, of the ringdown?

Kind regards

niklz
Thanks for the offer, a link to the .pxp file can be found at:

http://rcpt.yousendit.com/847000391/1cf992188d15426876e366ef91132e23

The data I'm analysing is the onres2 in the editing and analysis folder.

If you FFT this, you should find a sharp component at 1.8MHz, and a smaller one at a slightly lower frequency.

Not that offres, is the background (equipment ring down \& no NMR), is it possible to perhaps Zero both of these to the y, and subtract one from the other?

Also please try out the curve fitting, see if you get the same results


Kind regards

niklz

April 2, 2010 at 03:08 am - Permalink

niklz-

I've spent some time playing with your data. You have a pretty challenging problem there.

1) You need to change your equation to account for the fact that sin() takes radians:
f(time) =D + A*sin(B*time)*exp(-time/C)
Once you do that, and use the correct frequency (about 1.853e6) the sin wave matches pretty well.

2) I think you did your FFT of the entire data set, but the part you're interested in is a small part of that set. Put cursors on the graph to show that part of it that you want to work on (presumably the part after that saturated bit, before the ringing gets to be the dominant factor). I was looking at about 0.00043 to about 0.0013.

3) It is extremely difficult to fit a sin function to so many periods. The slightest little error in the frequency will cause parts of the fit curve to be 180 degrees out of phase, and parts to be in phase. The parts that are out of phase will contribute maximally to chi-square; the easiest and most stable thing for a fit to do about that is to reduce the amplitude! So you wind up with a fit that has zero amplitude and has the y0 offset adjusted so that the resulting flat line passes through the middle of your data.

I think you're going to have to figure out how to get the envelope of your data and fit that. I will leave others to suggest how you might get the envelope.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

April 2, 2010 at 10:10 am - Permalink

Thanks for your reply.

Yes this was my original plan.

So new question, how would I find this envelope?

Also, is there away of subtracting one data set from the other, the on/off res. This should leave only the pure FID singal that I'm after.

Sorry if I seem to be abusing the forums with stupid questions, but I'm truly grateful for any time anyone spends! :)

Regards,

niklz

April 2, 2010 at 12:06 pm - Permalink

niklz:

Subtracting one data set from another is one of the easiest things to do in Igor. To avoid destroying your data, first make a duplicate:
Duplicate onRes, onResSubtracted

Then just subtract:
onResSubtracted -= offRes


You can do the duplication using a dialog if that is easier for you: Data->Duplicate Waves

You would benefit from looking at the Getting Started manual: Help->Getting Started. This includes a Guided Tour, step-by-step instructions to do various common things in Igor. It takes a couple hours but will save you much more than that pretty quickly.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

April 2, 2010 at 01:15 pm - Permalink

Thanks for the comments.

I will try and read the getting started bit, in-between writing up :).

Couple of things, I tried to model the function on a smaller range, thinking that It should make only minor statistical difference, but even in radians I'm still wildy out in terms of frequency.. Orders of magnitude =/. For frequency I have put "2*pi*1.8e6".

What I need to obtain the positive envelope, is a routine that finds local maxima, and deletes the rest, I could do it manually, but it would take hours and hours.

Any suggestions?

Kind regards,

niklz

April 3, 2010 at 04:26 am - Permalink

I am wondering if this is the part of your data you are concerned about:

Window Graph0() : Graph
    PauseUpdate; Silent 1       // building window...
    String fldrSav0= GetDataFolder(1)
    SetDataFolder root:'Editing & Analysis':
    Display /W=(298,227,1197,630) offres,onres
    SetDataFolder fldrSav0
    ModifyGraph rgb(onres)=(0,0,65535)
    SetAxis left -0.0731718800000001,0.0567968799999999
    SetAxis bottom 0.000416998776009792,0.000450998776009792
EndMacro


(copy that to the procedure window and then run it)

I think that you could easily do a fit to that portion of onres if that would get what you want.

But from the graph, I don't see how subtracting offres would be meaningful.

April 3, 2010 at 07:52 am - Permalink

Thank you for the reply,

I cant do what you ask though, sorry to sound like I don't want help. But I don't have the time to learn the language, so I cant interpret what you have written, I have tried doing my fit procedure on a smaller sample of the data, with no luck. Although it is entirely possible that I have done this wrong.

Is that procedure just taking a smaller sample of the data. It's telling me its expecting a wave name when I try to run it..

regards

niklz

EDIT:

Managed to find the range you have suggested, I think you may have misinterpreted what I want to do, I would like to extract the time constant for the exponential governing the decay of the sine wave. Am I right in thinking this sample is very close the the beginning of the data?

I need to make sure I average over those beats along the envelope, otherwise the result would be non-physical.

Is there an easy way of removing everything but the tips of all sine maxima?

again, regards,

niklz

April 3, 2010 at 08:08 am - Permalink

I was just wanting to know if the particular zoomed in portion of your data is the portion of interest.

You don't need to understand the macro (which is a recreation macro generated by Igor when you close a graph.)
Just copy the above text to the procedure window, and then choose Graph0 from the Windows->Graph Macros menu.

If the exponential time constant for the sine decay in that portion is what you are after, then I can help with a fit function. That would be more accurate and easier to do than devising an envelope extraction routine.

April 3, 2010 at 08:19 am - Permalink

OK yeah, how would you suggest I got about it?

Just use my regular fit function on that data sample?

Can anyone replicate the frequency issue I'm having, it's not important, but Im concerned that others seem to get other results.

thanks

niklz

EDIT:

Just noticed some strange issues, change the frequency by orders of magnitude does not correspond on the graph, that is, 1.8e8 appears higher in frequency to 1.8e9 =/

Also tried fitting, didnt work :(

April 3, 2010 at 08:52 am - Permalink

For the record, here is what I did to find the envelope. This required quite a bit of trial and error and if I were to do it over again, I would probably extract out a subset of the data, reset the origin to zero and work on that.

Starting with the provided experiment, copy the following to the procedure window:
Function myFit(w,x)
    Wave w
    Variable x
   
    Variable t0= w[0]                       // time offset to start of decay (fixed)
    Variable dc= w[1]                       // instrumental offset
    Variable sx= sin(w[2]*(x-t0) + w[3])        // our sine
   
    return dc + (w[5]*exp(-w[4]*(x-t0)) + w[6])*sx      // w[6] because the final amp is not zero
End


(That is an old-style fit function.)

Bring the command window to the front and execute the following as individual groups so you can follow along:
SetDataFolder root:'Editing & Analysis':
Duplicate/O offres,offresT0
Display offres,offresT0
SetAxis/A=2 left
SetAxis bottom 0.000417,0.000451
ModifyGraph rgb(offresT0)=(0,0,65535)

Make/D/O coef={0.000417, 0, 1.8e6*2*pi, 0, 0.1e6, 0.05,-0.001}
offresT0= myFit(coef,x)

FuncFit/H="1000000"/NTHR=0 myFit coef  offres(0.000417,0.000451) /D=offresT0 /R

Make/O envelope
SetScale x,0.000417,0.000451,"s", envelope
envelope= coef[5]*exp(-coef[4]*(x-coef[0])) + coef[1]
Append envelope


My resulting graph is attached.
Graph0_0.png Graph0_0.png Graph0_0.png

April 3, 2010 at 09:47 am - Permalink

OK, thanks very much. Think I'm starting to work out what'g going on here.

I shall attempt to do the same procedure on the onres data set. However the envelope was only necessary if it was impossible to get the FIT to work, however you did this.

Need to figure out how to reset the axis, but i'll get there, hopefully.

Thanks alot

niklz

April 3, 2010 at 10:38 am - Permalink

WOW! I think I have it working, might redo it, in an place without the offres ringdown.

Thank you all so much, I'll keep you posted on how it goes.

Kind Regards

niklz

April 3, 2010 at 11:06 am - Permalink

Managed to do the fit over a nice physical range, got my timeconst.

Thanks all so much, you've really saved me, am truly grateful.

Check graph if you're curious.

Kind regards and best wishes

niklz

EDIT: Ignore the x-axis units, mistake :)
FIDwithenv2.pdf

April 4, 2010 at 05:27 am - Permalink

niklz wrote:
Managed to do the fit over a nice physical range, got my timeconst.


Glad you got it figured out. I'm not sure how Larry got this fit to work where I didn't! My guess is that the initial guess of the frequency is pretty critical.

Something that's worth noting: I think your original problem with the frequency being off by an order of magnitude was most likely caused by aliasing. By default, the auto-destination wave from a curve fit (that's what you get if you select Auto Trace in the Output Options tab of the Curve Fitting dialog, or if you use QuickFit) has 200 points. Consequently, the high-frequency sinusoid is way undersampled, and it looks like it has way too-low frequency.

You can a more benign version of this in the original data and Larry's fit: the envelope of the data appears to have a bit of ringing over 4 to 5 periods of the basic high-frequency sinusoid. I was really surprised to see Larry's fit reproduced it perfectly, because there isn't anything in the fit function that does that. Then I realized that the low-amplitude peaks are simply sampled on either side of the peak value, and the high amplitude peaks are sampled almost at the peak value.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

April 5, 2010 at 04:33 pm - Permalink

What did you need the time constant for. If it was for a T2 relaxation time couldn't you do a specialised pulse sequence (CPMG) instead of analysing the FID decay?

April 6, 2010 at 12:26 am - Permalink