GRPDLY and remove_digital_filter

[sciguy]

I've been using nmrglue to read Bruker data and even when I use remove_digital_filter, I get a lot of wiggles throughout the spectrum. I don't have this issue when the acqus file has a GRPDLY value above zero.

Looking at the remove_digital_filter code, when GRPDLY is < 0, the DECIM and DSPFSV values are used to find a value in the bruker_dsp_table look up table.

I've found that if I manually set the GRPDLY parameter to the look up table integer value (e.g. table value is 72.75, I set GRPDLY to 72), the resulting spectra looks much better.

I've tested this on the Processing-1D-Bruker-Data example with good results. For the example, there was no GRPDLY, only DSPFSV [12] and DECIM [32]. From the look up table the return value is 72.125. I added a GRPDLY value of 72 to the acqus file and ran the code from the example page: process_and_plot_nmrlgue.py. The only exception was I had to change the phase correction to p0=-50.0, instead of p0=-88.0.

Here is the spectrum with original settings from the example page:

Here is the spectrum using GRPDLY=72:

I'm relatively new to NMR, so I don't know if this approach is valid, all I know is the spectrum looks better.

Is this something that should be changed in nmrglue?

[jjhelmus]

This does look like a much better results. I have not had the time to really look but when the GRPDLY value is an integer then the processing in the remove_digital_filter function is limited to removing points from the beginning of the FID. Although this might give a better results graphically it may not be correct from a signal processing standpoint. Currently the remove_digital_filter function is attempting to reproduce the processing that NMRPipe does to remove the digital filter in Bruker files which is more complicated than simply removing points from the beginning of the trace. Give me some time to investigate this further.

[jjhelmus]

@sciguy I added an option to the remove_digital_filter to truncate the phase shift before applying the Bruker digital filter and make this the default setting. This should give better looking correction of Bruker data. Thanks for reporting this!

[ekwan]

I am having some problems with the digital filter correction. It mostly works, but the baseline still has a "frown" in it. I attached the raw data. Here's the code I used:

dic, FID = ng.bruker.read(filename)
FID = ng.bruker.remove_digital_filter(dic,FID)
C = ng.convert.converter()
C.from_bruker(dic, FID, udic)
pdic, pdata = C.to_pipe()
pdic, pdata = ng.pipe_proc.ft(pdic, pdata)                          # fourier transform
pdic, pdata = ng.pipe_proc.ps(pdic, pdata, p0=-40, p1=0)          # phase correction

Any ideas as to what might be wrong? Thanks!

sample.zip

[kaustubhmote]

This is a known problem with bruker datasets. See https://github.com/jjhelmus/nmrglue/pull/67 and the associated example files for details. the fix is, however, available only in the 0.7-dev version of nmrglue, so you might need to install the current master branch from github.

Briefly, this works in nmrglue 0.7-dev:

dic, FID = ng.bruker.read(filename)
pdata = ng.proc_base.fft(FID)                          
pdata = ng.bruker.remove_digital_filter(dic, pdata, post_proc=True)
pdata = ng.proc_base.ps(data=pdata, p0=-40, p1=0)          

If you cannot update your nmrglue, you can achieve the same thing by adding a 1st order phase correction to the data corresponding to the GRPDLY, i.e. remove the call to remove_digital_filter in the above example and apply a 1st order phase correction as follows:

pdata = ng.proc_base.ps(pdata, p0=-40, p1=360*dic['acqus']['GRPDLY'])
# edited to have p1 in degrees rather than radians

[jjhelmus]

@kaustubhmote did some really great work improving the logic in remove_digital_filter in the current development version. @ekwan Can you try out the dev version and see if that gives better results?

I'm happy to make a 0.7 release so that these improvements are more easy to obtain.

[robbyblum]

I was looking into this as well, back around when I submitted a PR for the TopSpin 4.0 change. I had also noticed that the release version (not the dev version with the improvements) didn't quite do what I expected it to do. It also didn't quite do what Bruker's own convdta command does. This is maybe not very helpful because I don't have time to dig into it again right now, but that command is another thing to potentially compare it to...

EDIT: Actually, I was working on the version that had @kaustubhmote 's changes added in. I think the remaining problem is in the "post_proc=false" case: the part where it flips the "ringup" and adds it to the beginning of the FID is the opposite sign from what Bruker's convdta command does (i.e., convdta subtracts the "ringup" from the beginning of the FID), and I don't think it makes sense to remove GRPDLY+2 points from the end, as opposed to just GRPDLY (this may be the correct thing to do in the case where GRPDLY doesn't exist, but I don't have data for that). Ultimately, if the "post_proc=false" procedure stopped after the IFFT, it might be more "correct," though you'd then have the "ringup" at the end of the FID and would have to worry about apodization commands.

[ekwan]

Thank you for your work everyone. I used the development version and it did indeed fix the problem. I suggest that this fix be placed in the next release, as it is not at all obvious from the examples that this kind of correction is necessary.

Am I right in saying that the artifact is caused by a delay between the last pulse and the start of acquisition, and the purpose of the remove_digital_filter method is to add a phase shift to counteract this? From the degree of ringing, it would seem that the delay is much longer than what one would see in Agilent spectra. Does that mean that correcting the time-domain data isn't always completely possible, because the data are discrete, and one can only shift over the data by an integer number of points?

[kaustubhmote]

At this point, I am considering the GRPDLY as just an initial guess for the number of points by which the FID has been delayed. If you go one integer lower or higher, you can still apply an additional first order phase correction to the data. In the figure below, I used bruker.remove_digital_filter with GRPDLYs of 75, 76 and 77 (76 is the correct one) and then manually + automatically phase corrected each resulting dataset (using proc_autophase.autops). The results are an exact match and the total first order phase correction is the same for all datasets (and very close to a delay of 76 points). I think the first order phase correction should take care of all of this even if the actual time shift is not an integer multiple. You will only end up seeing an additional non-zero first order phase correction.

Since the implementation details for the digital filters are trade secrets, the exact source of this delays, and more importantly, why is it so large for Bruker FIDs and not so much for Varian/Agilent ones, has been a matter of much debate and speculation. See here and here (2nd entry on the page) for some entertaining blog posts about this. Be sure to check out the discussions in the comments section

[ekwan]

Hi everyone, thanks for working to improve this feature. Here is a new test case illustrating the group delay correction that still seems to have a bit of a frown to it. I used nmrglue-0.8-dev. Any thoughts?

3.zip

[kfritzsc]

I thought @kaustubhmote had solved most of these problems. What was your "DIGMOD" set to during your experiment?

If the Bruker acquisition parameter "DIGMOD "is set to "baseopt "and "acqt0" is set correctly, the baseline problem is eliminated. If all goes well, the first order phasing (phc1) should also be zero.

More on the "acqt0" parameter: In many Bruker pulse sequences, there is a line "acqt0=x", where x is a float. The value is often negative half the excitation pulse. Sometimes there is a Pi in there. The parameter is supposed to indicate where the FID is thought to start relative to the start of acquisition. I guess if the baseline "smile" is a frown, the "acqt0" parameter needs to be more negative, and if the "smile" is happy, the setting needs to be more positive.

Using "DIGMOD"=basopt automatically changes the DSP filter setting from sharp(standard) to rectangular. I think that the main disadvantage is that the maximum spectral width is changed (on Avance III HD form ~ 625 kHz, as opposed to ~ 4 MHz). I don't believe this is a problem for most users. I wish there were more information about these parameters in the manual.

[kaustubhmote]

@ekwan, I think the additional 1st order phase correction you are doing using th ps command is the issue. If you apply only a 0th order phase correction as @kfritzsc mentions, the resulting baseline will be flat (or, as flat as what Topspin itself gives).

So, this works (with the data in Topspin processed with PHC0=-99.05 PHC1=0):

# pdata = ng.proc_base.ps(pdata, p0=97, p1=5)
pdata = ng.proc_base(pdata, p0=99.5)

[ekwan]

Wow! Thank you both for getting back to me so quickly! I really appreciate all the hard work you have put into this. This is an absolutely critical tool for my research and this is incredibly nice of you.

@kritzsc: The short answer is yes, all those things are set: DIGMOD is set to baseopt, acqt0 is set to -p1*2/pi as usual, and I am using the "rectangle" option. The one thing that is not optimized (because I don't know how) is DE, the pre-acquisition delay. I imagine this is similar to alfa or ddrtc on Varian systems. What's the best way to optimize DE?

Can either of you tell me why acqt0 has a factor of pi in it, given that it's already supposed to be in units of time?

@kaustubhmote: Your fix works! With a p0 of 99.0 only:

But watch this:

This apparently very equivalent command gives this:

Thank you so much for your help!

[kfritzsc]

@ekwan I don't take DE into account when setting acqt0. I think it must be done automatically.

If I interpret your question a little broader: similar to other systems, you should set the pre-scan delay as short as your equipment and experiment allows. (There is usually a lower limit based on the probe "ring down." There are other restrictions but modern powered transmit/receive switches, and other delay times are relatively fast. If you set DE too short, you will likely saturate the receiver and see an overflow error, even with otherwise correct receiver gain settings. I have a habit of setting my DE to be DW or 0.5 DW. With digital oversampling, I am not sure if my old pattern is necessary. )

I also don't have a reel answer to the factor of pi. I, too, did not expect it. Sheepishly, I admit that when writing pulse sequences, I usually guess and check this parameter.

As for your example, I think you are nicely showing what @kaustubhmote fixed in #67. Thepost_proc=True option seems to work great. As of yet, the algorithm for the removal of the digital filter in the time domain is not perfect.

[kangaroocry]

Thepost_proc=True option seems to work great. As of yet, the algorithm for the removal of the digital filter in the time domain is not perfect.

Yes, here I also do not understand the implementation. Why should the (in principle) same operation give different results depending on the type of domain you are working with. Applying remove_digital_filter with post_proc=True on the frequency domain data should give exactly the same result as applying remove_digital_filter with post_proc=False on the time domain data and FFT, shouldn't it?

[kangaroocry]

Ok, I examined this issue a bit in the source code: Mathematically, the post_proc=True method and the post_proc=False method are doing the same, with the one difference that for post_proc=False there is this strange fiddling with the shifting and skipping of data points:

    if post_proc:
        s = data.shape[-1]
        pdata = data * np.exp(2.j * np.pi * phase * np.arange(s) / s)
        pdata = pdata.astype(data.dtype)
        return pdata

    else:
        # frequency shift
        pdata = proc_base.fsh2(data, phase)

        # add points at the end of the specta to beginning
        pdata[..., :add] = pdata[..., :add] + pdata[..., :-(add + 1):-1]
        # remove points at end of spectra
        return pdata[..., :-skip]

If you comment out these two lines and add just return pdata both methods are the same (in their respective domains):

    if post_proc:
        s = data.shape[-1]
        pdata = data * np.exp(2.j * np.pi * phase * np.arange(s) / s)
        pdata = pdata.astype(data.dtype)
        return pdata

    else:
        # frequency shift
        pdata = proc_base.fsh2(data, phase)

        # add points at the end of the specta to beginning
        #pdata[..., :add] = pdata[..., :add] + pdata[..., :-(add + 1):-1]
        # remove points at end of spectra
        #return pdata[..., :-skip]
        return pdata

Here are the respective plots:

Original:

Improved (blue and red on top of each other):

[kaustubhmote]

@kangaroocry, I took a look at the post_proc=False case, and it does seem to be doing the same thing as post_proc=True in the initial steps using the fsh2 function. It takes a FT and then applies a 1st order phase shift and then takes a inverse-FT to give back the FID. So far its OK, but then there is the shift of points that I dont understand as well. As written in the comments after the docstring, it seems this was done to match the output of nmrpipe. This is also what convdta function in TopSpin, which converts raw FIDs into the AMX format seems to do (maybe that is the origin of all of this?)

Looking at all of this, I now feel that post_proc=True should be the default way to go about this. Unless your really, really want the FID to be corrected, there is no reason to take a FT, then phase-correct, then IFT and then FT again, instead of a simple FT and phase correction. And on second thoughts, maybe post_ft would have been a better argument to use for this function.

[ekwan]

@kaustubhmote: I agree that going back and forth between the time and frequency domains is wasteful. However, if you use post_proc=True and zero filling together, the algorithm no longer works. Any suggestions?

Despite the nmrpipe restrictions, I think the current behavior of the function (i.e., different with and without post-proc=True) is not intuitive and should probably be changed.

[kaustubhmote]

@ekwan, I am not sure why zerofilling should break this. The following works as well as before for the same dataset:

dic, fid = ng.bruker.read('.')
fid = ng.proc_base.zf_size(fid, 16384)
ft = ng.proc_base.fft(fid)
ft = ng.bruker.remove_digital_filter(dic, ft, post_proc=True)
ft = ng.proc_base.ps(ft, p0=99.05)

[ekwan]

@kaustubhmote: I was mistaken. Thank you for clarifying.

[kangaroocry]

After performing some tests and brain work I start to understand all this a bit better. The post_proc=True method is working fine for the absorptive (real) part of the spectrum. Still, the digital filter distortion is present in the dispersive part. It can be removed without affecting the (already nice) real part of the spectrum by removing the digital filter from the FID in a certain way. Being familiar with FFT you know that having a spectrum (here: distortion spectrum) with real part =0 means you have time domain data that are symmetric in a certain way. This is exactly why there are these steps to add part of the digital filter (now at the end of the FID) mirrored to the beginning of the FID. Unfortunately it is implemented wrong as you all could see. I figured out already how to do this in a correct way but I have to do more tests to be sure and I will report soon.

[ekwan]

Thank you for working so hard on this! I really appreciate your efforts.

[ajb204]

Hi chaps, I know this post is closed: the above discussion/examples are very interesting and helpful! Nmrglue is excellent! The 2/pi subtraction: you expect from a 90^o rectangular pulse to have spins evolve for a duration of (2/pi)*pw, which will show up as a p1 phase correction. In multi-dimensional sequences that have been well setup, you'll see subtraction of this time of this amount in places like indirect evolution periods to get the phasing right. So subtracting this is in principle a good idea for the direct dimension if you're aiming to minimise the need for a first order phase correction.