{ "cells": [ { "cell_type": "markdown", "id": "5089c228", "metadata": {}, "source": [ "# Circulant Matrices" ] }, { "cell_type": "markdown", "id": "2fef3a02", "metadata": {}, "source": [ "## Overview\n", "\n", "This lecture describes circulant matrices and some of their properties.\n", "\n", "Circulant matrices have a special structure that connects them to useful concepts\n", "including\n", "\n", "- convolution \n", "- Fourier transforms \n", "- permutation matrices \n", "\n", "\n", "Because of these connections, circulant matrices are widely used in machine learning, for example, in image processing.\n", "\n", "We begin by importing some Python packages" ] }, { "cell_type": "code", "execution_count": null, "id": "d9d0ad42", "metadata": { "hide-output": false }, "outputs": [], "source": [ "import numpy as np\n", "from numba import njit\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "id": "7974275a", "metadata": { "hide-output": false }, "outputs": [], "source": [ "np.set_printoptions(precision=3, suppress=True)" ] }, { "cell_type": "markdown", "id": "68b25f5c", "metadata": {}, "source": [ "## Constructing a Circulant Matrix\n", "\n", "To construct an $ N \\times N $ circulant matrix, we need only the first row, say,\n", "\n", "$$\n", "\\begin{bmatrix} c_{0} & c_{1} & c_{2} & c_{3} & c_{4} & \\cdots & c_{N-1} \\end{bmatrix} .\n", "$$\n", "\n", "After setting entries in the first row, the remaining rows of a circulant matrix are determined as\n", "follows:\n", "\n", "\n", "\n", "$$\n", "C=\\left[\\begin{array}{ccccccc}\n", "c_{0} & c_{1} & c_{2} & c_{3} & c_{4} & \\cdots & c_{N-1}\\\\\n", "c_{N-1} & c_{0} & c_{1} & c_{2} & c_{3} & \\cdots & c_{N-2}\\\\\n", "c_{N-2} & c_{N-1} & c_{0} & c_{1} & c_{2} & \\cdots & c_{N-3}\\\\\n", "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n", "c_{3} & c_{4} & c_{5} & c_{6} & c_{7} & \\cdots & c_{2}\\\\\n", "c_{2} & c_{3} & c_{4} & c_{5} & c_{6} & \\cdots & c_{1}\\\\\n", "c_{1} & c_{2} & c_{3} & c_{4} & c_{5} & \\cdots & c_{0}\n", "\\end{array}\\right] \\tag{4.1}\n", "$$\n", "\n", "It is also possible to construct a circulant matrix by creating the transpose of the above matrix, in which case only the\n", "first column needs to be specified.\n", "\n", "Let’s write some Python code to generate a circulant matrix." ] }, { "cell_type": "code", "execution_count": null, "id": "4403cab8", "metadata": { "hide-output": false }, "outputs": [], "source": [ "@njit\n", "def construct_cirlulant(row):\n", "\n", " N = row.size\n", "\n", " C = np.empty((N, N))\n", "\n", " for i in range(N):\n", "\n", " C[i, i:] = row[:N-i]\n", " C[i, :i] = row[N-i:]\n", "\n", " return C" ] }, { "cell_type": "code", "execution_count": null, "id": "a8bac9e9", "metadata": { "hide-output": false }, "outputs": [], "source": [ "# a simple case when N = 3\n", "construct_cirlulant(np.array([1., 2., 3.]))" ] }, { "cell_type": "markdown", "id": "2d4326e5", "metadata": {}, "source": [ "### Some Properties of Circulant Matrices\n", "\n", "Here are some useful properties:\n", "\n", "Suppose that $ A $ and $ B $ are both circulant matrices. Then it can be verified that\n", "\n", "- The transpose of a circulant matrix is a circulant matrix. \n", "- $ A + B $ is a circulant matrix \n", "- $ A B $ is a circulant matrix \n", "- $ A B = B A $ \n", "\n", "\n", "Now consider a circulant matrix with first row\n", "\n", "$$\n", "c = \\begin{bmatrix} c_0 & c_1 & \\cdots & c_{N-1} \\end{bmatrix}\n", "$$\n", "\n", "and consider a vector\n", "\n", "$$\n", "a = \\begin{bmatrix} a_0 & a_1 & \\cdots & a_{N-1} \\end{bmatrix}\n", "$$\n", "\n", "The **convolution** of vectors $ c $ and $ a $ is defined as the vector $ b = c * a $ with components\n", "\n", "\n", "\n", "$$\n", "b_k = \\sum_{i=0}^{n-1} c_{k-i} a_i \\tag{4.2}\n", "$$\n", "\n", "We use $ * $ to denote **convolution** via the calculation described in equation [(4.2)](#equation-eqn-conv).\n", "\n", "It can be verified that the vector $ b $ satisfies\n", "\n", "$$\n", "b = C^T a\n", "$$\n", "\n", "where $ C^T $ is the transpose of the circulant matrix defined in equation [(4.1)](#equation-eqn-circulant)." ] }, { "cell_type": "markdown", "id": "5c16feb7", "metadata": {}, "source": [ "## Connection to Permutation Matrix\n", "\n", "A good way to construct a circulant matrix is to use a **permutation matrix**.\n", "\n", "Before defining a permutation **matrix**, we’ll define a **permutation**.\n", "\n", "A **permutation** of a set of the set of non-negative integers $ \\{0, 1, 2, \\ldots \\} $ is a one-to-one mapping of the set into itself.\n", "\n", "A permutation of a set $ \\{1, 2, \\ldots, n\\} $ rearranges the $ n $ integers in the set.\n", "\n", "A [permutation matrix](https://mathworld.wolfram.com/PermutationMatrix.html) is obtained by permuting the rows of an $ n \\times n $ identity matrix according to a permutation of the numbers $ 1 $ to $ n $.\n", "\n", "Thus, every row and every column contain precisely a single $ 1 $ with $ 0 $ everywhere else.\n", "\n", "Every permutation corresponds to a unique permutation matrix.\n", "\n", "For example, the $ N \\times N $ matrix\n", "\n", "\n", "\n", "$$\n", "P=\\left[\\begin{array}{cccccc}\n", "0 & 1 & 0 & 0 & \\cdots & 0\\\\\n", "0 & 0 & 1 & 0 & \\cdots & 0\\\\\n", "0 & 0 & 0 & 1 & \\cdots & 0\\\\\n", "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n", "0 & 0 & 0 & 0 & \\cdots & 1\\\\\n", "1 & 0 & 0 & 0 & \\cdots & 0\n", "\\end{array}\\right] \\tag{4.3}\n", "$$\n", "\n", "serves as a **cyclic shift** operator that, when applied to an $ N \\times 1 $ vector $ h $, shifts entries in rows $ 2 $ through $ N $ up one row and shifts the entry in row $ 1 $ to row $ N $.\n", "\n", "Eigenvalues of the cyclic shift permutation matrix $ P $ defined in equation [(4.3)](#equation-eqn-examplep) can be computed by constructing\n", "\n", "$$\n", "P-\\lambda I=\\left[\\begin{array}{cccccc}\n", "-\\lambda & 1 & 0 & 0 & \\cdots & 0\\\\\n", "0 & -\\lambda & 1 & 0 & \\cdots & 0\\\\\n", "0 & 0 & -\\lambda & 1 & \\cdots & 0\\\\\n", "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n", "0 & 0 & 0 & 0 & \\cdots & 1\\\\\n", "1 & 0 & 0 & 0 & \\cdots & -\\lambda\n", "\\end{array}\\right]\n", "$$\n", "\n", "and solving\n", "\n", "$$\n", "\\textrm{det}(P - \\lambda I) = (-1)^N \\lambda^{N}-1=0\n", "$$\n", "\n", "Eigenvalues $ \\lambda_i $ can be complex.\n", "\n", "Magnitudes $ \\mid \\lambda_i \\mid $ of these eigenvalues $ \\lambda_i $ all equal $ 1 $.\n", "\n", "Thus, **singular values** of the permutation matrix $ P $ defined in equation [(4.3)](#equation-eqn-examplep) all equal $ 1 $.\n", "\n", "It can be verified that permutation matrices are orthogonal matrices:\n", "\n", "$$\n", "P P' = I\n", "$$" ] }, { "cell_type": "markdown", "id": "cbaeb3cd", "metadata": {}, "source": [ "## Examples with Python\n", "\n", "Let’s write some Python code to illustrate these ideas." ] }, { "cell_type": "code", "execution_count": null, "id": "358057be", "metadata": { "hide-output": false }, "outputs": [], "source": [ "@njit\n", "def construct_P(N):\n", "\n", " P = np.zeros((N, N))\n", "\n", " for i in range(N-1):\n", " P[i, i+1] = 1\n", " P[-1, 0] = 1\n", "\n", " return P" ] }, { "cell_type": "code", "execution_count": null, "id": "a3607c02", "metadata": { "hide-output": false }, "outputs": [], "source": [ "P4 = construct_P(4)\n", "P4" ] }, { "cell_type": "code", "execution_count": null, "id": "b0930b17", "metadata": { "hide-output": false }, "outputs": [], "source": [ "# compute the eigenvalues and eigenvectors\n", "πœ†, Q = np.linalg.eig(P4)" ] }, { "cell_type": "code", "execution_count": null, "id": "3ad2a8ce", "metadata": { "hide-output": false }, "outputs": [], "source": [ "for i in range(4):\n", " print(f'πœ†{i} = {πœ†[i]:.1f} \\nvec{i} = {Q[i, :]}\\n')" ] }, { "cell_type": "markdown", "id": "bce77e41", "metadata": {}, "source": [ "In graphs below, we shall portray eigenvalues of a shift permutation matrix in the complex plane.\n", "\n", "These eigenvalues are uniformly distributed along the unit circle.\n", "\n", "They are the **$ n $ roots of unity**, meaning they are the $ n $ numbers $ z $ that solve $ z^n =1 $, where $ z $ is a complex number.\n", "\n", "In particular, the $ n $ roots of unity are\n", "\n", "$$\n", "z = \\exp\\left(\\frac{2 \\pi j k }{N} \\right) , \\quad k = 0, \\ldots, N-1\n", "$$\n", "\n", "where $ j $ denotes the purely imaginary unit number." ] }, { "cell_type": "code", "execution_count": null, "id": "b48b9766", "metadata": { "hide-output": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(2, 2, figsize=(10, 10))\n", "\n", "for i, N in enumerate([3, 4, 6, 8]):\n", "\n", " row_i = i // 2\n", " col_i = i % 2\n", "\n", " P = construct_P(N)\n", " πœ†, Q = np.linalg.eig(P)\n", "\n", " circ = plt.Circle((0, 0), radius=1, edgecolor='b', facecolor='None')\n", " ax[row_i, col_i].add_patch(circ)\n", "\n", " for j in range(N):\n", " ax[row_i, col_i].scatter(πœ†[j].real, πœ†[j].imag, c='b')\n", "\n", " ax[row_i, col_i].set_title(f'N = {N}')\n", " ax[row_i, col_i].set_xlabel('real')\n", " ax[row_i, col_i].set_ylabel('imaginary')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "173fc239", "metadata": {}, "source": [ "For a vector of coefficients $ \\{c_i\\}_{i=0}^{n-1} $, eigenvectors of $ P $ are also eigenvectors of\n", "\n", "$$\n", "C = c_{0} I + c_{1} P + c_{2} P^{2} +\\cdots + c_{N-1} P^{N-1}.\n", "$$\n", "\n", "Consider an example in which $ N=8 $ and let $ w = e^{-2 \\pi j / N} $.\n", "\n", "It can be verified that the matrix $ F_8 $ of eigenvectors of $ P_{8} $ is\n", "\n", "$$\n", "F_{8}=\\left[\\begin{array}{ccccc}\n", "1 & 1 & 1 & \\cdots & 1\\\\\n", "1 & w & w^{2} & \\cdots & w^{7}\\\\\n", "1 & w^{2} & w^{4} & \\cdots & w^{14}\\\\\n", "1 & w^{3} & w^{6} & \\cdots & w^{21}\\\\\n", "1 & w^{4} & w^{8} & \\cdots & w^{28}\\\\\n", "1 & w^{5} & w^{10} & \\cdots & w^{35}\\\\\n", "1 & w^{6} & w^{12} & \\cdots & w^{42}\\\\\n", "1 & w^{7} & w^{14} & \\cdots & w^{49}\n", "\\end{array}\\right]\n", "$$\n", "\n", "The matrix $ F_8 $ defines a [Discete Fourier Transform](https://en.wikipedia.org/wiki/Discrete_Fourier_transform).\n", "\n", "To convert it into an orthogonal eigenvector matrix, we can simply normalize it by dividing every entry by $ \\sqrt{8} $.\n", "\n", "- stare at the first column of $ F_8 $ above to convince yourself of this fact \n", "\n", "\n", "The eigenvalues corresponding to each eigenvector are $ \\{w^{j}\\}_{j=0}^{7} $ in order." ] }, { "cell_type": "code", "execution_count": null, "id": "d10b4df1", "metadata": { "hide-output": false }, "outputs": [], "source": [ "def construct_F(N):\n", "\n", " w = np.e ** (-complex(0, 2*np.pi/N))\n", "\n", " F = np.ones((N, N), dtype=complex)\n", " for i in range(1, N):\n", " F[i, 1:] = w ** (i * np.arange(1, N))\n", "\n", " return F, w" ] }, { "cell_type": "code", "execution_count": null, "id": "60c6065a", "metadata": { "hide-output": false }, "outputs": [], "source": [ "F8, w = construct_F(8)" ] }, { "cell_type": "code", "execution_count": null, "id": "e849b2e6", "metadata": { "hide-output": false }, "outputs": [], "source": [ "w" ] }, { "cell_type": "code", "execution_count": null, "id": "9216c725", "metadata": { "hide-output": false }, "outputs": [], "source": [ "F8" ] }, { "cell_type": "code", "execution_count": null, "id": "a440dae8", "metadata": { "hide-output": false }, "outputs": [], "source": [ "# normalize\n", "Q8 = F8 / np.sqrt(8)" ] }, { "cell_type": "code", "execution_count": null, "id": "d5ebdd0b", "metadata": { "hide-output": false }, "outputs": [], "source": [ "# verify the orthogonality (unitarity)\n", "Q8 @ np.conjugate(Q8)" ] }, { "cell_type": "markdown", "id": "eb7b2d1a", "metadata": {}, "source": [ "Let’s verify that $ k $th column of $ Q_{8} $ is an eigenvector of $ P_{8} $ with an eigenvalue $ w^{k} $." ] }, { "cell_type": "code", "execution_count": null, "id": "7410b602", "metadata": { "hide-output": false }, "outputs": [], "source": [ "P8 = construct_P(8)" ] }, { "cell_type": "code", "execution_count": null, "id": "e8677611", "metadata": { "hide-output": false }, "outputs": [], "source": [ "diff_arr = np.empty(8, dtype=complex)\n", "for j in range(8):\n", " diff = P8 @ Q8[:, j] - w ** j * Q8[:, j]\n", " diff_arr[j] = diff @ diff.T" ] }, { "cell_type": "code", "execution_count": null, "id": "fbce2122", "metadata": { "hide-output": false }, "outputs": [], "source": [ "diff_arr" ] }, { "cell_type": "markdown", "id": "e6ad3575", "metadata": {}, "source": [ "## Associated Permutation Matrix\n", "\n", "Next, we execute calculations to verify that the circulant matrix $ C $ defined in equation [(4.1)](#equation-eqn-circulant) can be written as\n", "\n", "$$\n", "C = c_{0} I + c_{1} P + \\cdots + c_{n-1} P^{n-1}\n", "$$\n", "\n", "and that every eigenvector of $ P $ is also an eigenvector of $ C $.\n", "\n", "We illustrate this for $ N=8 $ case." ] }, { "cell_type": "code", "execution_count": null, "id": "c4e42377", "metadata": { "hide-output": false }, "outputs": [], "source": [ "c = np.random.random(8)" ] }, { "cell_type": "code", "execution_count": null, "id": "b365d889", "metadata": { "hide-output": false }, "outputs": [], "source": [ "c" ] }, { "cell_type": "code", "execution_count": null, "id": "28bb0f26", "metadata": { "hide-output": false }, "outputs": [], "source": [ "C8 = construct_cirlulant(c)" ] }, { "cell_type": "markdown", "id": "410884e5", "metadata": {}, "source": [ "Compute $ c_{0} I + c_{1} P + \\cdots + c_{n-1} P^{n-1} $." ] }, { "cell_type": "code", "execution_count": null, "id": "f84e32c1", "metadata": { "hide-output": false }, "outputs": [], "source": [ "N = 8\n", "\n", "C = np.zeros((N, N))\n", "P = np.eye(N)\n", "\n", "for i in range(N):\n", " C += c[i] * P\n", " P = P8 @ P" ] }, { "cell_type": "code", "execution_count": null, "id": "daaed7c4", "metadata": { "hide-output": false }, "outputs": [], "source": [ "C" ] }, { "cell_type": "code", "execution_count": null, "id": "e8987b58", "metadata": { "hide-output": false }, "outputs": [], "source": [ "C8" ] }, { "cell_type": "markdown", "id": "92f74844", "metadata": {}, "source": [ "Now let’s compute the difference between two circulant matrices that we have constructed in two different ways." ] }, { "cell_type": "code", "execution_count": null, "id": "f1285f08", "metadata": { "hide-output": false }, "outputs": [], "source": [ "np.abs(C - C8).max()" ] }, { "cell_type": "markdown", "id": "f27396d4", "metadata": {}, "source": [ "The $ k $th column of $ P_{8} $ associated with eigenvalue $ w^{k-1} $ is an eigenvector of $ C_{8} $ associated with an eigenvalue $ \\sum_{h=0}^{7} c_{j} w^{h k} $." ] }, { "cell_type": "code", "execution_count": null, "id": "c5cb146b", "metadata": { "hide-output": false }, "outputs": [], "source": [ "πœ†_C8 = np.zeros(8, dtype=complex)\n", "\n", "for j in range(8):\n", " for k in range(8):\n", " πœ†_C8[j] += c[k] * w ** (j * k)" ] }, { "cell_type": "code", "execution_count": null, "id": "f2601cb7", "metadata": { "hide-output": false }, "outputs": [], "source": [ "πœ†_C8" ] }, { "cell_type": "markdown", "id": "aeea85ca", "metadata": {}, "source": [ "We can verify this by comparing `C8 @ Q8[:, j]` with `πœ†_C8[j] * Q8[:, j]`." ] }, { "cell_type": "code", "execution_count": null, "id": "e181090e", "metadata": { "hide-output": false }, "outputs": [], "source": [ "# verify\n", "for j in range(8):\n", " diff = C8 @ Q8[:, j] - πœ†_C8[j] * Q8[:, j]\n", " print(diff)" ] }, { "cell_type": "markdown", "id": "6c7de95c", "metadata": {}, "source": [ "## Discrete Fourier Transform\n", "\n", "The **Discrete Fourier Transform** (DFT) allows us to represent a discrete time sequence as a weighted sum of complex sinusoids.\n", "\n", "Consider a sequence of $ N $ real number $ \\{x_j\\}_{j=0}^{N-1} $.\n", "\n", "The **Discrete Fourier Transform** maps $ \\{x_j\\}_{j=0}^{N-1} $ into a sequence of complex numbers $ \\{X_k\\}_{k=0}^{N-1} $\n", "\n", "where\n", "\n", "$$\n", "X_{k}=\\sum_{n=0}^{N-1}x_{n}e^{-2\\pi\\frac{kn}{N}i}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "91630230", "metadata": { "hide-output": false }, "outputs": [], "source": [ "def DFT(x):\n", " \"The discrete Fourier transform.\"\n", "\n", " N = len(x)\n", " w = np.e ** (-complex(0, 2*np.pi/N))\n", "\n", " X = np.zeros(N, dtype=complex)\n", " for k in range(N):\n", " for n in range(N):\n", " X[k] += x[n] * w ** (k * n)\n", "\n", " return X" ] }, { "cell_type": "markdown", "id": "6deb6600", "metadata": {}, "source": [ "Consider the following example.\n", "\n", "$$\n", "x_{n}=\\begin{cases}\n", "1/2 & n=0,1\\\\\n", "0 & \\text{otherwise}\n", "\\end{cases}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "d67ab81c", "metadata": { "hide-output": false }, "outputs": [], "source": [ "x = np.zeros(10)\n", "x[0:2] = 1/2" ] }, { "cell_type": "code", "execution_count": null, "id": "1fc30ac3", "metadata": { "hide-output": false }, "outputs": [], "source": [ "x" ] }, { "cell_type": "markdown", "id": "140f8dbf", "metadata": {}, "source": [ "Apply a discrete Fourier transform." ] }, { "cell_type": "code", "execution_count": null, "id": "e35b24e1", "metadata": { "hide-output": false }, "outputs": [], "source": [ "X = DFT(x)" ] }, { "cell_type": "code", "execution_count": null, "id": "4dce16bc", "metadata": { "hide-output": false }, "outputs": [], "source": [ "X" ] }, { "cell_type": "markdown", "id": "16232907", "metadata": {}, "source": [ "We can plot magnitudes of a sequence of numbers and the associated discrete Fourier transform." ] }, { "cell_type": "code", "execution_count": null, "id": "df5b9194", "metadata": { "hide-output": false }, "outputs": [], "source": [ "def plot_magnitude(x=None, X=None):\n", "\n", " data = []\n", " names = []\n", " xs = []\n", " if (x is not None):\n", " data.append(x)\n", " names.append('x')\n", " xs.append('n')\n", " if (X is not None):\n", " data.append(X)\n", " names.append('X')\n", " xs.append('j')\n", "\n", " num = len(data)\n", " for i in range(num):\n", " n = data[i].size\n", " plt.figure(figsize=(8, 3))\n", " plt.scatter(range(n), np.abs(data[i]))\n", " plt.vlines(range(n), 0, np.abs(data[i]), color='b')\n", "\n", " plt.xlabel(xs[i])\n", " plt.ylabel('magnitude')\n", " plt.title(names[i])\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "f1833bef", "metadata": { "hide-output": false }, "outputs": [], "source": [ "plot_magnitude(x=x, X=X)" ] }, { "cell_type": "markdown", "id": "34b4c9dd", "metadata": {}, "source": [ "The **inverse Fourier transform** transforms a Fourier transform $ X $ of $ x $ back to $ x $.\n", "\n", "The inverse Fourier transform is defined as\n", "\n", "$$\n", "x_{n} = \\sum_{k=0}^{N-1} \\frac{1}{N} X_{k} e^{2\\pi\\left(\\frac{kn}{N}\\right)i}, \\quad n=0, 1, \\ldots, N-1\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "9e6cc7ae", "metadata": { "hide-output": false }, "outputs": [], "source": [ "def inverse_transform(X):\n", "\n", " N = len(X)\n", " w = np.e ** (complex(0, 2*np.pi/N))\n", "\n", " x = np.zeros(N, dtype=complex)\n", " for n in range(N):\n", " for k in range(N):\n", " x[n] += X[k] * w ** (k * n) / N\n", "\n", " return x" ] }, { "cell_type": "code", "execution_count": null, "id": "ec9cc880", "metadata": { "hide-output": false }, "outputs": [], "source": [ "inverse_transform(X)" ] }, { "cell_type": "markdown", "id": "4796b16b", "metadata": {}, "source": [ "Another example is\n", "\n", "$$\n", "x_{n}=2\\cos\\left(2\\pi\\frac{11}{40}n\\right),\\ n=0,1,2,\\cdots19\n", "$$\n", "\n", "Since $ N=20 $, we cannot use an integer multiple of $ \\frac{1}{20} $ to represent a frequency $ \\frac{11}{40} $.\n", "\n", "To handle this, we shall end up using all $ N $ of the availble frequencies in the DFT.\n", "\n", "Since $ \\frac{11}{40} $ is in between $ \\frac{10}{40} $ and $ \\frac{12}{40} $ (each of which is an integer multiple of $ \\frac{1}{20} $), the complex coefficients in the DFT have their largest magnitudes at $ k=5,6,15,16 $, not just at a single frequency." ] }, { "cell_type": "code", "execution_count": null, "id": "fe56fa01", "metadata": { "hide-output": false }, "outputs": [], "source": [ "N = 20\n", "x = np.empty(N)\n", "\n", "for j in range(N):\n", " x[j] = 2 * np.cos(2 * np.pi * 11 * j / 40)" ] }, { "cell_type": "code", "execution_count": null, "id": "0b54961b", "metadata": { "hide-output": false }, "outputs": [], "source": [ "X = DFT(x)" ] }, { "cell_type": "code", "execution_count": null, "id": "33472233", "metadata": { "hide-output": false }, "outputs": [], "source": [ "plot_magnitude(x=x, X=X)" ] }, { "cell_type": "markdown", "id": "53a45bbf", "metadata": {}, "source": [ "What happens if we change the last example to $ x_{n}=2\\cos\\left(2\\pi\\frac{10}{40}n\\right) $?\n", "\n", "Note that $ \\frac{10}{40} $ is an integer multiple of $ \\frac{1}{20} $." ] }, { "cell_type": "code", "execution_count": null, "id": "7e81f52f", "metadata": { "hide-output": false }, "outputs": [], "source": [ "N = 20\n", "x = np.empty(N)\n", "\n", "for j in range(N):\n", " x[j] = 2 * np.cos(2 * np.pi * 10 * j / 40)" ] }, { "cell_type": "code", "execution_count": null, "id": "0269775b", "metadata": { "hide-output": false }, "outputs": [], "source": [ "X = DFT(x)" ] }, { "cell_type": "code", "execution_count": null, "id": "a8a47f39", "metadata": { "hide-output": false }, "outputs": [], "source": [ "plot_magnitude(x=x, X=X)" ] }, { "cell_type": "markdown", "id": "781b0f18", "metadata": {}, "source": [ "If we represent the discrete Fourier transform as a matrix, we discover that it equals the matrix $ F_{N} $ of eigenvectors of the permutation matrix $ P_{N} $.\n", "\n", "We can use the example where $ x_{n}=2\\cos\\left(2\\pi\\frac{11}{40}n\\right),\\ n=0,1,2,\\cdots19 $ to illustrate this." ] }, { "cell_type": "code", "execution_count": null, "id": "7cbefef5", "metadata": { "hide-output": false }, "outputs": [], "source": [ "N = 20\n", "x = np.empty(N)\n", "\n", "for j in range(N):\n", " x[j] = 2 * np.cos(2 * np.pi * 11 * j / 40)" ] }, { "cell_type": "code", "execution_count": null, "id": "b3933930", "metadata": { "hide-output": false }, "outputs": [], "source": [ "x" ] }, { "cell_type": "markdown", "id": "f014dfe5", "metadata": {}, "source": [ "First use the summation formula to transform $ x $ to $ X $." ] }, { "cell_type": "code", "execution_count": null, "id": "91de3263", "metadata": { "hide-output": false }, "outputs": [], "source": [ "X = DFT(x)\n", "X" ] }, { "cell_type": "markdown", "id": "3625d311", "metadata": {}, "source": [ "Now let’s evaluate the outcome of postmultiplying the eigenvector matrix $ F_{20} $ by the vector $ x $, a product that we claim should equal the Fourier tranform of the sequence $ \\{x_n\\}_{n=0}^{N-1} $." ] }, { "cell_type": "code", "execution_count": null, "id": "3f40002e", "metadata": { "hide-output": false }, "outputs": [], "source": [ "F20, _ = construct_F(20)" ] }, { "cell_type": "code", "execution_count": null, "id": "5eb9ff7f", "metadata": { "hide-output": false }, "outputs": [], "source": [ "F20 @ x" ] }, { "cell_type": "markdown", "id": "b7e37b90", "metadata": {}, "source": [ "Similarly, the inverse DFT can be expressed as a inverse DFT matrix $ F^{-1}_{20} $." ] }, { "cell_type": "code", "execution_count": null, "id": "6118178c", "metadata": { "hide-output": false }, "outputs": [], "source": [ "F20_inv = np.linalg.inv(F20)\n", "F20_inv @ X" ] } ], "metadata": { "date": 1714442504.8435717, "filename": "eig_circulant.md", "kernelspec": { "display_name": "Python", "language": "python3", "name": "python3" }, "title": "Circulant Matrices" }, "nbformat": 4, "nbformat_minor": 5 }