Vectorize 3D discrete Fourier transform
By Sophia Hammond
I'm trying to express the $N$-dimensional discrete Fourier transform (DFT) of an $N$-dimensional array as the product between a matrix and a vector.
If $N=2$ the problem is quite simple: given a $n\times m$ matrix $X$, since $DFT(X) = F_nXF_m^T = F_nXF_m$ (where $F_h$ is the $h \times h$ Fourier matrix) and since $vec(AXB)=(A^T \otimes B)vec(X)$, the desired matrix is $F_n\otimes F_m$ and the vector is $vec(X)$, so that
$$vec(DFT(X)) = (F_n \otimes F_m) vec(X)$$
What can be done if $N\geq 3$, especially if $N=3$?
EDIT: let me better specify the notation, by using MATLAB conventions. An $N$-dimensional array $X$ is a function from $\{1,\dots,n_1\} \times \{1, \dots, n_2\} \times \dots \times \{1, \dots, n_N\}$ to $\mathbb C$. In MATLAB notation, one can access its elements by $X(i_1, i_2, \dots, i_n)$ with suitable indices $i_j$.
The $N$-dimensional DFT of an $N$-dimensional array $X$ is exactly what MATLAB does by typing fftn(X). More details here.
Instead of calculating the $N$-dimensional DFT in $N$ steps by applying the 1-D DFT along each dimension, I would like to find a "one shot" matrix-vector product formulation of this operation. A procedure like
- Unwrap vector $X$ to a vectorized form;
- Multiply $X$ by a suitable matrix;
- Wrap back $X$ to an $N$ dimensional array form;
- The output is the same as
fftn(X).
This need comes from the study of a 3-D imaging problem where the 3-D image has to be reconstructed from an underdetermined set of measurements using a particular solving procedure. The matrix I'm trying to obtain with your help is (almost) the matrix I will use in the solving procedure.
$\endgroup$ 31 Answer
$\begingroup$Suppose $\mathcal{X}\in \mathbb{C}^{I_{1} \times I_{2} \times\dots \times I_{n}} $ is an $\text{n}^\text{th}$ order complex-valued tensor which is your original signal. If you want to get its n-dimensional discrete Fourier transform $\mathcal{Y}\in \mathbb{C}^{I_{1} \times I_{2} \times\dots \times I_{n}} $, you can write this as follows $$ \mathcal{Y}(j_{1},j_{2},\dots,j_{n}) = \sum_{i_{n}=1}^{I_{n}} F_{I_{n}}(j_{n},i_{n}) \dots \sum_{i_{2}=1}^{I_{2}} F_{I_{2}}(j_{2},i_{2}) \sum_{i_{1}=1}^{I_{1}} F_{I_{1}}(j_{1},i_{1})\mathcal{X}(i_{1},i_{2},\dots,i_{n})$$ where $$F_{I_{k}}(j_{k},i_{k}) = \exp \left (-\frac{2\pi i}{I_{k}} (j_{k}-1)(i_{k}-1)\right) \quad \forall k \in \{1,2,\dots,n\}$$ with $i$ being the imaginary unit. In tensor notation, this summation can be expressed with successive k-mode products which is defined as $$ \mathcal{B} = \mathcal{A} \times_{k} \mathbf{U} \iff \mathcal{B}(i_1,\dots,i_{k-1},j,i_{k+1},\dots,i_{n}) = \sum_{i_{k}=1}^{I_{k}}U(j,i_k)\mathcal{A}(i_1,\dots,i_{k-1},i_{k},i_{k+1},\dots,i_{n})$$ for any $\mathcal{A}\in \mathbb{C}^{I_{1} \times\dots \times I_{k-1} \times I_{k} \times I_{k+1} \times \dots \times I_{n}}$ and $\mathbf{U}\in\mathbb{C}^{J\times I_{k}} $ yielding $\mathcal{B}\in \mathbb{C}^{I_{1} \times\dots \times I_{k-1} \times J \times I_{k+1} \times \dots \times I_{n}}$
The final closed form n-dimensional discrete Fourier transform can be written as $$ \mathcal{Y} = \mathcal{X}\times_{1}\mathbf{F}_{I_1}\times_{2}\mathbf{F}_{I_1}\times_{3}\dots\times_{n}\mathbf{F}_{I_n}$$ Vectorization of each side results in $$ \mathbf{y} =\text{vec}(\mathcal{Y}),\mathbf{x} =\text{vec}(\mathcal{X}) \quad \text{then} \quad \mathbf{y}=( \mathbf{F}_{I_n} \otimes \dots \otimes \mathbf{F}_{I_2} \otimes \mathbf{F}_{I_1})\mathbf{x}$$ where $\otimes$ is Kronecker product. So, your original question for 3D DFT can be solved as $$ \mathbf{x}^{*} = \text{3D-DFT}(\mathbf{x}) =(\mathbf{F}_{M} \otimes \mathbf{F}_{L} \otimes \mathbf{F}_{K})\mathbf{x}$$ for $\mathbf{x} =\text{vec}(\mathcal{X})$ and $\mathcal{X} \in \mathbb{C}^{K\times L\times M}$.
$\endgroup$
