CHAPTER 6 A whirlwind tour of Discrete Fourier and Haar Analysis

In this chapter we capitalize on the knowledge acquired while studying Fourier series to develop the simpler Fourier theory in a finite dimensional setting.

In this context, the discrete trigonometric basis in ℂ^{N} is a collection of N orthonormal vectors, hence an orthonormal basis. The discrete Fourier transform is the linear transformation that gives us the coefficients in the trigonometric basis. Many of the features of the Fourier series theory are still present, but without the nuisances and challenges of being in an infinite dimensional setting (L^{2}(𝕋)), where one has to deal with integration and infinite sums. The tools required in the finite setting are the tools of linear algebra.

For practical purposes, this finite theory is what is needed. Computers deal with finite vectors. The Discrete Fourier Transform (DFT) is calculated by applying an invertible N ×N matrix to a given vector in ℂ^{N}, and the Discrete Inverse Fourier Transform by applying the inverse matrix to the transformed vector. The surprising fact is that one can perform these matrix multiplications in order N log_{2} N operations, as opposed to the expected N^{2} operations, using the celebrated Fast Fourier Transform (FFT) and algorithm which we describe in Section 6.5. This gain in the number of operations is of invaluable importance for numerical applications, especially when dealing with large problems.

In this chapter we introduce another orthonormal basis, the discrete Haar basis, in preparation for later chapters when we will introduce the Haar basis in L^{2}(ℝ) and in L^{2}[0, 1]). We highlight the similarities and the differences with the discrete trigonometric basis. There is a Fast Haar Transform (FHT), of order N operations, and we describe it in Section 6.8. The fast Haar transform is a precursor of the Fast Wavelet Transform (FWT) that we will study in more depth in later chapters.

6.1. Fourier Series: a summary

Let us begin by summarizing the ingredients of the theory of Fourier series, as we have developed it in Chapters 1–5.

(1) To each function f : 𝕋 → ℂ we associate a sequence {a_{n}} = {f̂(n)}_{n∈ℤ} of complex numbers. These numbers are known as the Fourier coefficients of f.

(2) The Fourier coefficients are given by the formula

a_{n} = f̂(n) := \frac{1}{2π} ∫_{−π}^{π} f(θ)e^{−inθ} dθ, for n ∈ ℤ.

(3) The building blocks of the Fourier series are the trigonometric functions e_{n} : 𝕋 → ℂ given by

e_{n}(θ) := e^{inθ}, where n ∈ ℤ and θ ∈ 𝕋.

(They are called trigonometric rather than exponential to emphasize the connection with cosines and sines given by Euler’s formula e^{inθ} = cos nθ+ i sin nθ.) Thus the n^{th} Fourier coefficient of f can be written in terms of the L^{2} inner product as

a_{n} = f̂(n) = ⟨f, e_{n}⟩.

(4) The Fourier series associated to f is the doubly infinite series formed by summing the products of the n^{th} Fourier coefficients a_{n} with the n^{th} trigonometric functions e^{inθ}:

f(θ) ∼ ∑_{n=−∞}^{∞} a_{n}e^{inθ}.

We use ∼ rather than = in this formula in order to emphasize that for a given point θ, the Fourier series may not sum to the value f(θ), and indeed may not converge at all.

(5) The set {e^{inθ}}_{n∈ℤ} of trigonometric functions forms an orthonormal basis for the vector space L^{2}(𝕋) of square-integrable functions on 𝕋. Here orthonormal means with respect to the L^{2}(𝕋) inner product. In particular this means that f equals its Fourier series in L^{2}(𝕋), or equivalently the partial Fourier sums converge in L^{2}(𝕋),

∥f − ∑_{|n|≤N} f̂(n)e^{inθ}∥_{L^{2}(𝕋)} → 0, N → ∞,

and the energy is preserved, that is Parseval’s identity holds,

∥f∥_{L^{2}(𝕋)}^{2} =∼ ∑_{n=−∞}^{∞} |f̂(n)|^{2}.

In the so-called continuous setting we have just described, the signal is a function, which takes as input the continuously varying quantity θ ∈ 𝕋. By contrast, in the discrete setting of finite Fourier analysis, the signal is a vector in ℂ^{N}, which can be thought of as a function taking as input the values in the discrete list 0, 1, …, N − 1. In the continuous setting the Fourier transform produces the doubly-infinite sequence of Fourier coefficients. In the discrete setting the discrete Fourier transform produces another vector in ℂ^{N}.

6.2. The discrete Fourier basis

We present the analogues in the discrete setting of the five ingredients of Fourier analysis listed above. Fix a positive integer N, and let ω := e^{2πi/N}. Recall that ℤ/(Nℤ) := {0, 1, …, N − 1}.

(1) To each vector v = (z_{0}, …, z_{N−1}) ∈ ℂ^{N} we associate a second vector v̂ = (a_{0}, …, a_{N−1}) ∈ ℂ^{N}. The vector v̂ is known as the discrete Fourier transform of v.

(2) The entries of the discrete Fourier transform are given by the formula

a_{n} = v̂(n) := \frac{1}{\sqrt{N}} ∑_{k=0}^{N−1} z_{k}ω^{−kn}, for n ∈ \{0, 1, …, N − 1\}.

(3) The building blocks of the discrete Fourier transform are the N discrete trigonometric functions e_{l} : ℤ/(Nℤ) → ℤ/(Nℤ) given by

e_{l}(k) := \frac{1}{\sqrt{N}} ω^{kl}, k ∈ \{0, 1, …, N − 1\},

for l ∈ {0, 1, …, N − 1}. Thus the n^{th} Fourier coefficient of v can be written in terms of the ℂ^{N} inner product as

a_{n} = v̂(n) = ⟨v, e_{n}⟩ := ∑_{k=0}^{N−1} v(k)\overline{e_{n}(k)}.

(4) The original vector v can be exactly determined from its discrete Fourier transform. The k^{th} entry of v is given by

v(k) := z_{k} = ∑_{k=0}^{N−1} a_{n}e_{n}(k).

(5) The set {e_{0}, …, e_{N−1}} forms an orthonormal basis for the vector space ℂ^{N}. Here orthonormal means with respect to the standard inner product or dot product in ℂ^{N}.

Exercise 6.1. Verify that the set {e_{0}, e_{1}, …, e_{N−1}} is an orthonormal set in ℂ^{N}.

We make four remarks. First, the discrete setting has a certain symmetry in that both the original signal and the transformed signal are of the same form: vectors in ℂ^{N}. In the continuous Fourier series setting, the original signal is a square-integrable function f : 𝕋 → ℂ, while the transformed signal is a doubly infinite sequence {a_{n}} ∈ l^{2}(ℤ). In Chapter 7 we will develop the Fourier theory in yet another setting, that of the Fourier transform. There the original signal and the transformed signal can again be of the same form. For instance, they can both be functions in the Schwartz class S(ℝ), or they can both be functions in L^{2}(ℝ). We will also discuss generalizations of the Fourier transform to much larger classes of input signals, such as functions in L^{p}(ℝ), and, more generally, tempered distributions.

Second, our signal vectors v = (z_{0}, …, z_{N−1}) ∈ ℂ^{N} in the discrete setting can be thought of as functions v : ℤ/(Nℤ) → ℂ, in the sense that for each number n ∈ ℤ/(Nℤ) := {0, 1, …, N − 1}, the function v picks out a complex number v(n) = z_{n}. So

v = (v(0), v(1), …, v(N − 1)).

This formulation of the input signal is closer to the way we presented the Fourier series setting. With this notation the discrete Fourier transform reads

(6.1)  v̂(m) = ∑_{n=0}^{N−1} v(n)e^{−2πinm/N}.

Third, in the discrete setting there are no difficulties with convergence of series, since all the sums are finite. Thus the subtleties discussed in Chapters 3, 4 and 5 are absent from the setting of the finite Fourier transform. We denote ℂ^{N} with the standard inner-product as the space ℓ^{2}(ℤ/Nℤ), to emphasize the similarities with the Hilbert space L^{2}(𝕋).

Fourth, in the discrete setting we can establish Parseval’s Identity, which is also known as Plancherel’s Formula:

∑_{n=0}^{N−1} |v(n)|^{2} = ∑_{n=0}^{N−1} |v̂(n)|^{2}.

This identity is exactly the Pythagorean Theorem for vectors with N entries. Rewritten in a form closer to the way we saw it in the continuous setting, Parseval’s Identity reads

∥v∥_{ℓ^{2}(ℤ/Nℤ)} = ∥v̂∥_{ℓ^{2}(ℤ/Nℤ)},

Here the norm is the Euclidean norm induced by the dot product; we also call it the ℓ^{2}-norm of a vector.

As in the continuous setting, Parseval’s Identity says that the energy of the transformed signal is the same as that of the original signal. In the continuous setting, the energy means the L^{2}-norm of a function or the ℓ^{2}-norm of a doubly infinite series. In the discrete setting, the energy means the ℓ^{2}-norm of a vector. Parseval’s Identity also holds in the Fourier transform setting. For example when the original and transformed signals, f and f̂, are both functions in L^{2}(ℝ), the energy is represented by the L^{2}-norm, and

∥f∥_{L^{2}(ℝ)} = ∥f̂∥_{L^{2}(ℝ)}.

Aside 6.2. Historically people have preferred to work with the orthogonal basis {f_{0}, f_{1}, …, f_{N−1}} in ℂ^{N} defined by

f_{l} = \sqrt{N} e_{l},

instead of the orthonormal basis {e_{l}}, the prize being the appearance of a factor 1/N in the definition of the coefficients and also in Parseval’s formula (but not in the reconstruction formula). To mirror the Fourier series theory, one could work with the vectors {f_{0}, f_{1}, …, f_{N−1}} and consider them to be normalized, by defining the following inner-product for vectors v, w ∈ ℂ^{N}:

(6.2)  ⟨v, w⟩_{nor} := \frac{1}{N} ∑_{n=0}^{N−1} v(n)\overline{w(n)}.

We made an analogous convention in L^{2}(𝕋) when we chose to include the factor 1/(2π) in the definition of the inner product. However, in our formulation of the discrete theory we prefer to use the standard inner product or dot product on ℂ^{N}, not the normalized one. Reader beware: it is wise to check which normalization a given author is using.

Exercise 6.3. Show that formula (6.2) defines an inner product in ℂ^{N}.

Exercise 6.4. Consider the vectors {f_{0}, f_{1}, …, f_{N−1}}. Verify that they are orthonormal with respect to the normalized inner product in ℂ^{N}.

Exercise 6.5. Define the normalized discrete Fourier coefficients of v ∈ ℂ^{N} to be

v^{nor}(n) = ⟨v, f_{n}⟩_{nor}.

6.3. A PARENTHESIS ON BASES AND DUAL BASES IN ℂ^{N} 97

Use the theory of orthonormal bases to justify the reconstruction formula

v = ∑_{n=0}^{N−1} v^{nor}(n)f_{n}.

Verify Parseval’s identity in this setting, namely

∑_{n=0}^{N−1} |v^{nor}(n)|^{2} = ∥v^{nor}∥_{ℓ^{2}(ℤ/Nℤ)}^{2} = ∥v∥_{nor}^{2} = \frac{1}{N} ∑_{n=0}^{N−1} |v(n)|^{2}.

We give two examples explicitly showing the Fourier orthonormal bases obtained using the scheme described above, for ℂ^{5} and ℂ^{8}.

Example 6.6. (Fourier Orthonormal Basis for ℂ^{5}) Let N = 5, so that

ω = e^{2πi/5}.

Let

f_{0} := [\begin{matrix}1 \\ 1 \\ 1 \\ 1 \\ 1 \end{matrix}] , f_{1} := [\begin{matrix}1 \\ ω \\ ω^{2} \\ ω^{3} \\ ω^{4} \end{matrix}] , f_{2} := [\begin{matrix}1 \\ ω^{2} \\ ω^{4} \\ ω \\ ω^{3} \end{matrix}] , f_{3} := [\begin{matrix}1 \\ ω^{3} \\ ω \\ ω^{4} \\ ω^{2} \end{matrix}] , f_{4} := [\begin{matrix}1 \\ ω^{4} \\ ω^{3} \\ ω^{2} \\ ω \end{matrix}] .

Notice the interesting arrangement of the powers of ω in each vector f_{l}.

The vectors f_{0}, f_{1}, f_{2}, f_{3}, f_{4} are pairwise orthogonal. Dividing by the length \sqrt{5} of each f_{l}, we obtain an orthonormal basis for ℂ^{5}, namely e_{o} = \{(1/\sqrt{5})f_{0}, e_{1} = (1/\sqrt{5})f_{1}, e_{2} = (1/\sqrt{5})f_{2}, e_{3} = (1/\sqrt{5})f_{3}, e_{4} = (1/\sqrt{5})f_{4}\}. ♢

Example 6.7. (Fourier Orthonormal Basis for ℂ^{8}) Let N = 8, so that

ω = e^{2πi/8} = e^{πi/4} = \frac{(1 + i)}{\sqrt{2}} .

We have pairwise orthogonal vectors f_{0}, f_{1}, …, f_{7} given by, respectively,

[\begin{matrix}1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \end{matrix}] , [\begin{matrix}1 \\ ω \\ ω^{2} \\ ω^{3} \\ ω^{4} \\ ω^{5} \\ ω^{6} \\ ω^{7} \end{matrix}] , [\begin{matrix}1 \\ ω^{2} \\ ω^{4} \\ ω^{6} \\ 1 \\ ω^{2} \\ ω^{4} \\ ω^{6} \end{matrix}] , [\begin{matrix}1 \\ ω^{3} \\ ω^{6} \\ ω \\ ω^{4} \\ ω^{7} \\ ω^{2} \\ ω^{5} \end{matrix}] , [\begin{matrix}1 \\ ω^{4} \\ 1 \\ ω^{4} \\ 1 \\ ω^{4} \\ 1 \\ ω^{4} \end{matrix}] , [\begin{matrix}1 \\ ω^{5} \\ ω^{2} \\ ω^{7} \\ ω^{4} \\ ω \\ ω^{6} \\ ω^{3} \end{matrix}] , [\begin{matrix}1 \\ ω^{6} \\ ω^{4} \\ ω^{2} \\ 1 \\ ω^{6} \\ ω^{4} \\ ω^{2} \end{matrix}] , [\begin{matrix}1 \\ ω^{7} \\ ω^{6} \\ ω^{5} \\ ω^{4} \\ ω^{3} \\ ω^{2} \\ ω \end{matrix}] .

Notice the patterns that arise because N is not prime.

Dividing by the length \sqrt{8} of each f_{l}, we obtain an orthonormal basis for ℂ^{8}, namely e_{0} = (1/\sqrt{8})f_{0}, e_{1} = (1/\sqrt{8})f_{1}, …, e_{7} = (1/\sqrt{8})f_{7}. ♢

6.3. A parenthesis on bases and dual bases in ℂ^{N}

Before discussing in more detail the discrete Fourier transform, let us review some fundamental linear algebra facts about bases.

Definition 6.8. A basis in ℂ^{N} is a set of N vectors {v_{1}, v_{2}, …, v_{N}} such that any vector v ∈ ℂ^{N} can be written as a linear combination of the elements of the basis, that is there exists complex numbers a_{1}, a_{2}, …, a_{N} and

(6.3)  v = a_{1}v_{1} + a_{2}v_{2} + ⋯ + a_{N}v_{N}.

The complex numbers a_{1}, a_{2}, …, a_{N} are uniquely determined by the vector v.

Denote by B the N × N matrix whose j^{th}-column is the vector v_{j}. That is,

B = [\begin{matrix}| & & | \\ v_{1} & … & v_{N} \\ | & & | \end{matrix}] .

With this notation, we have the following matrix representation of equation (6.3),

v = Ba, where a = [a_{1}, a_{2}, …, a_{N}]^{t}.

Theorem 6.9. The following are equivalent statements

(1) A set of N vectors {v_{1}, v_{2}, …, v_{N}} is a basis in ℂ^{N}

(2) The vectors {v_{1}, v_{2}, …, v_{N}} are linearly independent. That is, if 0 = a_{1}v_{1} + a_{2}v_{2} + ⋯ + a_{N}v_{N}, then a_{1} = ⋯ = a_{N} = 0.

(3) The N × N matrix B whose columns are the vectors v_{j} is invertible, or non-singular. Therefore we can find the coefficients applying the inverse matrix B^{−1} to v,

a = B^{−1}v.

Exercise 6.10. Prove the equivalences in Theorem 6.9.

Given N linearly independent vectors {v_{1}, v_{2}, …, v_{N}} in ℂ^{N}, that is a basis, let B be the invertible matrix with columns v_{j}, and denote by w_{j} the complex conjugate of the j^{th}-row vector of its inverse B^{−1}, namely,

B^{−1} = [− \overline{w_{1}} − ; ⋮ ; − \overline{w_{N}} −] .

With this notation, performing the matrix multiplication, we see that,

B^{−1}v = [− \overline{w_{1}} − ; ⋮ ; − \overline{w_{N}} −] [ v(1) ; ⋮ ; v(N)] = [ ⟨v, w_{1}⟩ ; ⋮ ; ⟨v, w_{N}⟩] .

Exercise 6.11. Verify that the j^{th} entry of the vector B^{−1}v is the inner product between the vector v and the vector w_{j}.

The coefficients of v in the basis {v_{1}, v_{2}, …, v_{N}} are given exactly by B^{−1}v, that means that not only they are uniquely determined, but we have an algorithm to find them,

a_{j} = ⟨v, w_{j}⟩,

where the conjugates of the vectors w_{j} are the rows of B^{−1}. That is there exists a set of dual vectors (or dual basis) {w_{1}, w_{2}, …, w_{N}} to the basis {v_{1}, v_{2}, …, v_{N}}, so that to compute the coefficients of a given vector v in the basis all we need to do is compute the inner products with the dual vectors.

Exercise 6.12. Verify that given a basis {v_{1}, v_{2}, …, v_{N}} and its dual vectors {w_{1}, w_{2}, …, w_{N}} satisfy the following “orthonormality” condition,

(6.4)  ⟨v_{k}, w_{j}⟩ = δ_{k,j}.

Furthermore, given a basis {v_{1}, v_{2}, …, v_{N}}, if you can locate vectors {w_{1}, w_{2}, …, w_{N}} with the “orthonormality” property (6.4), then the matrix whose rows are the conjugates of the vectors {w_{1}, w_{2}, …, w_{N}} must be the inverse matrix to the matrix whose columns are the vectors in the basis. That is the vectors {w_{1}, w_{2}, …, w_{N}} are the dual basis of {v_{1}, v_{2}, …, v_{N}}.

We will encounter dual basis in an infinite dimensional context when we discuss biorthogonal wavelets.

In particular if the basis {v_{1}, v_{2}, …, v_{N}} is orthonormal, the dual vectors are the vectors themselves. On the matrix side, this is equivalent to observing that

B^{−1} = \overline{B^{t}}.

Definition 6.13. An N × N matrix whose columns are orthonormal vectors in ℂ^{N} is a unitary matrix.

Unitary matrices are always invertible and their inverses are calculated by taking the conjugate of its transpose.

6.4. The discrete Fourier transform and its inverse

It is not hard to see that the transformation that maps the vector v ∈ ℂ^{N} into the vector v̂ ∈ ℂ^{N} is a linear transformation. We call this transformation the discrete Fourier transform or the Fourier transform in ℓ_{2}(ℤ/(Nℤ)). The notation means that we are viewing ℂ^{N} as an inner-product vector space whose elements are functions defined on the discrete set ℤ/(Nℤ) = {0, 1, …, N −1}, with the standard Euclidean inner-product.

The finite Fourier transform is completely determined by the N × N complexvalued symmetric matrix \frac{1}{\sqrt{N}} \overline{F_{N}}, where F_{N} is the matrix whose columns are exactly given by the orthogonal Fourier vectors f_{l}. We call F_{N} the Fourier matrix. In other words, the following matrix formula holds for the discrete Fourier transform,

(6.5)  v̂ = \frac{1}{\sqrt{N}} \overline{F}_{N}v.

Example 6.14. (Fourier Matrix for ℂ^{8}) We display the Fourier matrix in the case when N = 8,

F_{8} = [\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & ω & ω^{2} & ω^{3} & ω^{4} & ω^{5} & ω^{6} & ω^{7} \\ 1 & ω^{2} & ω^{4} & ω^{6} & 1 & ω^{2} & ω^{4} & ω^{6} \\ 1 & ω^{3} & ω^{6} & ω & ω^{4} & ω^{7} & ω^{2} & ω^{5} \\ 1 & ω^{4} & 1 & ω^{4} & 1 & ω^{4} & 1 & ω^{4} \\ 1 & ω^{5} & ω^{2} & ω^{7} & ω^{4} & ω & ω^{6} & ω^{3} \\ 1 & ω^{6} & ω^{4} & ω^{2} & 1 & ω^{6} & ω^{4} & ω^{2} \\ 1 & ω^{7} & ω^{6} & ω^{5} & ω^{4} & ω^{3} & ω^{2} & ω^{1} \end{matrix}] ,

Notice that the matrix \frac{1}{\sqrt{N}} F_{N} is a symmetric unitary matrix, because its columns are orthonormal vectors. In particular its inverse is its complex conjugate \frac{1}{\sqrt{N}} \overline{F_{N}}. We conclude that

(\frac{1}{\sqrt{N}} \overline{F_{N}})^{−1} = \frac{1}{\sqrt{N}} F_{N}.

Therefore, the discrete inverse Fourier transform in ℓ_{2}(ℤ_{N}/(Nℤ)) becomes,

(6.6)  v = \frac{1}{\sqrt{N}} F_{N} v̂.

In this finite dimensional context, all we are doing is an orthogonal change of basis, from the standard basis to a representation in the Fourier basis.

Exercise 6.15. Verify that

\frac{1}{N} F_{N} \overline{F_{N}} = \frac{1}{N} \overline{F_{N}} F_{N} = I_{N},

where I_{N} denotes the N × N identity matrix.

Both F_{N} and \overline{F_{N}} are full matrices, a priori, computing the discrete Fourier transform or inverse Fourier transform involves an operation count of the order of N^{2} operations. However there is a lot of hidden structure in these matrices, and one can dramatically improve on this operation count. The celebrated Fast Fourier Transform (FFT) is the algorithm that allows to perform the finite Fourier transform and the inverse Fourier transform in N log_{2} N operations.

6.5. The Fast Fourier Transform (FFT)

A key to the practical success of the discrete Fourier Transform (6.5) is a fast algorithm to accomplish it. The N × N matrix F_{N} is full, with entries F_{N}(n, m) = e^{−2πinm/N}. Applying it therefore involves N^{2} complex multiplications. In the 1960’s the American applied mathematician James Cooley and the American statistician John Wilder Tukey [1915-2000] rediscovered a faster algorithm, which the German mathematician Johan Carl Friedrich Gauss [1777-1855] had found in the early 1800s. This famous Fast Fourier Transform (FFT) algorithm reduces the number of multiplications to order N log_{2} N. This improvement revolutionized digital signal processing (and of course at Gauss’s time it was completely unforeseen!).

Mathematically, the FFT is based on a factorization of the Fourier matrix into a collection of sparse matrices. The process is easiest to explain, and implement, when N = 2^{n}. It begins with a factorization that reduces F_{N} to two copies of F_{N/2}. This reduction continues for n = log_{2} N steps, until the original matrix is written as a product of 2 log_{2} N sparse matrices, half of which can be collapsed to a single permutation matrix. The total number of multiplications required to apply F_{N} becomes N log_{2} N instead of the brute force N^{2}. To understand why this works, we rearrange (6.1) to express v̂(m) as a linear combination of two terms A_{m} and B_{m}, by

v̂(m) = \frac{1}{\sqrt{N}} ∑_{n=0}^{N−1} v(n)e^{\frac{−2πimn}{N}}

= \frac{1}{\sqrt{N}} [ ∑_{n=0}^{N/2−1} v(2n)e^{\frac{−2πimn}{N/2}} + e^{\frac{−2πim}{N}} ∑_{n=0}^{N/2−1} v(2n + 1)e^{\frac{−2πimn}{N/2}} ]

= A_{m} + e^{\frac{−πim}{N/2}} B_{m} .

The key is to note that v̂(m + \frac{N}{2} ) can be written as another linear combination of the same two terms, by

v̂(m + \frac{N}{2} ) = \frac{1}{\sqrt{N}} ∑_{n=0}^{N−1} v(n)e^{\frac{−2πi(m+N/2)n}{N}}

= \frac{1}{\sqrt{N}} [ ∑_{n=0}^{\frac{N}{2} −1} v(2n)e^{\frac{−2πimn}{N/2}} + e^{\frac{−2πi(m+N/2)}{N}} ∑_{n=0}^{\frac{N}{2} −1} v(2n + 1)e^{\frac{−2πimn}{N/2}} ]

= A_{m} − e^{\frac{−πim}{N/2}} B_{m} .

By clever ordering of operations we can reuse A_{m} and B_{m} and gain on the total number of operations. Furthermore, the sums A_{m} and B_{m} have the same structure as the sum in v̂(m), so we can do this process recursively.

In the language of matrices, what we are doing in the above calculations, is factorizing the conjugate Fourier matrix \overline{F_{N}}, for N even, as a product of three N × N matrices,

(6.7)  \overline{F_{N}} = [\begin{matrix}I_{N/2} & D_{N/2} \\ I_{N/2} & −D_{N/2} \end{matrix}] [\begin{matrix}\overline{F_{N/2}} & 0_{N/2} \\ 0_{N/2} & \overline{F_{N/2}} \end{matrix}] [\begin{matrix}Even & N/2 \\ Odd & N/2 \end{matrix}] ,

here I_{M} means the M × M identity matrix, F_{M} is the M × M Fourier matrix, 0_{M} denotes the M × M zero matrix. We will now describe precisely the diagonal M × M matrix D_{M}, and the M × 2M matrices Even_{M} and Odd_{M}.

The diagonal matrix D_{M} has diagonal entries e^{−πim/M} for m = 0, 1, …, M −1. Here is D_{4},

D_{4} = [\begin{matrix}1 & 0 & 0 & 0 \\ 0 & e^{\frac{−2πi}{4}} & 0 & 0 \\ 0 & 0 & e^{\frac{−2π2i}{4}} & 0 \\ 0 & 0 & 0 & e^{\frac{−2π3i}{4}} \end{matrix} ] = [\begin{matrix}1 & 0 & 0 & 0 \\ 0 & −i & 0 & 0 \\ 0 & 0 & −1 & 0 \\ 0 & 0 & 0 & i \end{matrix}]

The M × 2M matrices Even_{M} and Odd_{M} grab in order the even entries and the odd entries of a given 2M-vector respectively. Here are Even_{4} and Odd_{4} acting on a generic vector in ℂ^{8}, v = [v(0), v(1), …, v(7)],

[\begin{matrix}Even & 4 \\ Odd & 4 \end{matrix}] v = [\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}] [\begin{matrix}v(0) \\ v(1) \\ v(2) \\ v(3) \\ v(4) \\ v(5) \\ v(6) \\ v(7) \end{matrix}] = [\begin{matrix}v(0) \\ v(2) \\ v(4) \\ v(6) \\ v(1) \\ v(3) \\ v(5) \\ v(7) \end{matrix}] .

Notice that the first matrix in the decomposition (6.7) requires only N multiplications. The second requires 2(N/2)^{2} = N^{2}/2 multiplications, but it has the Fourier matrix F_{N/2} in the diagonal, and we can use the recursive equation (6.7), to replace the (N/2) × (N/2) matrix \overline{F_{N/2}} by

\overline{F_{N/2}} = [\begin{matrix}I_{N/4} & D_{N/4} \\ I_{N/4} & −D_{N/4} \end{matrix}] [\begin{matrix}\overline{F_{N/4}} & 0_{N/4} \\ 0_{N/4} & \overline{F_{N/4}} \end{matrix}] [\begin{matrix}Even & N/4 \\ Odd & N/4 \end{matrix}] ,

and keep replacing until one cannot subdivide further.

Exercise 6.16. Verify that after two steps in the recursion the middle matrix in (6.7) is decomposed as the product of three N × N matrices as follows,

[\begin{matrix}I_{N/4} & D_{N/4} \\ I_{N/4} & −D_{N/4} \\ & & I_{N/4} & D_{N/4} \\ & & I_{N/4} & −D_{N/4} \end{matrix}] [\begin{matrix}F_{N/4} \\ & F_{N/4} \\ & & F_{N/4} \\ & & & F_{N/4} \end{matrix}] [\begin{matrix}Even & N/4 \\ Odd & N/4 \\ Even & N/4 \\ Odd & N/4 \end{matrix}] .

After two iterations, the first two matrices in the decomposition of F_{N}, require only N multiplications each. If N = 2^{j}, then one can do this j − 1 = log_{2} N times (until the block size is 2 × 2). After j − 1 steps, the matrix \overline{F_{N}} will be decomposed into the product of log_{2}N − 1 matrices each one requiring only N multiplications for a total of N(log_{2} N −1) multiplications, followed by the product of j − 1 “scrambling” matrices (these are all the Odd/Even matrices). A miracle then occurs, the j −1-scrambling matrices, when applied in the right order, become a very simple operation which basically costs nothing in terms of operation counts. The miracle is observed when describing the numbers 0, 1, …, N − 1 in binary notation: “the succesive application of the Odd/Even matrices will simply put in place n the m^{th} entry, where m is found by reversing the digits in the binary decomposition of n”. This is called in computer sciences bit reversal.

Let us see the miracle in action with an example when N = 8 = 2^{3}. Denote the succesive scrambling matrices by

S_{4} = [\begin{matrix}Even & 4 \\ Odd & 4 \end{matrix}] , S_{2} = [\begin{matrix}Even & 2 \\ Odd & 2 \\ Even & 2 \\ Odd & 2 \end{matrix}] .

We start with the numbers 0, 1, …, 7 in binary notation and apply the matrix S_{4}, to the outcome we apply the matrix S_{2},

[\begin{matrix}000 \\ 001 \\ 010 \\ 011 \\ 100 \\ 101 \\ 110 \\ 111 \end{matrix}] →^{S_{4}} [\begin{matrix}000 \\ 010 \\ 100 \\ 110 \\ 001 \\ 011 \\ 101 \\ 111 \end{matrix}] →^{S_{2}} [\begin{matrix}000 \\ 100 \\ 010 \\ 110 \\ 001 \\ 101 \\ 011 \\ 111 \end{matrix}] .

Notice that the binary numbers in the right-hand-side column are obtained from the binary numbers in the left-hand-side column by reading backwards the binary digits.

Exercise 6.17. Can you justify the miracle for any N = 2^{j}? Try an argument by induction.

Even without a miracle, the scrambling N × N matrices S_{N/2^{k}} are very sparse, only N entries are non-zero, in fact their non-zero entries are all ones. Applying each of the scrambling matrices requires only N operations, and applying j −1 such matrices (N = 2^{j}), requires at most N(log_{2} N − 1) operations. This means that the operation count for the Fourier matrix based in this decomposition is no more than 2N(log_{2} N − 1) operations. This is still an algorithm of order N log_{2} N.

Exercise 6.18. Given vectors v, w ∈ ℂ^{N}, define their discrete convolution by

v ∗ w(n) := \frac{1}{\sqrt{N}} ∑_{k=0}^{N−1} v(k)w(n − k), n = 0, 1, …, N − 1.

Notice that if n < k then −N ≤ n − k < 0, in that case define w(n − k) := w(N − (n − k)) (or think as w as extended periodically over ℤ so that w(n) = w(N + n) for all n ∈ ℤ). Notice that in matrix notation, the linear transformation T_{w} that maps the vector v into the vector v ∗ w is given by the circulant matrix

[\begin{matrix}w(0) & w(N − 1) & w(N − 2) & … & w(1) \\ w(1) & w(0) & w(N − 1) & … & w(2) \\ w(2) & w(1) & w(0) & … & w(3) \\ ⋮ & ⋮ & ⋮ & & ⋮ \\ w(N − 1) & w(N − 2) & w(N − 3) & … & w(0) \end{matrix}] ,

This is a priori a full matrix and requires N^{2} operations to apply it.

(a) Show that

{ \atop v ∗ w}(m) = v̂(m) ŵ(m).

(b) Describe a fast algorithm to compute convolution in order N log_{2} N operations.

Convolution with a fixed vector w is a linear transformation given by a circulant matrix as described in Exercise 6.18. These are also the only translation invariant or shift invariant linear transformations. This means that if we first shift a vector and then convolve the output is the same as first convolving and then shifting. More precisely, denote by S_{k} the linear transformation that shifts a vector v ∈ ℂ^{N} k units, that is S_{k}v(n) = v(n − k) (we are viewing the vector as extended periodically so that v(n − k) is defined for all integers), then S_{k}T_{w} = T_{w}S_{k}.

Exercise 6.19. Verify that if T_{w}v = v ∗ w then S_{k}T_{w} = T_{w}S_{k}. Furthermore verify that if a linear transformation T is shift invariant, that is S_{k}T = T S_{k} then there exists a vector w ∈ ℂ^{N} such that T v = v ∗ w.

The book by Frazier [Fra] takes a linear algebra approach and discusses these an much more in the first three chapters, before moving on to infinite dimentional space. Chapter 1 is a thorough review of linear algebra with emphasis on changes of basis and translation invariant linear transformations. Chapter 2 discusses the discrete Fourier theory. Chapter 3 discusses the discrete wavelet theory.

Another book that has strong linear algebra flavor is the book by Strang and Nguyen [SN].

6.6. The discrete Haar basis

In this section we discuss the discrete Haar basis. Unlike the discrete Fourier basis, the Haar basis is more “localized”, a concept that we now describe.

In this context by localized we mean, in a not very precise way, the vector is zero except for a few entries. The discrete trigonometric vectors have all entries non-zero, in fact each entry has absolute value 1/\sqrt{N}. Certainly the discrete trigonometric vectors are not localized. In particular when computing Fourier coefficients, the inner product of a vector v against a trigonometric vector e_{n}, all entries of v are involved in the computation.

An example of a basis which is as localized as it can be, is the standard basis in ℂ^{N}. These vectors have all but one entry equal to zero, the non-zero entry is 1. Denote by s_{j} the vector in the standard basis whose j^{th}-entry is 1,

s_{1} = [1, 0, …, 0], s_{2} = [0, 1, 0, …, 0], …, s_{N} = [0, …, 0, 1].

Exercise 6.20. Verify that the discrete Fourier transform of e_{j} is s_{j}, that is ê_{j} = s_{j}.

Although the Fourier basis is not localized at all, its Fourier transform is as localized as it can be. We say the Fourier basis is localized in frequency, but not in space or time.

Exercise 6.21. Verify that the discrete Fourier transform of s_{j} is the complex conjugate of e_{j}, that is ŝ_{j} = \overline{e_{j}}.

Although the standard basis is as localized as it can be, its Fourier transform is not localized at all. We say the standard basis is localized in space or time but not in frequency.

Exercise 6.22. Let {v_{1}, v_{2}, …, v_{N}} be an orthonormal basis in ℂ^{N}. Show that their discrete Fourier transforms {v̂_{1}, v̂_{2}, …, v̂_{N}} is an orthonormal basis of ℂ^{N}.

These are incarnations of a general principle, the Uncertainty Principle, that we will discuss in more depth in subsequent chapters. The Uncertainty Principle says that it is imposible for a vector to be simultaneously localized in space and frequency.

The discrete Haar basis will be intermediate between the discrete Fourier basis and the standard basis, in terms of its space and frequency localization.

Let us first see an example.

Example 6.23. Let N = 2^{3} = 8. Define the following vectors h̃_{n} for n = 0, 1, …, 7,

h̃_{0} := [1, 1, 1, 1, 1, 1, 1, 1],
h̃_{1} := [−1, −1, −1, −1, 1, 1, 1, 1],
h̃_{2} := [−1, −1, 1, 1, 0, 0, 0, 0],
h̃_{3} := [0, 0, 0, 0, −1, −1, 1, 1],
h̃_{4} := [−1, 1, 0, 0, 0, 0, 0, 0],
h̃_{5} := [0, 0, −1, 1, 0, 0, 0, 0],
h̃_{6} := [0, 0, 0, 0, −1, 1, 0, 0],
h̃_{7} := [0, 0, 0, 0, 0, 0, −1, 1].

Exercise 6.24. Verify that the vectors {h̃_{0}, h̃_{1}, h̃_{2}, h̃_{3}, h̃_{4}, h̃_{5}, h̃_{6}, h̃_{7}} are orthogonal vectors in ℂ^{8}. Verify that the vectors

(6.8)  h_{0} = \frac{h̃_{0}}{\sqrt{8}}, h_{1} = \frac{h̃_{1}}{\sqrt{8}}, h_{2} = \frac{h̃_{2}}{\sqrt{4}}, h_{3} = \frac{h̃_{3}}{\sqrt{4}},

h_{4} = \frac{h̃_{4}}{\sqrt{2}}, h_{5} = \frac{h̃_{5}}{\sqrt{2}}, h_{6} = \frac{h̃_{6}}{\sqrt{2}}, h_{7} = \frac{h̃_{7}}{\sqrt{2}},

are orthonormal. They are the discrete Haar basis in ℂ^{8}.

Exercise 6.25. Compute the discrete Fourier transform of each of the vectors in the Haar basis in ℂ^{8}, displayed in (6.8). Either analytically or using Matlab, in the later case, show pictures of h_{n} and ĥ_{n} and notice that the more localized h_{n} is the less localized ĥ_{n} would be.

When N = 2^{j} then it should be clear how to define the N Haar vectors. The first not normalized vector (h̃_{0}) is a string of ones; the second one (h̃_{1}) has the first half of the entries equal to −1, the second half equal to 1; the third one (h̃_{2}) has the first quarter of the entries equal to −1, the second quarter 1, and the rest of the entries are zero; etc. The last not-normalized vector (h̃_{N−1}) consists of N − 2 zeros, and the last entries are −1 and 1 respectively.

Notice that half of the vectors are very localized, having exactly two consecutive entries non-zero, these are the 2^{j−1} vectors {h̃_{2^{j−1}}, h̃_{2^{j−1}+1}, …, h̃_{2^{j}−1}}. The 2^{j−2} vectors {h̃_{2^{j−2}}, h̃_{2^{j−2}+1}, …, h̃_{2^{j−1}−1}} have all entries except four consecutive entries non-zero. The 2^{j−k} vectors {h̃_{2^{j−k}}, h̃_{2^{j−k}+1}, …, h̃_{2^{j−k+1}−1}} have all entries except 2^{k} consecutive entries non-zero. The vectors {h̃_{2}, h̃_{3}} have half of the entries zero. Finally h̃_{0} and h̃_{1} have all entries non-zero. Furthermore the non-zero entries consists of a string of consecutive −1’s followed by the same number of 1’s, except for h̃_{0}.

We are now going to give a precise definition of the N^{th} Haar basis. We should really tag the N vectors with an N to specify at what dimension we are. However, to make the notation less cumbersome, we will not tag them, and we expect the reader to be aware of the underlying dimension.

Given 0 ≤ n < N = 2^{j} there exist a unique 0 ≤ k < j such that

2^{k} ≤ n < 2^{k+1}, moreover n = 2^{k} + m, 0 ≤ m < 2^{k}.

We now define h̃_{n} in terms of those unique numbers k and m determined by n.

h̃_{n}(l) = \{ \begin{matrix}−1 & if mN2^{−k} ≤ l < (2m + 1)N2^{−k−1}, \\ 1 & if (2m + 1)N2^{−k−1} ≤ l < (m + 1)N2^{−k}, \\ 0 & otherwise \end{matrix}

Exercise 6.26. Verify that using this formula when N = 8 we get the same vectors as in Example 6.23

Notice that h̃_{n}(l) is non-zero if and only if m2^{j−k} ≤ l < (m + 1)2^{j−k}, in fact those entries 2^{j−k} = N/2^{k} are ±1. The Haar vectors are defined to be,

h_{n} = \frac{2^{k/2}}{\sqrt{N}} h̃_{n}.

Exercise 6.27. Verify that the vectors {h_{n}}_{n=0}^{N−1} are orthonormal.

Notice that except for h_{0}, all the other vectors are rescaled contractions of h_{1} and translates of those.

The parameters k, m, uniquely determined by n so that n = 2^{k} +m, are scaling and translation parameters respectively. The scaling parameter k, 0 ≤ k < j, tells us that we will decompose the set {0, 1, 2, …, 2^{j−1}} into 2^{k} disjoint subsets of 2^{j−k} = N/2^{k} consecutive numbers. The translation parameter m, 0 ≤ m < 2^{j−k}, indicates which subinterval we are considering.

It would make sense to index the n^{th} Haar vector with the pair (k, m), and this is exactly what we will do when we study the Haar basis, and the wavelet bases on L^{2}(ℝ). These are orthonormal bases, {ψ_{k,m}}_{k,m∈ℤ}, indexed by two integer parameters (k, m) that are found by dilating and translating one function, ψ, the wavelet. More precisely,

ψ_{k,m}(x) = 2^{k/2}ψ(2^{k}x − k).

The functions ψ_{k,m} retain the shape of the wavelet ψ. If we assume that ψ is concentrated around zero on the interval [−1, 1], then the functions ψ_{k,m} is now localized around the point m2^{−k}, and concentrated on the interval [(m−1)2^{−k}, (m+ 1)2^{−k}] of length 2×2^{−k}. As k → ∞ the resolution increases (that is, the localization improves), as k → −∞ the localization gets worst and worst.

6.7. The Discrete Haar Transform

The discrete Haar transform is the linear transformation that applied to a vector v ∈ ℂ^{N}, N = 2^{j}, gives its coefficients in the Haar basis. The inverse discrete Haar transform is the linear transformation that reconstructs the vector v from its Haar coefficients {c_{n}}_{n=1}^{N−1}. Notice that

v = c_{0}h_{0} + c_{1}h_{1} + ⋯ + c_{N−1}h_{N−1}.

In matrix notation the following is the inverse discrete Haar transform,

v = H_{N} c, where c = [c_{0}, c_{1}, …, c_{N−1}]^{t},

and the Haar matrix H_{N} is the N × N matrix whose j^{th}-column is the Haar vector h_{j−1}. That is,

H_{N} = [\begin{matrix}| & & | \\ h_{0} & … & h_{N} \\ | & & | \end{matrix}] .

The Haar matrix H_{N} is a unitary matrix (orthonormal columns), hence its inverse is its conjugate transpose. H is also a real matrix (the entries are the real numbers 0, ±1), hence, H_{N}^{−1} = H_{N}^{t}. We can recover the coefficients c ∈ ℂ^{N} of the vector v ∈ ℂ^{N} by applying the transpose of H_{N} to the vector v. That is the discrete Haar transform,

c = H_{N}^{t} v.

Notice that either from the above matrix formula or from the fact that the Haar basis is an orthonormal basis, the j^{th} coefficient is calculated by taking the inner product of v with the j^{th} Haar vector,

c_{j} = ⟨v, h_{j}⟩.

Example 6.28. Here we display the matrix H_{8}^{t},

H_{8}^{t} = [\begin{matrix}\frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} \\ \frac{−1}{\sqrt{8}} & \frac{−1}{\sqrt{8}} & \frac{−1}{\sqrt{8}} & \frac{−1}{\sqrt{8}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} \\ \frac{−1}{\sqrt{4}} & \frac{−1}{\sqrt{4}} & \frac{1}{\sqrt{4}} & \frac{1}{\sqrt{4}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{−1}{\sqrt{4}} & \frac{−1}{\sqrt{4}} & \frac{1}{\sqrt{4}} & \frac{1}{\sqrt{4}} \\ \frac{−1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{−1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{−1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{−1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{matrix} ] .

Note that the Haar matrix H_{8}^{t} is sparse i.e it has lots of zero entries. If we count the nonzero entries we note that the first and second rows are both full. The third and fourth rows are half-full so they make one full row together. The fifth, sixth, seventh and eighth rows together make one full row of nonzero entries. Adding we get 2 + 1 + 1 = 4 full rows in H_{8}^{t} where n = 8 = 2^{3}. In H_{16}^{t} we get the equivalent of four full rows from the first 8 rowss as in H_{8}^{t} and the 9th to 16th would make up one more. i.e. 2 + 1 + 1 + 1 = 5 full rows in H_{16}^{t} where N = 16 = 2^{4}. The general formula is if N = 2^{j} we get j + 1 full columns, each of length N. So the total number of nonzero entries is N(j + 1) where j = log_{2}N. Hence, multiplying a vector of length N by H_{N} takes N(log_{2}N + 1) multiplications. This implies that

Discrete Haar Transform can be performed in order Nlog_{2}N operations. The brute force operation count is of the same order as the FFT, because of the sparcity of the Haar matrix. Here is an argument for seeking localized basis vectors. Can this count be improved by using some smart algorithm as we did for the discrete Fourier transform? Yes, the operation count can be brought down to order N operations. Such algorithm we will call the Fast Haar Transform, and it is the first example of the Fast Wavelet Transform.

6.8. The Fast Haar Transform

We will try to illustrate with some examples how one could apply the discrete Haar transform more efficiently. We choose here to argue, as we did for the FFT, in terms of matrix decomposition. We will revisit this algorithm when we discuss the multiresolution analysis for wavelets.

First consider N = 2. Let us explicitely compute the action of the matrix H_{2}^{t} on a vector v = [v(0), v(1)]^{t} ∈ ℂ^{2}.

H_{2}^{t}v = \frac{1}{\sqrt{2}} [\begin{matrix}1 & 1 \\ −1 & 1 \end{matrix}] [\begin{matrix}v(0) \\ v(1) \end{matrix}] = [\frac{v(1)+v(0)}{\sqrt{2}} ; \frac{v(1)−v(0)}{\sqrt{2}}] .

Notice that the output consist of scaled averages and differences of v(0) and v(1). In terms of multiplications, we only need two multiplications, namely: \frac{1}{\sqrt{2}}v(0) and \frac{1}{\sqrt{2}}v(1), then we add these two numbers and subtract them. For operation counts what is costly is the multiplications. Applying H_{2}^{t} requires only 2 multiplications.

Consider now the case N = 4. Let us explicitely compute the action of the matrix H_{4}^{t} on a vector v = [v(0), v(1), v(2), v(3)]^{t} ∈ ℂ^{4}.

H_{4}^{t}v = [\begin{matrix}\frac{1}{\sqrt{4}} & \frac{1}{\sqrt{4}} & \frac{1}{\sqrt{4}} & \frac{1}{\sqrt{4}} \\ \frac{−1}{\sqrt{4}} & \frac{−1}{\sqrt{4}} & \frac{1}{\sqrt{4}} & \frac{1}{\sqrt{4}} \\ \frac{−1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 & 0 \\ 0 & 0 & \frac{−1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{matrix} ] [\begin{matrix}v(0) \\ v(1) \\ v(2) \\ v(3) \end{matrix}] = [ \frac{v(3)+v(4)+v(1)+v(0)}{\sqrt{4}} ; \frac{(v(3)+v(2))−(v(1)+v(0))}{\sqrt{4}} ; \frac{v(1)−v(0)}{\sqrt{2}} ; \frac{v(3)−v(2)}{\sqrt{2}} ] .

Denote the scaled averages of pairs of consecutive entries a_{m}^{1}, and their differences by d_{m}^{1}, in our case j = 2.

a_{0}^{1} := \frac{v(1) + v(0)}{\sqrt{2}}, d_{0}^{1} := \frac{v(1) − v(0)}{\sqrt{2}}.

a_{1}^{1} := \frac{v(3) + v(2)}{\sqrt{2}}, d_{1}^{1} := \frac{v(3) − v(2)}{\sqrt{2}}.

Notice that they are the outputs of the 2 × 2 matrix H_{2}^{t} applied to the vectors [v(0), v(1)]^{t}, [v(2), v(3)]^{t} ∈ ℂ^{2}. With this notation, the output of the matrix H_{4}^{t}v is the vector

H_{4}^{t}v = [\frac{a_{1}^{1}+a_{0}^{1}}{\sqrt{2}} \frac{a_{1}^{1}−a_{0}^{1}}{\sqrt{2}} d_{0}^{1} d_{1}^{1}]^{t} .

The first two entries are the output of the matrix H_{2}^{t} applied to the average vector a^{1} = [a_{0}^{1}, a_{1}^{1}], the last two entries are the output of the 2×2 identity matrix applied to the difference vector d^{1} = [d_{0}^{1}, d_{1}^{1}].

Aside 6.29. The superscript 1 is at this stage unnecessary. However we are computing averages and differences of the averages. We can introduce a^{0} = a_{0}^{0} := \frac{a_{1}^{1}+a_{0}^{1}}{\sqrt{2}} , d^{0} := d_{0}^{0} := \frac{a_{1}^{1}−a_{0}^{1}}{\sqrt{2}} . With this notation we have that

H_{4}^{t}v = [a^{0} d^{0} d^{1}]^{t} .

Notice that the vectors d^{0} and a^{0} have dimension 1 = 2^{0}, and d^{1} has dimension 2 = 2^{1}, so the output has dimension 4.

It should now be clear that we can decompose H_{4}^{t} as the product of two 4 × 4 matrices,

H_{4}^{t} = [\begin{matrix}\frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} & 0 \\ \frac{−1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] [\begin{matrix}H_{2}^{t} & 0_{2} \\ 0_{2} & H_{2}^{t} \end{matrix}] .

The second matrix requires 2 × 2 multiplications (involves applying H_{2}^{t} twice), the first matrix involves only 2 multiplications (involves applying H_{2}^{t} once and the 2 × 2 identity matrix which involves no multiplications). The total multiplication count for H_{4}^{t} is then 2 × 2 + 2 = 6 = 2(4 − 1).

Notice that the frist matrix can be decomposed as the product of a block diagonal matrix with H_{2}^{t} and I_{2} (2 × 2 identity matrix) in the diagonal times a permutation matrix, more precisely,

[\begin{matrix}\frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} & 0 \\ \frac{−1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] = [\begin{matrix}H_{2}^{t} & 0_{2} \\ 0_{2} & I_{2} \end{matrix}] [\begin{matrix}1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] , .

Exercise 6.30. Verify that H_{8}^{t} can be decomposed as the product of three 8 × 8 matrices, as follows,

H_{8}^{t} = [\begin{matrix}H_{2}^{t} & 0_{2} \\ 0_{2} & I_{6} \end{matrix}] [\begin{matrix}1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}] [\begin{matrix}H_{4}^{t} & 0_{4} \\ 0_{4} & H_{4}^{t} \end{matrix}] .

Notice that the last matrix requires 2 × 6 = 12 multiplications (involves applying H_{4}^{t} twice), the first matrix involves only 2 multiplications (involves applying H_{2}^{t} once and I_{6} the 6 × 6 identity matrix which involves no multiplications). The middle matrix, is a permutation matrix that we denote P_{8}, and it involves no multiplications. The total multiplication count for H_{8}^{t} is then 2 × 6 + 2 = 14 = 2(8 − 1).

We can see a pattern emerging. The matrix H_{N}^{t}, N = 2^{j} can be decomposed as the product of three N × N matrices, the first and the last are block diagonal matrices, the middle matrix is a permutation matrix.

H_{N}^{t} = [\begin{matrix}H_{2}^{t} \\ & I_{N−2} \end{matrix}] P_{N} [\begin{matrix}H_{N/2}^{t} & 0_{N/2} \\ 0_{N/2} & H_{N/2}^{t} \end{matrix}] .

Where I_{N−2} is the (N − 2) × (N − 2) identity, and P_{N} is an N × N permutation matrix (that is a permutation of the columns of I_{N})

Notice that the last matrix in the decomposition requires 2 × k_{N/2} multiplications (involves applying H_{N/2}^{t} twice, and k_{N/2} is the number of multiplications required to compute H_{N/2}^{t}), the first matrix involves only 2 multiplications (involves applying H_{2}^{t} once and the (N −2)×(N −2) identity matrix which involves no multiplications), the permutation matrix costs nothing in terms of multiplications. The total multiplication count for H_{N}^{t} is then 2 × k_{N/2} + 2.

Exercise 6.31. Verify by induction that the if the above decomposition holds, the number of multiplications required to apply H_{N}^{t} is 2(N − 1). That is, k_{N} = 2(N − 1).

Returning to the averages and the differences but now in dimension N = 2^{j}. We start with the vector v ∈ ℂ^{N} and we create two vectors a^{j−1} and d^{j−1} of dimension N/2 = 2^{j−1} each by taking averages and diferences of pairs of consecutive entries. We now take the vector a^{j−1} ∈ ℂ^{N/2} and create two vectors of dimension N/4 = 2^{j−2} by the same procedure, and keep repeating until we reach dimension 1 = 2^{0}. This scheme can be repeated j times, and can be represented by the following tree or pyramid scheme,

\begin{aligned}v := a^{j} & → & a^{j−1} & → & a^{j−2} & → & ⋯ & a^{1} & a^{0} \\ & ↘ & & ↘ & & ↘ & & ↘ \\ & & d^{j−1} & & d^{j−2} & ⋯ & d^{1} & d^{0} \end{aligned} .

With this notation it can be seen that,

H_{N}^{t}v = [a^{0} d^{0} d^{1} … d^{j−1}]^{t} .

Notice that the vectors a^{j−k}, d^{j−k} have dimension 2^{j−k} = N2^{−k}, so the vector [a^{0} d^{0} d^{1} … d^{j−1}] has the correct dimension,

1 + 1 + 2 + 2^{2} + ⋯ + 2^{j−1} = 2^{j} = N.

We will encounter again averages and differences when discussing wavelets in later chapters.

Exercise 6.32. Can you be more precise in the description of the permutation matrix P_{N} in the decomposition of H_{N}^{t}?