Читать книгу The Sparse Fourier Transform - Haitham Hassanieh - Страница 11

2

Preliminaries

This chapter will introduce the notation and basic definitions that we will use throughout Part I of this book.

2.1 Notation

We use ω = e^−2πi/n as the n-th root of unity and as the -th root of unity. For any complex number a, we use ϕ(a) ∊ [0, 2π]to denote the phase of a. For a complex number a and a real positive number b, the expression a ± b denotes a complex number a′ such that |a − a′| ≤ b.

For a vector x ∊ ℂⁿ, its support is denoted by supp(x) ⊂ [n]. We use ||x||₀ to denote |supp(x)|, the number of non-zero coordinates of x. Its Fourier spectrum is denoted by x̂, with

For a vector of length n, indices should be interpreted modulo n, so x_−i = xn_−i. This allows us to define convolution

and the coordinate-wise product (x · y)_i = xiyi, so .

We use [n] to denote the set {1,…, n}. All operations on indices in are taken modulo n. Therefore we might refer to an n-dimensional vector as having coordinates {0, 1,…, n − 1} or {0, 1,…, n/2, − n/2 + 1, …, −1} interchangeably. When i ∈ ℤ is an index into an n-dimensional vector, sometimes we use |i| to denote min_{j≡i (mod n)} |j|. Finally, we assume that n is an integer power of 2.

For the case of 2D Fourier transforms which will appear in Chapter 5, we assume that is a power of 2. We use [m] × [m] = [m]² to denote the m × m grid {(i, j) : i ∊ [m], j ∊ [m]}. For a 2D matrix , its support is denoted by .We also use ||x||₀ to denote |supp(x)|. Its 2D Fourier spectrum is denoted by x̂, with

Finally, if y is in frequency-domain, its inverse is denoted by y̌.

2.2 Basics

2.2.1 Window Functions

In digital signal processing [Oppenheim et al. 1999] one defines window functions in the following manner.

Definition 2.1 We define a (∊, δ, w) standard window function to be a symmetric vector F ∈ ℝⁿ with supp(F) ⊆ [−w/2, w/2] such that for all i ∊ [−∊n, ∊n], and for all i ∉ [−∊n, ∊n].

Claim 2.1 For any ∊ and δ, there exists an standard window function.

Proof This is a well-known fact [Smith 2011]. For example, for any ∊ and δ, one can obtain a standard window by taking a Gaussian with standard deviation and truncating it at . The Dolph-Chebyshev window function also has the claimed property but with minimal big-Oh constant [Smith 2011] (in particular, half the constant of the truncated Gaussian). ■

The above definition shows that a standard window function acts like a filter, allowing us to focus on a subset of the Fourier coefficients. Ideally, however, we would like the pass region of our filter to be as flat as possible.

Definition 2.2 We define a (∊, ∊′, δ, ω) flat window function to be a symmetric vector F ∈ ℝⁿ with supp(F) ⊆ [−w/2, w/2] such that for all i ∈ [−∊′n, ∊′n] and for all i ∉ [−∊n, ∊n].

A flat window function (like the one in Figure 2.1) can be obtained from a standard window function by convolving it with a “box car” window function, i.e., an interval. Specifically, we have the following.

Claim 2.2 For any ∊, ∊′, and δ with ∊′ < ∊, there exists an flat window function.

Note that in our applications we have δ < 1/n^O((1) and ∊ = 2∊′. Thus the window lengths w of the flat window function and the standard window function are the same up to a constant factor.

Figure 2.1 An example flat window function for n = 256. This is the sum of 31 adjacent (1/22, 10⁻⁸, 133) Dolph-Chebyshev window functions, giving a (0.11, 0.06, 2 × 10⁻⁹, 133) flat window function (although our proof only guarantees a tolerance δ = 29 × 10⁻⁸, the actual tolerance is better). The top row shows G and in a linear scale, and the bottom row shows the same with a log scale.

Proof Let f = (∊ − ∊′)/2, and let F be an standard window function with minimal . We can assume ∊, ∊′ > 1/(2n) (because [−∊n, ∊n] = {0} otherwise), so . Let be the sum of 1 + 2(∊′ + f)n adjacent copies of , normalized to . That is, we define

so by the shift theorem, in the time domain

Since for |i| ≤ fn, the normalization factor is at least 1. For each i ∈ [−∊′n, ∊′n], the sum on top contains all the terms from the sum on bottom. The other 2∊′n terms in the top sum have magnitude at most δ/((∊′ + ∊)n) = δ/(2(∊′ + f)n), so . For |i| > ∊n, however, . Thus F′ is an (∊, ∊′, δ, w) flat window function, with the correct w. ■

2.2.2 Permutation of Spectra

Following Gilbert et al. [2005a], we can permute the Fourier spectrum as follows by permuting the time domain:

Definition 2.3 Suppose σ⁻¹ exists mod n. We define the permutation P_{σ, a, b} by

We also define πσ, _b(i) = σ(i − b) mod n.

Claim 2.3 .

Proof

Lemma 2.1 If j ≠ 0, n is a power of two, and σ is a uniformly random odd number in [n], then Pr[σj ∊ [−C, C]] ≤ 4C/n.

Proof If j = m2^l for some odd m, then the distribution of σj as σ varies is uniform over m′2^l for all odd m′. There are thus 2 · round(C/2^l+1) < 4C/2^l+1 possible values in [−C, C] out of n/2^l+1 such elements in the orbit, for a chance of at most 4C/n. ■

Note that for simplicity we will only analyze our algorithm when n is a power of two. For general n, the analog of Lemma 2.1 would lose an n/φ(n) = O(log log n) factor, where φ is Euler’s totient function. This will correspondingly increase the running time of the algorithm on general n.

Claim 2.3 allows us to change the set of coefficients binned to a bucket by changing the permutation; Lemma 2.1 bounds the probability of non-zero coefficients falling into the same bucket.

2.2.3 Subsampled FFT

Suppose we have a vector x ∈ ℂⁿ and a parameter B dividing n, and would like to compute ŷi = x̂i_(n/B) for i ∈ [B].

Claim 2.4 ŷ is the B-dimensional Fourier transform of . Therefore, ŷ can be computed in O(|supp(x)| + B log B) time.

Proof

2.2.4 2D Aliasing Filter

The aliasing filter presented in Section 1.1.2. The filter generalizes to two dimensions as follows:

Given a 2D matrix , and Br, Bc that divide , then for all (i, j) ∈ [Br] × [Bc] set

Then, compute the 2D DFT ŷ of y. Observe that ŷ is a folded version of x̂: