2.2 Methodical Background
2.2.1 Reconstruction Procedure for One Frequency
We begin by supposing that the function f possesses only a single energetic frequency ω∈
−N
2
+ 1, . . . ,N
2 with corresponding Fourier coefficientcω :=cω(f)∈C\ {0}.
We will later extend the reconstruction procedure to several energetic frequencies. Then f is of the form
f(x) =cω·eiωx.
Thus, f is completely determined if we can recover ω and cω. We will do this utilizing vectorsaM of equidistant samples of f, where we keep the slightly unintuitive notation ofaM used in [Iwe10, Iwe13, BZI19, Bit17c].
Definition 2.5 (Vector of Equidistant Samples)Let f ∈ C2π and M ∈ N. By aM ∈CM we denote the vector ofM equidistant samples of f, i.e.,
aM := aMj M−1 j=0 :=
f
2πj M
M−1 j=0
.
Recall that we learned in Section 1.2.1 that we can recover f by applying the FFT to the vector aN of N equidistant samples of f in O(NlogN) time. The central idea of Algorithm 2 in [Iwe10] is to obtain a method with shorter runtime than the FFT by applying the CDFT to several vectors of equidistant samples with lengths that are small compared to the bandwidthN. We choose the CDFT here so that we can consider func-tions with energetic frequencies contained in the interval
−N
2
+ 1, . . . ,N
2 centered around0. In order for this idea to work, the vector lengths have to be chosen in a very specific way, which we will discuss later on. First, we will investigate what happens if we apply the CDFT to the vectoraM of M N equidistant samples off. We find that
acMν = 1 M
M−1
X
j=0
e−2πijνM ·cωe2πijωM
= cω M
M−1
X
j=0
ωMj(ν−ω)
=
(cω if ν ≡ω mod M,
0 otherwise, ∀ν ∈
− M
2
+ 1, . . . , M
2
. (2.1)
Hence, acM has exactly one nonzero entry, namely the Fourier coefficient corresponding to ω, and its index is the residue of ω moduloM,
acMω mod M =cω = 1 2π
2π
Z
0
cωeiωx·e−iωxdx.
Consequently, (2.1) gives us the value of the Fourier coefficient cω and the residue of the frequencyω moduloM while actually knowing neither cω norω. If we compute the CDFTs of several such vectors with different lengths, we obtain a system of simultaneous congruencies. Under certain conditions on the occurring moduli, i.e., the lengths of the vectors of equidistant samples, such a system can be solved via the well-known Chinese Remainder Theorem. For a proof and more details, see [Lan05], Chapter I, §4.
2.2 Methodical Background
Theorem 2.6 (Chinese Remainder Theorem (CRT)) Let M1, . . . , ML be pairwise relatively prime integers and N ≤ QL
l=1Ml. Further, let rl ∈ n
−l
Ml
2
m
+ 1, . . . ,j
Ml
2
ko for all l ∈ {1, . . . , L}. Then there exists a unique solution modulo N of the system of simultaneous congruencies
x≡r1 mod M1, ...
x≡rL mod ML.
The unique solution of such a system can be computed using the following algorithm adapted from the implementationChineseReminGAP, an open source system for com-putational discrete algebra (see [The18a]).
Algorithm 1 CRT Reconstruction Algorithm
Input: Residuesr1, . . . , rL modulo pairwise relatively prime integersM1, . . . , ML. Output: Integer ν ∈ {0, . . . , M −1} such that ν ≡ rl mod Ml ∀l ∈ {1, . . . , L} and
M ≥QL
l=1Ml.
1: Setl= 1,ν =r1 andM =M1.
2: whilel < L do
3: Setl=l+ 1.
4: (g, u, v)←extended_gcd(M, Ml)
5: if g6= 1 and rl−ν mod g6= 0 then
6: Error .The residues have to be equal modulo g.
7: end if
8: ν =M
rl−ν g ·u
mod Ml
+ν
9: M = Mgl·M
10: end while
11: ν=ν mod M Output: ν,M.
The function extended_gcd in line 4 of Algorithm 1 computes the greatest common divisorg of two integersaand busing the extended Euclidean algorithm, as well as two integersu and v satisfying Bézout’s identity, see [Bos06], Chapter 2.4, Theorem 15.
Theorem 2.7 (Euclidean Algorithm) Let R be a Euclidean ring. For two elements a, b∈R\ {0} consider the sequence z0, z1, . . .∈R which is given inductively by
z0=a, z1=b, zk+1=
(zk−1 mod zk if zk 6= 0,
0 otherwise.
Then there exists a smallest index n∈Nsuch that zn+1 = 0. It satisfies zn= gcd(a, b).
Furthermore, there exists an explicit representation of the greatest common divisor of a and b in the form
g= gcd(a, b) =ua+vb for some u, v∈R, which is known as Bézout’s identity.
It can be shown that for a > b the extended Euclidean algorithm has a runtime of O(logb); see, e.g., [CLRS09], Chapter 31.2. Thus, Algorithm 1 requiresO
PL
l=1logMl arithmetical operations, as we will show in the proof of Theorem 2.22.
If acM is known for sufficiently many pairwise relatively prime moduli Ml N, l ∈ {1, . . . , L}, with N ≤ QL
l=1Ml, the CRT implies that we can uniquely recover the energetic frequencyω. The corresponding Fourier coefficientcω is already given by (2.1),
adMlωmod Ml=cω ∀l∈ {1, . . . , L}.
This means that, instead of computing the CDFT of length N of acN, it suffices to calculate L CDFTs of length M1, . . . , ML N and reconstruct ω from its residues moduloM1, . . . , ML.
Let us illustrate the reconstruction procedure from Algorithm 1 by an example.
Example 2.8 Letf ∈C2π,f(x) =ei·210x with bandwidthN = 1,000. According to the CRT, the single energetic frequencyω is uniquely determined modulo N by its residues moduloM1= 10, M2 = 11andM3 = 13, asM1, M2andM3 are pairwise relatively prime and their product is 1,430. Locating the nonzero entry ofac10,ac11 andac13, we find
ac100 = 1 ⇒ ω≡0 mod 10 ⇒ r1:= 0, ac111 = 1 ⇒ ω≡1 mod 11 ⇒ r2:= 1 and ac132 = 1 ⇒ ω≡2 mod 13 ⇒ r3:= 2.
Now we can recover ω from its residues r1, r2 and r3 using Algorithm 1. We begin by settingM :=M1 = 10andω :=r1= 0. Then, with
gcd(M, M2) = 1 =−1·10 + 1·11, it follows thatu=−1, and thus obtain
ω:=M ·(((r2−ω)·u) mod M2) +ω
= 10·(((1−0)·(−1)) mod 11) + 0
= 10·10
= 100.
The frequencyω= 100satisfies thatω≡0 mod 10andω≡1 mod 11. Now we update M :=M2·M = 110. Since
gcd(M, M3) = 1 =−2·110 + 17·13, we have thatu=−2, so we redefine
ω:= 110·(((2−100)·(−2)) mod 13) + 100
= 110·(196 mod 13) + 100
= 210.
Finally, we setM :=M3·M = 1,430. Asω ∈
−N
2
+ 1, . . . ,N
2 ={−499, . . . ,500}, the output of Algorithm 1 isω= 210and M = 1,430. The obtained frequency ω indeed satisfies the required congruencies ω ≡ 0 mod 10, ω ≡ 1 mod 11 and ω ≡ 2 mod 13.
2.2 Methodical Background
The Fourier coefficient cω is given by (2.1) as, e.g., cω =ac100 = 1.
For the above computation three CDFTs of lengths 10, 11 and 13 were necessary; hence, only 34 instead of N = 1,000 samples of f were used. Moreover, as we discussed in Section 1.1.1, there exist fast algorithms for the DFT of vectors of arbitrary length, and, by Remark 1.9, also for the CDFT, so the three CDFTs can be computed in O
P3
l=1MllogMl
time, instead of calculating one CDFT with complexityO(NlogN).
We will show in the proof of Theorem 2.22 that the frequency reconstruction needs O
P3
l=1logMl
arithmetic operations, which is insignificant compared to the
compu-tational costs of the CDFTs. ♦
However, as soon as the function we aim to recover has more than one energetic fre-quency, the residues of these frequencies can coincide modulo various integers. If we choose the moduliMl arbitrarily, it can happen thatω1≡ω2 mod M1, so it is impossi-ble to distinguish these two frequencies moduloM1. Furthermore, we cannot determine the values of the Fourier coefficients fromadM1. This means that in the case of aB-sparse function, we have to choose the moduliMlcarefully in order to avoid ambiguities. With-out a priori knowledge of the energetic frequencies, guaranteeing unique recoverability requires a more involved reconstruction procedure.