32  Fourier Analysis as Change of Basis

Fourier analysis provides a framework for representing functions and data in terms of frequency components. In engineering, this transformation is critical for signal processing, image compression, and solving differential equations where the Laplacian operator is diagonalized by the Fourier basis.

32.1 Discrete Fourier Transform (DFT)

NoteDefinition: Periodic Grid Signal

A discrete signal \(\bx\in\fC^n\) is interpreted as one period of a periodic sequence: \[ \begin{align} x_{j+n}=x_j. \end{align} \] The DFT represents this signal in the basis of discrete complex exponentials.

NoteDefinition: DFT Matrix

Unitary matrix \(\mathbf{F} \in \fC^{n \times n}\) with entries: \[ \begin{align} F_{jk} = \frac{1}{\sqrt{n}} \omega_n^{jk}, \quad \omega_n = e^{-2\pi i/n}. \end{align} \]

  • Unitary Property: \(\mathbf{F}^* \mathbf{F} = \bI\). Orthonormal columns.

  • Inverse DFT: \(\mathbf{F}^{-1} = \mathbf{F}^*\).

NoteTheorem: Parseval’s Theorem

The DFT is an isometry: \(\|\mathbf{F}\bx\|_2 = \|\bx\|_2\). Total energy is preserved between time and frequency domains.

Since \(\mathbf{F}\) is unitary, \(\mathbf{F}^*\mathbf{F}=\bI\). Hence \[ \begin{align} \|\mathbf{F}\bx\|_2^2 = (\mathbf{F}\bx)^*(\mathbf{F}\bx) = \bx^*\mathbf{F}^*\mathbf{F}\bx = \bx^*\bx = \|\bx\|_2^2. \end{align} \] Taking square roots gives the claim.

WarningExercise

  1. Use the result above to compute the \((j,\ell)\) entry of \(\mathbf{F}^*\mathbf{F}\).

  2. Evaluate the finite geometric sum that appears.

  3. Conclude \(\mathbf{F}^*\mathbf{F}=\bI\).

  4. Explain why this is the finite Fourier analog of orthonormal coordinates.

NoteDefinition: Frequency Bins

For samples on \([0,2\pi)\), the integer frequency attached to index \(k\) is usually ordered as \[ \begin{align} 0,1,2,...,\lfloor n/2\rfloor,-\lfloor(n-1)/2\rfloor,...,-1. \end{align} \] This is the ordering returned by np.fft.fftfreq after scaling by \(n\).

TipRemark

(Normalization convention) Some libraries put the factor \(1/n\) in the inverse transform instead of using the unitary factor \(1/\sqrt n\) in both directions. The formulas change by constants, but the frequency content is the same.

NoteDefinition: Nyquist Frequency and Aliasing

On \(n\) equally spaced samples of a \(2\pi\)-periodic signal, frequencies that differ by multiples of \(n\) are indistinguishable: \[ \begin{align} e^{i(k+n)x_j}=e^{ikx_j}. \end{align} \] The largest resolvable unsigned frequency is the Nyquist frequency, approximately \(n/2\) modes per period.

WarningExercise

  1. Let \(x_j=2\pi j/n\). Prove \(e^{i(k+n)x_j}=e^{ikx_j}\) for every grid point.

  2. For \(n=8\), list the frequency bin ordering from the result above.

  3. Sample \(\sin(7x)\) on \(n=8\) points and identify its aliased frequency.

32.2 Convolution and Circulant Matrices

NoteDefinition: Circulant Matrix

Matrix \(\mathbf{C}\) where each row is a cyclic shift of the first column \(\bc\).

  • Diagonalization: Every circulant matrix is diagonalized by the DFT: \(\mathbf{C} = \mathbf{F}^* \text{diag}(\mathbf{F}\bc) \mathbf{F}\).

  • Eigenvalues: The eigenvalues of \(\mathbf{C}\) are the DFT coefficients of its first column.

NoteTheorem: Convolution Theorem

Circular convolution \(\bz = \bc * \bx\) is equivalent to pointwise multiplication in frequency: \[ \begin{align} \mathbf{F}(\bc * \bx) = (\mathbf{F}\bc) \odot (\mathbf{F}\bx). \end{align} \] Complexity: Reduced from \(O(n^2)\) (direct) to \(O(n \log n)\) (FFT).

NoteExample

(Circular convolution for \(n=3\)) If \(\bc=(c_0,c_1,c_2)^T\), then \[ \begin{align} \bc * \bx = \begin{pmatrix} c_0&c_2&c_1\\ c_1&c_0&c_2\\ c_2&c_1&c_0 \end{pmatrix} \begin{pmatrix}x_0\\x_1\\x_2\end{pmatrix}. \end{align} \] The wraparound terms are a consequence of the periodic convention in the result above.

32.3 The Fast Fourier Transform (FFT)

NoteTheorem: Cooley-Tukey Algorithm

Computes the DFT in \(O(n \log n)\) flops for \(n = 2^p\).

  • Mechanism: Recursively splits \(n\)-point DFT into two \(n/2\)-point DFTs (even/odd indices).

  • NumPy: Use np.fft.fft and np.fft.ifft.

WarningExercise

  1. Split the DFT sum into even and odd input indices.

  2. Show that each part is an \(n/2\)-point DFT multiplied by phase factors.

  3. Write the recurrence \(T(n)=2T(n/2)+O(n)\).

  4. Solve the recurrence to obtain the cost in the result above.

TipRemark

Spectral Differentiation: For smooth periodic signals, differentiate by multiplying frequency components \(\hat{x}_k\) by \(ik\). Provides exponential accuracy compared to polynomial accuracy of finite differences.

NoteTheorem: Fourier Differentiation Rule

For a smooth \(2\pi\)-periodic function with Fourier series \[ \begin{align} f(x)=\sum_k \hat f_k e^{ikx}, \end{align} \] the derivative has coefficients \[ \begin{align} \widehat{f'}_k=ik\hat f_k. \end{align} \]

32.4 Exercises

WarningExercise

  1. Write out the \(4 \times 4\) DFT matrix explicitly. Verify it is unitary using the result above.

  2. Solve the periodic Poisson equation \(-u'' = f\) via FFT. (Hint: Laplacian is circulant).

  3. Compare the speed of np.fft.fft(x) vs. F @ x for \(n=4096\).

  4. Implement spectral differentiation for \(f(x) = \sin(x)\) using the result above and check error vs. \(n\).