Sparse phase retrieval of one-dimensional signals by Prony's method

0.32 MB·

by  Robert Beinert

Sparse phase retrieval of one-dimensional signals by Prony’s method
Robert Beinert1 and Gerlind Plonka2
1Institut fürMathematikundWissenschaftliches Rechnen
Karl-Franzens-UniversitätGraz
2InstitutfürNumerischeundAngewandte Mathematik
Georg-August-Universität Göttingen
7
1
Abstract: In this paper, we show that sparse signals f representable as a linear combination
0
of a finite number N of spikes at arbitrary real locations or as a finite linear combination of
2
B-splines of order m with arbitrary real knots can be almost surely recovered from O(N2)
n intensity measurements |F[f](ω)|2 up to trivial ambiguities. The constructive proof consists
a
J of two steps, where in the first step the Prony method is applied to recover all parameters of
1 the autocorrelation function and in the second step the parameters of f are derived. Moreover,
3 we present an algorithm to evaluate f from its Fourier intensities and illustrate it at different
numerical examples.

]
A Key words: Sparse phase retrieval; sparse signals, non-uniform spline functions; finite support;
N Prony’s method
.
h AMS Subject classifications: 42A05, 94A08, 94A12
t
a
m
1. Introduction
[

1
Phase retrieval problems occur in many scientific fields, particularly in optics and com-
v
7 munications. They have a long history with rich literature regarding uniqueness of solu-
4 tions and existence of reliable algorithms for signal reconstruction, see e.g. [SEC+15]
9
and references therein. Usually, the challenge in solving one-dimensional phase retrieval
8
0 problems is to overcome the strong ambiguousness by determining appropriate further
1. information on the solution signal. Previous literature on characterization of ambiguities
0 of the phase retrieval problem with given Fourier intensities is often concerned with the
7 discrete problem, where a signal x in RN or CN has to be recovered. For an overview
1
: on the complete characterization of nontrivial ambiguities is this discrete case as well as
v
on appropriate additional signal information we refer to our survey [BP15a] and further
i
X
recent results in [BP17, Bei17a, Bei17b].
r
a Contributionof thispaper. Inthispaper,weconsiderthecontinuousone-dimensional
sparse phase retrieval problem to reconstruct a complex-valued signal from the modulus
of its Fourier transform. Applications of this problem occur in electron microscopy,
wave front sensing, laser optics [SST04, SSD+06] as well as in X-ray crystallography and
speckle imaging [RCLV13]. For the posed problem, we will show that for sparse signals
the givenFourier intensities arealready sufficient foran almostsureunique recovery, and
we will give a construction algorithm to recover f.
2
RobertBeinertandGerlindPlonka
We assume that the sparse signal is either of the form
N
(0)
f(t)= c δ(t−T ) (1.1)
j j
j=1
X
or, for m > 0,
N
(m)
f(t)= c B (t) (1.2)
j j,m
j=1
X
with c(m) ∈ C, T ∈ R for j = 1,...,N, where δ denotes the Delta distribution, and B
j j j,m
is the B-spline of order m being determined by the (real) knots T < T < ... <T .
j j+1 j+m
We want to recover these signals from the Fourier intensities |f(ω)|2 and will show that
only O(N2) samples are needed to recover f, i.e. all coefficients c(m), j = 1,...,N
j
and knots T , j = 1,...,N +m, almost surely up to trivial ambbiguities. The proposed
j
procedure is constructive and consists in two steps. In a first step, we employ Prony’s
method in oder to recover all parameters of the (squared) Fourier intensity function
|F[f](ω)|2. In a second step, we recover the parameters T and the complex coefficients
j
c that determine the desired signal.
j
Related work on sparse phase retrieval. While the general phase retrieval problem
has been extensively studied for a long tome, the special case of sparse phase retrieval
grew to a strongly emerging field of research only recently, particularly often connected
with ideas from compressed sensing. Most of the papers consider a discrete setting,
where the N-dimensional real or complex k-sparse vector x has to be reconstructed from
measurements of the more general form |ha ,xi|2 with vectors a forming the rows of a
j j
measurement matrix A ∈ CM×N. The needed number M of measurements depends on
the sparsity k.
If A presents rows of a Fourier matrix, this setting is close to the sparse phase re-
trieval problem considered in optics, seee.g. [JOH13]. Here the problem is first rewritten
as (non-convex) rank minimization problem, then a tight convex relaxation is applied
and the optimization problem is solved by a re-weighted l -minimization method. The
1
related approach in [ESM+15] employs the magnitudes of the short-time Fourier trans-
form and applies the occurring redundancy for unique recovery of the desired signal. A
corresponding reconstruction algorithm is then based on an adaptation of the GESPAR
algorithm in [SBE14].
In [LV13], the measurement matrix A is taken with random rows and the PhaseLift
approach [CSV13] leads to a convex optimization problem that recovers the sparse solu-
tion with high probability. Employing a thresholded gradient descent algorithm to a
non-convex empirical risk minimization problem that is derived from the phase retrieval
problem, Cai et al. [CLM16] have established the minimax optimal rates of convergence
for noisy sparse phase retrieval under sub-exponential noise.
Other papers rely on the compressed sensing approach to construct special frame vec-
tors a to ensure uniqueness of the phase retrieval problem with high probability, where
j
the number of needed vectors is O(k), see e.g. [WX14, OE14, IVW17].
3
Sparsephaseretrievalofone-dimensionalsignalsbyProny’smethod
We would like to emphasize that all approaches employing general or random meas-
urement matrices in phase retrieval are quite different in nature from our phase retrieval
problem based on Fourier intensity measurements. In this paper, we want to stick on
considering Fourier intensity measurements because of their particular relevance in prac-
tice.
Early attempts to exploitsparsity of adiscrete signalforunique recovery usingFourier
intensities go back to unpublished manuscripts by Yagle [Yaga, Yagb], where a variation
of Prony’s method is applied in a non-iterative algorithm to sparse signal and image
reconstruction. Unfortunately, the algorithm proposed there not always determines the
signal support correctly.
The continuous one-dimensional phase retrieval problem has been rarely discussed in
the literature, see [Wal63, Hof64, RCLV13, Bei17b, BP15b]. In the preprint [RCLV13],
the authors also considered the recovery of sparse continuous signals of the form (1.1).
However, in that paper the sparse phase retrieval problem is in turn transferred into a
turnpike problem that is computationally expensive to solve. Moreover there exist cases,
where a unique solution cannot be found, see [Blo75]. Our method circumvents this
problem by proposing an iterative procedure to fix the signal support (resp. the knots of
the signal represented as a B-spline function) where the corresponding signal coefficients
are evaluated simultaneously.
Organization of this paper. InSection 2, we shortly recall the mathematical formula-
tionoftheconsidered sparsephaseretrievalproblemandthenotionoftrivialambiguities
of the phase retrieval problem that always occur.
Section 3 is devoted to the special case of phase retrieval for signals of the form (1.1).
Using Prony’s method, we give a constructive proof for the unique recovery of the N-
sparsesignal f upto trivialambiguities using 3/2N(N−1)+1Fourier intensity measure-
ments. Here we have to assume that the knot differences T −T are pairwise different.
j k
In Section 4, the ansatz is generalized to the unique recovery of spline functions of
the form (1.2) where we need to employ 3/2(N +m)(N +m−1)+1 Fourier intensity
measurements. In Section 5, we present an explicit algorithm for the considered sparse
phase retrieval problem and illustrate it at different examples.
2. Trivial ambiguities of the phase retrieval problem
We wish to recover an unknown complex-valued signal f: R → C of the form (1.1) or
(1.2) with compact support from its Fourier intensity |F[f]| given by
∞
|F[f](ω)| := f(ω) := f(t)e−iωtdt (ω ∈ R).
(cid:12) (cid:12) (cid:12)−Z∞ (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
b
(cid:12) (cid:12) (cid:12) (cid:12)
Unfortunately, the recovery of the signal f is complicated because of the well-known
ambiguousness of the phase retrieval problem. Transferring [BP15a, Proposition 2.1] to
our setting, we can recover f only up to the following ambiguities.
4
RobertBeinertandGerlindPlonka
Proposition 2.1. Let f be of a signal of the form (1.1) or a non-uniform spline function
of the form (1.2). Then
(i) the rotated signal eiαf for α∈ R,
(ii) the time shifted signal f(·−t ) for t ∈R,
0 0
(iii) and the conjugated and reflected signal f(−·)
have the same Fourier intensity |F[f]|.
Proof. Applying the properties of the Fourier transform, we have
(i) F[eiαf]= eiα F[f];
(ii) F[f(·−t )] = e−iωt0 F[f];
0
(iii) F[f[−·]]= F[f].
Considering the absolute value of each equation yields the assertion. (cid:3)
Although the ambiguities in Proposition 2.1 always occur, they are of minor interest
because of their close relation to the original signal. For this reason, we call ambiguities
caused by rotation, time shift, conjugation and reflection, or by combinations of these
trivial. Inthefollowing, we willshowthatfortheconsidered sparsesignals theremaining
non-trivial ambiguities only occur in rare cases.
3. Phase retrieval for distributions with discrete support
Initially, we restrict ourselves to the recovery of signals f of the form (1.1) with complex-
(0)
valued coefficients c and spike locations T < ··· < T .
j 1 N
N
f(ω) = c(0)e−iωTj (ω ∈ R),
j
j=1
X
b
and the known squared Fourier intensity |F[f]|2 can be represented by
N N
2
f(ω) = c(0)c(0)e−iω(Tj−Tk). (3.1)
j k
(cid:12) (cid:12) Xj=1Xk=1
(cid:12) (cid:12)
b
(cid:12) (cid:12)
Thus, in order to recover f being determined by the coefficients c(0) ∈ C and the knots
j
T ∈ R, j = 1,...,N, we will recover all parameters of the exponential sum in (3.1) in a
j
first step and then derive the desired parameters of f in a second step.
5
Sparsephaseretrievalofone-dimensionalsignalsbyProny’smethod
3.1. First step: Parameter recovery by Prony’s method
Assuming that the non-zero knot differences T − T with j 6= k are pairwise dif-
j k
ferent, and denoting the distinct frequencies T − T in increasing order by τ with
j k ℓ
ℓ = −N(N−1)/2,...,N(N−1)/2, we can rewrite (3.1) as
N(N−1)/2 N(N−1)/2
2
P(ω):= f(ω) = γ e−iωτℓ = γ + γ e−iωτℓ +γ eiωτℓ (3.2)
ℓ 0 ℓ ℓ
(cid:12)(cid:12) (cid:12)(cid:12) ℓ=−NX(N−1)/2 Xℓ=1 (cid:0) (cid:1)
b
(cid:12) (cid:12)
(0) (0)
with the related coefficients γ := c c for the non-zero frequencies τ = T −T and
ℓ j k ℓ j k
γ := N |c(0)|2 for the zero frequency. Since τ = −τ , the coefficients in (3.2) fulfill
0 j=1 j −ℓ ℓ
the conjugated symmetry γ = γ .
P −ℓ ℓ
In order to recover the parameters τ and the unknown coefficients γ from the expo-
ℓ ℓ
nential sum (3.2) we employ Prony’s method [Hil87, PT14]. Let h > 0 be chosen such
that hτℓ < π for all ℓ = 1,...,N(N−1)/2.
2
Using the intensity values P(hk) = f(hk) , k = 0,...,2N(N −1)+1, the unknown
parameters γ and τ in (3.2) can be d(cid:12)etermin(cid:12)ed by exploiting the algebraic Prony poly-
ℓ ℓ
(cid:12) (cid:12)
b
nomial Λ(z) defined by (cid:12) (cid:12)
N(N−1)/2 N(N−1)+1
Λ(z) := z−e−ihτℓ = λ zk, (3.3)
k
ℓ=−NY(N−1)/2(cid:16) (cid:17) Xk=0
where λ denote thecoefficients inthe monomial representation ofΛ(z). Obviously, Λ(z)
k
is always a monic polynomial, which means that λ = 1.
N(N−1)+1
Using the definition of the Prony polynomial Λ(z) in (3.3), we observe that
N(N−1)+1 N(N−1)+1 N(N−1)/2
λ P(h(k+m)) = λ γ e−ih(k+m)τℓ
k k ℓ
Xk=0 Xk=0 ℓ=−NX(N−1)/2
N(N−1)/2
= γ e−ihmτℓΛ e−ihτℓ = 0
ℓ
ℓ=−NX(N−1)/2 (cid:16) (cid:17)
for m = 0,...,N(N − 1). Consequently, the vector of remaining coefficients λ :=
(λ ,...,λ )T oftheProny polynomial Λ(z)canbedeterminedbysolvingthelinear
0 N(N−1)
equation system
Hλ = −h (3.4)
with H := (P(h(k+m)))N(N−1) and h := (P(h(N(N −1)+1+m)))N(N−1). Since the
m,k=0 m=0
Hankel matrix H can be written as
H = VTdiag γ ,...,γ V
−N(N−1)/2 N(N−1)/2
(cid:0) (cid:1)
6
RobertBeinertandGerlindPlonka
with the Vandemonde matrix V := (e−hkτℓ)N(N−1)/2,N(N−1)+1, the linear equation system
ℓ=−N(N−1)/2,k=0
(3.4)possessesauniquesolutionifandonlyiftheunimodularvaluese−ihτℓ differpairwise
for ℓ = −N(N−1)/2,...,N(N−1)/2. This assumption has been ensured by choosing an h
such that hτ ∈ (−π,π), since the τ had been supposed to be pairwise different.
ℓ ℓ
Knowing the coefficients λ of Λ(z), we can determine the unknown frequencies τ
k ℓ
by evaluating the roots of the Prony polynomial (3.3). The coefficients γ can now be
ℓ
computed by solving the over-determined equation system
N(N−1)/2
γ e−ihkτℓ = P(hk) (k = 0,...,2N(N −1)+1) (3.5)
ℓ
ℓ=−NX(N−1)/2
with a Vandermonde-type system matrix.
The procedure summarized above is the usual Prony method, adapted to the non-
negative exponential sum P(ω) in (3.2). In the numerical experiments in Section 5,
we will apply the approximate Prony method (APM) in [PT10]. APM is based on the
above considerations but it is numerically more stable and exploits the special properties
γ−ℓ = γℓ and τ−ℓ = −τℓ for ℓ = 0,...,N(N−1)/2.
Let us now investigate the question, how many intensity values are at least necessary
for the recovery of P(ω) in (3.2). Counting the number of unknowns of P(ω) in (3.2),
we only need to recover the 3/2N(N −1)+1 real values γ0 and Reγℓ, Imγℓ, τℓ, for ℓ =
1,...N(N−1)/2. We will show now that using the special structure of the real polynomial
P(ω) in (3.2) and of the Prony polynomial Λ(z) in (3.3), we indeed need only 3/2N(N −
1)+1 exact equidistant real measurements P(kh), k = 0,...,3/2N(N −1) to recover all
parameters determining P(ω). This can be seen as follows.
Reconsidering Λ(z) in (3.3) with τ = 0 and τ = −τ , we obtain
0 ℓ −ℓ
N(N−1)/2
Λ(z) = (z−1) z−eihτℓ z−e−ihτℓ
ℓY=1 (cid:16) (cid:17)(cid:16) (cid:17)
N(N−1)/2 N(N−1)+1
= (z−1) z2−2zcos(hτ )+1 = λ zk,
ℓ k
ℓ=1 k=0
Y (cid:0) (cid:1) X
where all occurring coefficients λ are real. Moreover, since
k
N(N−1)/2
z−(N(N−1)+1)/2Λ(z) = (z1/2−z−1/2) z−2cos(hτ )+z−1
ℓ
ℓ=1
Y (cid:0) (cid:1)
is antisymmetric, it follows that
λN(N−1)+1−k = −λk (k =0,...,N(N−1)/2),
and particularly λ = −λ = 1. In order to determine the unknown coefficients
N(N−1)+1 0
7
Sparsephaseretrievalofone-dimensionalsignalsbyProny’smethod
λk, k = 1,...,N(N−1)/2 of
N(N−1)/2
Λ(z) = λ zk −zN(N−1)+1−k ,
k
Xk=0 (cid:16) (cid:17)
we employ (3.2) and observe that for m = 0,...,N(N−1)/2−1,
N(N−1)/2
λ [P(h(k+m))−P(h(N(N −1)+1+m−k))]
k
k=0
X
N(N−1)/2 N(N−1)/2
= λ γ e−ih(k+m)τℓ −e−ih(N(N−1)+1+m−k)τℓ
k ℓ

Xk=0 Xℓ=1 (cid:16) (cid:17)
N(N−1)/2 
+ γ eih(k+m)τℓ −eih(N(N−1)+1+m−k)τℓ
ℓ 
Xℓ=1 (cid:16) (cid:17)
N(N−1)/2 N(N−1)/2 
= γ e−ihmτℓ λ e−ihkτℓ −e−ih(N(N−1)+1−k)τℓ
ℓ k
Xℓ=1 Xk=0 (cid:16) (cid:17)
N(N−1)/2 N(N−1)/2
+ γ eihmτℓ λ eihkτℓ −eih(N(N−1)+1−k)τℓ
ℓ k
Xℓ=1 Xk=0 (cid:16) (cid:17)
N(N−1)/2 N(N−1)/2
= γ e−ihmτℓΛ(e−ihτℓ)+ γ eihmτℓΛ(eihτℓ)= 0.
ℓ ℓ
ℓ=1 ℓ=1
X X
Therefore, the vector of unknown coefficients λ := (λ ,...,λ )T can be already
1 N(N−1)/2
evaluated from the system
N(N−1)/2
λ [P(h(k+m))−P(h(N(N −1)+1+m−k))]
k
k=1
X
= [P(hm)−P(h(N(N −1)+1+m))] (m =0,...,N(N−1)/2−1).
The parameters τ are then extracted from the zeros of Λ(z), and the coefficients γ ,
ℓ ℓ
ℓ = 0,...,N(N−1)/2, are computed as in (3.5) but with k = 0,...,3/2N(N −1).
3.2. Second step: Unique signal recovery
Having determined the parameters τ as well as the corresponding coefficients γ of (3.2),
ℓ ℓ
(0)
we want to reconstruct the parameters T and c , j = 1,...,N, of f in (1.1)in a second
j j
step.
Theorem 3.1. Let f be a signal of the form (1.1), whose knot differences T −T differ
j k
(0) (0)
pairwise for j,k ∈ {1,...,N} with j 6= k, and whose coefficients satisfy |c |=6 |c |.
1 N
8
RobertBeinertandGerlindPlonka
Further, let h be a step size such that h(T −T ) ∈ (−π,π) for all j,k. Then f can be
j k
uniquely recovered from its Fourier intensities |F[f](hℓ)| with ℓ = 0,...,3/2N(N−1) up
to trivial ambiguities.
Proof. Applying Prony’s method to the given data |F[f](hℓ)|, we can compute the
frequenciesτ andtherelatedcoefficients γ ofthesquaredFourier intensity (3.2). Again,
ℓ ℓ
we assume that the frequencies τ occur in increasing order and, further, denote the list
ℓ
of positive frequencies by T := {τ }N(N−1)/2.
ℓ ℓ=1
Obviously, the maximal distance τ now corresponds to the length T −T of
N(N−1)/2 N 1
the unknown f in (1.1). Due to the trivial shift ambiguity, we can assume without
loss of generality that T = 0 and T = τ . Further, the second largest distance
1 N N(N−1)/2
τ corresponds either to T −T or to T −T . Due to the trivial reflection
(N(N−1)/2)−1 N−1 1 N 2
andconjugationambiguity,wecanassumethatT = τ . Bydefinition,there
N−1 (N(N−1)/2)−1
exists a τℓ∗ > 0 in our sequence of parameters T such that τℓ∗+τ(N(N−1)/2)−1 = τN(N−1)/2,
and τℓ∗ hence corresponds to the knot difference TN −TN−1. Thus, we obtain
(0) (0) (0) (0) (0) (0)
cN c1 = γN(N−1)/2, cN−1c1 = γ(N(N−1)/2)−1, and cN cN−1 = γℓ∗.
These equations lead us to
c(0) = γN(N−1)/2, c(0) = γ(N(N−1)/2)−1,
N c(0) N−1 c(0)
1 1
and thus to
c(0) 2 = γN(N−1)/2γ(N(N−1)/2)−1.
1 γ ∗
ℓ
(cid:12) (cid:12)
(cid:12) (cid:12)
Since f can only be recovere(cid:12)d up(cid:12) to a global rotation, we can assume that c(0) is real
1
(0) (0) (0)
and non-negative, which allows us to determine the coefficients c , c , and c in a
1 N N−1
unique way.
To determine the remaining coefficients and support knots of f, we notice that the
third largest distance τ corresponds either to T −T or to T −T . As
(N(N−1)/2)−2 N 2 N−2 1
before, we always find a frequency τℓ∗ such that τ(N(N−1)/2)−2+τℓ∗ = τN(N−1)/2.
Case 1: If τ = T −T , then we have
(N(N−1)/2)−2 N 2
τℓ∗ = τN(N−1)/2−τ(N(N−1)/2)−2 = (TN −T1)−(TN −T2)= T2−T1
(0) (0) (0) (0)
with the related coefficient γℓ∗ = c2 c1 . Moreover, we have γ(N(N−1)/2)−2 = cN c2 such
that
c(0) = γℓ∗ = γ(N(N−1)/2)−2. (3.6)
2 (0) (0)
c c
1 N
Case 2: If τ = T −T , then we have
(N(N−1)/2)−2 N−2 1
τℓ∗ = τN(N−1)/2−τ(N(N−1)/2)−2 = (TN −T1)−(TN−2−T1)= TN −TN−2
9
Sparsephaseretrievalofone-dimensionalsignalsbyProny’smethod
(0) (0) (0) (0)
with the related coefficient γℓ∗ = cN cN−2 and γ(N(N−1)/2)−2 = cN−2c1 . Thus,
c(0) = γℓ∗ = γ(N(N−1)/2)−2. (3.7)
N−2 (0) (0)
c c
N 1
However, only one of the two equalities in (3.6) and (3.7) can be true, since if both were
(0) (0) (0) (0)
true then γℓ∗cN = c1 γ(N(N−1)/2)−2 and c1 γℓ∗ = cN γ(N(N−1)/2)−2 lead to
c(N0) = γ(N(N−1)/2)−2 = c(10)
(cid:12)(cid:12)c(10)(cid:12)(cid:12) (cid:12) γℓ∗ (cid:12) (cid:12)(cid:12)c(N0)(cid:12)(cid:12)
(cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
contradicting the assumption(cid:12)tha(cid:12)t |c0(cid:12)| =6 |c0|. Con(cid:12)sequ(cid:12)ently(cid:12), either the equation in (3.6)
N 1
(0)
or theequation in(3.7)holds true and we can either determine T with c or T with
2 2 N−2
(0)
c . Removing all parameters τ from the sequence of distances T that correspond to
N−2 ℓ
the differences T −T of the recovered knots, we can repeat this approach to find the
j k
remaining coefficients and knots of f inductively. (cid:3)
If we identify the space of complex-valued signals of the form (1.1) with the real space
R3N, the condition that two knot differences T −T and T −T are equal for fixed
j1 k1 j2 k2
indicesj ,j ,k ,andk definesahyperplanewithLebesguemeasurezero. Ananalogous
1 2 1 2
(0) (0)
observation follows for the condition |c |= |c |. The signals excluded in Theorem 3.1
1 N
hence form a negligible null set.
Corollary 3.2. Almost all signals f in (1.1) can be uniquely recovered from their Fourier
intensities |F[f]| up to trivial ambiguities.
Remark 3.3. 1. Since the proof of Theorem 3.1 is constructive, it can be used to recover
an unknown signal (1.1) analytically and numerically. If the number N of spikes is
known beforehand then the assumption of Theorem 3.1 can be simply checked during
the computation. If the assumption regarding pairwise different distances T − T is
j k
not satisfied, then the application of Prony’s method in the first step yields less than
N(N −1)+1 pairwise distinct parameters τ . The second assumption |c0 | =6 |c0| can be
ℓ N 1
(0) (0) (0)
verified in the second step, where c , c , and c are evaluated.
1 N−1 N
2. A similar phase retrieval problem had been transferred to a turnpike problem
in [RCLV13]. The turnpike problem deals with the recovery of the knots T from an
j
unlabeled set of distances. Although this problem is solvable under certain conditions, a
backtracing algorithm can have exponential complexity in the worst case, see [LSS03].
4. Retrieval of spline functions with arbitrary knots
Inthissection, wegeneralizeourfindings tosplinefunctions oforderm ≥ 1. Letusrecall
that the B-splines B in (1.2) being generated by the knot sequence T < ··· < T
j,m 1 N+m
are recursively defined by
B (t):= t−Tj B (t)+ Tj+m−t B (t)
j,m Tj+m−1−Tj j,m−1 Tj+m−Tj+1 j+1,m−1
10
RobertBeinertandGerlindPlonka
with
1 t ∈ [T ,T ),
B (t):= 1 (t):= j j+1
j,1 [Tj,Tj+1) (0 else,
see for instance [Boo78, p. 131]. Further, we notice that for 0 ≤ k ≤ m − 2 the kth
derivative of the spline f in (1.2) is given by
dk N+k
(m−k)
f(t)= c B (t), (4.1)
dtk j j,m−k
j=1
X
(m−k)
where the coefficients c are recursively defined by
j
(m−k+1) (m−k+1)
c −c
(m−k) j j−1
c := (m−k) (j = 1,...,N +k),
j T −T
j+m−k j
(m−k+1) (m−k+1)
with the convention that c = c = 0, see [Boo78, p. 139]. For k = m−1,
0 N+k
equation (4.1) coincides with a step function, i.e., with the right derivative of the linear
spline f(m−2). Further, in a distributional manner, the mth derivative of f is given by
dm N+m
(0)
f(t)= c δ(t−T ) (4.2)
dtm j j
j=1
X
with the coefficients
(0) (1) (0) (1) (0) (1) (1)
c := c , c := −c , c := c −c (j = 2,...,N +m−1),
1 1 N+m N+m−1 j j j−1
and the Dirac delta distribution δ.
Applying the Fourier transform to (4.2), we now obtain
N+m
f(m)(ω) = (iω)mf(ω)= c(0)e−iωTj. (4.3)
j
j=1
X
d b
and thus
N+mN+m
2
ω2m f(ω) = c(0)c(0)e−iω(Tj−Tk). (4.4)
j k
(cid:12) (cid:12) Xj=1 Xk=1
(cid:12) (cid:12)
b
(cid:12) (cid:12)
Since the exponential sum on the right-hand side of (4.4) has exactly the same struc-
ture as the exponential sum in (3.2), we can immediately generalize Theorem 3.1 by
considering
(N+m)(N+m−1)/2
2
P(ω) := ω2m f(ω) = γ e−iωτℓ. (4.5)
ℓ
(cid:12) (cid:12) ℓ=−(N+mX)(N+m−1)/2
(cid:12) (cid:12)
b
(cid:12) (cid:12)