{
"cells": [
{
"cell_type": "markdown",
"id": "8949524d",
"metadata": {},
"source": [
"# Circulant Matrices"
]
},
{
"cell_type": "markdown",
"id": "cedcbd92",
"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": "74a54642",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"from numba import jit\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fda729d0",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"np.set_printoptions(precision=3, suppress=True)"
]
},
{
"cell_type": "markdown",
"id": "92e1d300",
"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",
"<a id='equation-eqn-circulant'></a>\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": "7bc58f88",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"@jit\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": "82a406f2",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# a simple case when N = 3\n",
"construct_cirlulant(np.array([1., 2., 3.]))"
]
},
{
"cell_type": "markdown",
"id": "b7adcf84",
"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",
"<a id='equation-eqn-conv'></a>\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": "55d7d794",
"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",
"<a id='equation-eqn-examplep'></a>\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": "a42a801a",
"metadata": {},
"source": [
"## Examples with Python\n",
"\n",
"Letβs write some Python code to illustrate these ideas."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f963ffae",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"@jit\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": "8199ae1a",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"P4 = construct_P(4)\n",
"P4"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "aad9c1c3",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# compute the eigenvalues and eigenvectors\n",
"π, Q = np.linalg.eig(P4)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "41a1eea3",
"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": "c6025ded",
"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": "94fb4e0e",
"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": "25f48f99",
"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": "cde20594",
"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": "6fc4b73b",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"F8, w = construct_F(8)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0afe018a",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"w"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "44c4ae5e",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"F8"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b657c1e4",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# normalize\n",
"Q8 = F8 / np.sqrt(8)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e14757ca",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# verify the orthogonality (unitarity)\n",
"Q8 @ np.conjugate(Q8)"
]
},
{
"cell_type": "markdown",
"id": "b5e0ac50",
"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": "967c91a0",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"P8 = construct_P(8)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "720ef88b",
"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": "213877bc",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"diff_arr"
]
},
{
"cell_type": "markdown",
"id": "40d1e986",
"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": "78ef6c05",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"c = np.random.random(8)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "56512a8b",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"c"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4637ff62",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"C8 = construct_cirlulant(c)"
]
},
{
"cell_type": "markdown",
"id": "4ee497b1",
"metadata": {},
"source": [
"Compute $ c_{0} I + c_{1} P + \\cdots + c_{n-1} P^{n-1} $."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ef09cdbe",
"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": "af08dbfa",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"C"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f0e28675",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"C8"
]
},
{
"cell_type": "markdown",
"id": "c8787d93",
"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": "d65b4f24",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"np.abs(C - C8).max()"
]
},
{
"cell_type": "markdown",
"id": "417fb134",
"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": "edd4ecde",
"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": "74cc4aa7",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"π_C8"
]
},
{
"cell_type": "markdown",
"id": "f4ea26e4",
"metadata": {},
"source": [
"We can verify this by comparing `C8 @ Q8[:, j]` with `π_C8[j] * Q8[:, j]`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d460d839",
"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": "8950cf88",
"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": "7671e861",
"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": "8f1a3a02",
"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": "0f849e9a",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"x = np.zeros(10)\n",
"x[0:2] = 1/2"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fdb71970",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"x"
]
},
{
"cell_type": "markdown",
"id": "1dd47b60",
"metadata": {},
"source": [
"Apply a discrete Fourier transform."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d4d0ee6a",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"X = DFT(x)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "41556c4c",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"X"
]
},
{
"cell_type": "markdown",
"id": "70c77fd0",
"metadata": {},
"source": [
"We can plot magnitudes of a sequence of numbers and the associated discrete Fourier transform."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4b0b43f3",
"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": "45be5db9",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_magnitude(x=x, X=X)"
]
},
{
"cell_type": "markdown",
"id": "df8c907a",
"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": "e516c5ab",
"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": "93eb9950",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"inverse_transform(X)"
]
},
{
"cell_type": "markdown",
"id": "5cb92dfc",
"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": "adba59a8",
"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": "b66d404b",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"X = DFT(x)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a3069e46",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_magnitude(x=x, X=X)"
]
},
{
"cell_type": "markdown",
"id": "e8de18e6",
"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": "a7795a27",
"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": "70d3a27d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"X = DFT(x)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6756e5a0",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_magnitude(x=x, X=X)"
]
},
{
"cell_type": "markdown",
"id": "e1ea2de4",
"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": "d6262973",
"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": "cf370a28",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"x"
]
},
{
"cell_type": "markdown",
"id": "d961167e",
"metadata": {},
"source": [
"First use the summation formula to transform $ x $ to $ X $."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "28c94b78",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"X = DFT(x)\n",
"X"
]
},
{
"cell_type": "markdown",
"id": "bbce688c",
"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": "c9d52c83",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"F20, _ = construct_F(20)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "19a3342a",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"F20 @ x"
]
},
{
"cell_type": "markdown",
"id": "65f5be78",
"metadata": {},
"source": [
"Similarly, the inverse DFT can be expressed as a inverse DFT matrix $ F^{-1}_{20} $."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d2bfdbb8",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"F20_inv = np.linalg.inv(F20)\n",
"F20_inv @ X"
]
}
],
"metadata": {
"date": 1751441149.6289139,
"filename": "eig_circulant.md",
"kernelspec": {
"display_name": "Python",
"language": "python3",
"name": "python3"
},
"title": "Circulant Matrices"
},
"nbformat": 4,
"nbformat_minor": 5
}