**ECSE 343: Numerical Methods in Engineering** *Assignment 4: _Discrete Fourier Transform_* Due: Thursday, March 31^th April 7^th, 2022 at 11:59pm EST on [myCourses](https://mycourses2.mcgill.ca/) Final weight: __20%__ We suggest you take time to read the entire handout _before_ starting your work. # Assignment Policies Download and modify the standalone `Python` script, renaming the file according to your student ID as `YourStudentID`.py $\longrightarrow$ e.g., if your id is **234567890**, your submitted filename should be **234567890.py** !!! Tip: Every time you submit a new file, your previous submission will be overwritten. We will only grade the **final submitted file**, so feel free to submit often as you progress through the assignment. ## Late policy and plagiarism All the assignments are to completed individually, unless stated otherwise. You are expected to respect the [late policy](https://www.cim.mcgill.ca/~derek/ecse343.html#AssignmentLatePolicy) and [collaboration/plagiarism](https://www.cim.mcgill.ca/~derek//ecse343.html#collaboration&plagiarism) polices. If you have questions or doubts regarding these policies, please **contact the TA** or speak to the Professor *in class*. ## Library `import`s and writing tests The code we provide includes a superset of all the `import`s you are allowed to use when implementing your solution. !!! ERROR: You must not use any additional `import`s for your solution, other than the ones provided by us. Do not edit the `import`s at the top of the `.py` file. Doing so will result in an assignment grade of **0%**. One exception to this policy is that, in the `__main__` block, you are allowed to `import` any additional routines you would like in order to test your solution code. Do __not__ use `import` in any of your solution routines __nor__ at the top of the Python file; only optionally in `__main__`. !!! Tip: It is _very difficult_ to successfully complete this assignment without writing your own tests in `__main__`. We strongly recommend against relying exclusively on the example tests we provide. Think about what kind of input/output behaviours you can write and test against. # Discrete Fourier Transform Given an input continuous function $f$ sampled at $N$ equally-spaced discrete points on an integer lattice $x \in \{0, \ldots, {N-1}\}$ in the primal domain (e.g., time, space, etc.), we denote this discretized signal as $f(x) \equiv f[x]$. Note that, in general, $f$ can be a _complex-valued_ function (and, so too the sampled values $f[\cdot]$); however, our examples will only consider real-valued input functions $f : \mathbb{R} \rightarrow \mathbb{R}$. The __discrete Fourier transform__ (DFT) $\hat{f}[\omega]$ of the sampled signal is then given by: \begin{equation} \label{eq:dft} \hat{f}[\omega] = \sum\limits_{k=0}^{N-1} f[k] \exp{\left(\frac{-2\pi \imath k \omega}{N}\right)}, \end{equation} which can be equivalently expressed as a matrix-vector operation (see lecture slides) as $$ \mathbf{\widehat{f}} = \mathbf{M}_N\, \mathbf{f}~, $$ where $\mathbf{f} = (f[0], f[1], \ldots, f[{N-1}])$ is a vector of all the sampled function values, $\mathbf{M}_N$ encodes the linear operator form of the DFT for a length-$N$ input, and $\mathbf{\widehat{f}} = (\hat{f}[0], \hat{f}[1], \ldots, \hat{f}[N-1])$ is a vector of the sampled DFT of $f$. The (discrete) signal $\hat{f}[x]$ can be interpreted as a discrete sampling of the continuous Fourier transform $\mathcal{F}\{f_N(x)\}$ of an $N$-periodic summation of $f(x)$ -- i.e., $f(x)$ repeated periodically. Note that, regardless of whether the input signal is real- or complex-valued, the DFT $\hat{f}[x]$ **always** has complex-valued elements: $\operatorname{Re}(\hat{f})$ encodes the frequency amplitudes and $\operatorname{Im}(\hat{f})$ their phases. The __inverse discrete Fourier transform__ (iDFT) allows us to transform back from this (discrete, complex-valued) sampled frequency spectrum to the original (possibly only real-valued) sampled signal in the primal domain, as \begin{equation} \label{eq:idft} f[x] = \frac{1}{N}\sum\limits_{k=0}^{N-1} \hat{f}[k] \exp{\left(\frac{2\pi \imath k x}{N}\right)}. \end{equation} Note the similarities in the mathematical definitions of the DFT and iDFT in Equations \ref{eq:dft} and \ref{eq:idft}, comprising of a sign change in the complex exponential and different normalization constants. You will leverage these similarities to avoid code duplication: your DFT and iDFT implementations will rely on an implementation of a _generalized_ (i.e., bidirectional) DFT. Specifically, you will implement this generalized DFT formula for a sampled input $\Upsilon[\cdot]$, $$ \Psi[a;\mbox{s}] = \sum\limits_{b=0}^{N-1} \Upsilon[b] \exp{\left(\mbox{s}\frac{2 \pi \imath b a}{N}\right)}, $$ where we obtain the DFT by setting $\Upsilon = f$ and $s = -1$ (and $a \equiv \omega$ and $b \equiv x$), and the iDFT with $\Upsilon = \tfrac{1}{N} \hat{f}$ and $s = 1$ (and $a \equiv x$ and $b \equiv \omega$). Your implementations will act on, and return, an **entire vector** of sampled primal- (i.e., $\mathbf{f}$) or frequency-domain (i.e., $\widehat{\mathbf{f}}$) values. !!! warning: **Deliverable 1 [17.5%]**
Implement the Generalized DFT Routine: Complete the implementation of the generalized DFT function `DFT` as specified by $\Psi$ and with the default parameter setting of $\mbox{s} = -1$. As discussed above, the generalized DFT equation $\Upsilon$ can be used to realize both a forward and inverse DFT. !!! warning: **Deliverable 2 [2.5%]**
Implement the inverse DFT Routine: Complete the implementation of the inverse DFT function `iDFT` that relies on calling your generalized `DFT` routine, above. !!! note: A careful `numpy` implementation should not rely on __any__ `for` loops or related `Python` operations (e.g., loop comprehensions). You may also wish to leverage the fact that $\exp{(\imath x)} = \cos{x} + \imath \sin{x}$. # 2D Discrete Fourier Transform Let's now consider a generalization to 2D, where we sample a continuous function $f : \mathbb{R}^2 \rightarrow \mathbb{R}$ on the (square) 2D integer lattice leading to a discretized 2D signal $f[x,y]$ at $(x,y) \in \{0, \ldots, N-1\}^2$. We can similarly define the sampled 2D DFT $\hat{f}[\omega_x, \omega_y]$ of $f[x,y]$ as \begin{equation} \label{eq:dft2d} \hat{f}[\omega_x, \omega_y] = \sum\limits_{k_x=0}^{N-1} \sum\limits_{k_y=0}^{N-1} f[k_x, k_y] \exp{\left(\frac{-2\pi \imath (k_x \omega_x + k_y \omega_y)}{N^2}\right)}, \end{equation} which admits a similarly compact matrix-vector expression as $$ \mathbf{\widehat{F}} = \mathbf{M}_{N,N}\, \mathbf{F}~, $$ where we can arrange the sampled primal function values in a matrix $\mathbf{F}$ with elements $(\mathbf{F})_{i,j} = f[i,j]$, with $\mathbf{M}_{N,N}$ encoding the linear operator form of the 2D DFT for a square length-$N \times N$ sampled input, and $\mathbf{\widehat{F}}$ is the matrix of the sampled 2D DFT coefficients with $(\mathbf{\widehat{F}})_{i,j} = \hat{f}[i,j]$. Again, similarly to the 1D setting, DFT elements $\hat{f}[\cdot, \cdot]$ are always complex-valued with $\operatorname{Re}(\hat{f})$ encoding frequency amplitudes and $\operatorname{Im}(\hat{f})$ their phases. The 2D inverse DFT is now \begin{equation} \label{eq:idft2d} f[x,y] = \frac{1}{N^2} \sum\limits_{k_x=0}^{N-1} \sum\limits_{k_y=0}^{N-1} \hat{f}[k_x, k_y] \exp{\left(\frac{2\pi \imath (k_x x + k_y y)}{N^2}\right)} \end{equation} and a similar extension of the _generalized DFT_ routine applies in the 2D setting. ## Leveraging Factorization After little manipulation, one can observe that the 2D DFT and iDFT in Equations \ref{eq:dft2d} and \ref{eq:idft2d} can be _factorized_ into several applications of their 1D counterparts. By applying a 1D DFT independently to each of the _rows_ of $\mathbf{F}$, before applying the 1D DFT independently to the _columns_ of the (partially-)transformed $\mathbf{F}$, we arrive at the 2D DFT transform of $\mathbf{F}$. In matrix form, this amounts to \begin{equation} \label{eq:dft2dmatrix} \mathbf{\widehat{F}} = \mathbf{M}_{N,N}\, \mathbf{F} = \mathbf{M}_{N} \left(\mathbf{M}_{N} \mathbf{F}^{\text{T}}\right)^{\text{T}}~, \end{equation} You should leverage the factorization to accelerate your implementation. !!! warning: **Deliverable 3 [22.5%]**
Implement a Generalized 2D DFT Routine: Complete the implementation of the generalized 2D DFT function `DFT2D`, which extends the 1D generalized routine to 2D. !!! note: Keep in mind that you __don't__ have to explicitly form these matrices at any point in order to take advantage of the aforementioned factorization property. !!! warning: **Deliverable 4 [2.5%]**
Implement the inverse 2D DFT Routine: Complete the implementation of the inverse 2D DFT function `iDFT2D` that relies on calling your generalized `DFT2D` routine. # Fourier for Image Convolution & Deconvolution In Assignments 1 and 3, we explored convolution and deconvolution using linear operators that act on sampled signals _in the primal (i.e., pixel, spatial) domain_. We constructed the convolution operators as a function of the underlying convolution kernel, and we implemented a _wrap_-based boundary condition in instances during convolution where kernels only partially overlapped the input domain. Manually treating the indexing, which additionally required "_flattening_" 2D images into 1D vectors, required care. Frequency-space representations are particularly well-suited to certain applications, among them convolution. We will use the DFT to perform image convolution and (regularized) deconvolution, leveraging these advantages: 1. convolution in the primal domain amounts to multiplication in the frequency domain, and 2. in the frequency domain, we can perform computation _directly_ with native 2D signals (e.g., the image and kernel), instead of having to "flatten" them for the sake of formulating the forward (and backward) problems as systems of linear equations. ## 2D Fourier Convolution The (forward) convolution of a 2D function $I(x,y)$ with a 2D kernel $k(x,y)$ yields a 2D output function $B(x,y) = I(x,y) \circledast k(x,y)$. In the image convolution settings we have explored, $I(x,y) \equiv I[x,y]$ is a discretely-sampled input image, $k(x,y) \equiv k[x,y]$ is a (typically square, odd edge-length) discrete kernel, and $B(x,y) \equiv B[x,y]$ is the blurred output image. Given a matrix $\mathbf{I}$ of the 2D sampled image pixel values, i.e., with $(\mathbf{I})_{i,j} = I(i,j)$, we obtain its DFT $\mathbf{\widehat{I}}$ as per Equation \ref{eq:dft2dmatrix}. For the kernel, we can also form a matrix $\mathbf{K}$ of its sampled values before computing its DFT $\mathbf{\widehat{K}}$; here, however, we need take into account two additional points: 1. the kernel matrix $\mathbf{K}$ needs to be 0-padded to match the size of the image matrix $\mathbf{I}$ -- this will be needed to allow for an element-wise product of their DFTs, later on; and, 2. the "location"/placement of the non-zero kernel elements in the matrix $\mathbf{K}$ needs to be chosen carefully so to respect the same "wrap" boundary condition behaviour that we maintained in the primal domain[^kernel]. !!! warning: **Deliverable 5 [25%]**
Fourier Kernel Matrix: Complete the implementation of the `FourierKernelMatrix` that takes as input an $n \times n$ input blur kernel (with odd $n$), a 2-tuple of the shape (i.e., dimensions) of the input/output image (`image_shape`), and (optionally) a forward 2D DFT function (`DFT2D_func`) that you will call in your `FourierKernelMatrix` solution. The function returns the $\mathbf{\widehat{K}}$ matrix corresponding to the DFT of $\mathbf{K}$. The non-zero kernel elements in $\mathbf{K}$ should be placed in the matrix so to admit the same post-convolution wrap-based blurring behaviour as we saw in Assignments 1 and 3. !!! Tip: When implementing and debugging `FourierKernelMatrix`, you don't necessarily have to rely on your `DFT2D` routine (e.g., by passing something other than `DFT2D` as the `DFT2D_func` parameter). [^kernel]: The shifting property of convolutions can give you a hint of exactly how to structure and place the non-zero kernel matrix elements in $\mathbf{K}$, before taking its DFT. Given the appropriately constructed $\mathbf{K}$, and its DFT $\mathbf{\widehat{K}}$, we can obtain the DFT $\mathbf{\widehat{B}}$ of the sampled output image $\mathbf{B}$ (i.e., with $(\mathbf{B})_{i,j} = B(i,j)$) as $$ \mathbf{\widehat{B}} = \mathbf{\widehat{I}} \star \mathbf{\widehat{K}}~, $$ where $\star$ is an _element-wise_ product. Taking the real component of the inverse DFT of $\mathbf{\widehat{B}}$, we arrive at the final output blurred image matrix $\mathbf{B}$. ## 2D Fourier Deconvolution In an unregularized setting, we can readily express the deconvolution problem using the DFT as \begin{equation} \label{eq:dftdeblur} \mathbf{\widehat{I}} = \mathbf{\widehat{B}} \star^{-1} \mathbf{\widehat{K}}~, \end{equation} where $\star^{-1}$ is an _element-wise_ division. Taking the real component of the inverse DFT of $\mathbf{\widehat{I}}$ yields the unblurred image matrix $\mathbf{I}$. Since the conditioning of the DFT and iDFT are 1, the underlying conditioning of the convolution/deconvolution system remain as (un)stable as they were in the primal domain. !!! note: Unlike Assignments 1 and 3, we will not be providing you with the ground truth deblurred image for your reference, below. You can code up example tests for forward convolution and non-noisy/unregularized deblurring to sanity check your DFT/iDFT implementations; your code from A3 can come in handy, here! ### Laplacian Regularization We provide you with a blurred image polluted by mild noise, as well as an $n \times n$ blur kernel (which you will need to extend appropriately into $\mathbf{K}$, as above). As mentioned briefly in Assignment 3, _Laplacian regularization_ is well-suited to image deblurring problems: it adds a regularization term that penalizes large second derivatives in image-space, a statistical property of many natural images. In the primal domain, the $3 \times 3$ _Laplacian regularization kernel_ is \begin{align} l(x,y) \equiv l[x,y] = \begin{bmatrix} 0 & 1 & 0\\ 1 & -4& 1\\ 0 & 1 & 0\\ \end{bmatrix} \end{align} which you can extend and arrange into a matrix $\mathbf{L}$ with dimensions that match the image resolution and with structure that respects the _wrapping_ boundary condition treatment (i.e., similar to the blur kernel's matrix $\mathbf{K}$ structure). We can solve the regularized least squares deconvolution problem, expressed in the primal domain, \begin{equation} \min_{I} \|B(x,y) - I(x,y) \circledast k(x,y) \|_2^2 + \lambda \|I(x,y) \circledast l(x,y) \|_2^2 \end{equation} in the frequency domain as \begin{equation} \label{eq:dftlaplacian} \mathbf{\widehat{I}} = \mathbf{\widehat{B}} \star \underbrace{\left(\frac{\mathbf{\widehat{K}}}{\left|\mathbf{\widehat{K}}\right|^2 + \lambda \left|\mathbf{\widehat{L}}\right|^2} \right)}_{\mathbf{\widehat{K}}^{-1}_{\text{reg}}}~, \end{equation} where $\mathbf{\widehat{L}}$ is the DFT of $\mathbf{L}$ and the absolute value $| \cdot |$ of a complex value $a + b \imath$ is $\sqrt{a^2 + b^2}$. Taking the real component of the inverse DFT of $\mathbf{\widehat{I}}$ yields our regularized and deblurred image matrix $\mathbf{I}$. Note that in Equation \ref{eq:dftlaplacian} the regularized kernel $\mathbf{\widehat{K}}^{-1}_{\text{reg}}$ is already in an inverted form, which is why an element-wise multiplication ($\star$) is employed instead of an element-wise division ($\star^{-1}$, as in Equation \ref{eq:dftdeblur}); i.e., with $\lambda = 0$, Equations \ref{eq:dftdeblur} and \ref{eq:dftlaplacian} are equivalent. !!! warning: **Deliverable 6 [30%]**
Laplacian-regularized Fourier Deconvolution Kernel: Complete the implementation of `LaplacianInverseFourierKernel` to form and return $\mathbf{\widehat{K}}^{-1}_{\text{reg}}$. The function parameters match those of `FourierKernelMatrix`. Choose and hardcode a $\lambda$ regularization coefficient value that yields a qualitatively good deblurred image (i.e., when applied to our `blurred_noisy_image` test data using Equation \ref{eq:dftlaplacian}; see our test code in `__main__`: _part of your grade for this deliverable will be based on whether we can discern important details from the deblurred image_. # You're Done! Congratulations, you've completed the 4^th assignment. Review the submission procedures and guidelines before submitting your `Python` script file assignment solution. Recall that any test code you incude in your mainline (`__main__`) will not be graded, however it must still respect the [collaboration/plagiarism](https://www.cim.mcgill.ca/~derek//ecse343.html#collaboration&plagiarism) polices.