5.6 Sparse Fast IDFT and Sparse Fast IDCT-II
5.6.1 Detecting the Support Sets
cu02k0+1−y\(j+1)2k0+1 <
cu02k0+1+y\(j+1)2k0+1 ,
and y(j+1) =u1 otherwise. The required entry of uc0 can be computed from y(j) using O(m) operations,
cu02k0+1 =
2j+1−1
X
l=0
ω2j+1(2k0+1)lu0l
=
µ(j)+m−1
X
l=µ(j)
ω2j+1(2k0+1)lyl(j)+
2j+1−1−µ(j)
X
l=2j+1−m−µ(j)
ω2j+1(2k0+1)lyl−2(j)j,
since the support of y(j) and thus ofu0 is known from the previous iteration step.
The recovery of y(j+1) fromy(j)and an oddly indexed nonzero entry ofy\(j+1) via the procedure for case B requires at most2msamples ofybandO(m)arithmetical operations, so this approach also contributes sufficiently to obtaining an overall runtime of our sparse IDFT method that is sublinear in the vector length 2N. Note that, as for the method introduced in Section 5.5.1, a priori knowledge of the block lengthmofyis not necessary.
5.6 Sparse Fast IDFT and Sparse Fast IDCT-II
In Section 5.5 we presented the procedures necessary to derive the new sparse fast IDFT for vectors y ∈ R2N, N = 2J−1, that have a reflected block support and satisfy (5.4).
Using Lemma 5.4, we obtain at the same time a new sparse fast IDCT-II algorithm for vectors with one-block support. Note that neither of the procedures for reconstructing y(j+1) fromy(j)andybintroduced in Section 5.5 requires a priori knowledge of the length of the blocks in y. However, what these procedures do require is the support structure of the vector y(j) from the previous iteration step. More precisely, we need to detect whether y(j) has a one-block or a two-block support and determine the corresponding first and last support indices µ(j) and ν(j), as well as the support length m(j) or block length n(j).
5.6.1 Detecting the Support Sets
In Sections 5.5.1 and 5.5.2 we showed how to compute y(j+1) from y(j) and by. Both reconstruction methods introduced above rely heavily on the fact that the support struc-ture of y(j) is already known from the previous step. However, if y(j) has a one-block support, the reconstruction method from Theorem 5.20 does not provide us with both
the block or support length n(j+1) or m(j+1) and the first support index µ(j+1), even though we need both to recover y(j+2) from y(j+1) in the next iteration step. In this section we introduce methods for the stable and efficient detection of the first support index and the block or support lengthn(j+1) or m(j+1). These methods are designed for noisy data, as input data is usually noisy in practical applications, but can of course also be applied to exact data. Note that for noisy data the found block lengths n(j+1) for y(j+1) with two-block support might not be the same as the exact block length m ofx.
For detecting the support sets efficiently, we choose a threshold ε >0 depending on the noise level of the data.
A) One-block support:
A1) y(j) has the one-block support S(j) = Iµ(j)(j),2j−1−µ(j) of length m(j) < 2j centered around the middle of the vector.
Theorem 5.19, caseA1 implies thaty(j+1) has the two-block support S(j+1)=Iµ(j+1)(j+1), ν(j+1)∪I2(j+1)j+1−1−ν(j+1),2j+1−1−µ(j+1)
of lengthn(j+1)≤m(j). Further, we know that the support intervalIµ(j+1)(j+1), ν(j+1) of the first block of y(j+1) is a subset of S(j). By the proof of Theorem 5.19 and the symmetry guaranteed by Lemma 5.9 (iii), the indices corresponding to the signifi-cantly large entries of the first block ofy(j+1) have to be contained in the set
T(0)(j+1):=n k∈n
µ(j), µ(j)+ 1, . . . ,2j −1−µ(j)o :
yk(j+1)
> εo
=:{t1, . . . , tK}, wheret1<· · ·< tK. Thus, we define
µ(j+1):=t1 and n(j+1) :=tK−t1+ 1.
Hence, in order to find the first support index and the support length, we have to inspect the m(j) entries of y(j+1) corresponding to indices in T(0)(j+1), which has a runtime ofO m(j)
.
A2) y(j) has the one-block supportS(j)=I0,2(j)j−1 of lengthm(j) = 2j.
If j < J −1, Theorem 5.19, case A2 yields that y(j+1) has a one-block support whose location and length are unknown. As y(j+1) is symmetric, the indices of its significantly large entries are
T(j+1) :=n
k∈I0,(j+1)2j+1−1 : y(j+1)k
> ε
o
=:
t1, . . . , tK,2j+1−1−tK, . . . ,2j+1−1−t1 =:{t1, . . . , t2K}, where tK ≤ 2j −1 and tK+1 ≥ 2j. If t2K −t1 + 1 > 2j, we set µ(j+1) := 0 and m(j+1) := 2j+1, as we are not able to correctly determine the first support index for support lengths that are greater than half of the vector length. Otherwise, we need to detect whether the support block is centered around the middle or the boundary ofy(j+1). We define
d0 :=tK+1−tK and d1 := (t1−t2K) mod 2j+1,
5.6 Sparse Fast IDFT and Sparse Fast IDCT-II
since not all entries belonging to the support block have to be significantly large.
If d0 < d1, i.e., if the distance between the first significantly large entry to the left and right of the middle of the vector is smaller than the periodic distance between the first and last significantly large entry of the vector, then y(j+1) has a one-block support centered around 2j−1 and 2j. Analogously, if d0 > d1, the support block is centered around 0 and2j+1−1. Ifd0=d1, then y(j+1) must have full support of length m(j+1) = 2j+1, so we can conclude that
µ(j+1) :=
t1 if d0 < d1, tK+1 if d0 > d1, 0 if d0 =d1, and
m(j+1):=
t2K−t1+ 1 if d0 < d1, 2j+1−tK+1+tK+ 1 if d0 > d1, 2j+1 if d0 =d1.
By symmetry of y(j+1) it suffices to findt1, . . . , tK; hence, we only have to inspect 2j =m(j) entries of y(j+1).
In the case that j = J −1, there are some additional possibilities. Due to its symmetry, it suffices to check
T(0)(J) :=
k∈
0, . . . ,2J−2−1 :|yk|> ε =:{t1, . . . , tK} and T(1)(J) :=
k∈
2J−2, . . . ,2J−1−1 :|yk|> ε =:{u1, . . . , uL}.
If T(0)(J) =∅, theny has a one-block support centered around the middle with µ(J):=u1 and m(J):= 2 2J−1−u1
= 2m.
If T(1)(J) =∅, theny has a one-block support centered around the boundary with µ(J):= 2J−1−tK and m(J):= 2 (tK+ 1) = 2m.
Otherwise, there are three possibilities for the support of y, where we do not know the correct one a priori:
(i) S(J):=It(J)
1,2J−1−t1, (ii) S(J):=I2(J)J−1−uL, uL,
(iii) y has a two-block support with two blocks of possibly different lengths and unknown positions.
However, if case A2 occurs for j = J −1, then y(J−1) has a one-block support of full length m(J−1) = 2J−1 ≤ 2m. Consequently, the vector x is not really sparse and the first support index of x might not be uniquely determined. Additionally, the iteration stops with y, so not being able to detect its support uniquely does not pose an algorithmic problem. As we have to check 2J−1 = m(J−1) entries of y(J), this yields an arithmetical complexity of O m(J−1)
.
A3) y(j) has the one-block support S(j) = Iµ(j)(j),2j−1−µ(j) of length m(j) < 2j centered around the boundary of the vector andj < J−1.
It follows from Theorem 5.19, caseA3 thaty(j+1) has a one-block support of length m(j+1) :=m(j) withµ(j+1) =µ(j) or µ(j+1) = 2j +µ(j). We compare the entries at the possible locations of the support block and set
e0 :=
µ(j)+m(j)−1
X
k=µ(j)
yk(j+1)
and e1 :=
(2j+µ(j)+m(j)−1) mod 2j+1 X
k=2j+µ(j)
yk(j+1)
. Since for exact data one of the sums has only vanishing summands, see Figure 5.8, we choose
µ(j+1) :=
(µ(j) if e0 > e1, 2j+µ(j) if e0 < e1. This detection procedure has a computational effort of O m(j)
.
A4) y(J−1) has the one-block support S(J−1)=Iµ(J−1)(J−1),2J−1−1−µ(J−1) of lengthm(J−1) <
2J−1 centered around the boundary.
Theorem 5.19, caseA4 implies thaty(J) =y has the support
S(J)=Iµ(J(J)),2J−1−µ(J)∪Iη(J)(J),2J−1−η(J) with µ(J)<2J−1≤η(J). (5.17) This is either a two-block support with two separated blocks of possibly different lengths or, as a boundary case, a one-block support, where one of the two blocks in (5.17) vanishes. In the case of two blocks, one is centered around the mid-dle of the vector, i.e., around 2J−1 −1 and 2J−1, and its support is a subset of Iµ(J)(J−1),2J−1−µ(J−1), and the other one is centered around the boundary of the vec-tor, i.e., around2J−1and0, and its support is a subset ofI(J)
2J−1+µ(J−1),2J−1−1−µ(J−1). Ifyhas a one-block support, which can only happen if the first support indexµx of x is0 or its last support indexνx is2J−1−1, then one of the two blocks in (5.17) is empty. Since the support blocks have even length by Lemma 5.4 (i), set
T(0)(J):=n k∈n
µ(J−1), . . . ,2J−1−1o
:|yk|> εo
=:{t1, . . . , tK} and T(1)(J):=n
k∈n
2J−1+µ(J−1), . . . ,2J −1 o
:|yk|> ε o
=:{u1, . . . , uL},
see also Figure 5.9. IfT(0)(J)=∅, thenyhas a one-block support centered around the boundary. In this case we set
µ(J):=u1 and m(J):= 2· 2J−u1
= 2m
to obtain the support interval S(J):= Iµ(J)(J),2J−1−µ(J) of y. IfT(1)(J) =∅, then y has a one-block support centered around the middle of the vector, and we set
µ(J):=t1 and m(J):= 2· 2J−1−t1
= 2m, implying that the support interval of y is S(J) :=I(J)
µ(J), µ(J)+m(J)−1. If neitherT(0)(J) norT(1)(J) is empty,yhas a two-block support with two separated blocks of possibly
5.6 Sparse Fast IDFT and Sparse Fast IDCT-II
different lengths. Recall that we denote the first index of the block centered around the middle by µ(J) and the first index of the block centered around the boundary by η(J) and set
µ(J):=t1 and η(J):=u1.
We obtain the support setS(J) :=Iµ(J)(J),2J−1−µ(J)∪Iη(J)(J),2J−1−η(J), where the block centered around the middle of the vector has length n(J)(0) := 2 2J−1−t1
and the block centered around the boundary has length n(J)(1) := 2 2J−u1
. Thus, the support detection requires us to inspect 2·m(J−1)/2 entries of y(J), which has an arithmetical complexity of O m(J−1)
. B) Two-block support:
y(j)has the two-block supportS(j)=Iµ(j)(j), ν(j)∪I(j)
2j−1−ν(j),2j−1−µ(j) with block length n(j)=m.
We know from the proof of Theorem 5.19, case B, the proof of Theorem 5.23 and Figure 5.10 thaty(j+1) has a two-block support of block length n(j+1):=n(j)=m.
Further, the first support indexµ(j+1) of y(j+1) is eitherµ(j) or 2j+µ(j) with µ(j+1):=
µ(j) if
cu02k0+1−y\(j+1)2k0+1 <
cu02k0+1+y\(j+1)2k0+1 , 2j+µ(j) else,
where u0 and u1 denote the two possibilities for y(j+1). In this case the support detection has a computational effort ofO(1).
The support detection methods described in this section cover all possible cases for the support of y(j+1) for a given y(j). Additionally, as we will prove in Section 5.6.3, they are stable and efficient, as they require at most O m(j)
operations.
5.6.2 Sparse Fast IDFT for Vectors with Reflected Block Support