{
"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
}