## Abstract

We analyze the problem of optical superresolution (OSR) of a one-dimensional (1D) incoherent spatial signal from undersampled data when the support of the signal is known in advance. The present paper corrects and extends our previous work on the calculation of Fisher information (FI) and the associated Cramer-Rao lower bound (CRB) on the minimum error for estimating the signal intensity distribution and its Fourier components at spatial frequencies lying beyond the optical band edge. The faint-signal and bright-signal limits emerge from a unified noise analysis in which we include both additive noise of detection and shot noise of photon counting via an approximate Gaussian statistical distribution. For a large space-bandwidth product, we derive analytical approximations to the exact expressions for FI and CRB in the faint-signal limit and use them to argue why achieving any significant amount of unbiased bandwidth extension in the presence of noise is a uniquely challenging proposition. Unlike previous theoretical work on the subject of support-assisted bandwidth extension, our approach is not restricted to specific forms of the system transfer functions, and provides a unified analysis of both digital and optical superresolution of undersampled data.

© 2009 Optical Society of America

## 1. Introduction

Beating the detector sampling limit to achieve improved resolution from undersampled image data has been shown to be eminently possible, particularly if the detector SNR is not too small. The sought high-resolution information can be encoded, in full agreement with the generalized sampling theorem (GST), into low-resolution image data, as typically done via progressively larger sub-pixel image shifts in a sequence of images [1]. However, such digital superresolution (DSR) does not yield any extension of the spatial bandwidth for the reconstructed image beyond the optical diffraction limit of the imaging system. Indeed, the maximum achievable DSR is strictly upper-bounded by this diffraction limited bandwidth, as one of us argued in [2] with the help of the GST.

It is essential to have some meaningful *a priori* knowledge in addition to the image data in order to perform true bandwidth extension, or optical superresolution (OSR), by means of data inversion. In fact, most algorithms that claim to achieve OSR incorporate prior knowledge of one form or another, often implicitly, into their data inversion protocol. Enforcing a finite-support constraint in the physical domain is the essential basis of certain OSR algorithms [3, 4, 5]. Examples of other forms of prior knowledge include positivity, as in maximum-likelihood (ML) and other positive-constraint-type algorithms [6, 7, 8, 9, 10, 11, 12], smoothness, as in the maximum entropy method (MEM) [13, 14], and morphology, as in CLEAN [15, 16] and sampling-adjusted deconvolution [17, 18, 19]. A good discussion of the important astronomical image-deconvolution algorithms that achieve OSR is given in [20].

In a conference paper [2], one of us presented preliminary results of a calculation for one dimensional (1D) OSR when the support of the original signal, which we shall henceforth call the *object*, is known in advance. This work was an attempt to generalize a large body of previous theoretical literature [21, 22, 23, 24, 25, 26] over the past four decades that was typically devoted to 1D OSR under a support constraint in a clear-aperture coherent imaging system for which the so-called prolate spheroidal wave functions (PSWF) provide an excellent eigenbasis for formulating the OSR problem. By contrast, our work is not tied to any specific imaging system, and makes use of sampling-type sinc functions to represent the finite-support object in the spatial-frequency domain.

In the present paper, we correct and greatly expand upon our previous work [2], with carefully obtained new results. The noise model we employ allows for both the additive Gaussian noise of CCD detection and the Possion noise of photon counting via an approximate Gaussian probability density function (PDF). We start with a sequence of sub-pixel-shifted low-resolution (LR) 1D signals, which we shall henceforth loosely call *images*, in which information about the object out to the diffraction limit of the imager can be encoded. We can recover still higher object frequencies by exploiting prior knowledge of the object support, as we shall demonstrate here. A finite object support couples these latter, superresolving (SR) frequencies lying beyond the imager’s optical bandwidth into that bandwidth so the image data, in principle, carry information about the SR frequencies. We validate this idea quantitatively by means a Fisher-information (FI) based statistical analysis of the maximum fidelity with which SR frequencies can be estimated. Previous authors [26, 27] have discussed OSR and DSR from the perspective of FI and the associated Cramér-Rao lower bound (CRB), but our work is the first to provide a unified description of the two aspects of superresolution via this statistical approach.

We present in Sec. 2 a brief discussion of the imaging model in which we relate the image data to a spatial-frequency-domain integral of the object spectrum and the imager’s optical and digital transfer functions. In Sec. 3 we introduce the essential ingredient of our approach, namely a sampling-type expansion in the spatial-frequency domain involving “critical” samples of the spatial spectrum of a finite-spatial-support object brightness function. The next section introduces the notion of a complex Hermitian FI matrix with respect to a set of complex paramaters to be estimated and its inverse relation to the CRB on the sum of variances of the real and imaginary parts of the parameters. Here we also derive an expression for the FI matrix for the superresolution problem in the presence of a mixed Poisson-Gaussian noise of CCD detection. In Sec. 5, we evaluate this expression and its matrix inverse approximately for the case of a large space-bandwidth product (SBP), and demonstrate an approximate permanence of the CRB, independent of the degree of sought OSR, for the data-resolved (DR) frequencies lying well within the diffraction limit. Here we also establish the mathematical basis for the extreme sensitivity of estimation of the SR frequencies lying outside the diffraction limit to detection noise. In Sec. 6, we present the results of a numerical evaluation of our CRB calculations and discuss them by relating them to the physical properties of the OSR problem in the spatial-frequency domain. Similar considerations in the image domain form the content of the next section, where we shall confirm a previously observed [5, 26] space-variant character of OSR with a higher fidelity of OSR near the support boundaries than well in the interior. Some concluding remarks appear in Sec. 8.

## 2. Imaging model

For a space-invariant imaging system with optical PSF given by *h*
_{0}(*x*), a continuously defined image, *g*(*x*), may be expressed in terms of the object intensity function, *f* (*x*), as the convolution

A pixelated sensor array on which an image is typically formed integrates and samples the continuous image to generate an array of discrete numbers, one for each pixel, that may be expressed in the form

$$=\int dx\int dx\prime p(x-{x}_{k}){h}_{0}(x-x\prime )f\left(x\prime \right)+{n}_{k},$$

where *p*(*x*) is the pixel sensitivity function, *x _{k}*=

*kδw*is the position coordinate of the center of the

*k*th pixel in terms of the pixel pitch,

*δw*, and

*n*is the additive read-out noise at the

_{k}*k*th pixel. A complete noise model for the sensor must also include the Poisson noise of counting for which the first term on the right-hand side (RHS) must be regarded as the mean count rate at the

*k*th pixel.

For undersampled sensor data, any improvement of resolution amounts first to a DSR out to the optical diffraction limit followed by an OSR that extends the resolution still further. We treat DSR and OSR in an integrated fashion as part of a single program of realizing improved resolution. We work with a sequence of LR image datasets with sub-pixel shifts that are chosen to be regularly spaced to accord maximum noise robustness [1, 27] to DSR. In terms of the product of the optical transfer function (OTF), *H*
_{0}(*u*), and the detector-pixel transfer function (PTF), *P*(*u*), we may express the intensity at the *k*th pixel of the *m*th image as [2]

where *t _{m}* is the sub-pixel shift of that image and

*n*

^{(m)}

*is the additive noise at its*

_{k}*k*th pixel.

For idealized sensor arrays with no dead space between successive pixels, the PTF has the form, *P*(*u*)=*Cδ w*sinc(*uδw*), where the sinc function is defined as sinc(*x*)≡sin*πx*/(*πx*) and *C* is the quantum efficiency of the sensor pixels. We shall assume *C*=1 in the rest of the paper. For a clear aperture of optical bandwidth *B*
_{0}, the incoherent OTF vanishes outside the spatial-frequency range |*u*|<*B*
_{0} where its functional form is (1-|*u*|/*B*
_{0}). Expression (3) may thus be rewritten as the following integral over the normalized frequency *u*̄=*u*/*B*
_{0}:

where the spatial variables are now expressed in units of the critical sampling interval, *δx _{c}*=1/(2

*B*

_{0}), specifically

*x*̄

*=*

_{k}*x*/

_{k}*δx*and

_{c}*t*̄

*=*

_{m}*t*/

_{m}*δx*. The symbol

_{c}*χ*denotes the ratio,

*δw*/

*δxc*, of detector-sampling and critical optical-sampling intervals. For

*χ*<1, the detector oversamples the image while for

*χ*>1 it undersamples the image.

## 3. Finite object support and the object spectrum

When *f* (*x*) is supported compactly on the interval [-*L*,*L*], it has a Fourier-series expansion,

where

The Fourier transform (FT) of this expansion yields the full object spectrum, *F*(*u*), in terms of its spatial-frequency samples, *F*
* _{ℓ}*=

*F*(

*u*=

*ℓ/2L*), as

in a manner analogous to the usual Shannon-Nyquist expansion of a band-limited function in terms of its critical spatial samples. The implication of this expansion for the imaging equation (4) is that the information in the spatial-frequency samples outside the measurement band, [-*B*
_{0},*B*
_{0}], is aliased into the measurement band and thus in principle contained in the acquired data. The finiteness of the support is directly responsible for this aliasing. Note, however, that because the width of the sinc function in Eq. (7) is of order 1/(2*L*), the maximum bandwidth extension that can be expected is also of that order. This is consistent with previous results on unbiased, primary OSR [28, 25].

From the present perspective, of DSR and OSR, only OSR is technically a de-aliasing problem. In DSR, high-frequency information out to the optical passband is already present without having to be aliased into the PTF, while in support-based OSR, information at still higher frequencies aliased into the optical passband due to a compact object support must ultimately be de-aliased.

Use of expression (7) in Eq. (4) yields the following image data value in terms of the object spatial-frequency samples:

*Q*=2*LB*
_{0} is the space-bandwidth product (SBP) and *u*̄ has been simply relabeled as *u*.

Because of the decaying tail of the sinc aliasing function in Eq. (9), the inverse problem of recovering such frequencies from noisy image data becomes increasingly unstable as frequencies increasingly beyond the band edge are considered. One must therefore truncate the sum (8) to a maximum magnitude, *ℓ _{max}*, of the index,

*i.e*., to the range -

*ℓ*≤

_{max}*ℓ*≤

*ℓ*. The ratio

_{max}*S*=

*ℓ*/

_{max}*Q*, when it exceeds 1, is the sought OSR factor. In the absence of noise, the (2

*SQ*+1) object spectral components contained in the truncated sum can be solved for by acquiring at least (2

*SQ*+1)/N independent LR images, each with

*N*pixels.

## 4. Fisher information and estimation error for fourier samples

We now consider the error in the unbiased estimation of the Fourier samples of a 1D object supported on the interval, [-*L*,*L*]. A simply calculable lower bound on this error is the CRB, determined by the diagonal elements of the inverse of the FI matrix [29].

A lower bound on the sum of variances for estimating the real and imaginary parts of a complex parameter like *Fℓ*, sometimes called its pseudo-variance, is given by the corresponding diagonal element of the inverse of a complex Hermitian FI matrix. The *ℓℓ*′ element of the FI matrix is defined as a statistical self-average in two equivalent ways,

where *P*=*P*({*g*
^{(m)}
* _{k}*}|{

*F*,

_{ℓ}*F**

*}) denotes the PDF of the image data set, {*

_{ℓ}*g*

^{(m)}

*k*}, corresponding to the given object. The reality of the object implies

*F**

*=*

_{ℓ}*F*-

*ℓ*. It is therefore sufficient to consider only non-negative

*ℓ*,

*ℓ*′ values when considering the matrix elements of

**J**.

Most superresolution algorithms make in effect a *biased* estimation of SR frequencies, so our analysis based on using simply the inverse of the FI matrix to provide a lower bound on the variance of any *unbiased* estimation must be regarded as merely suggestive of the actual sensitivity limits on the performance of most practical OSR algorithms. Limits on a biased estimation of the SR frequencies may also be calculated in the FI formalism, but such limits depend critically on the gradient of the estimator bias relative to the parameters being estimated by the specific algorithm and, unlike unbiased-estimation CRBs, are thus algorithm-dependent. It is worth noting that the absence or presence of estimator bias in OSR has been utilized recently to sharpen the discussion of OSR in terms of primary and secondary OSR [26].

If the noise is limited to additive sensor noise alone, then the data PDF consists of a product of Gaussian factors, one for each pixel of each image, with the mean value for the *k*th pixel of the *m*th image given by the first term on the RHS of Eq. (8) and variance *σ*
^{2}
* _{D}*. In this case, as we argue in Appendix A, the average (10) is evaluated easily in terms of the coefficients,

*c*

^{(m)}

*k*ℓ,

In view of the factorized form of the terms inside the sum (11), the FI matrix assumes a compact form,

where *c*
^{(m)}-* _{k}* is the column vector (

*c*

^{(m)}

_{k0},

*c*

^{(m)}

_{k1},…)

*and*

^{T}*c*

^{(m)†}

_{-k}its Hermitian transpose.

When the Poisson noise of counting cannot be neglected, then we must combine it with the additive sensor noise for a more complete noise model. Whenever either *σ _{D}*≫1 or the mean count rate 〈

*g*

^{(m)}

*〉≫1, we may do this approximately via a Gaussian PDF with mean 〈*

_{k}*g*

^{(m)}

*〉 and variance (*

_{k}*σ*

^{2}

*+〈*

_{D}*g*

^{(m)}

*〉). Under such a Gaussian approximation to the exact PDF for the mixed-noise model, which we assume to be valid throughout this paper, the FI (12) is modified by bringing*

_{k}*σ*

^{2}

*inside the*

_{D}*m*,

*k*double sum and replacing it by the sum noise variance,

*σ*

^{2}

*+〈*

_{D}*g*

^{(m)}

*〉, pixel-by-pixel,*

_{k}A detailed justification of these approximations is presented in Appendix A.

Because of the sum form (13) of the FI matrix, its maximum rank is equal to the number of terms in the sum, namely the product of the number of images and the number of detector pixels per image. This is to be expected since the data-based FI matrix contains only as much information as has been acquired through independent measurements.

## 5. The large-Q limit

Let us regard the FI matrix as consisting of four submatrices, two diagonal and two off-diagonal, corresponding to the four possible domains in the (*ℓ*, *ℓ*′) plane of matrix-element labels, namely 0≤*ℓ*, *ℓ*′<*Q*; *Q*≤*ℓ*, *ℓ*′≤*ℓ _{max}*; 0≤

*ℓ*<

*Q*≤

*ℓ*′≤

*ℓ*; and 0≤

_{max}*ℓ*′<

*Q*≤′≤

*ℓ*. The off-diagonal submatrices, which are Hermitian adjoints of each other, represent the cross-sensitivity of the estimation of the SR frequencies on the DR frequencies and vice-versa.

_{max}In the limit of large SBP, *Q*≫1, the aliasing sinc function inside the integral (9) is sharply peaked at *u*=*ℓ/Q* with a width of order 1/*Q*. Hence for the spatial frequencies well inside the passband for which |*Q*-*ℓ*|≫1, we may approximate the integral by replacing the relatively slowly varying factor in the integrand, namely (1-|*u*|)sinc(*χu*/2), by its value at *u*=*ℓ/Q* and extending the limits of integration about *v*≡(*u*-*ℓ*/*Q*)=0 to ±∞. Thus, we have

Using the integral representation,

in Eq. (14), interchanging the order of the two resulting integrations, and then noting that the inner integral may be evaluated in terms of the Dirac *δ* function, we may express the integral in Eq. (14) as

where *t*̃* _{m}*=

*t*̄

*/*

_{m}*χ*is the sub-pixel shift of the mth image as a fraction of an image-pixel width. The integral (16) evaluates as (1/

*Q*)Θ(1-

*χ*|

*k*-

*t*̃

*|/*

_{m}*Q*), where Θ denotes the Heaviside unit step function that takes the value 1 when its argument is positive and 0 otherwise. Putting all these simplifications together, we have the following approximate expression for the coefficient

*c*

^{(m)}

*valid for large*

_{k}*Q*:

Substituting expression (17) in Eq. (11) yields the following approximation for the FI matrix elements correponding to object spatial-frequency samples well inside the optical passband:

$$\times \sum _{m,k}^{\prime}\frac{\mathrm{exp}\left[i\pi \left(\ell -\ell \prime \right)\chi \left(k-{\tilde{t}}_{m}\right)\u2044Q\right]}{{\sigma}_{D}^{2}+\u3008{g}_{k}^{\left(m\right)}\u3009},$$

where the prime on the sum indicates the restriction imposed by the Θ function in Eq. (17) that |*k*-*t*̄* _{m}*|<

*Q*/

*χ*. This constrains the sum for each image m to be over only those pixels

*k*over which the image signal is supported.

#### 5.1. The large-additive-noise limit

The double sum in Eq. (18) may be performed exactly when the additive noise dominates, *i.e*. when *σ*
^{2}
* _{D}*≫〈

*g*

^{(m)}

*k*〉 for all

*m*,

*k*, so the denominator of the summand can be replaced simply by

*σ*

^{2}

*. The double sum then reduces to a simple geometric-series inner sum over*

_{D}*k*running from -

*Q*/

*χ*to

*Q/χ*, which equals (2

*Q*/

*χ*)

*ℓ′, followed by a trivial outer sum over*

_{δℓ}*m*, which equals

*M*, the total number of LR images. This leads to the following approximately diagonal form for the diagonal FI submatrix corresponding to frequency samples inside the passband:

For the off-diagonal FI submatrices, the derivation may once again be performed approximately. The details of the calculation for all of the FI submatrix elements are given in Appendix B, where it is shown that for large *Q*, the elements of the off-diagonal submatrices are essentially zero, at least well in their interior. As a result, regardless of the form of the remaining, diagonal FI submatrix that involves only the SR frequencies, the inverse of the overall FI matrix, when similarly sub-divided into four submatrices, consists of diagonal submatrices that are simply inverses of the corresponding diagonal submatrices of the overall FI matrix. In other words, because of the approximate vanishing of the cross-sensitivity-submatrix elements, the CRBs on the errors of an unbiased estimation of the DR frequencies well within the passband or of the SR frequencies lying well outside are independent of each other for large SBP.

Based on this analysis, the CRBs on the variances of any unbiased estimation of the object spatial frequencies lying within the system passband are simply the reciprocals of the diagonal matrix elements (19),

$$\approx \frac{2Q{\sigma}_{D}^{2}}{{\chi}^{M}{\left(1-\frac{\mid \ell \mid}{Q}\right)}^{2}{\mathrm{sinc}}^{2}\left(\chi \frac{\ell}{2Q}\right)},\mid \ell \mid ,\mid \ell \prime \mid <Q,$$

independent of any estimation of the superresolving frequencies.

Note that the approximate form (20) of the CRB is erroneously divergent for those DR frequencies, *u _{ℓ}*=

*ℓ/Q*, where the detector PTF (

*χ*/2)sinc(

*χuℓ*/2) vanishes, which can happen for

*χ*>2. This divergence may be traced to expression (14) which fails to provide an accurate approximation of the exact expression (9) near such frequencies. The finite width of the support-based aliasing sinc function cannot be neglected in Eq. (9) even in the large-

*Q*limit, so that the exact CRB is in fact finite for such frequencies. The prior knowledge of the object support injects information at the zero crossings of the PTF within the optical passband, thus enabling DSR throughout that band. Without prior information, such object frequencies would be filtered out perfectly by the sensor PTF and could not be estimated.

As we show in Appendix B, the FI submatrix corresponding to the SR frequencies alone has elements given by the spatial-frequency integral,

This submatrix, while not diagonal, has a number of characteristics that distinguish it from the previous submatrix (19). First, all its elements scale as 1/*Q*
^{2} unlike the 1/*Q* scaling of expression (19). Thus, since *Q*≫1, it has much smaller elements, and its inverse, which determines the CRBs, has much larger elements than those given by Eq. (20). The CRBs for these SR frequencies are, in fact, considerably larger than the simple scaling argument suggests, since when *ℓ*, *ℓ*′≫*Q*, the *ℓ*, *ℓ*′ dependence of expression (21) may be approximated as being of a factorized form, (-1)*ℓ*-*ℓ*′/(*ℓ*
*ℓ*′). Such matrices are highly ill-conditioned, since their rank is small (if the factorization were exact for all *ℓ*, *ℓ*′, then the rank would be 1, independent of the dimension of the submatrix). Physically, this is due to the fact that the image data contained within the diffraction passband become increasingly insensitive to – and thus increasingly lacking in information about – increasingly larger SR frequencies. We expect the rise of the CRB for these SR frequencies to become dramatically steep as *ℓ _{max}* is increased well beyond

*Q*.

## 6. Numerical results and discussion of OSR in the spatial-frequency plane

In our numerical evaluations, we considered a support-constrained Gaussian-shaped object of form,

where *w* is a fractional width parameter and *K* may be interpreted as the integrated source flux whenever *w* is small compared to 1. For this object, we now illustrate several essential features of the variation of the CRB for a joint, unbiased estimation of its spatial frequencies for four different values of *Q*, namely 5, 10, 15, and 20. On the plots the spatial-frequency samples are labeled by their index *ℓ*, which takes integer values between 0 and *ℓ _{max}*. In our numerical calculations of the exact expression (13) for the FI and its inverse, the set of jointly estimated samples was enlarged incrementally by adding the next higher-frequency sample pair at a time,

*i.e*., by incrementing

*ℓ*by 1. Three values of

_{max}*ℓ*, namely

_{max}*Q*,

*Q*+1, and

*Q*+2, were considered corresponding to the inclusion of zero, one, and two pairs of samples outside the optical bandwidth, respectively. Incrementing

*ℓ*by one unit corresponds to a fractional bandwidth extension of amount 1/

_{max}*Q*, which is a 20% extension for

*Q*=5 but only a 5% extension for

*Q*=20.

In Figs. 1, we plot the CRB for unbiased estimation of the various object spatial frequency components of a broad Gaussian object (*w*=1) from a single (*M*=1) critically sampled (*χ*=1) LR image for the four values of *Q*. The peak SNR, defined as the ratio of the maximum object pixel intensity to the standard deviation, *σ _{D}*, of the additive noise, was taken to be very small at 0.001, so the additive noise was dominant and the CRBs, according to Eq. (13), essentially independent of the object brightness distribution. Here the CRBs were normalized by dividing out a factor of

*σ*

^{2}

*before plotting them. A number of observations may be made. First, the CRBs for estimating the DR frequencies remain essentially unchanged for moderate to large values of*

_{D}*Q*(Figs. 1(c) and (d)) even as more and more of the SR frequencies are jointly estimated. This is in consistency with our earlier result that at least for large

*Q*the cross-sensitivity of the DR and SR frequencies in their joint estimation is expected to be small. Furthermore, the approximate result (20), indicated by the dashed line segments on each plot, is quite accurate even for

*Q*=15, although for smaller

*Q*values the greater cross-sensitivity of DR and SR frequencies causes DR-frequency estimation to become more noisy as more SR frequencies are jointly estimated. Second, the normalized CRB, which may be regarded as the inverse of the maximum SNR with which the object spatial power spectrum may be estimated, rises dramatically, in accordance with the approximate analysis of Sec. 5, as the optical band edge is approached from below and then exceeded. Finally, with increasing

*Q*the CRB for estimating the DR frequencies increases linearly with

*Q*, when expressed as a function of the scaled sample index,

*ℓ/Q*, but the SR-frequency-estimation variance has dramatically larger CRBs.

The next set of figures correspond to the same set of parameters as Figs. 1 except that the peak SNR was now chosen to be 10^{3}. In this bright-source limit, the estimation problem is dominated by the shot (Poisson) noise of counting, and the CRBs were normalized by dividing them by the square of the mean source intensity (which is proportional to the strength of the dc component of the source). The source width parameter was sufficiently large so that the source intensity had only a weak variation over its support, and the general behavior of the CRB in the present bright-source limit is, as expected, much the same as that we saw for the uniform-additive-noise considered in the previous figures. However, the normalized CRB value at each spatial frequency is about 3–4 orders smaller because of the large source brightness.

The situation changes when the source width is decreased in relation to the support size. Figures 3 (a)–(c) explore the dependence of CRB on a changing Gaussian-source-width parameter value in the bright-source limit (peak SNR=10^{3}). In all of our figures, *L*=1 and *Q*=20, so the optical system bandwidth, *B*
_{0}=*Q*/(2*L*), equals 10. The spatial spectrum of the Gaussian source (22) has a characteristic width of order 2/(*πwQ*), when expressed as a fraction of the optical system bandwidth. The support-based sinc sampling function, on the other hand, has a chacteristic frequency width of order 2/*Q*. For *w*=1, the source spectral width is then less than the sampling-function bandwidth, which itself is a small fraction of the optical system bandwidth for *Q*≫1. Thus, there is little signal at the SR frequencies that can get coupled into the optical bandwidth, leading to essentially no support-induced gain in information at the DR or SR frequencies for such wide Gaussian sources. By contrast, for an extremely narrow source, *w*=0.05, the source spectral profile has a significant value at the optical band edge. In spite of the narrow bandwidth of the spectral sampling function (for large *Q*), the relatively slow variation of the object spectrum near the band edge implies a relatively strong aliasing of the SR frequencies of the source into the optical band. However, this makes it difficult to disambiguate and thus estimate the SR frequencies individually. The optimal condition holds when the source spectral profile is neither too small nor too slowly varying at the band edge so the support-based aliasing yields good sensitivity of the data just below the band edge to the individual SR frequencies lying just above. This condition is met for *w*=0.1 for which we find that the CRBs for an unbiased estimation of the frequency samples on the whole are smaller by a factor of 2–3 compared to those for other source widths. The two limits of very narrow and very wide sources lead to little source-based information gain when *Q* is not small.

The form of the object is of course irrelevant for the additive-noise-dominated estimation, but when the noise is neither additive nor shot noise but a comparably weighted combination of the two, as for peak SNR of 10, then the results change dramatically. Indeed, when the source support is narrow compared to the width of the Gaussian (22), we gain much information from the support, as seen in Figs. 4 in the sharply lower CRBs for the wider Gaussian sources.

Next we consider the dependence of the CRBs on the number, *M*, of sub-pixel-shifted LR images used to estimate the SR frequencies. An important motivation for our present work has been to advance an integrated treatment of DSR and OSR when images are in general under-sampled by the detector pixels. For the case of critically sampled images, the use of more than one LR image can affect only the relative variance of estimation noise by reducing it by a factor of *M* corresponding to the standard noise reduction due to multiple independent measurements of a specific quantity. This is seen in Figs. 5 for the same additive-noise-dominated estimation CRBs as covered in Fig. 1, except that *M*=5. We verified the same factor-*M* reduction factor for critically sampled images in the bright-source limit too.

However, when these images are undersampled (*χ*>1), the availability of multiple sub-pixel-shifted LR images should provide for nontrivial DSR as well as noise reduction. We show this in the next set of figures where *χ* equals 1 (Fig. 6(a)), 2 (Fig. 6(b)), and 5 (Fig. 6(c)), while keeping *M*=5, in the bright-source limit. For *M* larger than *χ*, not only is DSR possible but the normalized CRBs are also reduced by the noise-averaging factor *M/χ*. This is readily seen for *χ*=1 by comparing Fig. 6(a) with Fig. 2(d), particularly at the lower spatial frequencies. By contrast, as the last figure (Fig. 6(d)) in this set shows, having fewer images than the under-sampling factor, *i.e*., *M*<*χ*, in a joint estimation of the DR and SR frequencies dramatically increases the noise of estimating even those DR frequencies that would be normally resolved by the detector. The same, although not shown here, is true in the faint-source limit as well. These considerations underscore the fact that the number of independent spatial degrees of freedom for a support-limited, band-limited imaging environment is given by the total number of critical samples, spaced consecutively by distance 1/(2*B*) over the spatial support, 2*L*, of the object, namely by 2*Q*, as has been widely recognized before [21, 22, 25].

The cusps seen in the curves on Fig. 6(c) at *ℓ*=8 and 16 signal reduced information at spatial frequencies that are multiples of 2/*χ* (times *B*) where the detector PTF vanishes identically. Without any support-based convolutions in the spatial-frequency domain, these frequencies would have infinite CRBs (i.e., no data sensitivity) and these cusps would be replaced by true divergences, as we noted earlier in Sec. 5.1. This is an instance of information at certain spatial frequencies, even those within the optical passband, that without any prior knowledge would be simply inaccessible to measurements [1]. Whenever *χ*>2, the CRB plots exhibit cusp-like behavior at those spatial frequencies whose index, *ℓ*, is an integral multiple of *ℓ*
_{0}≡2*Q*/*χ*. For *χ*=2, as in Fig. 6(b), a PTF zero occurs precisely at the optical band edge. The vanishing of both the PTF and OTF there yields dramatically higher CRBs for frequency samples near the band edge, as a comparison with Fig. 6(a) easily shows.

In all of our figures containing CRB plots it is clear that regardless of the value of *Q* the minimum noise variance of adding a single pair of SR spatial frequencies to the reconstructed image, *i.e*., increasing *ℓ _{max}* by 1, is essentially a constant factor, of order 10

^{1.5}to 10

^{2}, larger than before. This represents a logarithmic scaling of the number of estimated SR frequencies with image SNR, whether in the additive or Poisson-noise-dominated regime, a conclusion that has also been reached by previous researchers [28, 25, 26].

We repeated our study for different brightness profiles, including those that contain only a few purely sinusoidal waveforms along with the dc term, corresponding to the presence of only a handful of Fourier spatial frequencies including some just beyond the diffraction limit of the imager. Without finite support boundaries, such sources would have infinite extension, so the knowledge of the support does provide useful information, particularly in the moderate-SNR regime, when one or more of the source frequencies lie just outside the band edge. We found many of our CRB results to be qualitatively similar to those obtained for the Gaussian source.

## 7. OSR in the image domain

We now shift attention to the physical plane and evaluate the CRBs for estimating the image-domain pixels. The physical-domain FI matrix elements are related to those of the frequency-domain FI matrix considered so far via a pair of forward and inverse FTs, namely

where *ℓ _{max}*/(2

*L*) is the largest spatial frequency considered in the object. Note that

*ℓ*=

_{max}*Q*refers to frequency at the optical band edge of the imager. It is assumed that the object pixel indices

*k*,

*k*′ run between -

*ℓ*and

_{max}*ℓ*-1 as well, corresponding to successive pixels that are

_{max}*L*/

*ℓ*apart. Note that this is consistent with the Shannon-Nyquist sampling theorem, since for a signal of bandwidth

_{max}*ℓ*/(2

_{max}*L*) there are in all 2

*ℓ*independent, critical, physical-domain samples successively separated by distance

_{max}*L/ℓ*over the [-

_{max}*L,L*) support of the signal. Thus, as increasingly higher degree of bandwidth extension, or OSR, is sought and achieved, it is possible to estimate physical-domain pixels that are increasingly finer.

The physical-domain CRBs for unbiased estimation are given by inverting the corresponding FI matrix (23), which because of the unitarity of the Fourier basis functions, namely

may be expressed as the following forward-inverse pair FT of the inverse of the Fourier-domain FI matrix:

In the next set of figures, we plot the physical-domain CRB (25) for the following three-frequency sinusoidal object intensity distribution:

in which the dc term *f*
_{0} is taken to be a sufficiently large positive number so the object intensity is non-negative everywhere over the support, [-*L*,*L*]. The spatial frequencies, *ν*
_{1},*ν*
_{2}, and *ν*
_{3}, were so chosen that at least one of them was just above the optical band edge, *B*
_{0}=*Q*/(2*L*). The object and its critically sampled mean image are shown in Fig. 7. Note the large loss of spatial resolution in the image.

The peak SNR was chosen at 10^{3} corresponding to the bright-source limit, while *Q* was fixed at 20. For *L*=1, *i.e*., for a support size of 2 units, the largest of the three sinusoidal frequencies, 8, 9, and 11, we chose to illustrate our results lies just above the band edge at 10. A variety of values for the detector undersampling factor *χ* and the number *M* of LR images were considered to give us a fuller picture of how the theoretical prediction of estimation fidelity fares against the physically expected results. From the figures, we arrive at the following important conclusions. First, the CRB is 2–4 orders of magnitude lower near the support boundaries than deep inside the support interior, at least when some OSR (*ℓ _{max}*>

*Q*) is sought. This is expected since the knowledge of support provides the most information about the object near its boundaries, rather than in its interior, as the object is required, with no error, to have zero pixel intensity everywhere outside the support. This space-variant character has also been noted by other investigators in both one and two dimensional studies [5, 26].

Second, even for the undersampled images (Figs. 8(b)–(d)), it is possible to achieve a reasonable fidelity of estimation, provided sufficiently many sub-pixel-shifted images carrying high-frequency information are gathered first. This is most dramatically seen when we compare Fig. 8(b) to 8(c): A subpixel shift of 1/2 corresponds, for *χ*=2, *M*=2, to critically sampling the object, but that misses its spatial frequency 11 that lies just beyond the optical band width. On the other hand, for *M*=3, the subpixel shift is a third of a pixel from one image to the next, which carries significant information about the SR frequency 11, unlike the coarser sub-pixel shift of 1/2. This is seen in dramatically reduced CRBs when *M* is changed from 2 to 3. The improvement continues, but a far slower pace as *M* is increased further, in a manner consistent with the 1/*M* reduction of the squared noise-to-signal ratio with multiple, independently performed measurements (Fig. 8 (d)), as noted earlier.

Third, by comparing Figs. 8(a) and (d) we see that it is possible to reach quantitatively similar estimation error lower bounds even when images are ×2 undersampled, provided *M* is sufficiently large. This quantifies the cost overhead of DSR, which requires a sufficiently large LR image sequence before the fidelity of OSR can begin to approach that for a single critically sampled image.

Finally, when *ℓ _{max}*<

*Q*, as true for the bottom two curves on each sub-plot, the CRBs are essentially uniform over the entire image except close to the support boundaries. Mathematically, this is so because

**J**

^{-1}is a diagonal matrix in this case, as seen from Eq. (20), and thus the CRB of the

*k*image pixel,

*C*

^{(ph)}

*, as given by Eq. (25), becomes independent of the pixel index*

_{k}*k*.

## 8. Conclusions

This paper has explored estimation-theoretic limits on the extent of bandwidth extrapolation, or OSR, from a set of sub-pixel-shifted undersampled images based on a prior knowledge of object support in one dimension. Our approach makes use of a special sampling-type basis set in the spatial-frequency domain to expand the spatial-frequency amplitude of the object intensity distribution in terms of its spatial-frequency samples that are successively separated by the critical sampling interval. Since such an expansion of the object spatial spectrum does not depend on the imager transfer functions, our approach may be employed without any difficulty to analyze DSR and OSR in arbitrary imaging systems.

While this problem has been studied extensively over the past four decades, the present work further clarifies the essential principles behind OSR and theoretical upper bounds on it when the support of the object brightness distribution is known in advance. Like the work of Matson and Tyler [26], we have exploited the concept of FI and associated error lower bounds to characterize a fundamental estimation-theoretic limit on the fidelity of performing any bandwidth extension. Such error bounds, which can be large near and above the optical band edge, were computed for unbiased estimation and graphically illustrated both in the spatial-frequency and image domains. Most practical image reconstructions from statistical image data are biased estimations with typically far lower errors, however, so the present results give us confidence that OSR from image data may be practical even when the size of the object is its only additional property known *a priori*.

The present work on 1D signals can be extended to two-dimensional (2D) images by generalizing the expansion (7) by means of appropriate basis functions defined over the 2D support of the images. The 1D sinc basis functions are replaced by the Fourier-Bessel basis functions for 2D circular supports and by a direct-product basis of the sinc functions for 2D rectangular supports. For general support geometries, the basis functions may be further generalized, e.g., via a variational formulation, but the relation of the spatial-frequency spectrum to the coefficients of the expansion, which are the discrete samples analogous to *F _{ℓ}* of the present paper, becomes more involved. These generalizations will be discussed in future papers.

#### A. Gaussian Approximation to the Mixed Gaussian-Poisson Noise PDF and Approximate FI

Consider a sum of two independent random variates,

where *X* is distributed according to a Poisson distribution with mean 〈*X*〉 and *N* according to a zero-mean Gaussian PDF with variance *σ*
^{2}. Its characteristic function, *Q _{Y}* (λ), defined as the expected value of exp(

*iλY*), is the product of the characteristic functions of

*X*and

*N*evaluated at λ, namely exp{〈

*X*〉[exp(

*iλ*)-1]} and exp(-

*λ*

^{2}

*σ*

^{2}/2), respectively. Expanding the exponent of the product of these two functions in powers of λ yields the following expression for

*QY*(

*λ*):

Note that the first exponential factor in Eq. (28) is significant only when |*λ*|≤*O*(1)/(*σ*
^{2}+〈*X*〉)^{1/2}. For either *σ* or 〈*X*〉 large compared to 1, this implies that only values of |*λ* | small compared to 1 are important. For such small |*λ* |, the exponent of the second exponential factor in Eq. (28) is of order 〈*X*〉*λ*
^{3}≤*O*(1)〈*X*〉/(*σ*
^{2}+〈*X*〉)^{3/2}, which is small compared to 1 provided either 〈*X*〉≫1 or *σ*≫1. In either case, the second exponential factor may be set approximately equal to 1 over the range of values of *λ* over which the first exponential factor is significant. For a given 〈*X*〉, the PDF of *Y*, which is the inverse Fourier transform of *QY* (*λ*), may thus be approximated as

$$\approx \frac{1}{2\pi}{\int}_{-\infty}^{\infty}{e}^{-i\lambda \left(y-\u3008X\u3009\right)-{\lambda}^{2}\left({\sigma}^{2}+\u3008X\u3009\right)/2}d\lambda $$

$$=\frac{1}{\sqrt{2\pi \left({\sigma}^{2}+\u3008X\u3009\right)}}{e}^{-{\left(y-\u3008X\u3009\right)}^{2}/\left[2\left({\sigma}^{2}+\u3008X\u3009\right)\right],}$$

which is the sought Gaussian form that is valid whenever either *σ*≫1 or 〈*X*〉≫1.

Taking the second derivative of the negative logarithm of *P _{Y}* (

*y*|〈

*X*〉) with respect to 〈

*X*〉 generates a quantity whose expected value is the FI relative to 〈

*X*〉. With approximation (29), which is valid if

*σ*

^{2}+〈

*X*〉≫1, this set of steps yields the following FI:

$$=\frac{1}{{\sigma}^{2}+\u3008X\u3009}+\frac{1}{{2\left({\sigma}^{2}+\u3008X\u3009\right)}^{2}}\approx \frac{1}{{\sigma}^{2}+\u3008X\u3009},$$

where repeated use was made of the fact that 〈(*y*-〈*X*〉)〉=0.

If 〈*X*〉 is a differentiable function of a vector of parameters, *λ*̱=(*λ*
_{1},…,*λ _{n}*), then the elements of the FI matrix with respect to those parameters are given by a generalization of Eq. (30) that follows immediately from the chain rule of differentiation,

$$=\frac{\partial \u3008X\u3009}{\partial {\lambda}_{\ell}}\xb7J\left(Y\mid \u3008X\u3009\right)\xb7\frac{\partial \u3008X\u3009}{\partial {\lambda}_{\ell \prime}}\text{}$$

$$\approx \frac{\partial \u3008X\u3009}{\partial {\lambda}_{\ell}}\xb7\frac{1}{{\sigma}^{2}+\u3008X\u3009}\xb7\frac{\u3008\partial X\u3009}{\partial {\lambda}_{\ell \prime}}.$$

This expression may be easily generalized further for the case of *M* independent measurement channels for which 〈*X*〉 is a vector with each component depending functionally on the parameter vector *λ*̱. In this case we need simply append a label to the mean value, say 〈*X _{m}*〉, to indicate a specific measurement channel and then sum expression (31) over all channels,

*m*=1,…,

*M*. For the case of a linear functional dependence of 〈

*X*〉 on

_{m}*λ*̱, the desired expression (13) is obtained in this way. By setting the mean values equal to 0 in that expression, we recover the case of additive noise alone, namely Eq. (11).

#### B. FI and CRB in the Large-Q, Additive-Noise-Dominated Limit

For *σ*
^{2}
* _{D}*≫max(

*g*

^{(m)}

*), the FI matrix elements are given by Eq. (11). Substituting expression (9) into expression (11), first for*

_{k}*ℓ*and then in its complex-conjugated form for

*ℓ*′ using a different integration variable, say

*u*′, yields a double frequency integral for

*Jℓℓ*′,

$$\times \mathrm{sinc}\left(Qu-\ell \right)\mathrm{sinc}\left(Qu\prime -\ell \prime \right)S\left(u-\ell \prime \right),$$

where we have exchanged the order of the double integral and the double sum, denoting the latter by *S*(*u*-*u*′),

in which as before *t*̃* _{m}*=

*m/M*. By separating out the summand into a product of an

*m*-dependent factor and a

*k*-dependent factor, we may separate the double sum into a product of two single geometric-series sums, which may be evaluated simply. The following expression for the double sum (33) thus results:

$$=\mathrm{exp}\left[-\mathit{i\pi}\left(u-u\prime \right)\chi \left(1-\frac{1}{2M}\right)\right]\frac{\mathrm{sin}\left[\pi Q\left(u-u\prime \right)\right]}{\mathrm{sin}\left[\pi \frac{\chi}{2M}\left(u-u\prime \right)\right]},$$

where repeated use was made of the identity, 1-exp(-*ix*)=2*i*exp(-*ix*/2) sin(*x*/2) to obtain the second equality from the first.

In the limit *Q*→∞, the ratio of the two sines in the final expression in Eq. (34) is highly peaked at *u*=*u*′ within a width of order 1/*Q*, so for the purposes of evaluating the frequency integrals, it may be replaced in terms of the Dirac delta function by (2*M/χ*)*δ* (*u*-*u*′). Because of the vanishing of the denominator of this ratio at |*u*-*u*′|=2*Mn*/*χ, n*=1,2,…, however, this replacement is not quite accurate, but whenever *M*≥*χ*, this is of no consequence since |*u*-*u*′| cannot exceed 2 inside the bandlimited double-frequency integral (32). Physically speaking, assuming *M*≥*χ* is sensible since we must acquire at least as many pixels of information as there are parameters to estimate. For a down-sampling factor *χ*, each LR image has only of order 2*Q/χ* meaningful pixels over which the image is supported, so we must acquire at least *M*=*χ* images to estimate (2*Q*+1) or more (for finite OSR) spatial-frequency samples of the object brightness distribution.

In view of this *δ* function, we may perform one of the two integrals in Eq. (32), say the *u*′ integral, trivially so the following expression for the FI matrix elements is obtained in the limit of large *Q*:

Let us now consider this expression for the different ranges of values of *ℓ* and *ℓ*ℓ′.

*B.1.* |*ℓ*|,|*ℓ*′|<*Q*

For these matrix elements, note that the latter two sinc functions in Eq. (35) are both highly peaked at *u*=*ℓ/Q* and *u*=*ℓ*′/*Q* within a width of order 1/*Q*. But since two dissimilar values of *ℓ* and *ℓ*′ differ by at least 1, this implies that for *Q*≫1 the integral (35) will be vanishingly small unless *ℓ*=*ℓ*′. The narrow width of these sinc functions enables us to perform the integral approximately by replacing all slowly varying functions of u in the integrand by their values at *u*=*ℓ/Q* and extending the integration limits to ±∞. Use of the well known integral,

then confirms the previously obtained result (19) for *Jℓℓ*′.

*B.2*. |*ℓ*|,|*ℓ*′|>*Q*

For the SR frequency samples, *ℓ*,*ℓ*′>*Q*, we may simplify the product of the latter two sinc functions in Eq. (35) as

$$=\frac{\mathrm{cos}\left[\pi (\ell -\ell \prime )\right]-\mathrm{cos}\left[2\pi Qu-(\ell +\ell \prime )\pi \right]}{2{\pi}^{2}\left(Qu-\ell \right)(Qu-\ell \prime )}$$

Because of the rapid oscillation of the second term in the numerator of Eq. (37) for large *Q*, its contribution to the integral (35) is of order 1/*Q* when compared to the contribution of the first term. We may thus ignore the rapidly oscillating term, and the expression (35) reduces to the form,

where we used the fact that cos(*ℓ*-*ℓ*′)*π*=(-1)* ^{ℓ}*-

*′.*

^{ℓ}*B.3*. |*ℓ*|<*Q*, |*ℓ*′|>*Q*
*or* |*ℓ*|>*Q*, |ℓ′|<*Q*

When |*ℓ*|<*Q* but |*ℓ*′|>*Q*, then the sinc(*Qu*-*ℓ*) factor inside the integral (35) is peaked within the range of integral (35) and may be replaced by (1/*Q*)*δ*(*u*-*ℓ*/*Q*) in the large *Q* limit. This enables us to evaluate the integral as

since sinc(*ℓ*-*ℓ*′)=0 for all integral *ℓ*,*ℓ*′. The other sub-case, |*ℓ*|>*Q*,|*ℓ*′|<*Q*, is symmetrical, and the same result (39) holds once again.

All of the key results obtained in this Appendix are approximate but they become increasingly more accurate as *Q* becomes increasingly large and *ℓ*,*ℓ*′ take values well in the interior of their ranges considered here.

## Acknowledgments

This work was inspired by numerous discussions on the subject of optical superresolution with Professor M. Fiddy in a number of program-review workshops held under the PERIODIC project. The research was supported in part by IARPA and executed under contract with DMEA, and by AFOSR under grant nos. FA9550-08-1-0151 and FA9550-09-1-0495.

## References and links

**1. **S. Prasad, “Digital superresolution and the generalized sampling theorem,” J. Opt. Soc. Am. A **24**, 311–325 (
2007) [CrossRef]

**2. **S. Prasad, “Digital and optical superresolution of low-resolution image sequences,” Unconventional Imaging Conference, 2007 Annual SPIE Meeting, San Diego, CA, Aug 25–29, 2007: Proc. SPIE **6712**, 67120E 1–11 (
2007).

**3. **R. Gerchberg, “Superresolution through error energy reduction,” Opt. Acta **21**, 709–721 (
1974). [CrossRef]

**4. **A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits Syst. **CAS-22**, 735–742 (
1975). [CrossRef]

**5. **S. Plevritis and A. Macovski, “Spectral extrapolation of spatially bounded images,” IEEE Trans. Medical Imaging **14**, 487–497 (
1995). [CrossRef]

**6. **W. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. **62**, 55–59 (
1972). [CrossRef]

**7. **L. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. **79**, 745–754 (
1974) [CrossRef]

**8. **L. Shepp and Y. Vardi, “Maximum likelihood reconstruction in positron emission tomography,” IEEE Trans. Medical Imaging **1**, 113–122 (
1982). [CrossRef]

**9. **T. Holmes, “Maximum-likelihood image restoration adapted for noncoherent optical imaging,” J. Opt. Soc. Am. A **5**, 666–673 (
1988). [CrossRef]

**10. **T. Holmes and Y.-H. Liu, “Richardson-Lucy/maximum likelihood image restoration algorithm for fluorescence microscopy: further testing,” Appl. Opt. **28**, 4930–4938 (
1989). [CrossRef] [PubMed]

**11. **J.-A. Conchello, “Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images,” J. Opt. Soc. Am. A **15**, 2609–2619 (
1998). [CrossRef]

**12. **E. Boukouvala and A. Lettington, “Restoration of astronomical images by an iterative superresolving algorithm,” Astron. Astrophys. **399**, 807–811 (
2003). [CrossRef]

**13. **B. R. Frieden, “Image enhancement and restoration,” in *Picture Processing and Digital Filtering*,
T. Huang, ed. (Springer-Verlag, NY,
1975), 177–248.

**14. **R. Narayan and R. Nityananda,“Maximum entropy image restoration in astronomy,” Ann. Rev. Astron. Astrophys. **24**, 127–170 (
1986). [CrossRef]

**15. **J. Hogboom, “Aperture synthesis with a non-regular distribution of interferometer baselines,” Astron. Astrophys. Supp. **15**, 417–426 (
1974).

**16. **D. Fried, “Analysis of the CLEAN algorithm and implications for superresolution,” J. Opt. Soc. Am. A **12**, 853–860 (
1995). [CrossRef]

**17. **P. Magain, F. Courbin, and S. Sohy, “Deconvolution with correct sampling,” Astrophys. J. **494**, 472–477 (
1998). [CrossRef]

**18. **F. Pijpers, “Unbiased image reconstruction as an inverse problem,” Mon. Not. Roy. Astron. Soc. **307**, 659–668 (
1999). [CrossRef]

**19. **R. Puetter and R. Hier, “Pixon sub-diffraction space imaging,” in Proc. SPIE **7094**, 709405-709405-12 (
2008).

**20. **J. Starck, E. Pantin, and F. Murtagh, “Deconvolution in astronomy: A review,” Publ. Astron. Soc. Pacific **114**, 1051–1069 (
2002). [CrossRef]

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

**22. **G. Toraldo Di Francia, “Degrees of freedom of an image,” J. Opt. Soc. Am. **59**pp. 799–805 (
1969).

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

**24. **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]

**25. **M. Bertero and C. De Mol, “Superresolution by data inversion,” Progress in Optics **XXXVI**, 129–178 (
1996). [CrossRef]

**26. **C. Matson and D. Tyler, “Primary and secondary superresolution by data inversion,” Opt. Express **14**, 456–473 (
2006). [CrossRef] [PubMed]

**27. **D. Robinson and P. Milanfar, “Statistical performance analysis of superresolution,” IEEE Trans. Image Process. **15**, 1413–1428 (
2006). [CrossRef] [PubMed]

**28. **P. Sementilli, B. Hunt, and M. Nader, “Analysis of of the limit of superresolution in incoherent imaging,” J. Opt. Soc. Am. A **10**, 2265–2276 (
1993). [CrossRef]

**29. **H. Van Trees, *Detection, Estimation, and Modulation Theory* (Wiley, New York,
1968).