Frequency-Constrained QR: Signal and Image Reconstruction

Abstract

Because a finite set of measurements is limited in the amount of spectral content it can represent, the reconstruction process from discrete samples is inherently band-limited. In the case of 1D sampling using ideal measurements, the maximum bandwidth of regular and irregular sampling is well known using Nyquist and Gröchenig sampling theorems and lemmas, respectively. However, determining the appropriate reconstruction bandwidth becomes difficult when considering 2D sampling geometries, samples with variable apertures, or signal to noise ratio limitations. Instead of determining the maximum bandwidth a priori, we derive an inverse method to simultaneously reconstruct a signal and determine its effective bandwidth. This inverse method is equivalent to incrementally computing a band-limited inverse using a frequency-constrained QR decomposition (FQR). Comparisons between reconstruction results using FQR and QR decompositions illustrate how FQR is less sensitive to noisy measurement errors, but it is more sensitive to high-frequency components. These methods are particularly useful in the reconstruction of remote sensing images from such as microwave radiometers and scatterometers.

1. Introduction

Signal reconstruction from a finite set of measurements is a common problem found in signal and image processing. When measurements are modeled as idealized delta-like samples, the Nyquist sampling theorem [1] and Gröchenig’s sampling lemma [2] can be used to determine a reconstruction’s maximum spectral band limit in the 1D sampling case. However, when dealing with measurements with variable apertures, SNR limitations, or 2D sampling considerations, finding the reconstructable bandwidth can be difficult [3].

To address these reconstruction issues, we review the concept of a band-limited sampling matrix, its inverse, and an iterative method to determine a reconstruction’s reconstructable bandwidth as treated in [4]. In this paper, we present an inverse method for simultaneously reconstructing a signal and determining its effective bandwidth, i.e., the maximum amount of the recoverable signal bandwidth without significant noise amplification. This work is motivated by the need for finer resolution from microwave remote sensors such as radiometers and scatterometers, which have irregular sampling with variable apertures [5,6,7,8].

We first consider the need for a band-limited inverse by reviewing how a discrete sampling matrix is inherently band-limited in order to find a unique, nonaliased solution. We also review how noise is amplified and how model aliasing occurs in the computation of a sampling matrix inverse. We then show how to incrementally invert a sampling matrix’s band-limited content and derive a method for computing an incremental band-limited inverse using a frequency-constrained QR decomposition (FQR). We then compare the difference in reconstruction using FQR and QR decompositions.

2. Materials and Methods

2.1. Discrete Sampling

To motivate the need for a bandlimited inverse, we consider how finite measuring systems are inherently band-limited. We start by defining an ideal finite measurement as a sample of a continuous linear system,

$$z_i=\int m_i\left(r\right) a\left(r\right) dr,$$

(1)

where $z_i$ is the ith discrete measurement, $m_i\left(r\right)$ is the measurement’s response function (MRF), $a\left(r\right)$ is the truth signal, and r is a positional index. Given a finite number of measurements N, inverting the measurement system in Equation (1) is an ill-posed problem, since the existence, uniqueness, and stability of a solution depend on the configuration, spacing, and number of $m_i\left(r\right)$, as well as assumptions about $a\left(r\right)$.

To illustrate the difficulty in estimating $a\left(r\right)$ from a set of $z_i$, imagine if each $m_i\left(r\right)$ is a delta function centered at position $r_i$ (i.e., $m_i\left(r\right)=\delta (r-r_i)$) and then $z_i=a\left(r_i\right)$. Although each $z_i$ may accurately represent the value of the truth signal at the point $r_i$, it is impossible to uniquely infer the entire continuous signal ($a\left(r\right)$) without further information, since an infinite number of signals may also cross at the same discrete values. However, if $a\left(r\right)$ is assumed to be band-limited, sufficient sampling configurations can uniquely determine $a\left(r\right)$. Under band-limited sampling theory, only a periodic band-limited signal $a\left(r\right)$ can be perfectly reconstructed from a finite set of samples [3].

We define a band-limited signal $a\left(r\right)$ to be a signal which only posses frequency content below a maximum cutoff frequency, or band limit B (i.e., the Fourier transform of $a\left(r\right)$ has a finite region of support). Let $A\left(\omega \right)$ be the Fourier transform of $a\left(r\right)$, $A\left(\omega \right)$ that is band-limited if $\left|A\right(\omega \left)\right|=0\forall \omega >B$. According to Nyquist, a uniformly sampled band-limited signal can be perfectly reconstructed from a set of ideal discrete samples $a\left[n\right]$ as

$$a\left(r\right)=\sum _{n}a\left[n\right] d(r-r_n),$$

(2)

where $a\left[n\right]=a\left(r_n\right)$, $d\left(r\right)$ is a Dirichlet kernel (i.e., band-limited identity function) corresponding to the band limit of $a\left(r\right)$, and n represents an arbitrary sample index corresponding to position $r_n$ when the sample distances between $r_n$ is at least twice the highest frequency. Since $d\left(r\right)$ is a periodic band-limited identity, $a\left(r\right)$ is periodic. Note that Equation (2) is equivalent to bandlimited interpolation. Even though multiple sampling configurations (i.e., the particular spacing and number of sample positions $r_n$) can satisfy the reconstruction in Equation (2), all such configurations are equivalent under the band-limited assumption [3]. For convenience, we can assume that $r_n$ is spaced on a uniform grid, where the grid spacing is $\mathrm{\Delta }r<=1/\left(2f_c\right)$ with $f_c$ as the highest frequency within the band limit of $a\left(r\right)$ (i.e., $1/\left(2f_c\right)$ is the Nyquist sampling distance).

The relationship between $z_i$ and $a\left[n\right]$ can be found by combining Equations (1) and (2),

$$\begin{array}{cc}\hfill z_i& =\int m_i\left(r\right) \left[\sum _{n}a\left[n\right] d(r-r_n)\right] dr,\hfill \end{array}$$

$$\begin{array}{cc}\hfill z_i& =\sum _n\left[\int m_i\left(r\right) d(r-r_n) dr\right] a\left[n\right],\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill z_i& =\sum _n{m}_i\left[n\right] a\left[n\right],\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill m_i\left[n\right]& =\int m_i\left(r\right) d(r-r_n) dr,\hfill \end{array}$$

(5)

where $m_i\left[n\right]$ is a discrete band-limited sample of $m_i\left(r\right)$, and $d(r-r_n)$ is the Dirichlet kernel function centered at $r_n$. Note that $m_i\left[n\right]\ne m_i\left(r_n\right)$, unless $m_i\left(r\right)$ is already band-limited within $d\left(r\right)$’s corresponding bandlimit.

Equation (3) can be expressed in vector and matrix form as

$$\begin{array}{cc}\hfill z_i& ={\mathbf{m}}_i^T\mathbf{a},\hfill \\ \hfill \mathbf{z}& =\mathbf{M}\mathbf{a},\hfill \end{array}$$

(6)

where ${\mathbf{m}}_i$ and a are length ‘N’ vectors of $m_i\left[n\right]$ and $a\left[n\right]$, respectively, z is a length M vector of $z_i$, and M is a matrix with rows of ${\mathbf{m}}_i^T$ and is $\left[MxN\right]$, where M is the total number of $z_i$, and N is the total number of $r_n$. From Equation (6), it can be seen that the band-limited reconstruction process is equivalent to finding the inverse of M; however, the invertibility of M depends on the chosen band-limited discretization and spacing of each MRF.

In the case of uniformly spaced samples and constant MRF apertures, i.e., $m_i=m_j$, M is square and invertible (assuming each $m_i\left[n\right]$ is defined over the assumed band limit). When measurements are dependent, i.e., they do not increase the rank of the inverse of M, they do not improve the resolution of the estimate but do provide redundancy for noise reduction. Irregular sampling and variable apertures tend to have more dependent measurements. A unique band-limited solution only exists if a band-limited inverse of M exists.

In order to define the band-limited inverse of M, it is helpful to study the discretization of ideal and band-limited measurements. Under band-limited assumptions, the discrete equivalent to an ideal measurement is found by substituting $m_i\left(r\right)=d(r-r_i)$ into Equation (5):

$$\begin{array}{cc}\hfill m_i\left[n\right]& =\int d(r-r_i) d(r-r_n) dr=d(r_n-r_i)=d_i\left[n\right],\hfill \end{array}$$

(7)

where $d_i\left[n\right]$ is a Dirichlet kernel centered at $r_i$ that is discretely sampled over the grid defined by $r_n$.

To express Equation (2) in matrix form, we define a band-limtied identity matrix D constructed with rows $d(r-r_n)$ sampled over the grid defined by $r_n$. The multiplication of a by D is an identity operation, as a is unchanged by the band limit of D. This band-limited identity operation, $\mathbf{a}=\mathbf{D}\mathbf{a}$, is the discrete band-limtied reconstruction equivalent to Equation (2). Using this result, we can express a band-limited inverse of M:

$$\begin{array}{cc}\hfill \mathbf{D}& =\mathbf{W}\mathbf{M},\hfill \\ \hfill \mathbf{a}& =\mathbf{W}\mathbf{z},\hfill \end{array}$$

(8)

where W is a band-limited inverse of M. The existence of W depends on M being full rank over a band-limited subspace. When this criteria is met, W is the SVD pseudo-inverse of M. When M is full row rank, W is the Moore–Penrose pseudo-inverse, and D becomes the identity matrix. The existence of W depends on the sampling positions $r_i$, the choice of discrete positions $r_n$, and the band limit. For a more rigorous study on how sampling conditions affect the existence of this band-limited inverse, see [1,4,9,10].

In summary, an ideal band-limited reconstruction can be accomplished from a finite number of measurements if a band-limited inverse of M exists, where M is the discrete band-limited sampling matrix of all measurement response functions.

2.2. Background

2.2.1. Sampling Matrix Inversion

The process of estimating a discrete signal a from a discrete set of measurements z is equivalent to estimating a particular inverse of the sampling matrix M. The estimated signal can be found by using the sampling inverse of $\mathbf{M}$,

$$\begin{array}{cc}\hfill \widehat{\mathbf{a}}& ={\mathbf{M}}^{†}\mathbf{z},\hfill \end{array}$$

(9)

where $\widehat{\mathbf{a}}$ is the estimated signal, and ${\mathbf{M}}^{†}$ is a unique inverse of the sampling matrix M. While there may be a plethora of matrix inversion methods that could be used to estimate ${\mathbf{M}}^{†}$, their various constraints, assumptions, and approximations can affect the quality of the reconstructed signal $\widehat{\mathbf{a}}$.

2.2.2. Inverse Error Amplification

In order to better understand noise amplification and reconstruction artifacts, we consider the additive noise model:

$$\begin{array}{cc}\hfill \tilde{\mathbf{z}}& =\mathbf{M}\mathbf{a}+\mathit{\nu },\hfill \end{array}$$

(10)

where $\tilde{\mathbf{z}}$ is a set of noisy measurements, and $\mathit{\nu }$ is measurement noise. Assuming ${\mathbf{M}}^{†}$ exists, let

$$\begin{array}{cc}\hfill {\mathbf{a}}_{\nu }& ={\mathbf{M}}^{†}\mathit{\nu },\hfill \end{array}$$

(11)

where ${\mathbf{a}}_{\nu }$ is the inverse noise. Let

$$\begin{array}{cc}\hfill {\mathbf{a}}_{\mu }& =(\mathbf{I}-{\mathbf{M}}^{†}\mathbf{M})\mathbf{a},\hfill \end{array}$$

(12)

where ${\mathbf{a}}_{\mu }$ is the aliasing of inverse model artifacts,

$$\begin{array}{cc}\hfill \tilde{\mathbf{a}}& =\mathbf{a}+{\mathbf{a}}_{\mu }+{\mathbf{a}}_{\nu },\hfill \end{array}$$

(13)

and $\tilde{\mathbf{a}}$ is the estimated signal corrupted by amplified noise and model artifacts.

In the ideal case, ${\mathbf{M}}^{†}$ is a perfect inverse of M, i.e., ${\mathbf{M}}^{†}\mathbf{M}=\mathbf{I}$ and ${\mathbf{a}}_{\mu }=0$, such that there are no model errors or reconstruction artifacts. However, when sampling model errors are introduced into the computation of ${\mathbf{M}}^{†}$, ${\mathbf{M}}^{†}\mathbf{M}\ne \mathbf{I}$, the content from a can be aliased in ${\mathbf{a}}_{\mu }$. Additionally, since ${\mathbf{M}}^{†}$ is sensitive to small errors, computational errors can be amplified such that features in ${\mathbf{a}}_{\nu }$ can overpower recoverable signal features.

2.2.3. Deconvolution Example

To illustrate how errors are amplified and aliased in the inversion of a sampling matrix, we examine the case of discrete deconvolution using a shift invariant sampling function m. The measurement vector is

$$\begin{array}{cc}\hfill \mathbf{z}& =\mathbf{m}⊛\mathbf{a}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =\mathbf{M}\mathbf{a},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \mathbf{M}& =\left[\begin{array}{cccc}{m}_0& m_1& \cdots & m_{N-1}\\ m_{N-1}& m_0& ⋮& m_{N-2}\\ m_{N-2}& m_{N-1}& \ddots & m_{N-3}\\ ⋮& ⋮& \ddots & ⋮\\ m_1& m_2& \cdots & m_0\end{array}\right]\hfill \end{array}$$

(15)

where z is the convolution of vectors m and a, $(⊛)$ is the circular convolution operator, and $\mathbf{M}$ is a circulant matrix of $\mathbf{m}$. Because $\mathbf{M}$ is a circulant matrix, it can be factored using the Discrete Fourier Transform (DFT) matrix F [11]:

$$\begin{array}{cc}\hfill \mathbf{M}& ={\mathbf{F}}^{\ast }(\mathbf{F}\mathbf{m}·\mathbf{F})\hfill \\ \hfill & ={\mathbf{F}}^{\ast }\mathrm{diag}\left(\mathbf{F}\mathbf{m}\right)\mathbf{F}\hfill \end{array}$$

(16)

where ∗ denotes the conjugate transpose, $(·)$ is an element-wise multiplier, $\mathrm{diag}\left(\mathbf{x}\right)$ denotes forming a diagonal matrix with the vector x, and F is normalized such that ${\mathbf{F}}^{\ast }\mathbf{F}=\mathbf{I}$.

The deconvolution of m can likewise be represented in matrix form. In the ideal case, the deconvolution of m is w, where

$$\mathbf{w}={\mathbf{F}}^{\ast }{\left(\mathbf{F}\mathbf{m}\right)}^{-1}.$$

(17)

The deconvolution then is

$$\begin{array}{cc}\hfill \tilde{\mathbf{a}}& =\mathbf{w}⊛\mathbf{z}\hfill \\ \hfill & =\mathbf{W}\mathbf{z},\hfill \\ \hfill & =\mathbf{W}\mathbf{M}\mathbf{a},\hfill \end{array}$$

(18)

where $\mathbf{W}$ is a circulant matrix of $\mathbf{w}$, and $\mathbf{W}$ is also the ideal inverse of $\mathbf{M}$. Then, we have the following:

$$\begin{array}{cc}\hfill \mathbf{W}& ={\mathbf{F}}^{\ast }\mathrm{diag}{\left(\mathbf{F}\mathbf{m}\right)}^{-1}\mathbf{F},\hfill \\ \hfill & ={\mathbf{M}}^{-1}.\hfill \end{array}$$

(19)

However, the deconvolution of m is sensitive to small spectral values within $\mathbf{F}\mathbf{m}$. The inversion of small spectral values is a common cause of noise amplification. When the deconvolution does not match or equalize the spectrum of $\mathbf{F}\mathbf{m}$, aliasing of spectral features can occur as follows:

$$\begin{array}{cc}\hfill {\mathbf{a}}_{\mu }& ={\mathbf{F}}^{\ast }\mathrm{diag}(1-(\mathbf{F}\mathbf{w}·\mathbf{F}\mathbf{m}\left)\right)\mathbf{F}\mathbf{a}\hfill \end{array}$$

(20)

where ${\mathbf{a}}_{\mu }$ defines how spectral features are aliased by deconvolution.

To demonstrate noise amplification and model error aliasing, a simple Gaussian MRF was used to blur a 300 sample truth signal with 3 Dirichlet kernels of different degrees. The deconvolution of this synthetically blurred scene was performed twice. The first deconvolution was done after noise $\mathit{\nu }$ added to the scene in order to achieve a 30 dB SNR. The second deconvolution was done after adding corruption to the Gaussian MRF to achieve a 10 dB SNR of corruption. Both deconvolutions are shown in Figure 1a. The FFT of the ideal Gaussian MRF, i.e., $\mathbf{F}\mathbf{w}·\mathbf{F}\mathbf{m}$ alongside its equalization at an expected 30 dB SNR level and equalization using the corrupted MRF, are shown in Figure 1b. Note that even a small amount of error, without any conditioning, significantly distorted the equalization such that $1-(\mathbf{F}\mathbf{w}·\mathbf{F}\mathbf{m})$ amplified the noise floor to a visible level.

Regularization techniques can be used to avoid inverting small errors in the frequency domain, lessening the sensitivity of deconvolution to small errors. However, regularization inherently changes the assumed measurement response spectrum and distorts the reconstructed signal. To preserve the recoverable measurement response spectrum, the sampling matrix can instead be band-limited such that deconvolution only inverts a subset of the recoverable measurement response spectrum, minimizing the effects of small errors.

2.2.4. Band-Limited Inverse

To mitigate the amplification of noise and reconstruction artifacts, band-limited filtering can be used to remove frequency content within the model before it is aliased in the calculation of the matrix inverse. One way to accomplish this is by band-limiting each row of M (i.e., taking its DFT, applying a band-limiting filter by multiplication, then taking the inverse DFT) and then calculating the pseudo-inverse of M. Assuming that the same band-limiting filter is applied to each row, we can express this process in linear algebra terms by recognizing that a band-limiting frequency filter h can also be applied by a band-limited matrix $\mathbf{D}=\mathbf{F}\mathrm{diag}\left(\mathbf{h}\right)\mathbf{F}$. We define the right-hand band-limited inverse of M bandlimited by D as

$$\begin{array}{cc}\hfill {\mathbf{M}}^{\beta }& ={\left(\mathbf{M}\mathbf{D}\right)}^{†},\hfill \end{array}$$

(21)

where ${\mathbf{M}}^{\beta }$ is the band-limited inverse, and D is a band-limiting matrix that applies the frequency filter h to each row of M.

While there may not be a computational advantage to computing the band-limited inverse using Equation (21), since band-limiting with a matrix multiplication is $\mathcal{O}\left(n^3\right)$, and contrastingly, the DFT is $\mathcal{O}\left(n^2\log n\right)$, Equation (21) can help us understand how a band-limited inverse is related to the ideal pseudo-inverse, for instance, ${\mathbf{M}}^{\beta }={\mathbf{M}}^{-1}$ when M is full rank and D is an identity matrix (maximum discrete bandwidth). ${\mathbf{M}}^{\beta }={\mathbf{M}}^{†}$ when M is not full rank and D is a band-limited identity matrix (i.e., when D is a Toeplitz Dirichlet kernel matrix) with a higher bandwidth than each row of M. However, when D is lower rank than M, i.e., a lower bandwidth than the rows of M, then ${\mathbf{M}}^{\beta }\ne {\mathbf{M}}^{†}$.

Ideally, to preserve model content and remove noise, a band-limiting filter matching the recoverable model spectrum should be used. However, since the recoverable model spectrum depends on measurement sampling locations and MRF bandwidths, it can be hard to determine for variable aperture measurements a priori. One approach to teasing out the recoverable spectrum of variable aperture measurements is suggested in [3], which involves incrementally bandlimiting a sampling matrix and repetitively calculating its rank.

Instead of treating the bandlimiting and inverse steps sequentially, it becomes advantageous to compute a band-limited matrix and matrix inverse simultaneously. In this paper, we derive a method for calculating an iterative band-limited inverse, which we later show is equivalent to a partial QR decomposition using a frequency basis (see Appendix A).

2.3. Frequency QR

Here, we define a method for computing a band-limited inverse using single-rank updates. When computed in full, this method is equivalent to performing QR decomposition on a frequency basis (FQR). When computed sequentially, it can progressively construct a band-limited inverse from low to high frequencies, stopping when the magnitude of each frequency component falls below the assumed SNR noise floor in the frequency domain. This stops the reconstruction at the point where measurement noise theoretically overtakes the magnitude of reconstructable content. In the noise-free case, these criteria stop the reconstruction at the inherent sampling band limit without going over. An additional example comparing the result of FQR to that of a QR decomposition is provided in the following section.

2.3.1. Partitioned Band-Limited Matrix

Let ${\mathbf{f}}_n$ represent an orthonormal basis vector of the frequency domain such that ${\mathbf{f}}_i^{\ast }{\mathbf{f}}_j=0,\forall i\ne j$.

$$\begin{array}{cc}\hfill \mathbf{F}& =\left[\begin{array}{c}{\mathbf{f}}_0,{\mathbf{f}}_1,\dots ,{\mathbf{f}}_{N-1}\end{array}\right],\hfill \\ \hfill {\mathbf{F}}^{\ast }\mathbf{F}& =\mathbf{I},\hfill \end{array}$$

(22)

where F is a matrix of the basis vectors ${\mathbf{f}}_n$ and I is the identity matrix. Typically, each ${\mathbf{f}}_n$ corresponds to a DFT kernel within the calculation of the DFT. However, other frequency bases may be used when advantageous. More examples of frequency bases and their advantages are discussed in a later section.

A band-limiting frequency filter can be applied in the frequency domain by applying the filter weights $h_n$ to each corresponding basis vector ${\mathbf{f}}_n$. Taking advantage of the orthogonal basis, $\mathbf{M}\mathbf{D}$ can be partitioned into a set of singular matrices:

$$\begin{array}{cc}\hfill \mathbf{M}\mathbf{D}& =\mathbf{M}\mathbf{F}\mathrm{diag}\left(\mathbf{h}\right){\mathbf{F}}^{\ast },\hfill \\ \hfill \mathbf{M}\mathbf{D}& =({\mathbf{B}}_0+{\mathbf{B}}_1+\cdots +{\mathbf{B}}_{N-1}),\hfill \end{array}$$

(23)

where $\mathrm{diag}\left(\mathbf{x}\right)$ denotes the formation of a diagonal matrix with x along its diagonal, h is a frequency filter, and each ${\mathbf{B}}_n$ is a singular matrix defined by the outer product:

$$\begin{array}{cc}\hfill {\mathbf{B}}_n& =\mathbf{M}{\mathbf{f}}_n{h}_n{\mathbf{f}}_n^{\ast },\hfill \\ \hfill {\mathbf{B}}_n& ={\mathbf{b}}_n{\mathbf{f}}_n^{\ast },\hfill \end{array}$$

(24)

where $h_n$ are the values of h corresponding to each basis vector ${\mathbf{f}}_n$, and ${\mathbf{b}}_n$ is a column vector of $\mathbf{M}\mathbf{F}$ weighted by $h_n$:

$$\begin{array}{cc}\hfill {\mathbf{b}}_n& =\mathbf{M}{\mathbf{f}}_n{h}_n.\hfill \end{array}$$

(25)

Note that the matrix product ${\mathbf{B}}_i{\mathbf{B}}_j^{\ast }=0,\forall i\ne j$, because the row space of each ${\mathbf{B}}_n$ is orthogonal. However, the column spaces of ${\mathbf{B}}_n$ are not orthogonal, meaning that ${\mathbf{B}}_i^{\ast }{\mathbf{B}}_j\ne 0,\forall i\ne j$.

Using the partitioned representation of $\mathbf{M}\mathbf{D}$, we can express the progression of band-limited sampling matrices ${\left(\mathbf{M}\mathbf{D}\right)}_n$ as

$$\begin{array}{cc}\hfill {\left(\mathbf{M}\mathbf{D}\right)}_n& =\sum _{k=0}^n{\mathbf{b}}_k{\mathbf{f}}_k^{\ast },\hfill \end{array}$$

(26)

where each ${\left(\mathbf{M}\mathbf{D}\right)}_n$ progressively increases in bandwidth as n increases. The direct inversion of each ${\left(\mathbf{M}\mathbf{D}\right)}_n$ leads us to the existence of a progressive band-limited inverse:

$$\begin{array}{cc}\hfill {\mathbf{M}}_n^{\beta }& ={\left(\mathbf{M}\mathbf{D}\right)}_n^{†}.\hfill \\ \hfill {\mathbf{M}}_n^{\beta }& ={({\mathbf{B}}_0+{\mathbf{B}}_1+\cdots +{\mathbf{B}}_n)}^{†}.\hfill \end{array}$$

(27)

Regrettably, while each ${\mathbf{B}}_n$ may have orthogonal row spaces, the progressive band-limited inverse cannot be formed by the sum of ${\mathbf{B}}_n^{†}$, i.e.,

$$\begin{array}{cc}\hfill {\mathbf{M}}_n^{\beta }& \ne \sum _{k=0}^n{\mathbf{B}}_n^{†}, \hfill \end{array}$$

(28)

since each ${\mathbf{B}}_n$ does not have orthogonal column spaces. However, given that sequential ${\left(\mathbf{M}\mathbf{D}\right)}_n$ only differ by singular matrix updates, it stands to reason that ${\mathbf{M}}_n^{\beta }$ can also be partitioned into singular band-limited matrices.

2.3.2. Derivation

Our derivation of the progressive band-limited inverse ${\mathbf{M}}_n^{\beta }$ begins with the following theorem, which relates the generalized inverse of the sum of two matrices to the generalized inverse of each individual matrix [12]:

$$\begin{array}{cc}\hfill \mathrm{if} \mathrm{rank}(\mathbf{X}+\mathbf{Y})& =\mathrm{rank}\left(\mathbf{X}\right)+\mathrm{rank}\left(\mathbf{Y}\right),\hfill \\ \hfill \mathrm{then} {(\mathbf{X}+\mathbf{Y})}^{†}& =(\mathbf{I}-\mathbf{S}){\mathbf{X}}^{†}(\mathbf{I}-\mathbf{T})+\mathbf{S}{\mathbf{Y}}^{†}\mathbf{T},\hfill \\ \hfill \mathbf{S}& ={\left({\mathbf{P}}_{Y^{*}}\left[\mathbf{I}-{\mathbf{P}}_{X^{*}}\right]\right)}^{†},\hfill \\ \hfill \mathbf{T}& ={\left(\left[\mathbf{I}-{\mathbf{P}}_X\right]{\mathbf{P}}_Y\right)}^{†},\hfill \end{array}$$

(29)

where X and Y are two arbitrary matrices, ${\mathbf{P}}_X$ denotes the projection onto the column space of X, S represents the disjoint portion of the row space of Y from the row space of X, and T represents the disjoint portion of the column space of Y from the column space of X.

Recursively applying the theorem in Equation (29), we can find each progressive band-limited inverse. Applying the theorem to ${\mathbf{B}}_n$ and the zero matrix, we trivially find the initial bandlimited inverse:

$${\mathbf{M}}_0^{\beta }={\mathbf{B}}_0^{†}={\displaystyle \frac{{\mathbf{f}}_0{\mathbf{b}}_0^{\ast }}{{\mathbf{b}}_0^{\ast }{\mathbf{b}}_0}}.$$

(30)

The next bandlimited inverse can be found by performing a rank update to ${\mathbf{B}}_0^{†}$ using ${\mathbf{B}}_1$:

$$\begin{array}{cc}\hfill \mathrm{given} \mathrm{rank}({\mathbf{B}}_0+{\mathbf{B}}_1)& =\mathrm{rank}\left({\mathbf{B}}_0\right)+\mathrm{rank}\left({\mathbf{B}}_1\right),\hfill \\ \hfill \mathrm{then}{({\mathbf{B}}_0+{\mathbf{B}}_1)}^{†}& =(\mathbf{I}-{\mathbf{S}}_1){\mathbf{B}}_0^{†}(\mathbf{I}-{\mathbf{T}}_1)+{\mathbf{S}}_1{\mathbf{B}}_1^{†}{\mathbf{T}}_1,\hfill \\ \hfill {\mathbf{S}}_1& ={\left({\mathbf{f}}_1{\mathbf{f}}_1^{\ast }(\mathbf{I}-{\mathbf{f}}_0{\mathbf{f}}_0^{\ast })\right)}^{†}={\mathbf{f}}_1{\mathbf{f}}_1^{\ast },\hfill \\ \hfill {\mathbf{T}}_1& ={\left(\left[\mathbf{I}-{\displaystyle \frac{{\mathbf{b}}_0{\mathbf{b}}_0^{\ast }}{{\mathbf{b}}_0^{\ast }{\mathbf{b}}_0}}\right]{\displaystyle \frac{{\mathbf{b}}_1{\mathbf{b}}_1^{\ast }}{{\mathbf{b}}_1^{\ast }{\mathbf{b}}_1}}\right)}^{†},\hfill \end{array}$$

(31)

where ${\mathbf{T}}_1$ is the disjoint portion of the column space of ${\mathbf{B}}_1$ from ${\mathbf{B}}_0$, and ${\mathbf{S}}_1$ is the row space projection of ${\mathbf{B}}_1$. As the row space of each ${\mathbf{B}}_1$ is already disjoint ${\mathbf{S}}_1{\mathbf{B}}_1^{†}={\mathbf{B}}_1^{†}$ and $(\mathbf{I}-{\mathbf{S}}_1){\mathbf{B}}_0^{†}={\mathbf{B}}_0^{†}$, upon simplifying using this result, we find the following:

$${\mathbf{M}}_1^{\beta }={({\mathbf{B}}_0+{\mathbf{B}}_1)}^{†}={\mathbf{M}}_0^{\beta }(\mathbf{I}-{\mathbf{T}}_1)+{\mathbf{B}}_1^{†}{\mathbf{T}}_1,$$

(32)

where ${\mathbf{M}}_1^{\beta }$ can be found using ${\mathbf{M}}_0^{\beta }$, ${\mathbf{T}}_1$, and ${\mathbf{B}}_1^{†}$. Applying the theorem to the nth rank update, we can form a recursive relationship, which iteratively constructs a band-limited inverse:

$$\begin{array}{cc}\hfill {\mathbf{M}}_n^{\beta }& ={\mathbf{M}}_{n-1}^{\beta }(\mathbf{I}-{\mathbf{T}}_n)+{\mathbf{B}}_n^{†}{\mathbf{T}}_n,\hfill \\ \hfill {\mathbf{T}}_n& ={\left(\left[\mathbf{I}-{\mathbf{P}}_{{\left(\mathbf{M}\mathbf{D}\right)}_{n-1}}\right]{\displaystyle \frac{{\mathbf{b}}_n{\mathbf{b}}_n^{\ast }}{{\mathbf{b}}_n^{\ast }{\mathbf{b}}_n}}\right)}^{†},\hfill \end{array}$$

(33)

where ${\mathbf{M}}_n^{\beta }$ can be found by applying a rank update to ${\mathbf{M}}_{n-1}^{\beta }$ with the disjoint column space projection of each ${\mathbf{B}}_n^{†}$.

To find ${\mathbf{T}}_n$ recursively, we use the fact that each ${\mathbf{T}}_n$’s column space is disjoint and $\{{\mathbf{T}}_0,{\mathbf{T}}_1,\dots ,{\mathbf{T}}_n\}$ spans the column subspace of ${\left(\mathbf{M}\mathbf{D}\right)}_n$ to re-express $(\mathbf{I}-{\mathbf{P}}_{{\left(\mathbf{M}\mathbf{D}\right)}_{n-1}})$ with the sum of the disjoint column spaces of ${\mathbf{T}}_n$:

$$\begin{array}{cc}\hfill (\mathbf{I}-{\mathbf{P}}_{{\left(\mathbf{M}\mathbf{D}\right)}_{n-1}})& =(\mathbf{I}-{\mathbf{P}}_{{\mathbf{T}}_{n-1}})(\mathbf{I}-{\mathbf{P}}_{{\mathbf{T}}_{n-2}})\dots (\mathbf{I}-{\mathbf{P}}_{{\mathbf{T}}_0}),\hfill \\ \hfill (\mathbf{I}-{\mathbf{P}}_{{\left(\mathbf{M}\mathbf{D}\right)}_{n-1}})& =\left(\mathbf{I}-\sum _{k=0}^{n-1}{\mathbf{P}}_{{\mathbf{T}}_k}\right).\hfill \end{array}$$

(34)

By substitution, we find a recursive relationship to calculate each singular matrix ${\mathbf{T}}_n$:

$$\begin{array}{cc}\hfill {\mathbf{T}}_n& ={\left(\left[\mathbf{I}-\sum _{k=0}^{n-1}{\mathbf{P}}_{{\mathbf{T}}_k}\right]{\displaystyle \frac{{\mathbf{b}}_n{\mathbf{b}}_n^{\ast }}{{\mathbf{b}}_n^{\ast }{\mathbf{b}}_n}}\right)}^{†},\hfill \\ \hfill {\mathbf{T}}_n& ={\mathbf{b}}_n{\left({\mathbf{b}}_n-\sum _{k=0}^{n-1}{\mathbf{P}}_{{\mathbf{T}}_k}{\mathbf{b}}_n\right)}^{†},\hfill \\ \hfill {\mathbf{T}}_n& ={\mathbf{b}}_n{\mathbf{t}}_n^{\ast },\hfill \end{array}$$

(35)

where ${\mathbf{t}}_n^{\ast }$ is the vector inverse of ${\mathbf{q}}_n$

$$\begin{array}{cc}\hfill {\mathbf{q}}_n& ={\mathbf{b}}_n-\sum _{k=0}^{n-1}{\displaystyle \frac{{\mathbf{t}}_k{\mathbf{t}}_k^{\ast }}{{\mathbf{t}}_k^{\ast }{\mathbf{t}}_k}}{\mathbf{b}}_n,\hfill \end{array}$$

(36)

and ${\mathbf{q}}_n$ is a column space basis vector of ${\mathbf{M}}_n^{\beta }$. Note that this recursive relationship is equivalent to the Gram–Schmidt process used in QR factorization [13,14,15].

Substituting ${\mathbf{T}}_n={\mathbf{b}}_n{\mathbf{t}}_n^{\ast }$ into the recursive relationship for ${\mathbf{M}}_n^{\beta }$ in Equation (33), we have the following:

$$\begin{array}{cc}\hfill {\mathbf{M}}_n^{\beta }& ={\mathbf{M}}_{n-1}^{\beta }(\mathbf{I}-{\mathbf{b}}_n{\mathbf{t}}_n^{\ast })+{\mathbf{B}}_n^{†}{\mathbf{b}}_n{\mathbf{t}}_n^{\ast },\hfill \end{array}$$

(37)

which reveals that the row space of ${\mathbf{M}}_n^{\beta }$ is spanned by orthogonal ${\mathbf{t}}_n$. As each band-limited update is singular, we can slove for the column space vector ${\mathbf{m}}_n^{†}$, which updates each ${\mathbf{M}}_n^{\beta }$. By multiplying both sides of Equation (37) by ${\mathbf{q}}_n$ and noting that ${\mathbf{B}}_n^{†}{\mathbf{b}}_n={\mathbf{f}}_n$, we find an expression for ${\mathbf{m}}_n^{†}$:

$$\begin{array}{cc}\hfill {\mathbf{m}}_n^{†}& ={\mathbf{M}}_n^{\beta }{\mathbf{q}}_n,\hfill \\ \hfill {\mathbf{m}}_n^{†}& ={\mathbf{f}}_n+{\mathbf{M}}_n^{\beta }\left[{\mathbf{q}}_n-{\mathbf{b}}_n\right],\hfill \\ \hfill {\mathbf{m}}_n^{†}& ={\mathbf{f}}_n-{\mathbf{M}}_n^{\beta }\left[\sum _{k=0}^{n-1}{\mathbf{q}}_k{\mathbf{t}}_k^{\ast }{\mathbf{b}}_n\right],\hfill \\ \hfill {\mathbf{m}}_n^{†}& ={\mathbf{f}}_n-\sum _{k=0}^{n-1}{\mathbf{m}}_k^{†}{\mathbf{t}}_k^{\ast }{\mathbf{b}}_n.\hfill \end{array}$$

(38)

Note that Equation (38) is equivalent to the process of forward substitution used to solve a matrix inverse using a QR decomposition. Each ${\mathbf{m}}_n^{†}{\mathbf{t}}_n^{\ast }$ forms a band-limited singular matrix partition of ${\mathbf{M}}_n^{\beta }$, which—when summed together—form progressive band-limited inverses:

$$\begin{array}{cc}\hfill {\mathbf{M}}_n^{\beta }& =\sum _{k=0}^n{\mathbf{m}}_k^{†}{\mathbf{t}}_k^{\ast }.\hfill \end{array}$$

(39)

It is important to note that although the row space of each ${\mathbf{m}}_n^{†}{\mathbf{t}}_n^{\ast }$ is orthogonal, and the column space (spanned by ${\mathbf{m}}_n^{†}$) is not orthogonal, which is similar to the inequality in Equation (28). Specifically, each ${\mathbf{m}}_n^{†}$ may contain content from previous ${\mathbf{m}}_k^{†}$, which are removed through forward substitution each band-limited update. This means that each band-limited inverse ${\mathbf{M}}_n^{\beta }$ is not necessary equivalent to the pseudo-inverse of M band-limited by each previous ${\mathbf{B}}_n$,

$$\begin{array}{cc}\hfill {\mathbf{M}}_n^{\beta }& \ne \left[\sum _{k=0}^n{\mathbf{B}}_k\right]{\mathbf{M}}^{†},\hfill \end{array}$$

(40)

since M^† may contain content associated with higher-frequency features than contained in $\left[{\sum }_{k=0}^n{\mathbf{B}}_k\right]$. Therefore, ${\mathbf{M}}_n^{\beta }=\left[{\sum }_{k=0}^n{\mathbf{B}}_k\right]{\mathbf{M}}^{†}$ only when ${\mathbf{B}}_n$ has a higher bandwidth than the maximum bandlimit of M.

As pointed out in the derivation, the vectors ${\mathbf{t}}_n$ and ${\mathbf{m}}_n^{†}$ can be solved for using Gram–Schmidt orthogonalization and forward substitution methods, respectively. These are the basic steps for QR decomposition and inverse solutions. Appendix A shows how FQR is equivalent to applying a frequency rotation to a matrix and then solving for its partial inverse using QR iterations. A summary of the algorithm steps matching the notation given in this paper is shown in Algorithm 1. Note that many algorithmic improvements can be made to make the orthogonalization and forward substitution methods more numerically stable and efficient [16,17]

Algorithm 1: Frequency-constrained QR.

${\mathbf{f}}_i^{\ast }{\mathbf{f}}_j=0, {\mathbf{f}}_i^{\ast }{\mathbf{f}}_i=1, \forall i\ne j$

${\mathbf{b}}_n\leftarrow \mathbf{M}{\mathbf{f}}_n{h}_n$

${\mathbf{q}}_n\leftarrow {\mathbf{b}}_n-{\sum }_{k=0}^{n-1}{\mathbf{q}}_k{\mathbf{t}}_k^{\ast }{\mathbf{b}}_n$

${\mathbf{t}}_n\leftarrow {\mathbf{q}}_n/\left({\mathbf{q}}_n^{\ast }{\mathbf{q}}_n\right)$

${\mathbf{m}}_n^{†}\leftarrow {\mathbf{f}}_n-{\sum }_{k=0}^{n-1}{\mathbf{m}}_k^{†}{\mathbf{t}}_k^{\ast }{\mathbf{b}}_n$

${\mathbf{M}}_n^{\beta }\leftarrow {\mathbf{M}}_{n-1}^{\beta }+{\mathbf{m}}_n^{†}{\mathbf{t}}_n^{\ast }$

2.3.3. Observations

When computed over all N nonzero singular updates, this method is identical to taking a full inverse. In theory, this progressive band-limited inverse can be used as an inversion method for any arbitrary matrix X. However, the computational complexity of this method is no better than using QR decomposition and forward substitution methods using $\mathcal{O}\left(n^3\right)$.

The main difference between this inversion method and other inversion methods is the order in which mutually independent information is inverted. For comparison, we describe how the SVD and QR decompositions break up a matrix inverse. The computation of the inverse using SVD decomposition breaks a matrix down starting with its largest eigenvectors, cutting off the inversion when eigenvalues become too small. The computation of the inverse using QR decomposition breaks a matrix down starting with its first column, cutting off the inversion when new column’s contributions become too small. The computation of the inverse using FQR breaks a matrix down starting with its lowest frequency projection, cutting off the inversion when higher frequency’s contributions become too small. All three methods (SVD, QR, and FQR) are essentially rotations of each other. In the noise-free case, all three of these methods are equivalent when run to completion, but in the presence of computation noise or model errors, the methods posses different sensitivities.

As a sampling matrix is inherently band-limited, the FQR inverse is advantageous to limit the inversion of noise content outside the sampling bandlimit.

2.3.4. SNR Filtering

When the inherent bandlimit of M is unknown, but the expected noise floor or signal-to-noise ratio (SNR) is known, FQR can be modified to stop inversion when the magnitude of new ${\mathbf{t}}_n$ drops below the expected noise level. This is very similar to the process of Wiener filtering [18], except that no regularization is added to the inversion process. As the expected SNR is decreased, the recoverable bandwidth shrinks and FQR iterations reach the raised expected noise level earlier. Further SNR filtering can be implemented by adding a weighting or taper on the magnitude of each ${\mathbf{t}}_n$. Such a taper can be used to mitigate the Gibbs phenomenon present at the edge of band-limited signals.

3. Results

To demonstrate the difference between QR and FQR decompositions, we iteratively performed reconstruction on two example 1D truth signals. To contrast both decompositions, two simple synthetic signals of 101 samples were generated and used as truth for each simulation. The first signal consisted of three low-bandwidth peaks, and the second signal resembles a rect function with higher frequency trends. The average magnitude of the second signal was set to 100× that of the first signal.

In order to form a variable aperture sample matrix M, 101 MRFs were synthetically generated selecting by shifting three MRFs apertures centered at random indices. The three MRF apertures were one Gaussian centered with a standard deviation of 25 indices, one left skewed Gaussian approximately 13 indices wide, and one skewed right that mirrored the previous MRF. Each MRF was normalized to sum to 1. The same sampling matrix, a 101 × 101 matrix, was used on each input truth signal to create two corresponding sets of noise free measurements. To both of these measurement sets, additive Gaussian noise was added to achieve a 30 dB SNR. Example MRFs are shown in Figure 2 to visualize the centered and skewed apertures.

In general the main difference between the reconstruction using QR vs. FQR is that QR reconstructs onto a basis of orthogonal MRFs, while FQR reconstructs onto an orthogonal frequency basis. To illustrate this, three selected MRFs are shown in Figure 3 along with the low-order truth signal for reference. Next to this plot in Figure 3, an example QR reconstruction using the same three MRFs is shown. Additionally, a FQR reconstruction using the first three discrete cosine transform (DCT) bins of all MRFs is shown. Not surprisingly, the FQR took fewer iterations than the QR to represent a low-bandwidth signal, since FQR uses all MRFs to reconstruct each frequency basis, while QR only uses one MRF for each iteration.

The full reconstruction for the low bandwidth input signal was performed using all measurements and frequency iterations. A few sample QR reconstructions using only the first 5, 11, and 21 MRFs are shown in Figure 4a. Similarly, sample FQR reconstructions using only the first 5, 11, and 21 DCT bins are also shown in Figure 4a. The same process was repeated for the higher bandwidth input signal in Figure 4b. Note that in the QR decomposition, a single MRF could greatly bias the reconstruction due to the mismatch between MRF shape and truth signal features. Theoretically, if all iteration steps are taken, the bias caused by a single measurement is removed in future iteration steps. Likewise, in FQR, a single frequency basis can bias the reconstruction due to sampling density sensitivities.

The total error with respect to each truth signal was calculated for each iteration for both the QR and FQR reconstruction methods. The total error per pixel for these four simulation are all shown on the same scale in Figure 5. Note that signal 2 is higher than signal 1 due to the intentional 100× magnitude to plot both simulations on the same total error curves. In both cases, the FQR method was able to resolve more low-bandwidth content before the amplification of noise overtook the reconstruction process.

4. Discussion

4.1. Equivalence

Given that FQR is equivalent to performing QR decomposition after a frequency-basis rotation (see Appendix A), it stands to reason that other matrix inversion techniques can also be used to compute a band-limited inverse. For example, an iterative SVD decomposition could be used to compute band-limited inverses (see Appendix B). If inverting a noise-free and well-conditioned matrix, the frequency-basis rotation becomes redundant, as an identity matrix could be theoretically used as a trivial rotation. However, in the presence of noise or poor conditioning, the basis rotation can reorder and prioritize reconstructing signal content to minimize amplifying noise.

The prioritizing of reconstructed content in a QR inverse computation has been a wide field of study. Many rigorous and specialized algorithms have been developed to solve very large or sparse matrices [19,20]. To extend the discussion of how the choice of frequency basis can prioritize signal content differently, we review several frequency basis choices and their benefits in 1D and 2D sampling.

4.2. Choice of Frequency Basis

Theoretically, any basis that is chosen results in the same inverse if the sampling matrix is full rank and all basis vector iterations are taken. However, not all basis vectors are needed to compute a band-limited inverse. Thus, the choice of frequency basis can affect when and how much noise amplification occurs during the inversion process. Depending on the assumed structure of each inverse, a recombination or reordering of frequency components can be useful in computing progressively band-limited inverses.

4.2.1. DFT vs. DCT

The DFT is typically preferred in 1D sampling, as it works for both real and complex signals. To progressively reconstruct lowpass signals with an increasing bandwidth, reordering the DFT frequencies from low to high is useful.

If we define the DFT basis vector as

$$\begin{array}{cc}\hfill f_{nk}& =\sqrt{{\displaystyle \frac{1}{N}}}exp\left({\displaystyle \frac{-2\pi ikn}{N}}\right),n,k\in [0,N-1]\hfill \end{array}$$

(41)

we find that the first half $n\in [1,\mathrm{floor}(N/2\left)\right]$ increases from low- to high- frequency components (often called the positive frequencies), while the second half $n\in \left[\mathrm{ceil}\right(N/2),N-1]$ decreases from high- to low- frequencies (often called negative frequencies). To rearrange these basis vectors in order of increasing bandwidth, we need to alternate positive and negative frequencies. The ordering is $[0,1,N-1,2,N-2\dots ]$. For real signals, the positive and negative frequencies of the DFT are conjugates of each other, making one redundant in computing the inverse. Thus, for real signals, only the first half of the complex DFT basis vectors is needed to compute each band-limited inverse $[0,1,2,\mathrm{floor}(N/2\left)\right]$. An equivalent real-only basis can be formed by using the discrete cosine transform (DCT), which takes advantage of conjugate symmetric pairs. The most common forms of the DCT are already ordered low to high frequencies (DCT-I, DCT-II, DCT-III, DCT-IV) by taking half-frequency steps [21].

4.2.2. 2D Frequency Basis

For the case of 2D sampling, similar bases can be used such as the 2D DFT for complex signals and 2D DCT for real signals. However, the extra dimensionality of the 2D basis introduces directional frequencies. While directional frequencies may not be a concern when inverting a uniform sampling scheme, some irregular 2D sampling schemes may be sensitive to directional sampling errors.

One method to mitigate these directional sensitivities is to combine the symmetry of the four quadrants of the 2D DFT. Just as the DCT in one dimension takes advantage of the positive and negative symmetry of DFT bins offset by $\pi /2$, it can also be used to pair symmetric quadrants in 2D. To take advantage of this symmetry for complex signals, the discrete sine transform (DST) can be used to fill in the other $\pi /4$ shifted quadrants, i.e., using the DCT and DST to respectively compute the real and imaginary portions of a symmetric DFT basis. Note that basis vectors still need to be ordered by their radial distance to perform progressively bandlimited reconstructions.

It is possible to further combine radially symmetric DFT bins by reducing the frequency step size and applying the appropriate 2D rotations to preserve the orthogonality of each unique radial DFT bin. Example 2D basis vectors for the DFT, DCT, and a shifted even DCT basis are shown in Figure 6.

4.3. Advantages

While there may be many practical uses for the incremental computation of a bandlimited inverse using FQR, we highlight three main advantages for its use in discrete reconstruction problems. Firstly, since it can be hard to determine the maximum reconstruction bandwidth a priori, FQR can be used to simultaneously reconstruct and find the maximum reconstruction bandwidth. Secondly, since FQR uses a frequency basis, modern filtering techniques can be leveraged to further constrain, condition, and regularize the FQR band-limited inverse. Thirdly, since each singular update to the FQR inverse uses all measurements, noisy measurement errors are mitigated when terminating the inversion process early.

4.4. Applications

Using a 2D frequency basis, the FQR or similar frequency-constrained reconstruction methods can be applied to enhance the spatial resolution of remote sensing systems, without the need to determine the maximum reconstruction bandwidth a priori. This is particularly useful for reconstructing irregular sampling schemes and measurements with variable apertures, often found in spaceborne radiometers and scatterometers. However, the requirement of FQR to reconstruct each individual frequency basis can be as computationally taxing as other pseudo-inverse methods. Iterative methods could be used to reduce the computational complexity by approximating a band-limited inverse, but may amplify more noise than an individual frequency-basis inverse. Furthermore, since a frequency basis implicitly assumes a spatial periodicity or band limit, additional spatial tapering or windowing may be required to preserve the temporal or spatial relevance between measurements.

5. Conclusions

A frequency-constrained QR decomposition improves reconstruction performance for bandlimited signals, particularly for measurement with variable apertures, SNR limitations, or 2D sampling considerations. Additionally, a frequency-constrained QR reconstruction can simultaneously reconstruct a signal and determine its effective bandwidth. While this paper focused on applying frequency constraints to QR reconstruction, in order to highlight the theoretical difference between standard matrix inverse techniques and a bandlimited inverse, similar frequency constraints can be applied to more practical and efficient inverse methods. Further work is needed to demonstrate the benefits of a bandlimited inverse, i.e., applying frequency constraints to existing inversion methods and remote sensing datasets.