Physically constrained Fourier transform deconvolution method

Francis A. Flaherty

doi:10.1364/JOSAA.26.001191

1. INTRODUCTION

Limitations on the resolution of physical measurements due to convolution with an instrumental response function, diffraction, or other causes are ubiquitous in the physical sciences. A considerable amount of research has been devoted to methods for deconvolving experimental data to obtain a more accurate representation of the original data. Deconvolution allows more information to be extracted from the data or the same amount of information to be obtained at lower resolution, allowing the use of simpler instrumentation. The approach employed here makes use of the a priori information that in many cases of interest the measured quantities are real and positive, which is employed in some of the better known methods [1, 2, 3, 4, 5].

2. THEORY

Convolution can be performed by multiplying the Fourier transform (FT) of the data by the FT of the instrument response convolving function and then applying the inverse Fourier transform (IFT). Let $D (τ)$ represent the original data consisting of N data points and $\tilde{D} (v)$ be the FT coefficients of the data. The tilde symbol will be used throughout to make it easier to distinguish the FT coefficients from the corresponding data series. $C (τ)$ are these data after convolution with an instrumental response function $R (τ)$ ; $C (τ)$ are the data actually recorded by the instrument. The convolved data $C (τ)$ are mathematically equivalent to the IFT of the product of the $\tilde{D} (v)$ FT coefficients and the $\tilde{R} (v)$ coefficients:

\tilde{D} (v) = N^{- 1} \sum_{τ = 0}^{N - 1} D (τ) \exp (- i 2 π (v ∕ N) τ),

\tilde{R} (v) = N^{- 1} \sum_{τ = 0}^{N - 1} R (τ) \exp (- i 2 π (v ∕ N) τ),

\tilde{C} (v) = \tilde{D} (v) \tilde{R} (v),

C (τ) = \sum_{v = 0}^{N - 1} \tilde{C} (v) \exp (i 2 π (τ ∕ N) v) .

The inverse process of convolution is deconvolution; it can be accomplished by inverse filtering, which involves dividing the

\tilde{C} (v)

coefficients of the convolved data by the

\tilde{R} (v)

coefficients of the convolving function and then applying the IFT. The problem with this simple approach is that in almost all cases of interest some of the

\tilde{R} (v)

coefficients are extremely small or equal to zero; frequently this occurs beyond some cutoff frequency. Information about these frequencies is lost and cannot be restored by this simple inverse filtering process, since it entails dividing by coefficients that are close to or equal to zero. If the coefficients of these frequencies are set equal to zero for lack of a better choice, the IFT yields a resolution that is higher than that of the convolved data but much lower than that of the original data. In addition, the deconvolved data will display spurious oscillations and may take on negative values even though the physical quantity can assume only positive values. The method presented here uses the a priori knowledge that the original data do not have negative values to restore missing frequency coefficients. Let

A (τ)

denote the inverse filtered data and M be the minimum magnitude of

\tilde{R} (v)

determined by the signal-to-noise ratio beyond which meaningful inverse filtering is not feasible:

\tilde{A} (v) = \tilde{C} (v) ∕ \tilde{R} (v) if ∣ R (v) ∣ ⩾ M,

This choice for

\tilde{A} (v)

avoids infinite or unreasonably large values when

\tilde{R} (v)

is close to or equal to zero:

A (τ) = \sum_{v = 0}^{N - 1} \tilde{A} (v) \exp (i 2 π (τ ∕ N) v) .

Let

B (τ)

represent the contribution to the data points that would be obtained by taking the IFT of the missing coefficients designated as

\tilde{B} (v) :

B (τ) = \sum_{v = 0}^{N - 1} \tilde{B} (v) \exp (i 2 π (τ ∕ N) v).

Since the Fourier coefficients of

A (τ)

and

B (τ)

are mutually exclusive, a restriction on the coefficients is

\tilde{A} (v) \tilde{B} (v) = 0.

Note that for conceptual simplicity the $\tilde{A} (v)$ coefficients can be regarded as nonzero below some cutoff frequency and the $\tilde{B} (v)$ as nonzero above the cutoff frequency. This method is not restricted to this case however, and the convolving-function Fourier coefficients could equal zero at various frequencies or ranges of frequencies that are not necessarily sequential as long as Eq. (8) holds. Ideally the restoration process should result in $D (τ) = A (τ) + B (τ)$ . The FT coefficients of $A (τ)$ are considered to be known from Eq. (5), while the FT coefficients of $B (τ)$ will be estimated. $A (τ) + B (τ)$ may yield negative values as a result of errors in the estimation of the $\tilde{B} (v)$ coefficients. This error will be used to improve the estimate of the $\tilde{B} (v)$ coefficients. The basic idea behind this method is to try to restore as much of the missing information as possible by making use of a priori information about the data, in particular the positive value of many physically measurable quantities such as, for example, the intensity of spectral lines. This can be accomplished by minimizing the sum of the squares of the negative portion of the deconvolved data. Proceeding in a manner similar to the approach of Howard [6], the Fourier expansion of a representative coefficient is

\tilde{B} (v) = N^{- 1} \sum_{τ = 0}^{N - 1} B (τ) \exp (- i 2 π (v ∕ N) τ) .

Using the definition of the Heaviside step function H(y),

H (y) = 0 y < 0, H (y) = 1 ∕ 2 y = 0, H (y) = 1 y > 0,

the sum of the squares of the negative values can be written as

\sum_{τ = 0}^{N - 1} [H (- {A (τ) + B (τ)})] {[A (τ) + B (τ)]}^{2} .

The minimization procedure cannot be applied directly to the Heaviside step function because it is not continuous. However, an alternative approximation to this function can be written that is continuous, making it possible to take derivatives. K in this expression is an arbitrary parameter that in the limit that

K \to \infty

reproduces the behavior of the Heaviside step function:

H (y) = \lim_{K \to \infty} {[1 + \exp (- K y)]}^{- 1} .

The sum of the squares of the negative values in this approximation is

\sum_{τ = 0}^{N - 1} {[1 + \exp (K {A (τ) + B (τ)})]}^{- 1} {[A (τ) + B (τ)]}^{2} .

The minimization condition can be satisfied by taking the derivatives with respect to the unknown coefficients and setting them equal to zero:

\frac{\partial}{\partial \tilde{B} (v)} \sum_{τ = 0}^{N - 1} {[1 + \exp (K {A (τ) + B (τ)}]}^{- 1} {[A (τ) + B (τ)]}^{2} = 0 .

Considering in detail the derivative of this expression with respect to

\tilde{B} (v)

combined with the definition of

B (τ)

from Eq. (7) and setting it equal to zero gives

\sum_{τ = 0}^{N - 1} 2 [A (τ) + B (τ)] {[1 + \exp (K {A (τ) + B (τ)})]}^{- 1} [\exp (i 2 π v τ ∕ N)] - {[A (τ) + B (τ)]}^{2} {[1 + \exp (K {A (τ) + B (τ)})]}^{- 2} [K \exp (K {A (τ) + B (τ)})] [\exp (i 2 π v τ ∕ N)] = 0 .

In the limit that

K \to \infty

for positive values of

[A (τ) + B (τ)]

, this expression approaches zero. In the limit that

K \to \infty

for negative values of

[A (τ) + B (τ)]

, the complex conjugate of this expression becomes

\sum_{τ = 0}^{N - 1} {[A (τ) + B (τ)]}_{-} \exp (- i 2 π v τ ∕ N) = 0,

where the notation

{[A (τ) + B (τ)]}_{-}

represents only the negative values of

[A (τ) + B (τ)]

and is otherwise zero. Similarly

{[A (τ) + B (τ)]}_{+}

will represent only the positive values of

[A (τ) + B (τ)]

thus:

[A (τ) + B (τ)] = {[A (τ) + B (τ)]}_{+} + {[A (τ) + B (τ)]}_{-} .

Taking the FT of this expression and subtracting Eq. (15) divided by N for a particular unknown frequency

v_{m}

, where the m subscript denotes a v missing from the measured data, the FT expansion gives

N^{- 1} \sum_{τ = 0}^{N - 1} \exp (- i 2 π v_{m} τ ∕ N) [A (τ) + B (τ)] = N^{- 1} \sum_{τ = 0}^{N - 1} \exp (- i 2 π v_{m} τ ∕ N) {[A (τ) + B (τ)]}_{+} .

Since the subscript m refers to the coefficients of the missing frequencies, Eq. (8) and the orthogonality of the summation eliminates $A (τ)$ from the left-hand side of Eq. (17), giving

N^{- 1} \sum_{τ = 0}^{N - 1} \exp (- i 2 π v_{m} τ ∕ N) B (τ) = N^{- 1} \sum_{τ = 0}^{N - 1} \exp (- i 2 π v_{m} τ ∕ N) {[A (τ) + B (τ)]}_{+} .

Note that only the values of

A (τ)

that result in a positive value when added to the estimate for

B (τ)

are summed on the right-hand side of Eq. (18); hence the orthogonality condition does not eliminate

A (τ)

. The left-hand side of Eq. (18) is precisely the

\tilde{B} (v)

coefficients from Eq. (9), giving a prescription for calculating these coefficients:

\tilde{B} (v) = N^{- 1} \sum_{τ = 0}^{N - 1} \exp (- i 2 π v τ ∕ N) {[A (τ) + B (τ)]}_{+} .

The basis of this method is the recognition that the

\tilde{B} (v)

coefficients can be obtained from the FT of

{[A (τ) + B (τ)]}_{+}

by employing an iterative method. Since initially

B (τ)

is unknown, the first approximation is obtained by taking the FT of the positive part of

A (τ)

and then performing the IFT, giving an estimate for

A (τ) + B (τ)

. The FT of

{[A (τ) + B (τ)]}_{+}

is calculated, and the

\tilde{A} (v)

coefficients are replaced with the values from the original inverse filtered data while the new

\tilde{B} (v)

coefficients are retained for the IFT. The iterative process is continued until the

\tilde{B} (v)

coefficients converge to the desired tolerance. Although the Heaviside function was approximated by the function involving an exponential in order to derive the equations, it is not necessary to use it in an actual computer algorithm to identify positive or negative coefficients.

An equivalent procedure that follows from Eq. (16) is to obtain the $\tilde{B} (v)$ coefficients by subtracting ${[A (τ) + B (τ)]}_{-}$ from the estimate of $[A (τ) + B (τ)]$ and utilizing these in the IFT. The advantage is that the term can be subtracted with a coefficient greater than unity, resulting in faster convergence, as will be discussed below.

3. DATA AND ANALYSIS

The method was tested by examining the isotope shift of the red line of a deuterium hydrogen gas mixture. The 512 point spectrum was collected on a Jobin Yvon/Spex Triax 550 monochromator with the slits set at $0.08 mm$ and is shown as the upper trace in Fig. 1 .

The data were processed by performing the FT of the spectrum and were inverse filtered by multiplying the frequency coefficients by a smoothly increasing function up to a cutoff limit at the 71st complex coefficient, beyond which all Fourier coefficients were set equal to zero. This function was obtained from the inverse of the convolving function of the spectrometer estimated from the profile of an isolated spectral line. The inverse filtered spectrum is shown in the lower trace of Fig. 1. Although the resolution is enhanced, there are significant negative excursions and spurious peaks.

The results of the method after various numbers of iterations are shown in Fig. 2 . The method gives a fairly constant improvement in resolution with each iteration up to some limit that depends on the details of the data and convolving function.

The results of the method after 100 iterations are shown as the upper trace in Fig. 3 , and the middle trace is the spectrum recorded with a $0.01 mm$ slit width for comparison. The width of the lines shows a fourfold improvement over the original data, and the method succeeded in greatly reducing the oscillations and negative excursions. The deconvolved spectrum compares very well with the data recorded with an eight times smaller slit width. The separation of the deconvolved peaks agrees with the values obtained from the higher-resolution data to within the experimental uncertainty. The area of the smaller peak relative to the larger was overestimated by approximately 40%; however, extensive experimentation with recorded and simulated data has shown that this error is not typical and is due to inaccuracy in estimating the convolving function for inverse filtering rather than to a defect in the method itself.

Since the missing frequencies are identified by their negative contribution to the spectrum, this underestimates their magnitude by a factor of 2 since only the negative area of each is taken into consideration. Experimentation bears this out, and multiplying ${[A (τ) + B (τ)]}_{-}$ by a factor of 2 before subtracting results in an algorithm that converges to the same solution in fewer iterations. Factors larger than 2 speed up the algorithm even more but can result in noise and spurious peaks, which ultimately can cause the method to diverge.

Another technique that reduced the number of iterations for a given increase in resolution was to subtract a constant offset from the estimate of $[A (τ) + B (τ)]$ . The price of this increase in speed was again the introduction of more noise and small spurious peaks. The most efficient approach was to use the offset for most of the iterations for a quicker increase in resolution and then perform a few more iterations without the offset, which eliminated much of the noise and spurious speaks. The optimum offset depends on details such as the accuracy of the inverse filtering and the signal-to-noise ratio. The progress of the method was monitored visually by plotting the new spectrum after each iteration and was stopped when little further improvement occurred. This could just as easily have been automated by terminating the process when the change after each iteration became less than the preestablished tolerance. However, it was found that visually monitoring the progress provided insight into how well the method was working and exposed potential problems if the deconvolution was pursued too aggressively.

This method was also evaluated by comparing it with a Landweber method modified to include a nonnegativity constraint as implemented by Tadrous [7]. This method gave similar results with somewhat higher resolution, although at the price of producing two spurious peaks. The height of the smaller peak was overestimated in both methods, but experimentation on simulated data with both methods showed that this was caused by difficulty in establishing the exact form of the convolving function rather than by the methods themselves.

4. CONCLUSIONS

The method detailed in this paper is computationally efficient especially for large data sets because the fast Fourier transform can be employed. The most critical factor in obtaining a good restoration is the choice of the initial convolving function. A poor choice for this function can introduce spurious peaks and noise. With a reasonable choice for this function, this method greatly reduces spurious peaks, negative excursions, and noise while providing a further increase in resolution. Even with some uncertainty in the choice of this function and with the presence of noise, the method proves robust and converges to a consistent and reasonably accurate representation of the original data. In applying the method, it may be necessary to establish a local baseline such as, for example, when a spectrum is superimposed on an interference pattern from a window. The same is true in applying this method to other areas such as image deconvolution, where, for example, text might be superimposed on a uniform background. As with any deconvolution method the signal-to-noise ratio is very important in determining the ultimate level of resolution enhancement that is achievable. An advantage of this method is the simplicity of implementing it and the ease with which other a priori information can be incorporated. Conditions such as finite extent [8] (also known as compact support) and bounded upper limits can be handled in exactly the same way.

ACKNOWLEDGMENTS

The author thanks Efstratios Manousakis, Cecilia Barnbaum, and Sean Barton for their suggestions for improving this manuscript.

Fig. 1 Upper trace, unprocessed data; lower trace, data after inverse filtering. The traces are not plotted to the same scale and have been offset for clarity.

Download Full Size | PDF

Fig. 2 Deconvolved data after 1 (top), 10 (middle), and 100 (bottom) iterations. The traces are not plotted to the same scale and have been offset for clarity.

Download Full Size | PDF

Fig. 3 Upper trace, data deconvolved by the method presented here from the $0.08 mm$ slit width data. Middle trace, $0.01 mm$ slit width data. Bottom trace, $0.08 mm$ slit width data deconvolved by the Landweber method with a nonnegativity constraint. The traces are not plotted to the same scale and have been offset for clarity.

Download Full Size | PDF

1. Y. Biraud, “A new approach for increasing the resolving power by data processing,” Astron. Astrophys. 1, 124–127 (1969).

2. P. A. Jansson, R. H. Hunt, and E. K. Plyler, “Resolution enhancement of spectra,” J. Opt. Soc. Am. 60, 596–599 (1970). [CrossRef]

3. B. R. Frieden, “Restoring with maximum likelihood and maximum entropy,” J. Opt. Soc. Am. 62, 511–517 (1972). [CrossRef] [PubMed]

4. B. R. Frieden and J. J. Burke, “Restoring with maximum entropy, II: superresolution of photographs of diffraction-blurred impulses,” J. Opt. Soc. Am. 62, 1202–1210 (1972). [CrossRef]

5. P. A. Jansson, Deconvolution of Images and Spectra (Academic, 1997).

6. S. J. Howard, “Continuation of discrete Fourier spectra using a minimum-negativity constraint,” J. Opt. Soc. Am. 71, 819–824 (1981). [CrossRef]

7. P. J. Tadrous, BiaQIm image processing software, Version 21.01, 2008, URL http://www.bialith.com.

8. S. J. Howard, “Method for continuing Fourier spectra given by the fast Fourier transform,” J. Opt. Soc. Am. 71, 95–98 (1981).

Physically constrained Fourier transform deconvolution method

Abstract

1. INTRODUCTION

2. THEORY

3. DATA AND ANALYSIS

4. CONCLUSIONS

ACKNOWLEDGMENTS

Cited By

Figures (3)

Equations (21)

Journal of the Optical Society of America A