Reference Deconvolution

It all began one month ago when Ed commented on this blog praising the Reference Deconvolution and the way it is implemented by Chenomx. I had heard about the subject, but never studied it, and apparently could not, because there is no scientific library in my town. Luckily, the Chenomx web site contains an article which not only introduces the basic concepts, but also contains some useful references. The latter of them:
KR Metz, MM Lam, and AG Webb. 2000. "Reference Deconvolution: A Simple and Effective Method for Resolution Enhancement in Nuclear Magnetic Resonance Spectroscopy". Concepts in Magnetic Resonance. 12:21-42.
is even freely accessible on the web and also makes for accessible reading. Even if you have never worked with complex numbers before, there you can learn how to perform reference deconvolution.
With this minimal theoretical background, without further reading, I have experimented with Reference Deconvolution and, actually, I have not found it to be a great improvement over the traditional lorentz-to-gaussian weighting. All my spectra are well shimmed and my lineshapes are quite near to the ideal lorentzian function, but even the pictures I have found on the review did not impress me so much. Here I will try to explain in even simpler terms what Reference Deconvolution is all about, assuming that the reader is familiar with the practice of weighting the FID in FT-NMR. You don't have to understand the equations, it's enough to understand the process graphically. iNMR users perform weighting by dragging a control with the mouse. The more they drag it, the more the shapes change. Almost everybody believes that you have to know what's happening at the math level, you need to know which functions you are using and with which parameters. I don't think so. It's enough that you open your eyes to have a complete understanding of the process. When you have your own musical tastes they will not change if you study the history of music. What matters is what you like: if you like a weighted spectrum, and if your software is well-written, you can soon learn weighting. You will continue to ignore the theory but your spectra will look much better. If you don't like weighted spectra, studying chapters of theory will not change neither your attitude nor the appearance of your spectra.
Everybody will agree with me on the following example, at least. The reason why the number of points of an NMR spectrum is a power of two is because the Cooley-Tukey algorithm requires so. Everybody knows it, but only a few spectroscopists actually know how the algorithm works and everybody agrees this is not a necessity!
Convolution is like dressing somebody who is already dressed: after the process he will have two clothes. A convoluted peak will have two shapes. A deconvoluted peak is like a naked person. Like the physician needs an undressed patient, the spectroscopist may need an undressed multiplet, to better measure its coupling constant or, more simply, to count the number of components. Convolution and Deconvolution in frequency domain are equivalent to multiplication or division in time domain. For the sake of simplicity, I will not mention the domains explicitly, from now on. For example: when I will say that I divide spectrum A by spectrum B, I mean that first I apply an inverse FT to both spectra, then I divide them, then apply a direct FT to the result. Division in the other domain is not only faster than direct deconvolution, but more feasible: we don't have the formula for the deconvolution...
When you perform the traditional lorentz-to-gauss transformation you are performing a deconvolution with a lorentzian and a convolution with a gaussian. The rationale is that you assume that the starting shape of your peaks is lorentzian. Let's suppose that this is true and you don't know it. Can you obtain the pure shape of your peak, in the format of a weighting function? The answer is yes. First step: you correct the phase of your peak. Second step: you move it to the position of the transmitter. This can't be done with common NMR software, yet it's trivial to write a program that shifts/rotate the spectrum; otherwise move the transmitter. Third step: you amplify or de-amplify the spectrum in order to set the area to one. At this point you have lost all the information about phase, frequency and intensity. All that remains of your original peak is pure shape. Perform an inverse FT and you obtain a real function (an exponential, to be more precise) that starts from 1 and decays towards zero. Once you know this function, you can divide by it any other spectrum and undress its peaks. Even if the other spectrum contains a different number of points, falling at different positions on the time axis, you can safely calculate the missing point by interpolation. A good thing of this approach is that, knowing that the result is a real function, you don't even need to correct the phase or change the frequency (steps 1 and 2). It's enough, if the spectrum contains a single peak, to calculate the absolute value of the FID.
You ask: if the sample is poorly shimmed, and the mathematical shape unknown, can I apply the same method to create a tailored weighting function? The answer is no. First reason: when the peak is not symmetric, the equivalent shape in time domain is not real but complex, and it means that you cannot eliminate the phase and frequency information. Second reason: even when the peak is symmetric, it goes from minus to plus infinite; if you are forced to truncate the spectrum (because of the noise, or to exclude other peaks) you are introducing unwanted components. The method outlined above does not work at all, because the distortions (caused by truncation) are heavy indeed. Reference Deconvolution works because it creates two sets of such distortions. You divide the first set by the second one, and the distortions will mutually cancel. The three players are: an experimental spectrum S, an isolated peak R (the reference peak) and its ideal shape I (the ideal peak). The deconvoluted spectrum D is simply obtained:
D = S (I / P).
Because I and P share the same frequency, phase and _truncation_effects_, all the distortions are mutually cancelled. You cannot, however, obtain the equivalent of I and P in time domain. They write that Reference Deconvolution is performed in time domain, but it's not the true, experimental, time domain. You cannot apply it directly on the FID. In all cases I have seen, you need to transform the spectrum to frequency domain and when, later in the process, you apply an inverse FT on it, you get something quite different from the FID. I think it's more correct to say that you are working into an artificial time domain.
From the point of view of the spectroscopist, it's a pity that any generic shape can't be described by a real function in time domain. Otherwise, reference deconvolution could be an extreme variant of weighting and inserted in the usual processing workflow. It would be possible to store the calculated shape and recycle it to weight different FIDs (through interpolation, as I has written). It would be possible to mix this shape (or, better, its inverse) with other weighting functions, etc...
In the real world, you are confined to the formula: D = S (I / P). You can also rearrange it into D = S / (P/I), as I am using to do in practice, but nothing really changes. The individual "I" and "P", in the artificial time domain, are of no use, only their ratio works.
There is some freedom in the definition of P and I. I have seen, for example, that the final result is very sensitive to the selection of P. The reference peak must be an isolated singlet with enough signal/noise. Once you have selected a region around it, you can verify that, enlarging or narrowing the region by a few points, changes the level of noise in D (the final result) in a non predictable way. In practice you need an adjustable control to test all the possible combinations in a couple of seconds, by trial and error (FT only take a fraction of second on recent computers). That part of the spectrum that is not selected is normally set to zero, but I wonder if there could be any advantage in smoothing it more slowly along the edges.
As for the shape of I, the review I have cited at the beginning suggests using a delta function, that is: all the area in concentrated into a single point, the rest is set to zero. In this case it is necessary to add a classic exponential or gaussian weighting (in the artificial time domain), either to I or to S, because the delta function, alone, amplifies the noise to intolerable levels. From what I have seen, Chenomx doesn't seem to adopt this recipe. It seems to me that they choose as "I" the lorentzian line graphically displayed by the program before the operation, which is probably truncated at the very same points as "P". In this implementation, the truncation effects on P and I are the same ones (first advantage) and there is much more control on the numerical value of the final linewidth (second advantage). This is why I wrote, in my previous post, that this implementation is optimized for line broadening. The delta function, instead, is better suited to achieve resolution enhancement. You can experiment by yourself, because the delta function alternative has been implemented by iNMR 2.1, and both iNMR and Chenomx's NMR suite can be downloaded and tried.
Theory says you can experiment with different ideal shapes, provided they have no doublet character (the FID of a doublet contains characteristic holes, in which you don't like to fall). Chenomx allows you to include the 29-Si satellites of DSS, which definitely possess a doublet character, and the algorithm still works great! The only possible explanation is that the noise has filled the holes. For a strange coincidence, the very first spectrum I experimented with, though sporting apparently suited singlets, didn't lend itself to reference deconvolution. I have tested the two programs and in both cases I've got wiggles everywhere (but at the reference peak). I can send this file to the interested reader.
It's like when they reprint a long-playing on CD. Quite often I like the new remix, but there are also cases when the original LP sounded better. My current opinion about Reference Deconvolution is that it's only a trick of Digital Signal Processing, not an essential part of experimental NMR. Trying it doesn't hurt, though. I have also found a reference to a drastically different application of the method.