jjhelmus / nmrglue

A module for working with NMR data in Python
BSD 3-Clause "New" or "Revised" License
208 stars 85 forks source link

Bug in unit_conversion().ppm_scale() #157

Closed LCageman closed 2 years ago

LCageman commented 2 years ago

I found that after creating a ppm_scale from a universal directory, the spectral width calculated from the first and last point of this ppm_scale does not match the original spectral width in the universal directory.

Running this code should make it clear

import nmrglue as ng

udic = ng.fileio.fileiobase.create_blank_udic(1)
udic[0]['sw'] = 5000.0
udic[0]['obs'] = 60.0
udic[0]['car'] = 200.0
udic[0]['size'] = 16000
uc = ng.fileio.fileiobase.uc_from_udic(udic)
ppm_scale = uc.ppm_scale()
sw_from_ppm_scale = (ppm_scale[0] - ppm_scale[-1]) * udic[0]['obs']
print(sw_from_ppm_scale)

>>> 4999.6875

I also found the issue, which is in the ppm_limits method in the unit_conversion class. Here the two outer limits of ppm_range are computed as '0' and 'size -1' (as you would normally define a list in python):

    def ppm_limits(self):
        """
        Return tuple of left and right edges in ppm
        """
        return self.ppm(0), self.ppm(self._size - 1)

However, these numbers are given to ppm(x), which calculates the ppm value using this formula, in which 'first' (most-left ppm value) and 'dx' (difference between points) are computed from the udic: ppm(x) = x * dx + first

Here, the issue lies in the fact that x = 0 will return the 1st value for the ppm_range and x = 'size -1' will return the 'last' - dx instead of just 'last', where 'last' is the most-right ppm value. x = 0 -> first x = size - 1 -> last - dx x = size -> last

Therefore, the ppm_limits method should be written like this:

    def ppm_limits(self):
        """
        Return tuple of left and right edges in ppm
        """
        return self.ppm(0), self.ppm(self._size)

I would write this in a pull request, but I see that other methods could also be affected (hz_limits, sec_limits, ms_limits, us_limits, etc.) and as this a central part of the module I think someone with more experience should look into this.

eisoab commented 2 years ago

Better watch out with this. I think it is correct as implemented, perhaps depending on the definition of "ppm_limits" and whether the "edge" in that function is meant to be the edge of the array (which it is now) or of the edge of the spectral region.

The spectral width (sw) should actually be 1 point broader than the difference between the first and last point of the spectrum, i.e. the points of the the array representng the spectrum does not cover exactly the complete spectral width, but falls one point short (usually on the rigtht side)

So if ppm_limits is defined as the left and right edge of the array representing the spectrum, the difference between these 2 should not be equal to the spectral width, and the line below should not give you the spectral width.

sw_from_ppm_scale = (ppm_scale[0] - ppm_scale[-1]) * udic[0]['obs'] To make it correct, you could change it to something like this:

size = udic[0]['size']
sw_from_ppm_scale = (ppm_scale[0] - ppm_scale[-1]) * udic[0]['obs']  * (size+1) / size

This is also related the following:

Any other definition of the ppm scale will lead to peaks having different ppm values depending on the amount of zero filling.

LCageman commented 2 years ago

Thanks, that makes a lot of sense. I already had the feeling I was overlooking something. I guess in other words, with an even number of points (which is always even) there is no center of the spectrum. Therefore the carrier frequency is at point 'size/2 +1' and the last point is ignored, I guess by convention and practical reasons.