## Abstract

Superresolution by data inversion is the extrapolation of measured Fourier data to regions outside the measurement bandwidth using postprocessing techniques. Here we characterize superresolution by data inversion for objects with finite support using the twin concepts of primary and secondary superresolution, where primary superresolution is the essentially unbiased portion of the superresolved spectra and secondary superresolution is the remainder. We show that this partition of superresolution into primary and secondary components can be used to explain why some researchers believe that meaningful superresolution is achievable with realistic signal-to-noise ratios, and other researchers do not.

© 2006 Optical Society of America

## 1. Introduction

Superresolution research has been conducted for many decades, both to develop ways to achieve superresolution and to characterize the properties of superresolution (see Refs.[1], [2], [3], and the references therein). In the research literature, the term “superresolution” is commonly used to describe the reduction in the point spread function (PSF) width in image-domain intensities as a result of applying some sort of processing methodology when compared to the PSF width in the corresponding unprocessed intensities. This definition results from the classical concept of resolution as proposed by Lord Rayleigh in which the resolution of a system is related to the minimum separation of two points that can be resolved. PSF width reduction can occur as the result of two different phenomena, either separately or in combination. The first phenomena is an increase in the Fourier amplitudes of the image-domain intensities within the measurement bandwidth. The Fourier amplitudes can be increased by postprocessing techniques or by including appropriate phase and/or amplitude pupil-plane filters in the optical system. This latter approach is often referred to as “optical superresolution” and is useful when a narrow pre-detection PSF is desired as is the case, for example, for optical disk readout [4] or for real-time resolution increase in adaptive optics data [5]. The second phenomena is the extrapolation of the measured Fourier data outside of the measurement bandwidth. By its very nature, only postprocessing can accomplish this second type of superresolution, and thus the term “superresolution by data inversion” has been coined to describe this type of superresolution [1]. The superresolution we explore in this paper is superresolution by data inversion.

There are two schools of thought on superresolution by data inversion, hereafter referred to simply as superresolution. The first is that meaningful superresolution (defined as accurate extrapolation of the measured data more than an incremental amount beyond the measurement bandwidth for measured data with space-bandwidth products (SBPs) ≳ 10 and realistic signal-to-noise ratios (SNRs)) is unachievable. The definition of SBP will be given in Section 2. This belief is based upon signal-to-noise analyses carried out under the assumption that the object of interest has finite support in the image domain [1], [6], where the support of an object is the region in the image domain outside of which the object’s intensities are zero. The second school of thought is that meaningful superresolution is possible. This belief is based on the appearance of structure in the Fourier spectra of postprocessed imagery outside of the measurement bandwidth that appears to be correlated with the Fourier structure of the true object [3], [7].

In this paper, we resolve this apparent paradox by analyzing the Fourier- and image-domain properties of superresolved spectra for finite-support objects. This analysis is carried out by first decomposing the superresolved spectrum into primary and secondary superresolution components. We define primary superresolution as that portion of the superresolved spectrum that can be estimated in an essentially unbiased manner given the data SNRs and define secondary superresolution as the remainder of the superresolved spectrum. The phrase “essentially unbiased” will be explained and quantified in Section 4. We then show that the inclusion of the primary superresolution component of the superresolved spectrum in an image reconstruction produces a space-invariant but minimal increase in spatial resolution. When both the primary and secondary components of the superresolved spectrum are included in the image reconstruction process, the increase in spatial resolution, averaged over the support region, is greater than that achieved using just the primary superresolution component, but is still minimal. However, we show that the increase in spatial resolution due to secondary superresolution is space variant. The space-variant property is the key to resolving the disagreement between the two schools of thought. In particular, we show that the first school of thought is correct when considering the average increase in resolution across the support region, while the second school of thought is correct when considering the space-variant increase in resolution. Specifically, we will show that there is a non-negligible increase in spatial resolution near the edges of the support region.

The analysis in this paper is carried out in the context of achieving superresolution in the estimate of a 1-D support-constrained object, where the measured data is contained in a low-pass region of the Fourier domain. This problem formulation is particularly convenient because the inverse problem can be expressed analytically. In Section 2, we look at the inverse problem and characterize it in terms of its eigenvectors and eigenvalues. In Section 3, we analyze superresolution using the algorithm-independent Cramér-Rao lower bound (CRB) approach. The CRB analysis indicates that any unbiased estimate of the object’s Fourier spectrum at any frequency outside the measurement bandwidth has an infinite variance. For this reason, in Section 4, the superresolution analysis is approached in terms of the eigenvector/eigenvalue expansion of the inverse problem where only a finite number of terms are included in order to keep the estimator variances finite. The bias properties of the superresolved spectrum are then explored as a function of the amount of extrapolation and the number of terms included in the expansion. From this analysis comes the decomposition of superresolved spectra into primary and secondary superresolution components. Additionally, in Section 4, properties of primary superresolution are described. In Section 5, properties of secondary superresolution are given. Finally, conclusions and future work are presented in Section 6.

## 2. Forward and inverse models

In this section, equations for the forward and inverse problems are given. Particular attention is placed upon an eigenvector/eigenvalue expansion of the inverse problem in the Fourier domain. Properties of these eigenfunctions and eigenvalues are presented in order to gain insight into the properties of the variances of unbiased estimates of the object’s Fourier spectrum. An expression is given for these variances in terms of an infinite summation involving the eigenfunctions and eigenvalues and it is shown that the variances of the superresolved spectrum appear to be infinite. Due to the lack of closed-form expressions for the eigenvectors, however, conclusive statements are not possible. This motivates the CRB approach developed in Section 3.

The forward model describing the functional dependence of the image i(x) on the noise-free object o(x) with support [-x_{o},x_{o}], the system PSF h(x), and the additive noise n(x), is given by

where h(x) is the inverse Fourier transform of an ideal low-pass filter with a cutoff frequency of f_{o} (in cycles/unit length) and n(x) is stationary zero-mean Gaussian noise with variance σ [6], [8]. It is straightforward to generalize this forward model to include other PSFs and noise models, but this simple and oft-used image model is chosen for purposes of clarity. The eigenfunctions of this forward model are the celebrated prolate spheroidal wave functions (PSWFs) of Slepian, Pollack, and Landau [9], [10], and are a function of the space-bandwidth product (SBP) 4x_{o}f_{o} [8]. These eigenfunctions and their associated eigenvalues can be used to rewrite Eq. (1) in a form that permits solving for o(x) in terms of the measured data i(x). Although o(x) can be estimated from ‘raw’ values of i(x), lower-noise estimates of o(x) can be obtained if only the portion of i(x) associated with data inside the measurement bandwidth is included in the estimation process [11]. We follow this latter approach for our results but carry out the calculations in the Fourier domain to facilitate analysis of the superresolved spectra.

To continue, we write the inverse problem in terms of an eigenvalue/eigenvector expansion of the Fourier transform I(f) of the measured data i(x). The expansion can be written in terms of the true object spectrum O(f) and the Fourier transform of the noise N(f) since I(f) = O(f) + N(f). The desired estimator expression is given by

where O_{e}(f) is the estimated (and superresolved) Fourier spectrum. By restricting the integration limits in Eq. (2) to ±f_{o}, only the signal-bearing portion of I(f) is included in the estimate of O(f). For this reason, the estimator in Eq. (2) has the same noise properties as the image-domain estimator proposed in [11]. The PSWFs {ϕ_{m}(x)} of Eq. (2) are scaled versions of the eigenfunctions of the forward model and their associated eigenvalues {λ_{m}} are their energies inside the measurement bandwidth [9]. Although, strictly speaking, the PSWFs in Eq. (2) are defined only on [-f_{o},f_{o}], they are analytic functions and can be extended to (-∞,∞). In the following analysis, this extension is assumed to have occurred. Because ϕ_{m} (f) is non-zero almost everywhere for each value of m, M in Eq. (2) must be set to infinity to obtain an unbiased estimate of O_{e}(f) at any frequency, in general. Although any practical application can include only a finite number of terms, analyzing the infinite-summation case is important because it provides insight into the SNR properties of unbiased estimates. From Eq. (2), it is easily seen that the expected value of O_{e}(f) for M=∞ is just the true object spectrum. To obtain the SNRs of O_{e}(f) at any frequency, an expression for the variances of O_{e}(f) are needed and can be obtained in a straightforward manner from Eq. (2), yielding

$$\phantom{\rule{4.5em}{0ex}}={\sigma}^{2}\sum _{m=0}^{M}{\lambda}_{m}^{-1}{\phi}_{m}^{2}\left(f\right).$$

The standard imaging definition of SNR is the expected value of the estimated quantity divided by the square root of its variance [12]. It can be seen from Eqs. (2) and (3) that the SNRs of O_{e}(f) depend on the properties of the PSWFs and their eigenvalues. For these reasons, these properties are explored next.

In general, PSWFs for any value of the SBP have similar properties when appropriately normalized. For example, any PSWF whose index satisfies m≲SBP has most of its energy inside the measurement bandwidth. Most of the energy of any PSWF whose index satisfies m≳SBP is outside the measurement bandwidth. We call every PSWF whose index satisfies m≳SBP a *superresolving PSWF* for this reason. To gain insight into additional properties of PSWFs, it is useful to look at a representative set of PSWFs as a function of frequency and index value. Using the method of Latham and Tilton [13], for any SBP we are able to generate PSWFs and their eigenvalues quickly and accurately down to λ_{m}≈10^{-30}.

Figure 1 is a movie displaying PSWFs as a function of their index for a SBP of 40. These PSWFs were generated in a 1024 element array of which 256 elements were in the measurement bandwidth. Their associated eigenvalues are plotted in Fig. 2. In addition to the energy properties described in the previous paragraph, notice that the superresolving PSWFs become essentially decaying sinusoids whose starting locations f_{m} in frequency space monotonically increases as a function of m. This implies that the energies of the superresolving PSWFs become increasingly concentrated at higher and higher frequencies. Recall, however, that the energy of any PSWF inside the measurement bandwidth is equal to its eigenvalue; thus, from Fig. 2, it is clear that every PSWF has non-zero energy inside the measurement bandwidth. As a result, noise-free knowledge of the true object’s Fourier spectrum inside the measurement bandwidth permits superresolving this knowledge to all frequencies.

The primary difference in the properties of PSWFs and their eigenvalues as the SBP is varied is the decay rate of the eigenvalues. When this difference impacts the properties of primary and secondary superresolution, it will be called out. Otherwise, results described here apply to PSWFs for any SBP value.

Returning now to the analysis of the SNRs of O_{e}(f) for all frequencies when M = ∞, recall that the expected value of O_{e}(f) is just the true object spectrum. Therefore, for the SNRs to remain finite, var{O_{e}(f)} must remain finite. From. Eq. (3), it can be seen that the question of whether or not var{O_{e}(f)} is finite for any frequency is determined by the index dependence of the ratio ϕ^{2}
_{m}(f)/λ_{m}. From Fig. 1, it can be seen that ϕ^{2}
_{m}(f)→0 as m → ∞ for any frequency f, and from Fig. 2 that λ_{m} → 0 as m → ∞. As a result, their relative rates of convergence to zero determine whether or not var{O_{e}(f)} remains finite. Because there does not exist an analytical expression for ϕ_{m} (f), Eq. (3) cannot be used to definitively answer this question. However, we can use the PSWFs and eigenvalues calculated for Figs. 1 and 2 to calculate var{O_{e}(f)} for values of λ_{m} down to ≈10^{-30} to see if it appears to be converging to a finite number for any frequency. A movie of var{O_{e}(f)} as a function of the upper limit of the summation in Eq.(3) is shown in Fig. 3. A key property seen in this movie is that var{O_{e}(f)} does not seem to be converging to a finite value for any frequency, even those inside the measurement bandwidth, although the rate of increase of var{O_{e}(f)} inside the measurement bandwidth is significantly less than outside. Because of this fact, it appears that the estimator given by Eq.(2) can provide an unbiased estimate of the object spectrum at any frequency only at the cost of an infinite variance.

The conclusion that an unbiased estimate of a support-constrained object’s Fourier spectrum at any frequency has infinite variance was based upon a specific estimator (i.e., Eq.(2)). It is natural to wonder if an alternate estimator exists that can produce an unbiased estimate of an object’s Fourier spectrum with finite variances for at least some frequencies. To generate an algorithm-independent answer, the next section uses a CRB approach to calculate the variances of any unbiased estimate of a support-constrained object’s Fourier spectrum.

## 3. Cramér-Rao bounds

In this section, an algorithm-independent approach for calculating a particular set of lower bounds to the variances of any unbiased estimate of O(f) is employed. These lower bounds, called CRBs [14], are the diagonal elements of the inverse of the Fisher information matrix (FIM) corresponding to O(f) and the conditional measurement probability density function (PDF) that is parameterized by O(f). Our approach to calculating these CRBs is to first generate the FIM corresponding to o(x), then take its inverse, then transform this inverse into the inverse of the FIM corresponding to O(f), and then extract its diagonal elements. Although the CRBs could have been calculated by creating the FIM for O(f) directly, it will be shown that the transformation process from the image-domain inverse FIM to the Fourier-domain inverse FIM provides important insight into the properties of the Fourier-domain CRBs.

The unbiased CRBs corresponding to a vector of parameters are lower bounds to the variances of any unbiased estimates of these parameters. CRBs for a function of these parameters can be calculated when the Jacobian of the function with respect to the parameters is known. Although CRBs are not necessarily achievable, in many cases they can be achieved and in other cases they can be approached closely. To generate the CRBs for a vector of parameters, the PDF of the uncertainty in the measurements, conditioned on the parameter values, must be known. For the imaging model of Eq. (1), this conditional PDF can be obtained from the noise PDF and the noise-free imaging model. However, because CRB theory applies to a vector of random variables, not a stochastic process, Eq. (1) must be rewritten in a sampled data format. For this reason, we replace the continuous variable x by a set of evenly-spaced locations covering the interval [-x_{o},x_{o}] for some arbitrary number of locations. Let **α** be a vector containing these locations and let **y**, **θ** and **η** be vectors that contain the values of i(**α**), o(**α**), and n(**α**), respectively. In addition, let **H** be the matrix associated with h(**α**) [15]. This permits rewriting Eq. (1) as a matrix-vector equation given by

where the PDF of **y** conditioned on **θ** is given by

and f_{η}(**η**-**Hθ**) is the PDF of **η** but with mean **Hθ**.

As mentioned earlier in this section, the first step in calculating the CRBs for any unbiased estimate of O(f) is to calculate the FIM **F**(**θ**) corresponding to **θ**. The element of **F**(**θ**) in the p^{th} row and the q^{th} column is given by

where ln denotes the natural logarithm and E[] is the expected value of the expression in the brackets. For the imaging model of Eq. (4), Eq. (6) becomes

$$=\frac{1}{{\sigma}^{2}}\sum _{k}h\left({\alpha}_{k}-{\alpha}_{p}\right)h\left({\alpha}_{k}-{\alpha}_{q}\right)$$

where the second equality in Eq. (7) is written in terms of H and the vector of x locations **α** rather than in terms of the elements of **H** to enhance the reader’s understanding of how we calculate **F**(**θ**). Notice that each element of **F**(**θ**) is the inverse of the noise variance multiplying a discrete approximation to the integral of the product of two shifted versions of the PSF. This discrete integration approximation must be accurate to greater than one part in 10^{30} to generate results comparable in accuracy to the noise variance calculations in Section 2. One approach to obtaining this accuracy is to densely sample the support of the object; however, that leads to large dimensions for **F**(**θ**) and prohibitively long calculations. A second approach, and the one we implemented, is to carry out the full integration implied by the summation in Eq. (7). We used a quad-precision Romberg integration routine based upon the Numerical Recipes’ qromb routine [16].

The second step is to calculate **F**(**θ**)^{-1}. Prior to the calculation, the eigenvalues of **F**(**θ**) were obtained to determine the stability of the inverse. A set of eigenvalues for the same set of parameters as the results in Section 2 (SBP of 40, vector size of 1024 elements, and an object support size of 161 elements) is plotted in Fig. 4 as a function of the FIM eigenvalue index. From Fig. 4, it can be seen that the FIM eigenvalue behavior as a function of the FIM eigenvalue index is quite similar to the PSWF eigenvalue behavior as a function of the PSWF index; specifically, the eigenvalues are close to one for values of the FIM eigenvalue index ≳ SBP and decay at approximately the same rate as do the PSWF eigenvalues for larger values of the index. This similarity can be explained by looking again at Eq. (7). It can be seen that this equation is a discrete approximation to the convolution of a sinc function with itself, which results in a sinc function. Since each row of the FIM is the previous row circularly-shifted by one pixel, the FIM is just a discrete version of the convolution operator in Eq. (1). Thus the eigenvectors and eigenvalues of the FIM are just PSWFs and their eigenvalues. This implies that, even though the FIM is invertible for any (finite) discretization of the interval [-x_{o},x_{o}], the condition number of the FIM approaches infinity as the number of elements in **α** approaches infinity. Thus unbiased estimates of o(x) with finite variances do not exist. However, this does not necessarily imply that unbiased estimators of O(f) do not exist, at least for some set of frequencies. The existence of unbiased estimators of O(f) depends upon the properties of the transformation of **F**(**θ**)^{-1} into the inverse of the FIM corresponding to O(f). This is discussed next.

The third (and final) step to calculating the desired CRBs associated with any unbiased estimate of O(f) is to transform **F**(**θ**)^{-1} into the inverse of the FIM corresponding to O(f) for any desired set of frequencies and extract the CRBs that are along its diagonal. This transformation is achieved by premultiplying **F**(**θ**)^{-1} by **G**, the Jacobian of the Fourier transform operation for the desired set of frequencies, and postmultiplying by **G**
^{T}. Because the eigenvalues of **F**(**θ**) converge to zero, and because computer calculations have finite precisions, the matrix product **GF**
_{n}(**θ**)^{†}
**G**
^{T} must be computed in a way to clearly illustrate at what frequencies, if any, unbiased estimates of O(f) can be obtained with finite variances. Our approach is to calculate the sequence of matrix products {**GF**
_{n}(**θ**)^{†}
**G**
^{T}} as a function of n=1,2,…, where **F**
_{n}(**θ**)^{†} is the pseudoinverse of **F**(**θ**) associated with the n largest eigenvalues of **F**(**θ**), and to examine the associated CRBs as a function of n. We believe that, if the sequence of CRBs at a given frequency are a non-converging and increasing function of n, it is plausible to assume that they approach infinity as n approaches infinity. Because of the finite precision of computer calculations, our results are limited to approximately thirty orders of magnitude; however, for all practical purposes, the CRBs are infinite at this point anyway.

Based upon the discussion in the previous paragraphs, the unbiased CRBs for estimating O(f) as a function of the number of FIM eigenvalues included in the calculation of **F**
_{n}(**θ**)^{†} are displayed in the movie in Fig. 5 for the same set of parameters as used for generating Fig. 4. Each frame of the movie corresponds to adding one more eigenvalue to the calculation of **F**
_{n}(**θ**)^{†} than in the previous frame. Two important observations can be made from this movie. The first is that none of the CRBs appear to be converging to a finite value. Clearly, the CRBs outside the measurement bandwidth are dramatically increasing in value as a function of n. Even the CRBs inside the measurement bandwidth do not appear to be converging. This slow increase in the CRBs inside the measurement bandwidth is undoubtedly due to the “leakage” of the large CRB values outside the measurement bandwidth into the measurement bandwidth brought about by the Fourier-domain correlations enforced by the support constraint [17]. The second observation is that the CRB plots in Fig. 5 are remarkably similar to the noise variance plots of Fig. 3, as expected, based upon the discussion of Fig. 4.

Because of the preceding CRB analysis as well as the results in Section 2, it seems clear that no unbiased estimate of O(f) for any superresolved frequency is possible with finite variance. Because values of O(f) inside the measurement bandwidth are available from the measurement directly, they can be estimated in an unbiased manner with finite variances merely by taking as their estimates the measurement values themselves. However, no unbiased superresolution with finite variance is possible. For this reason, biased estimators are explored in the next section, leading to the concepts of primary and secondary superresolution.

## 4. Primary superresolution

Although unbiased superresolution with finite variances is not achievable, in this section it is shown that essentially unbiased superresolution with finite variances is achievable. Because biases are a function of the estimator, first a PSWF-based estimator is proposed and its use justified. Next, a quantitative definition of the phrase “essentially unbiased” is derived, and the primary superresolution component of the superresolved spectrum is defined as that portion that is essentially unbiased. Finally, some properties of primary superresolution are presented.

To explore biased superresolution, an estimator must be specified. We choose to use a finite-summation version of the PSWF-based inverse problem formulation in Eq. (2) for three reasons. The first is that, for any finite upper limit M to the sum in Eq. (2), the PSWFs have the most energy in the measurement bandwidth of any basis set [18]. The benefit this brings to the inverse problem can be seen from Eq. (3), where it is shown that the noise variances are inversely proportional to the PSWF eigenvalues. Because these eigenvalues are the energies of the PSWFs inside the measurement bandwidth, the noise levels in the superresolved spectrum for any set of M basis functions are minimized by using PSWFs as the basis set. The second reason is that the similarities between the pseudo-inverse-based CRBs shown in Fig. 5 and the variances for this estimator shown in Fig. 3 imply that the PSWF-based estimator is in some sense qualitatively optimal. The third reason is that a closed-form solution to the inverse problem permits straightforward analysis of the bias properties of superresolution.

To begin the analysis, consider Eq. (2). Clearly the bias properties of this (and any other) estimator are a function of the object being estimated. However, if the biases in the superresolved spectrum produced by this estimator could be characterized in a way that is independent of the object being estimated, the result would be more general and thus more useful. Such a characterization is possible and is derived next based upon the frequency behavior of the PSWFs. Because the estimator in Eq. (2) generates an estimate of O(f) using linear combinations of PSWFs, it follows that PSWFs with minimal energy at a given frequency contribute little to O_{e}(f) at that frequency. In addition, the movie of the PSWFs in Fig. 1 shows that the energies of superresolving PSWFs are negligible in a connected region outside of but adjacent to the measurement bandwidth, and that the size of this region is an increasing function of the PSWF index. The combination of these two facts implies that the bias at a given superresolved frequency is a function of the energies at that frequency contained in PSWFs not included in the summation in Eq. (2) and that this bias is a decreasing function of the summation upper limit. At a given superresolved frequency, if the energies of the PSWFs not included in the estimator are negligible, the superresolved spectrum at this frequency can be said to be essentially unbiased.

To quantify this discussion, the sums of the energies of the PSWFs included in the inverse problem are calculated as a function of superresolved frequency and PSWF index upper limit. To present the results with increased clarity, these sums were divided by the total energy of the PSWFs at each superresolved frequency to show the fraction of energy included. The results of these calculations are shown in the movie in Fig. 6, where each frame of the movie is for a given superresolved frequency. The PSWFs used for this movie are the PSWFs used to generate the movie in Fig. 1, but the general shapes of the curves in Fig. 6 are similar for all SBPs. The impact resulting from the fact that the PSWFs are not identical for all SBPs will be discussed later in this section. Notice that the PSWF energy included in the inverse problem for a given superresolved frequency is essentially zero for low values of the PSWF index, then increases rapidly as a function of PSWF index and asymptotes to one. The plots clearly show a “knee in the curve” that is arrived at when approximately 98% of the PSWF energy is included. For this reason, we call the superresolved spectrum essentially unbiased when more than 98% of the total PSWF energy is included in the estimation of the object’s spectrum at that point. Although the 98% point is somewhat arbitrary, the conclusions to be drawn from this definition of primary superresolution are relatively insensitive to the exact percentage used in the definition.

Some properties of primary superresolution will now be discussed. First, as noted previously, the maximum spatial frequency where primary superresolution is achieved is an increasing function of the maximum PSWF index value allowed in the inverse problem. Because the maximum index value chosen is a function of the data SNRs due to the noise amplification process described by the summation in Eq. (3), this indicates that the amount of primary superresolution is a function of the data SNRs. Additional properties of primary superresolution can be seen by examining Fig. 7, where the amount of primary superresolution is plotted as a function of the inverse of the square root of the PSWF eigenvalues. As noted in the discussion following Eq. (3), the SNRs of the estimated Fourier spectrum are related to the inverses of the square roots of the PSWF eigenvalues, so Fig. 7 shows the amount of primary superresolution possible as a function of the recovery SNR. The vertical axis of Fig. 7 is the amount of primary superresolution achieved and is expressed in terms of the number of degrees of freedom added to the image reconstruction. In this context, a degree of freedom refers to the classical sampling degree of freedom that is defined as the inverse of two times the maximum frequency included in the image reconstruction [8].

Several additional properties of primary superresolution can be deduced from Fig. 7. The first is that, for SBPs ≳10, the amount of primary superresolution is essentially independent of the SBP. In other words, this means that primary superresolution is additive in the sense that, regardless of the object support size or the size of the Fourier-domain measurement bandwidth, as long as the SBP is not on the order of one, the amount of primary superresolution possible is independent of the SBP and adds degrees of freedom. The additive nature of superresolution has been discussed previously [19], [20]. The second observation is that the number of degrees of freedom added is a logarithmic function of the data SNR, a conclusion also reached previously by a number of researchers [1], [21]. The third observation is surprising, however. Notice that the amount of primary superresolution possible for poorly-resolved objects (i.e., for objects whose SBPs are on the order of one) is less than that possible for well-resolved objects. This conclusion apparently contradicts previous results [1]. We note that the results in [1] are in terms of the average total superresolution possible, not just primary superresolution; however, it will be shown in Section 5 that less superresolution is possible for poorly resolved objects than for well-resolved objects even when the total amount of superresolution is considered. The apparent contradiction between our results and the results in [1] is resolved by clarifying how superresolution is measured. We measure superresolution in terms of the number of degrees of freedom *added* to the reconstructed object. If, however, superresolution is viewed in terms of the *percentage increase* in the number of degrees of freedom, as is done in [1], the increase in resolution for poorly-resolved objects is much greater than for well-resolved objects because poorly-resolved objects start out with very little information content. We argue that the best way to view primary superresolution is in terms of adding degrees of freedom, rather than in terms of percentage increase in the number of degrees of freedom, because primary superresolution is an inherently additive phenomenon. In the next section, where secondary superresolution is explored, it will be shown that an additive degrees-of-freedom approach is also an appropriate way to view superresolution in the interior of the support region; however, the definition of degree of freedom changes. At the edges of the support region, secondary superresolution will be shown to be a function of the SBP and thus not additive.

## 5. Secondary superresolution

The concept of secondary superresolution is discussed in this section and its properties are compared and contrasted with the properties of primary superresolution. The amount of superresolution possible when including both the primary and secondary components of the superresolved spectrum is evaluated using a PSWF-based degrees-of-freedom approach. It is shown that secondary superresolution is spatially variant with significantly higher resolutions possible at the edges of the support region as compared to interior to the support region.

Secondary superresolution is defined as all of the superresolved spectrum that is not essentially unbiased. Another way to express this fact is to call secondary superresolution either biased or incomplete. For this reason, the classical degrees-of-freedom relationship relating the number of independent segments in the image domain to the maximum frequency included in the reconstruction process no longer holds. To illustrate this fact, consider the superresolved Fourier spectra of a simulated 1-D triple star shown in Fig. 8. The triple star was created in a 1024 element vector and consists of three equal-magnitude point sources located at elements 437, 512, and 552. The support region used for the calculations was [432,592]. The true Fourier amplitudes are plotted along with the reconstructed Fourier amplitudes using Eq. (2) with an upper limit to the summation of 46. The PSWFs from Fig. 1 were used for this reconstruction. Because these PSWFs were created using a SBP of 40, the upper summation limit of 46 means that six superresolving PSWFs were included in the reconstruction process. Although the data in Fig. 8 are noise free, including six superresolving PSWFs corresponds to assuming a data SNR of approximately100.

Notice that the region of primary superresolution is indeed essentially unbiased, while the region of secondary superresolution is strongly biased. Although the normalized frequency axis in Fig. 8 extends only to twice the measurement bandwidth, the superresolved Fourier
spectra is not bandlimited because the structure of superresolving PSWFs is that of decaying sinusoids. The classical degree-of-freedom model states that the number of degrees of freedom in the image is equal to 4x_{o}f_{s}, where f_{s} is the largest frequency included in the image reconstruction. Because f_{s}=∞ for any reconstruction where superresolving PSWFs are included, the classical degree-of-freedom approach states that there also are an infinite number of degrees of freedom in the reconstruction. Clearly this is in error. Thus the classical degree-of-freedom model is not valid to characterize the resolution properties of secondary superresolution. An alternate model must be used.

It has been shown previously that a PSWF of index m has m zeros in the support region used to generate the PSWFs [9]. Based upon this fact, it has been proposed that the resolution provided by the m^{th} PSWF can be characterized by the average distance between these zeros [1], [6]. Such a definition of resolution agrees with the classical Rayleigh resolution limit when m is equal to the SBP. Based upon this definition of resolution, the number of degrees of freedom in a reconstructed image using PSWFs is equal to one plus the number of PSWFs included in the reconstruction. This definition is consistent with the general concept of resolution in the classical sense and does not predict infinite resolution when including an infinite number of frequencies as is the case for the sampling degree of freedom definition.

Using the PSWF definition of degrees of freedom, the average increase in the resolution of a reconstruction brought about by including both the primary and secondary superresolution components in the reconstruction process is plotted in Fig. 9 as a function of the inverse of the square root of the PSWF eigenvalue and for several values of the SBP. There are several properties of the increased degrees of freedom brought about by including secondary superresolution in the reconstruction process worth noting. The first is that poorly-resolved objects still benefit the least from the superresolution process as measured by increased degrees of freedom. The second is that the amount of increase in degrees of freedom is no longer independent of the SBP for SBP ≳ 10. In fact, the number of degrees of freedom in a superresolved reconstruction of an object is an increasing function of the number of degrees of freedom in the measurement. Thus the inclusion of secondary superresolution makes the superresolution process more than additive in terms of our definition of additive superresolution given in Section 4. For example, primary superresolution produces one additional degree of freedom when the data SNR is 100 for all SBP values ≳ 10, while the combination of primary and secondary superresolution produces four additional degrees of freedom when the SBP is 10, and five additional degrees of freedom when the SBP is 40. This also indicates that secondary superresolution is a valuable component of the overall superresolution process. Finally, we have plotted in Fig. 9 the predicted increase in degrees of freedom found in [1], Table 1, for a SBP of 1.27 in order to compare our results to theirs. Notice the good agreement, especially considering that the results in [1] have a precision of only one significant digit.

However, defining the resolution of a given superresolving PSWF in terms of the average spacing of its zeros across the entire support region implies that the zeros are evenly spaced throughout the support region. As can be seen for a representative superresolving PSWF in Fig. 10, this is not the case. The zeros are relatively evenly spaced in the center of the support region, but their spacings decrease rapidly as one approaches the edge of the support. This implies that superresolving PSWFs produce higher resolutions near the edges of the support region and lower resolutions away from the edges of the support region than predicted by Fig. 9. To take the space-dependent resolution properties of superresolving PSWFs into account, the support region was divided into two regions, the center 75% and the remaining 25% at the edges. The average amount of total superresolution achieved in the center 75% of the support region is plotted in Fig. 11 and can be seen to be significantly less than predicted by Fig. 9 for the entire support region. To compare directly the secondary superresolution results to the original degrees of freedom in the image, the results plotted in Fig. 11 correspond to the increase in the number of degrees of freedom that would occur in the entire support region if the resolution calculated using the center 75% of the support region is the same throughout the entire support region. In addition, notice that the amount of total superresolution is essentially independent of the SBP for SBP t 10, just as is the case for primary superresolution, and thus is additive by our definition. As before, the amount of total superresolution achievable is less for SBP values on the order of one. A difference from the primary superresolution results is that the amount of superresolution is a little larger when including secondary superresolution in the reconstruction process.

The average amount of total superresolution achieved at the edges of the support region for the remaining 25% of the support region is shown in Fig. 12. Notice the large increase in the numbers of degrees of freedom at the edges of the support. As for the results in Fig. 11, the numbers in Fig. 12 correspond to the increase in the degrees of freedom that would occur in the entire support region if the resolution at the edges of the support region is the same throughout the support region.

To put the primary and secondary superresolution results in perspective, it is useful to consider the relative percentage increases in resolution for a specific case. To this end, consider an image that has a SBP of 40. If the data SNR allows the use of all the PSWFs whose eigenvalues are greater than 10^{-4} (that is, if the data SNR ≈ 100), the percent increase in the resolution throughout the entire support region due to primary superresolution would be 2.5% because primary superresolution only provides a single additional degree of freedom. The resolution increase including both primary superresolution and secondary superresolution in the center region of the support would be 5%, a factor of two increase over what primary superresolution provides, but still small. Finally, at the edges of the support region, the resolution increase including both primary and secondary superresolution would be 30%, a large and noticeable increase. To demonstrate visually the impact of these numbers, three noise-free reconstructions of the triple star whose Fourier amplitude spectrum is plotted in Fig. 8 are displayed in Fig. 13. The first reconstruction uses just the Fourier data inside the measurement bandwidth, the second reconstruction uses both the measured Fourier data and the primary component of the superresolved spectrum, and the third reconstruction uses all the superresolved spectrum along with the measured Fourier data. The two rightmost components of the triple star are contained inside the center 75% of the support; as a result, the amount of superresolution possible is predicted to be ≤ 5% for all three reconstructions and thus essentially negligible. The reconstructions in Fig. 13 confirm this prediction. However, the leftmost component of the triple star is near the edge of the support region and benefits from the increased resolution brought about by secondary superresolution. As predicted by the plots in Fig. 7, the primary superresolution reconstruction of this component differs very little from the no-superresolution reconstruction. However, as predicted by the plots in Fig. 12, the reconstruction that uses all the superresolved spectrum has an increased resolution of ~30%. The increase in resolution is greater on the side of the component near the edge of the support as compared to the side nearer the center of the support because of the decreasing zero spacings of PSWFs as one moves toward the support boundaries. Although the increased resolution is encouraging, it can be seen that the reconstruction using all the superresolved spectrum produces less accurate estimates of the relative magnitudes of the triple star components.

## 6. Discussion and future work

The results in the preceding sections provide insight into the properties of superresolved spectra. It was shown that these properties can be seen more clearly if the superresolved spectrum region is decomposed into two separate regions. The first region is connected to and outside of the measurement bandwidth and is the region where the essentially unbiased, or primary, component of the superresolved spectrum resides. The second region is outside of both of these regions and connected to the primary superresolution region and is the region where the biased, or secondary, component of the superresolved spectrum resides. Although the secondary superresolution region is theoretically infinite in size, the data SNR limits it to be finite.

Primary superresolution is additive in nature in that the amount of primary superresolution is independent of the SBP as long as the SBP ≳ 10 and depends only on the data SNR. For SBPs on the order of one, the amount of primary superresolution for a given data SNR is less than for larger values of the SBP. The increase in resolution in the image domain is constant across the support region. In addition, for objects with SBPs ≳ 10, the amount of primary superresolution achievable with realistic data SNRs is on the order of a few percent or less.

Secondary superresolution produces increases in resolution in the support region that are a function of the location in the support region. Across most of the support region away from its edges the resolution increase due to secondary superresolution is essentially constant and has many of the same properties as primary superresolution (i.e., the amount of resolution increase is independent of the SBP when the SBP ≳ 10, and for smaller SBPs, the amount of resolution increase is less). In addition, although secondary superresolution in the interior of the support region adds more resolution to the primary resolution results, it is still on the order of a few percent for data with SBPs ≳ 10 and data SNRs ≈ 100. These properties are due to the portion of the superresolved spectrum in the (small) region of the secondary-superresolution Fourier spectrum adjacent to the primary superresolution region. However, the amount of resolution increase at the edges of the support region due to secondary superresolution is much greater than in the interior of the support region. For example, for data with a SBP of 40 and an SNR ≈ 100, the resolution increase averaged over the edges of the support region is ~30%. These higher resolutions are associated with the long “tails” of the secondary superresolution spectrum in the Fourier domain.

Based upon these properties of primary and secondary superresolution, it is clear that both schools of thought regarding the achievability of meaningful superresolution can be correct depending upon the metric used to quantify the phrase “meaningful superresolution”. If the metric is the obtainment of more than a few percent resolution increase across the entire support region for data with realistic SNRs and SBPs ≳ 10, meaningful superresolution is not possible. If the metric is obtaining more than a few percent resolution increase at the edges of the support for any SBP, meaningful superresolution is possible.

This latter point is worth emphasizing. It is well known that the human visual system keys off of edges in an image. Therefore, if the object under consideration is an object with well-defined edges on a black background (such as is common in astronomical imaging), meaningful superresolution appears to be quite possible and worth pursuing. This possibility has been noted previously [7]. On the other hand, if the edges of interest are interior to the object support, little edge enhancement is possible.

The results presented in this paper have been independent of the actual object being imaged and superresolved. It is important to determine how effective primary and secondary superresolution are in improving the resolution of images in terms of looking at actual examples. Our initial results (Fig. 13) indicate that using secondary superresolution in the reconstruction process can increase resolution at the edges of the support region but at the cost of decreased radiometric accuracy. In addition, our results only included noise in the sense of limiting the number of PSWFs that are included in the reconstruction. We expect that noisy data will require additional constraints in the reconstruction process such as Fourier-domain regularization. For these reasons, we are continuing our research into the image domain effects of secondary superresolution and will report on the results a future publication.

## Acknowledgments

The authors thank the Air Force Office of Scientific Research for their financial support of this research; specifically, CLM is supported by grant LRIR02DE01COR and DWT is supported by grant F49620-02-1-0107 from the Directorate of Mathematics and Space Sciences. The authors also thank Dr. Latham and Mr. Tilton for providing technical assistance and the software used to generate the prolate spheroidal wave functions employed in this research project.

## References and links

**1. **M. Bertero and C. De Mol, “Super-resolution by data inversion,” in *Progress in Optics XXXVI*,
E. Wolf, ed. (Elsevier, Amsterdam, 1996), pp. 129–178.

**2. **S. Bhattacharjee and M. K. Sundareshan, “Mathematical extrapolation of image spectrum for constraint-set design and set-theoretic superresolution,” J. Opt. Soc. Am. A **20**, 1516–1527 (2003). [CrossRef]

**3. **B. R. Hunt, “Super-resolution of images: algorithms, principles, and performance,” Int. J. Imaging Syst. Technol. **6**, 297–304 (1995). [CrossRef]

**4. **H. Liu, Y. Yan, Q. Tan, and G. Jin, “Theories for the design of diffractive superresolution elements and limits of optical superresolution,” J. Opt. Soc. Am. A **19**, 2185–2193 (2002). [CrossRef]

**5. **V. F. Canales, D. M. de Juana, and M. P. Cagigal, “Superresolution in compensated telescopes,” Opt. Lett. **29**, 935–937 (2004). [CrossRef] [PubMed]

**6. **C. K. Rushforth and R. W. Harris, “Restoration, resolution, and noise,” J. Opt. Soc. Am. **58**, 539–545 (1968). [CrossRef]

**7. **J. J. Green and B. R. Hunt, “Improved restoration of space object imagery,” J. Opt. Soc. Am. A **16**, 2859–2865 (1999). [CrossRef]

**8. **B. R. Frieden, “Evaluation, design, and extrapolation methods for optical signals based on the use of prolate functions,” in *Progress in Optics IX*,
E. Wolf, ed. (North-Holland, Amsterdam, 1971), pp. 313–407.

**9. **D. Slepian and H. O. Pollak, “Prolate spheroidal wave functions, Fourier analysis and uncertainty – I,” Bell Syst. Tech. J. **40**, 43–63 (1961).

**10. **H. J. Landau and H. O. Pollak, “Prolate spheroidal wave functions, Fourier analysis and uncertainty – II,” Bell Syst. Tech. J. **40**, 65–84 (1961).

**11. **M. Bertero and E. R. Pike, “Resolution in diffraction-limited imaging, a singular value analysis I. The case of coherent illumination,” Opt. Acta **29**, 727–746 (1982). [CrossRef]

**12. **M. C. Roggemann and B. Welsh, *Imaging Through Turbulence* (CRC Press, Boca Raton, 1996), p. 49.

**13. **W. P. Latham and M. L. Tilton, “Calculation of prolate functions for optical analysis,” Appl. Opt. **26**, 2653–2658 (1987). [CrossRef] [PubMed]

**14. **B. Porat, *Digital Processing of Random Signals, Theory and Methods* (Prentice-Hall, Englewood Cliffs, 1994), pp. 65–67.

**15. **R. C. Gonzalez and R. E. Woods, *Digital Image Processing* (Addison-Wesley, Reading, 1992), ch. 5.

**16. **W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, *Numerical Recipes in FORTRAN, 2 ^{nd} ed.*, (Cambridge Press, Cambridge, 1996), pp. 134–135.

**17. **C. L. Matson, “Variance reduction in Fourier spectra and their corresponding images with the use of support constraints,” J. Opt. Soc. Am. A **11**, 97–106 (1994). [CrossRef]

**18. **H. J. Landau and H. O. Pollack, “Prolate spheroidal wave functions, Fourier analysis and uncertainty – III: the dimension of the space of essentially time-and band-limited signals,” Bell Syst. Tech. J. **41**, 1295–1336 (1962).

**19. **P. J. Sementilli, B. R. Hunt, and M. S. Nadar, “Analysis of the limit to superresolution in coherent imaging,” J. Opt. Soc. Am. A **10**, 2265–2276 (1993). [CrossRef]

**20. **C. L. Matson, “Fourier spectrum extrapolation and enhancement using support constraints,” IEEE Trans. Sig. Process. **42**, 156–163 (1994). [CrossRef]

**21. **Y. L. Kosarev, “On the superresolution limit in signal reconstruction,” Sov. J. Commun. Technol. Electron. **35**, 90–108 (1990).