2. Theory
The path through the origin, development, and reduction to practice of Burg’s deconvolution theory began with the mid-century goal of extracting weak harmonic signals in stationary time series, where the signals are buried in Gaussian noise. The procedure started with autocorrelation, an operation incompatible with spectroscopy. Applications evolved through multiple steps, including forward prediction [
24]. We do not attempt to trace the path here. We consider only the result, cast in a form suitable for the present application. This is Andersen’s approach [
22], later extended to complex coefficients by Kesler and Hayken [
25].
In spectroscopy, the procedure is based on the Fourier coefficients of the spectrum being analyzed. To define these procedures, let
be a positive-definite spectrum consisting of
real values
,
, projected onto the range
according to
The restriction to an odd number (2
N + 1) of data is required by the maximum-entropy derivation, and is also mathematically convenient. Next, let
,
be the set of digital Fourier coefficients associated with the
, determined as
Because the
are real,
. Then, the CME representation
of
of order
M is given by [
23]
where the
are solutions of
The Toeplitz (diagonal-constant) form originated with the need to perform autocorrelation in time-sequence applications, and is retained for spectroscopy because maximum-entropy theory requires it. The theory is also best suited to Lorentzian features, because as shown below, Equation (3) reduces to Lorentzian form under small-term conditions.
The solution of these equations is readily obtained either by Levinson recursion [
26,
27] or by inverting the Toeplitz matrix. The resulting uppermost equation is solved for
, after which the remaining coefficients are calculated. When the resulting spectrum
of Equation (3) is Fourier analyzed, the original coefficients
are found to be recovered exactly, while those for
continue the established trend. As previously shown in [
23], it follows that this approach possesses all the advantages of the brick-wall filter with none of the disadvantages, yielding the most accurate noise-free representation of spectral data with no apodization (filter-cutoff) errors.
Burg/Andersen (B/A) deconvolution [
22,
25] is entirely different, sharing only the starting Toeplitz matrix and the final pseudo-Lorentzian form, Equation (3). Here, the objective is to obtain
such that the denominator of Equation (3) approaches zero as closely as possible at the locations
of the
j features in the spectrum. The result is a sharpened or “whitened” version of the original. The calculation is a Levinson-like recursion, generating a spectrum
that when Fourier-analyzed, the original
are found to be replaced with new values having a smaller decay coefficient, as described below.
The specific B/A procedure is as follows. Assume that the elements
of the solution vector of the
Mth-order Toeplitz-matrix equation are known, but with the vector normalized to
. The term
in the right side vector of Equation (4) is replaced by a “power”
that is also determined from the
Mth-order solution. (The terminology “power” is relevant for stationary time series but has no meaning for spectroscopy). Next, increase the size of
by adding 1 row and 1 column, so
. Now evaluate
From a stationary-time-sequence perspective, Equation (5) is a filter running in the forward direction, i.e., predicting the output
for increasing
, whereas Equation (6) is the same filter operating in the reverse direction, predicting
for decreasing
. Next, evaluate
The
-order solutions
,
are then given by
The recursion calculation is initiated with the starting values
and
. The starting value of
is determined by minimizing
with respect to
. The results are
The calculation then proceeds with Equations (7)–(11). The bidirectional averaging implicit in Equations (5) and (6) minimizes odd-harmonic contributions, leading to deconvolution.
3. Analytic Investigation
We now consider analytic solutions, with the objectives of obtaining information about the deconvolution process, properties of the solutions, and useful estimates of the amount of sharpening that can be expected in specific situations. As with our previous work [
23], we base our analysis on the Fourier coefficients
of the continuum normalized Lorenzian
which for mathematical simplicity we place at the center of the range
, defined by Equation (1). The digital Fourier transform
of
is
In Equation (19), it is assumed that
is large enough that the apodization (“ringing”) term in Equation (18) can be discarded. The
solution
of Equations (3) and (4) is identical to Equation (19), except that it is exact. As shown in [
23], with
, all higher-order terms vanish. For small
and
, Equation (20) reduces to
which differs from Equation (16) only by the
normalizing factor of the continuum Fourier transform. Thus, the broadening is the square root of the closest approach of the denominator to zero. The full-width-half-maximum (FWHM) is
.
We now consider the equivalent B/A
M = 1 solution using Equations (3) and (14), providing detail as necessary. Substituting
in Equation (14) yields
where Equation (24) follows from Equation (23) by canceling the common sums, and Equation (25) from Equation (24) by multiplying the numerator and denominator by
. As a result of the special property of
, the sums in Equation (23) cancel regardless of the initial and final values of
n, so the resulting Equation (25) is independent of
N.
The remaining mathematics is more efficiently carried out in two parts, evaluating first the power prefactor
and next the denominator
. For the prefactor, we find
assuming that
N is large enough that the apodization term in the numerator can be ignored. Here, the sum depends on the starting value of
n. We can remove this ambiguity by noting that the digital transform Equation (2) uses all terms
, whereas Equation (26) is single-ended. Consistency is achieved by taking the average of sums starting at
and
, in which case each term is considered once and only once. The result is
The contribution of the denominator is
Combining the two expressions yields the
lineshape
, i.e., the
deconvolved version of
:
The structural similarity between Equations (19) and (32) is evident.
This equation is best understood by repeating the small-term expansion leading to Equation (21) in the CME case. Performing the series expansions of
,
, and
for small
and
, we find
Ignoring the small
term in the denominator, the lineshape is seen to be Lorentzian, with a FWHM of
This can be compared to the FWHM of 2Γ for the original Lorentzian. The relative narrowing is more relevant, so we divide Equation (35) by
, which yields
where in Equation (37) it is assumed that the second term in the radical can be ignored relative to the first. Since the dimensionless reference for
is the Fourier period 2π, we find the surprising result that the relative sharpening is also a function of the width of the structure relative to the width of the spectral segment being analyzed. For structures with values of Γ that are small compared to the intrinsic reference scale of 2π, the sharpening can be significant.
While the
B/A FWHM is obviously less than the original for
, a short calculation shows that it is always the case if
Because this inequality is satisfied for any , the B/A process for always results in a narrower Lorentzian line.
Because the width and peak values of the Lorentzian are related, narrowing is equivalent to moving the singularity in the denominator closer to zero. At
, we find
compared to
for the original normalized Lorentzian. The pole has thus moved 2 inverse powers of
closer to zero, consistent with the decreased width.
At this point we note the importance of averaging: had normalization been performed using only the first term in the denominator of Equation (23), the result would be , yielding the CME result Equation (19). Had Equation (23) been normalized by the second term, the result would be , violating a fundamental constraint of the theory. In essence, the B/A approach works because of averaging.
It is straightforward to take the calculation one step further, determining the result for
. Using the same strategy, we find
whence
This expression is significantly more complicated than that for
, although its properties are easily determined. We consider first the extrema, which can be found by setting the derivative of the denominator with respect to
equal to zero. The calculation yields
The factored term
shows that one extremum occurs at
, as expected by symmetry. However, Equation (47) also has second and third extrema at locations given by
Because for , it follows that the right side of Equation (48) is less than 1 under all conditions, hence by symmetry this solution is valid for any Γ. Thus for , the reconstruction of the single Lorentzian is split into two peaks.
We next consider the small-term expansion for
. Expanding Equation (47) to fourth order in
and
, we find
The denominator is the square of a parabola with singularities at
. This is consistent with the above analysis: the single peak of the original Lorentzian has split into two, with each new singularity located a distance
from the symmetry axis. The infinities at
in Equation (45) are eliminated in a 6th-order expansion of the denominator, which yields the more accurate small-term result
where
Because , the location is a local minimum.
Thus, for , the singularities are moved even closer to zero, although the price paid is that two peaks are present instead of one. While it is possible to obtain an expression yielding the FWHM of each peak, the splitting has made the concept of sharpening irrelevant. We note that in going from to , the peak went from being the reciprocal of a fourth-order quantity to the reciprocal to a sixth-order quantity.
Next, we note that small-term expansions may be expected to have limited utility. As an alternative, replacing
and setting the denominator equal to zero yields
Again, we see that the result contains two peaks, which are now separated by .
The above calculations are too crude to have significant quantitative value, but they illustrate a basic mechanism that we expect to be valid in other situations: when too many peaks (too large a value of M) are requested, the existing peaks are not only shifted, but extra features appear. We expect this to apply to orders of M beyond 2 as well.
4. Discussion
To better visualize the nature of the solutions, we turn to numerical methods.
Figure 1 compares the normalized deconvolved spectrum for
, calculated from Equation (32), to the original Lorentzian and the associated pseudo-Lorentzian calculated from Equations (16) and (19), respectively. The curves are shown on the expanded range
to better display differences. On this scale, the Lorentzian and the pseudo-Lorentzian are essentially identical, but the deconvolved result is larger than the others by a factor of about 2½. From Equation (35) the linewidth is reduced by a factor of 8.1, in agreement with
Figure 1. It can be noted from Equation (35) that the extreme sharpening seen here is a direct consequence of the relatively small value
relative to the
range of
θ. Had we chosen a larger value of
, the relative amount of sharpening would have been less.
As noted above, when
M exceeds the number of singularities, the system generates extra poles.
Figure 2 illustrates this in detail. Here, we show the results of the B/A procedure for
and 3, comparing the results to each other and to the original Lorentzian at the bottom. The enhancement of the amplitude of the
spectrum with respect to the original Lorentzian is clear. Another striking feature is that the curve for
has two peaks instead of one, with the peaks located approximately at
,as shown in the previous section. The curve for
recovers the single dominant peak, but its amplitude is reduced by about a factor of 2. In addition, spurious satellites appear at
. Thus, if
M exceeds the number of legitimate critical points in the spectrum, the result is not only extra features, but also a significantly reduced peak height when a feature coincides with a structure of the spectrum.
Because spectra usually contain more than one feature, the question arises as to how much of the simple single-line theory is transferrable to spectra containing multiple structures.
Figure 3 shows a model spectrum consisting of a pair of lines with
, symmetrically located about the origin with a separation of
. This is sufficient for the peaks to remain distinct while retaining some overlap. It is seen that the
lineshapes are much sharper than the original, so the B/A procedure is effective. However, an examination of the numerical values shows that the broadening is about 0.042, leading to a ratio of 0.17 compared to the predicted value of 0.123. The difference is about 30%. In addition, the singularities are found to lie at about
from the symmetry axis instead of
, a difference of 17%. Although no measures of broadening have appeared before, errors in the locations of singularities are well-known. We shall discuss this in more detail elsewhere.
Figure 4 compares the capabilities of the CME and B/A approaches to deal with singularities separated essentially by the broadening parameter. The synthetic spectrum consists of two poles, the first at
with an amplitude of 0.5, and the second at
with an amplitude of 1.0. The broadening parameter
is the same in both cases, and both calculations were performed to the order
. The CME is functioning as it should, i.e., generating a replica of the data with white-noise coefficients replaced with most-probable values. The B/A response is qualitatively different, with the sharpening evident and the two singularities clearly identified. However, the B/A procedure places these at
and +0.258, which is only qualitatively correct. The CME performs better, but only at higher
M. For
, two poles are also identified, located at
and +0.110, in better agreement with the correct values.