Search

We study the spectral properties of infinitely smooth multivariate kernel matrices when the nodes form a single cluster. We show that the geometry of the nodes plays an important role in the scaling of the eigenvalues of these kernel matrices. For the multivariate Dirichlet kernel matrix, we establish a criterion for the sampling set ensuring precise scaling of eigenvalues. Additionally, we identify specific sampling sets that satisfy this criterion. Finally, we discuss the implications of these results for the problem of super-resolution, i.e. stable recovery of sparse measures from bandlimited Fourier measurements., Comment: 24 pages, 12 figures

Given a reaction-diffusion equation with unknown right-hand side, we consider a nonlinear inverse problem of estimating the associated leading eigenvalues and initial condition modes from a finite number of non-local noisy measurements. We define a reconstruction criterion and, for a small enough noise, we prove the existence and uniqueness of the desired approximation and derive closed-form expressions for the first-order condition numbers, as well as bounds for their asymptotic behavior in a regime when the number of measured samples is fixed and the inter-sampling interval length tends to infinity. When computing the sought estimates numerically, our simulations show that the exponential fitting algorithm ESPRIT is first-order optimal, as its first-order condition numbers have the same asymptotic behavior as the analytic condition numbers in the considered regime.

In this paper we establish accuracy bounds of Prony's method (PM) for recovery of sparse measures from incomplete and noisy frequency measurements, or the so-called problem of super-resolution, when the minimal separation between the points in the support of the measure may be much smaller than the Rayleigh limit. In particular, we show that PM is optimal with respect to the previously established min-max bound for the problem, in the setting when the measurement bandwidth is constant, with the minimal separation going to zero. Our main technical contribution is an accurate analysis of the inter-relations between the different errors in each step of PM, resulting in previously unnoticed cancellations. We also prove that PM is numerically stable in finite-precision arithmetic. We believe our analysis will pave the way to providing accurate analysis of known algorithms for the super-resolution problem in full generality., Comment: 33 pages, 5 figures

We study recovery of amplitudes and nodes of a finite impulse train from noisy frequency samples. This problem is known as super-resolution under sparsity constraints and has numerous applications. An especially challenging scenario occurs when the separation between Dirac pulses is smaller than the Nyquist-Shannon-Rayleigh limit. Despite large volumes of research and well-established worst-case recovery bounds, there is currently no known computationally efficient method which achieves these bounds in practice. In this work we combine the well-known Prony's method for exponential fitting with a recently established decimation technique for analyzing the super-resolution problem in the above mentioned regime. We show that our approach attains optimal asymptotic stability in the presence of noise, and has lower computational complexity than the current state of the art methods., Comment: 5 pages, 10 figures

Inverse source problems are central to many applications in acoustics, geophysics, non-destructive testing, and more. Traditional imaging methods suffer from the resolution limit, preventing distinction of sources separated by less than the emitted wavelength. In this work we propose a method based on physically-informed neural-networks for solving the source refocusing problem, constructing a novel loss term which promotes super-resolving capabilities of the network and is based on the physics of wave propagation. We demonstrate the approach in the setup of imaging an a-priori unknown number of point sources in a two-dimensional rectangular waveguide from measurements of wavefield recordings along a vertical cross-section. The results show the ability of the method to approximate the locations of sources with high accuracy, even when placed close to each other.

We study the problem of super-resolution of a linear combination of Dirac distributions and their derivatives on a one-dimensional circle from noisy Fourier measurements. Following numerous recent works on the subject, we consider the geometric setting of "partial clustering", when some Diracs can be separated much below the Rayleigh limit. Under this assumption, we prove sharp asymptotic bounds for the smallest singular value of a corresponding rectangular confluent Vandermonde matrix with nodes on the unit circle. As a consequence, we derive matching lower and upper min-max error bounds for the above super-resolution problem, under the additional assumption of nodes belonging to a fixed grid., Comment: 28 pages, 5 figures

Following recent interest by the community, the scaling of the minimal singular value of a Vandermonde matrix with nodes forming clusters on the length scale of Rayleigh distance on the complex unit circle is studied. Using approximation theoretic properties of exponential sums, we show that the decay is only single exponential in the size of the largest cluster, and the bound holds for arbitrary small minimal separation distance. We also obtain a generalization of well-known bounds on the smallest eigenvalue of the generalized prolate matrix in the multi-cluster geometry. Finally, the results are extended to the entire spectrum., Comment: To appear in Applied and Computational Harmonic Analysis

We study rectangular Vandermonde matrices $\mathbf{V}$ with $N+1$ rows and $s$ irregularly spaced nodes on the unit circle, in cases where some of the nodes are "clustered" together -- the elements inside each cluster being separated by at most $h \lesssim {1\over N}$, and the clusters being separated from each other by at least $\theta \gtrsim {1\over N}$. We show that any pair of column subspaces corresponding to two different clusters are nearly orthogonal: the minimal principal angle between them is at most $$\frac{\pi}{2}-\frac{c_1}{N \theta}-c_2 N h,$$ for some constants $c_1,c_2$ depending only on the multiplicities of theclusters. As a result, spectral analysis of $\mathbf{V}_N$ is significantly simplified by reducing the problem to the analysis of each cluster individually. Consequently we derive accurate estimates for 1) all the singular values of $\mathbf{V}$, and 2) componentwise condition numbers for the linear least squares problem. Importantly, these estimates are exponential only in the local cluster multiplicities, while changing at most linearly with $s$.

We consider the problem of stable recovery of sparse signals of the form $$F(x)=\sum_{j=1}^d a_j\delta(x-x_j),\quad x_j\in\mathbb{R},\;a_j\in\mathbb{C}, $$ from their spectral measurements, known in a bandwidth $\Omega$ with absolute error not exceeding $\epsilon>0$. We consider the case when at most $p\le d$ nodes $\{x_j\}$ of $F$ form a cluster whose extent is smaller than the Rayleigh limit ${1\over\Omega}$, while the rest of the nodes are well separated. Provided that $\epsilon \lessapprox SRF^{-2p+1}$, where $SRF=(\Omega\Delta)^{-1}$ and $\Delta$ is the minimal separation between the nodes, we show that the minimax error rate for reconstruction of the cluster nodes is of order ${1\over\Omega}SRF^{2p-1}\epsilon$, while for recovering the corresponding amplitudes $\{a_j\}$ the rate is of the order $SRF^{2p-1}\epsilon$. Moreover, the corresponding minimax rates for the recovery of the non-clustered nodes and amplitudes are ${\epsilon\over\Omega}$ and $\epsilon$, respectively. These results suggest that stable super-resolution is possible in much more general situations than previously thought. Our numerical experiments show that the well-known Matrix Pencil method achieves the above accuracy bounds.

We prove sharp lower bounds for the smallest singular value of a partial Fourier matrix with arbitrary "off the grid" nodes (equivalently, a rectangular Vandermonde matrix with the nodes on the unit circle), in the case when some of the nodes are separated by less than the inverse bandwidth. The bound is polynomial in the reciprocal of the so-called "super-resolution factor", while the exponent is controlled by the maximal number of nodes which are clustered together. As a corollary, we obtain sharp minimax bounds for the problem of sparse super-resolution on a grid under the partial clustering assumptions.

Soft extrapolation refers to the problem of recovering a function from its samples, multiplied by a fast-decaying window and perturbed by an additive noise, over an interval which is potentially larger than the essential support of the window. A core theoretical question is to provide bounds on the possible amount of extrapolation, depending on the sample perturbation level and the function prior. In this paper we consider soft extrapolation of entire functions of finite order and type (containing the class of bandlimited functions as a special case), multiplied by a super-exponentially decaying window (such as a Gaussian). We consider a weighted least-squares polynomial approximation with judiciously chosen number of terms and a number of samples which scales linearly with the degree of approximation. It is shown that this simple procedure provides stable recovery with an extrapolation factor which scales logarithmically with the perturbation level and is inversely proportional to the characteristic lengthscale of the function. The pointwise extrapolation error exhibits a H\"{o}lder-type continuity with an exponent derived from weighted potential theory, which changes from 1 near the available samples, to 0 when the extrapolation distance reaches the characteristic smoothness length scale of the function. The algorithm is asymptotically minimax, in the sense that there is essentially no better algorithm yielding meaningfully lower error over the same smoothness class. When viewed in the dual domain, the above problem corresponds to (stable) simultaneous de-convolution and super-resolution for objects of small space/time extent. Our results then show that the amount of achievable super-resolution is inversely proportional to the object size, and therefore can be significant for small objects.

We provide an overview of some results on the "geometry of error amplification" in solving Prony system, in situations where the nodes near-collide. It turns out to be governed by the "Prony foliations" $S_q$, whose leaves are "equi-moment surfaces" in the parameter space. Next, we prove some new results concerning explicit parametrization of the Prony leaves.

The traditional sparse modeling approach, when applied to inverse problems with large data such as images, essentially assumes a sparse model for small overlapping data patches. While producing state-of-the-art results, this methodology is suboptimal, as it does not attempt to model the entire global signal in any meaningful way - a nontrivial task by itself. In this paper we propose a way to bridge this theoretical gap by constructing a global model from the bottom up. Given local sparsity assumptions in a dictionary, we show that the global signal representation must satisfy a constrained underdetermined system of linear equations, which can be solved efficiently by modern optimization methods such as Alternating Direction Method of Multipliers (ADMM). We investigate conditions for unique and stable recovery, and provide numerical evidence corroborating the theory.

Given two polynomials $P,q$ we consider the following question: "how large can the index of the first non-zero moment $\tilde{m}_k=\int_a^b P^k q$ be, assuming the sequence is not identically zero?". The answer $K$ to this question is known as the moment Bautin index, and we provide the first general upper bound: $K\leqslant 2+\mathrm{deg} q+3(\mathrm{deg} P-1)^2$. The proof is based on qualitative analysis of linear ODEs, applied to Cauchy-type integrals of certain algebraic functions. The moment Bautin index plays an important role in the study of bifurcations of periodic solution in the polynomial Abel equation $y'=py^2+\varepsilon qy^3$ for $p,q$ polynomials and $\varepsilon \ll 1$. In particular, our result implies that for $p$ satisfying a well-known generic condition, the number of periodic solutions near the zero solution does not exceed $5+\mathrm{deg} q+3\mathrm{deg}^2 p$. This is the first such bound depending solely on the degrees of the Abel equation.

We consider Fourier reconstruction problem for signals F, which are linear combinations of shifted delta-functions. We assume the Fourier transform of F to be known on the frequency interval [-N,N], with an absolute error not exceeding e > 0. We give an absolute lower bound (which is valid with any reconstruction method) for the "worst case" reconstruction error of F in situations where the nodes (i.e. the positions of the shifted delta-functions in F) are known to form an l elements cluster of a size h << 1. Using "decimation" reconstruction algorithm we provide an upper bound for the reconstruction error, essentially of the same form as the lower one. Roughly, our main result states that for N*h of order of (2l-1)-st root of e the worst case reconstruction error of the cluster nodes is of the same order as h, and hence the inside configuration of the cluster nodes (in the worst case scenario) cannot be reconstructed at all. On the other hand, decimation algorithm reconstructs F with the accuracy of order of 2l-st root of e.

We consider polynomial systems of Prony type, appearing in many areas of mathematics. Their robust numerical solution is considered to be difficult, especially in "near-colliding" situations. We consider a case when the structure of the system is a-priori fixed. We transform the nonlinear part of the Prony system into a Hankel-type polynomial system. Combining this representation with a recently discovered "decimation" technique, we present an algorithm which applies homotopy continuation to an appropriately chosen Hankel-type system as above. In this way, we are able to solve for the nonlinear variables of the original system with high accuracy when the data is perturbed.

We compare three approaches to studying the behavior of an analytic function $f(z)=\sum_{k=0}^\infty a_kz^k$ from its Taylor coefficients. The first is "Taylor domination" property for $f(z)$ in the complex disk $D_R$, which is an inequality of the form \[ |a_{k}|R^{k}\leq C\ \max_{i=0,\dots,N}\ |a_{i}|R^{i}, \ k \geq N+1. \] The second approach is based on a possibility to generate $a_k$ via recurrence relations. Specifically, we consider linear non-stationary recurrences of the form \[ a_{k}=\sum_{j=1}^{d}c_{j}(k)\cdot a_{k-j},\ \ k=d,d+1,\dots, \] with uniformly bounded coefficients. In the third approach we assume that $a_k=a_k(\lambda)$ are polynomials in a finite-dimensional parameter $\lambda \in {\mathbb C}^n.$ We study "Bautin ideals" $I_k$ generated by $a_{1}(\lambda),\ldots,a_{k}(\lambda)$ in the ring ${\mathbb C}[\lambda]$ of polynomials in $\lambda$. \smallskip These three approaches turn out to be closely related. We present some results and questions in this direction., Comment: arXiv admin note: substantial text overlap with arXiv:1301.6033

We consider the problem of recovering a linear combination of Dirac delta functions and derivatives from a finite number of Fourier samples corrupted by noise. This is a generalized version of the well-known spike recovery problem, which is receiving much attention recently. We analyze the numerical conditioning of this problem in two different settings depending on the order of magnitude of the quantity $N\eta$, where $N$ is the number of Fourier samples and $\eta$ is the minimal distance between the generalized spikes. In the "well-conditioned" regime $N\eta\gg1$, we provide upper bounds for first-order perturbation of the solution to the corresponding least-squares problem. In the near-colliding, or "super-resolution" regime $N\eta\to0$ with a single cluster, we propose a natural regularization scheme based on decimating the samples \textendash{} essentially increasing the separation $\eta$ \textendash{} and demonstrate the effectiveness and near-optimality of this scheme in practice.

We consider the problem of "algebraic reconstruction" of linear combinations of shifts of several known signals $f_1,\ldots,f_k$ from the Fourier samples. Following \cite{Bat.Sar.Yom2}, for each $j=1,\ldots,k$ we choose sampling set $S_j$ to be a subset of the common set of zeroes of the Fourier transforms ${\cal F}(f_\ell), \ \ell \ne j$, on which ${\cal F}(f_j)\ne 0$. It was shown in \cite{Bat.Sar.Yom2} that in this way the reconstruction system is "decoupled" into $k$ separate systems, each including only one of the signals $f_j$. The resulting systems are of a "generalized Prony" form. However, the sampling sets as above may be non-uniform/not "dense enough" to allow for a unique reconstruction of the shifts and amplitudes. In the present paper we study uniqueness and robustness of non-uniform Fourier sampling of signals as above, investigating sampling of exponential polynomials with purely imaginary exponents. As the main tool we apply a well-known result in Harmonic Analysis: the Tur\'an-Nazarov inequality (\cite{Naz}), and its generalization to discrete sets, obtained in \cite{Fri.Yom}. We illustrate our general approach with examples, and provide some simulation results.

We continue studying robustness of solving algebraic systems of Prony type (also known as the exponential fitting systems), which appear prominently in many areas of mathematics, in particular modern "sub-Nyquist" sampling theories. We show that by considering these systems at arithmetic progressions (or "decimating" them), one can achieve better performance in the presence of noise. We also show that the corresponding lower bounds are closely related to well-known estimates, obtained for similar problems but in different contexts.

Systems of Prony type appear in various signal reconstruction problems such as finite rate of innovation, superresolution and Fourier inversion of piecewise smooth functions. We propose a novel approach for solving Prony-type systems, which requires sampling the signal at arithmetic progressions. By keeping the number of equations small and fixed, we demonstrate that such "decimation" can lead to practical improvements in the reconstruction accuracy. As an application, we provide a solution to the so-called Eckhoff's conjecture, which asked for reconstructing jump positions and magnitudes of a piecewise-smooth function from its Fourier coefficients with maximal possible asymptotic accuracy -- thus eliminating the Gibbs phenomenon.

We consider the problem of ``algebraic reconstruction'' of linear combinations of shifts of several signals $f_1,\ldots,f_k$ from the Fourier samples. For each $r=1,\ldots,k$ we choose sampling set $S_r$ to be a subset of the common set of zeroes of the Fourier transforms ${\cal F}(f_\l), \ \l \ne r$, on which ${\cal F}(f_r)\ne 0$. We show that in this way the reconstruction system is reduced to $k$ separate systems, each including only one of the signals $f_r$. Each of the resulting systems is of a ``generalized Prony'' form. We discuss the problem of unique solvability of such systems, and provide some examples.

We consider the "moment vanishing problem" for a general class of piecewise-analytic functions which satisfy on each continuity interval a linear ODE with polynomial coefficients. This problem, which essentially asks how many zero first moments can such a (nonzero) function have, turns out to be related to several difficult questions in analytic theory of ODEs (Poincare's Center-Focus problem) as well as in Approximation Theory and Signal Processing ("Algebraic Sampling"). While the solution space of any particular ODE admits such a bound, it will in the most general situation depend on the coefficients of this ODE. We believe that a good understanding of this dependence may provide a clue for attacking the problems mentioned above. In this paper we undertake an approach to the moment vanishing problem which utilizes the fact that the moment sequences under consideration satisfy a recurrence relation of fixed length, whose coefficients are polynomials in the index. For any given operator, we prove a general bound for its moment vanishing index. We also provide uniform bounds for several operator families.

We consider "Taylor domination" property for an analytic function $f(z)=\sum_{k=0}^{\infty}a_{k}z^{k},$ in the complex disk $D_R$, which is an inequality of the form \[ |a_{k}|R^{k}\leq C\ \max_{i=0,\dots,N}\ |a_{i}|R^{i}, \ k \geq N+1. \] This property is closely related to the classical notion of "valency" of $f$ in $D_R$. For $f$ - rational function we show that Taylor domination is essentially equivalent to a well-known and widely used Tur\'an's inequality on the sums of powers. Next we consider linear recurrence relations of the Poincar\'e type \[ a_{k}=\sum_{j=1}^{d}[c_{j}+\psi_{j}(k)]a_{k-j},\ \ k=d,d+1,\dots,\quad\text{with }\lim_{k\rightarrow\infty}\psi_{j}(k)=0. \] We show that the generating functions of their solutions possess Taylor domination with explicitly specified parameters. As the main example we consider moment generating functions, i.e. the Stieltjes transforms \[ S_{g}\left(z\right)=\int\frac{g\left(x\right)dx}{1-zx}. \] We show Taylor domination property for such $S_{g}$ when $g$ is a piecewise D-finite function, satisfying on each continuity segment a linear ODE with polynomial coefficients.

Many reconstruction problems in signal processing require solution of a certain kind of nonlinear systems of algebraic equations, which we call Prony systems. We study these systems from a general perspective, addressing questions of global solvability and stable inversion. Of special interest are the so-called "near-singular" situations, such as a collision of two closely spaced nodes. We also discuss the problem of reconstructing piecewise-smooth functions from their Fourier coefficients, which is easily reduced by a well-known method of K.Eckhoff to solving a particular Prony system. As we show in the paper, it turns out that a modification of this highly nonlinear method can reconstruct the jump locations and magnitudes of such functions, as well as the pointwise values between the jumps, with the maximal possible accuracy., Comment: arXiv admin note: text overlap with arXiv:1211.0680

Prony mapping provides the global solution of the Prony system of equations \[ \Sigma_{i=1}^{n}A_{i}x_{i}^{k}=m_{k},\ k=0,1,...,2n-1. \] This system appears in numerous theoretical and applied problems arising in Signal Reconstruction. The simplest example is the problem of reconstruction of linear combination of $\delta$-functions of the form $g(x)=\sum_{i=1}^{n}a_{i}\delta(x-x_{i})$, with the unknown parameters $a_{i},\ x_{i},\ i=1,...,n,$ from the "moment measurements" $m_{k}=\int x^{k}g(x)dx.$ Global solution of the Prony system, i.e. inversion of the Prony mapping, encounters several types of singularities. One of the most important ones is a collision of some of the points $x_{i}.$ The investigation of this type of singularities has been started in \cite{yom2009Singularities} where the role of finite differences was demonstrated. In the present paper we study this and other types of singularities of the Prony mapping, and describe its global geometry. We show, in particular, close connections of the Prony mapping with the "Vieta mapping" expressing the coefficients of a polynomial through its roots, and with hyperbolic polynomials and "Vandermonde mapping" studied by V. Arnold., Comment: arXiv admin note: text overlap with arXiv:1301.1187

In this note we present a simple upper bound for the row-wise norm of the inverses of general confluent Vandermonde matrices.

In this paper we provide a reconstruction algorithm for piecewise-smooth functions with a-priori known smoothness and number of discontinuities, from their Fourier coefficients, posessing the maximal possible asymptotic rate of convergence -- including the positions of the discontinuities and the pointwise values of the function. This algorithm is a modification of our earlier method, which is in turn based on the algebraic method of K.Eckhoff proposed in the 1990s. The key ingredient of the new algorithm is to use a different set of Eckhoff's equations for reconstructing the location of each discontinuity. Instead of consecutive Fourier samples, we propose to use a "decimated" set which is evenly spread throughout the spectrum.

We consider the problem of reconstruction of planar domains from their moments. Specifically, we consider domains with boundary which can be represented by a union of a finite number of pieces whose graphs are solutions of a linear differential equation with polynomial coefficients. This includes domains with piecewise-algebraic and, in particular, piecewise-polynomial boundaries. Our approach is based on one-dimensional reconstruction method of [Bat]* and a kind of "separation of variables" which reduces the planar problem to two one-dimensional problems, one of them parametric. Several explicit examples of reconstruction are given. Another main topic of the paper concerns "invisible sets" for various types of incomplete moment measurements. We suggest a certain point of view which stresses remarkable similarity between several apparently unrelated problems. In particular, we discuss zero quadrature domains (invisible for harmonic polynomials), invisibility for powers of a given polynomial, and invisibility for complex moments (Wermer's theorem and further developments). The common property we would like to stress is a "rigidity" and symmetry of the invisible objects. * D.Batenkov, Moment inversion of piecewise D-finite functions, Inverse Problems 25 (2009) 105001, Comment: Proceedings of Complex Analysis and Dynamical Systems V, 2012

In this paper we consider several nonlinear systems of algebraic equations which can be called "Prony-type". These systems arise in various reconstruction problems in several branches of theoretical and applied mathematics, such as frequency estimation and nonlinear Fourier inversion. Consequently, the question of stability of solution with respect to errors in the right-hand side becomes critical for the success of any particular application. We investigate the question of "maximal possible accuracy" of solving Prony-type systems, putting stress on the "local" behavior which approximates situations with low absolute measurement error. The accuracy estimates are formulated in very simple geometric terms, shedding some light on the structure of the problem. Numerical tests suggest that "global" solution techniques such as Prony's algorithm and ESPRIT method are suboptimal when compared to this theoretical "best local" behavior.

This paper presents some results on a well-known problem in Algebraic Signal Sampling and in other areas of applied mathematics: reconstruction of piecewise-smooth functions from their integral measurements (like moments, Fourier coefficients, Radon transform, etc.). Our results concern reconstruction (from the moments or Fourier coefficients) of signals in two specific classes: linear combinations of shifts of a given function, and "piecewise $D$-finite functions" which satisfy on each continuity interval a linear differential equation with polynomial coefficients. In each case the problem is reduced to a solution of a certain type of non-linear algebraic system of equations ("Prony-type system"). We recall some known methods for explicitly solving such systems in one variable, and provide extensions to some multi-dimensional cases. Finally, we investigate the local stability of solving the Prony-type systems.

Accurate reconstruction of piecewise-smooth functions from a finite number of Fourier coefficients is an important problem in various applications. The inherent inaccuracy, in particular the Gibbs phenomenon, is being intensively investigated during the last decades. Several nonlinear reconstruction methods have been proposed, and it is by now well-established that the "classical" convergence order can be completely restored up to the discontinuities. Still, the maximal accuracy of determining the positions of these discontinuities remains an open question. In this paper we prove that the locations of the jumps (and subsequently the pointwise values of the function) can be reconstructed with at least "half the classical accuracy". In particular, we develop a constructive approximation procedure which, given the first $k$ Fourier coefficients of a piecewise-$C^{2d+1}$ function, recovers the locations of the jumps with accuracy $\sim k^{-(d+2)}$, and the values of the function between the jumps with accuracy $\sim k^{-(d+1)}$ (similar estimates are obtained for the associated jump magnitudes). A key ingredient of the algorithm is to start with the case of a single discontinuity, where a modified version of one of the existing algebraic methods (due to K.Eckhoff) may be applied. It turns out that the additional orders of smoothness produce a highly correlated error terms in the Fourier coefficients, which eventually cancel out in the corresponding algebraic equations. To handle more than one jump, we propose to apply a localization procedure via a convolution in the Fourier domain.

We consider the problem of exact reconstruction of univariate functions with jump discontinuities at unknown positions from their moments. These functions are assumed to satisfy an a priori unknown linear homogeneous differential equation with polynomial coefficients on each continuity interval. Therefore, they may be specified by a finite amount of information. This reconstruction problem has practical importance in Signal Processing and other applications. It is somewhat of a ``folklore'' that the sequence of the moments of such ``piecewise D-finite''functions satisfies a linear recurrence relation of bounded order and degree. We derive this recurrence relation explicitly. It turns out that the coefficients of the differential operator which annihilates every piece of the function, as well as the locations of the discontinuities, appear in this recurrence in a precisely controlled manner. This leads to the formulation of a generic algorithm for reconstructing a piecewise D-finite function from its moments. We investigate the conditions for solvability of the resulting linear systems in the general case, as well as analyze a few particular examples. We provide results of numerical simulations for several types of signals, which test the sensitivity of the proposed algorithm to noise.