window function

[BrunoVitorge]

Hello, I think I found a bug at least in the exponential window function.I did not check the others functions.
What I understood about exponential window function is that the last point of the fid should be multiply by the exp(-lb), so if the value of the last point is 1, and lb=5 after the application of exponential window function it value should be 0.006737946999085467
Then for the other points, the coefficient will depend on the number of point and on the point which is considered.
So the coefficient change for each points, for the point i, it will be i*lb/(number of points)
In the Bruker manual "procref" page 80 they also multiply by a "pi/2SWH" but I suppose it is to be consistent with the units.
Also the invert of the exponential should be more change also

Here is the way I programmed it :

def em_cor(data,lb = 0.0, inv = False):
    print 'shape', data.shape[-1]
    apod = np.exp(- (np.arange(data.shape[-1])+1) * (lb/data.shape[-1])).astype(data.dtype)
    if inv:
        apod = apod[::-1]
    return apod * data

I hope my explanation and the code above convinced you, that there is a problem with the window function.
I hope the other functions do not suffer such troubles
I found this bug because the signal to noise was decreasing along the spectrum (1D). And I was processing a fake fid with only 1 (a ones_like(fid)).
Hope this will help.
Thank a lot for the project.
Regards,
Bruno

[chatcannon]

The exponential window should depend on time, rather than number of points

def em_cor(data, lb=0.0, dt=1.0):
    apod = np.exp(-lb * (np.arange(data.shape[-1]) * dt))
    return apod * data

dt is the sampling period, often known as the dwell time

[jjhelmus]

Unfortunately there is no canonical definition of the window functions used in NMR (at least not in my opinion) and each program, spectrometer vendor and even researcher seems to have their own version which they use. For nmrglue I tried to make the the window function in proc_base be

1. As simple as possible, including a minimal number of parameters.
2. Have parameters expressed in units of points, or as much as possible not depend on experiment/nuclei specific values.

For the exponential window most everyone agrees that this should be of the form

exp(- something * some_index_of_the_data), what the something is and the some_index_of_the_data should be are where people disagree. Trying to make this as simple as possible I used:

apod = np.exp(-pi * np.arange(data.shape[-1]) * lb)

So the something is pi * lb and the some_index_of_the_data is just the index number (starting from 0). This matches the definition used by NMRPipe (http://spin.niddk.nih.gov/NMRPipe/ref/nmrpipe/em.html) if nmglue's lb parameter is set as NMRPipe's lb/sw, which makes the parameter be in units of points. It is also similar to the definition given in Hoch and Stern's NMR Data Processing book: exp(-pi * W * k * dt), so nmrglue's lb parameters is W * dt.

@chatcannon added a dt to the em window parameters, this would be more in like the version in NMRPipe and the definition in Hoch and Stern, but introduces an addition parameter and both parameters now have units that are not points (Hz and seconds). I think the single parameter version in nmrglue is simpler.

Also it should be noted that the mathematical form of all of the apodization function are given the in the nmrglue reference manual, e.g.
http://jjhelmus.github.io/nmrglue/current/reference/generated/nmrglue.process.proc_base.em.html

Finally it you obtain a similar version of an exponential window like @BrunoVitorge purposed by using the current em function in nmrglue as follows:

nmrglue.proc_base.em(data, lb / (data.shape[-1] * 3.1415926))

There is a index mismatch between the above and @BrunoVitorge function, but besides that they are identical.

[BrunoVitorge]

Hi guy I can see that you are very reactive,
sorry I was not clear enough, and I also did a mistake with the index.
The corrected version is this one
(For "pi" coefficient I don't have any opinion about it, so if you think it is needed added it):

def em_cor(data,lb = 0.0, inv = False):
    print 'shape', data.shape[-1]
    apod = np.exp(- np.arange(data.shape[-1]) * (lb/(data.shape[-1]-1))).astype(data.dtype)
    if inv:
        apod = apod[::-1]
    return apod * data

There is two main differences between nmrglue version and my version:
1- I DIVIDE LB BY THE NUMBER OF DATA POINT, the window function is then independent of the number of point in the fid (beside the "np.pi").
2- the inv condition, which seems more logic for me

Of course I can use nmrglue function with an modify lb_mod = lb/data.size
but I think it is less simple to use an lb of 0.001953125 if I have 1024 points, and 0.00390625 if I have 512 points in the fid, than 2 whatever is the number of point in the fid.

But now I understand why my lb coefficient was so low and changing in function of the number of points.
Please tell me if you change something.

[jjhelmus]

@BrunoVitorge I did not mean for my first comment to come across as reactive, I was trying to provide some background on why the current form of the em window was chosen. I now see more clearly the two differences in your em_cor function from the one currently in nmrglue. Adressing these.

1. An em function with or without a division by the number of points are both valid definitions of "the" em window, the question is which to include in nmrglue? My opinion is that the current version, without division is better, but this stems from my own use. I typically apply exponential line broadening in units of Hz, which can be done with the current version using ng.proc_base.em(data, lb_hz / sw) where sw is the spectral width in Hz and lb_hz is the broadening to apply in Hz. This would be much more complicated (and less elegant) with a em function that divides by the number of data points)
   There are also cases where the version including division by the number of data points would be simpler, ie your example of changing the size and wanting to use the same lb value, and so it may be appropriate to include both in nmrglue. Then the question is what to name the two functions (em and em2?) and if added functionality outweigh the confusion to users of having two similar but different em functions. I tend to lean towards thinking that the confusion outweights the added functionality and that there should only be one em function but a well written pull request with good documentation may change my mind.
2. Currently when the inv parameter is True the inverse of the window (1/window) is applied, this is consistent with a number of sources including NMRPipe's "-inv" flag. In this way you can recover the original signal (within a round-off error) by applying the inverse. For example:

data_0
data_1 = ng.proc_base.em(data_0, lb=2.0)
data_2 = ng.proc_base.emp(data_1, ln=2.0, inv=True)
# Now data_2 and data_0 should be nearly identical

What your are proposing is to apply the reversed window, which is different from the inverse. There are cases in which applying the reversed window would be useful, but it should be an additional parameter, maybe named rev, not a replacement of inv. I'd be open to a contribution to nmrglue that added a reverse window parameter to the apodization functions in proc_base, and may add such functionality in the future.

[BrunoVitorge]

@jjhelmus , thank you for your answer. I apology, if I was offensive in my previous message . I did not find how to make word in bold (I just wanted to emphases what was the difference between your and my version of em )

I was thinking that window function should not depend on the number of points in the fid, but It seems that I'm the only one :-) . After few test, Indeed nmrPipe EM function give the same result as nmrglue.proc_base.em function if you divide lb by sw (Hz), and this is what you do.

So I will probably do that now.

For the second point about the inverse exponential you are right again, and your suggestion about a "rev" option seems much better than what I proposed.
In fact I never use this kind of flag, Once I used a negative lb value to increase the resolution but that's all.
Sorry to have bother you, may I suggest to add a line in the comment of the nmrglue.proc_base.em explaining how to get the nmrPipe behaviour. By the way there is a typo, lb became lp in the comment at least in my nmrglue version.

Now I know how to apply an exponential function, next step, base line correction, :-)
Thank you for your useful comments, and for the nmrglue..

[jjhelmus]

Commit 19b5d90 adds a rev parameter to all the apodization functions in proc_base.py. Also fixed up the documentation to mention how to get results like NMRPipe and the typo.