"Make sure you have run all cells in your notebook in order before running the cell below, so that all images/graphs appear in the output. The cell below will generate a zip file for you to submit. **Please save before exporting!**"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7d10d7ed",
"metadata": {
"deletable": false,
"editable": false
},
"outputs": [],
"source": [
"# Save your notebook first, then run this cell to export your submission.\n",
"grader.export(pdf=False, run_tests=True)"
]
},
{
"cell_type": "markdown",
"id": "159a45fb",
"metadata": {},
"source": [
" "
]
}
],
"metadata": {
...
...
%% Cell type:code id:9692b307 tags:
%% Cell type:code id:0b9af790 tags:
``` python
# Initialize Otter
importotter
grader=otter.Notebook("Projections and signal restoration.ipynb")
Let us consider a signal with $N$ elements, i.e. a vector in $\mathbb{R}^N$.
Under our observations conditions, we can only recover partially the signal's values, the other remain unknown, i.e. we know that:
$x[k] = x_k$ for $k\in E = \{e_0, e_1, \dots,e_{M-1}\}$, with $E\in\mathbb{N}^M$ and $e_i\neq e_j \forall i\neq j$ (i.e. each known index is unique).
We also make the assumption that the signal is contained within **lower frequencies of the spectrum**. This can be expressed using the (normalized) Fourier matrix you have constructed last week $\hat{W}$.
In this notebook, we will try to reconstruct the signal by projecting its observation successively on the Fourier subspace defined above, then back to its original domain (with the constraint regarding its values), then on the Fourier domain, etc. This is a simplified version of a more general method called ["Projection onto convex sets" (POCS)](https://en.wikipedia.org/wiki/Projections_onto_convex_sets). The method is illustrated by the figure below (of course you do not project on lines but on spaces having larger dimension).
![POCS](images/pocs.png"POCS")
### 1. Low-pass filter
Let us first create example signals to validate our implementation of the filtering operation.
Let $X = \hat{W}x$, $X$ is the *discrete Fourier transform* of the input signal $x$. The frequency condition can then be rewritten as
$X[k] = 0$ for $w_c < k\leq N-w_c$.
2. Implement the function that returns a $N\times N$ matrix $P$, s.t. $PX$ satisfies the above condition for a given $w_c$. Make sure the input values are valid, if not raise a `ValueError` exception.
3. Now let us try the filtering on the test signals and make sure it behaves appropriately. Using the matrix $P$ defined above, choose the parameter $w_c$ approiately s.t. the filtered signals retain $w_1$ and $w_2$ but discard $w_3$.
In order to express the condition on $x[k]$ as a pure matrix operation we need to make use of an extension of the input signal and design a slightly different Fourier transform matrix to use properly those extended signals.
Let us denote by $x_E$ the vector from $\mathbb{R}^M$ containing the known values of $x$, i.e. the $x_k$ for $k\in E$.
For each vector $y$ in $\mathbb{R}^N$ we define its extension as $\tilde{y} = \begin{pmatrix}y \\ x_E\end{pmatrix}$.
**This notation will remain throughout the notebook**, i.e. a vector with a tilde denotes its extension made by adding the $x_E$ values at the end.
5. Let us define the matrix $B_E$ and $y=B_E x$, s.t. $y[k] = 0$ if $k\in E$ and $y[k] = x[k]$ otherwise. Write a function that returns its extended version $\tilde{B_E}$ s.t. $\tilde{y}=\tilde{B_E}\tilde{x}$, for any $x\in\mathbb{R}^N$.
- $\tilde{B_E}$ is a square matrix of size $N+M$.
- Check the validity of parameters and raise a `ValueError` in case of invalid inputs.
Let us know design $C_E$ as an operator from $\mathbb{R}^{N+M} \rightarrow \mathbb{R}^{N+M}$ such that when applied to any extended vector $\tilde{z}$ s.t. $\tilde{z}[k] = 0, \forall k\in E$ as input (i.e. for instance the output of $\tilde{B}_E$), it generates a vector $\tilde{z}_R$ s.t.:
$\tilde{z}_R[k] = \left\{\begin{matrix} x_k \mbox{ if } k\in E \\\tilde{z}[k] \mbox{ otherwise} \end{matrix}\right.$
8. Verify (numerically on an example) that $D_E$ is a projector. You can use the function `numpy.testing.assert_array_almost_equal` to check that arrays are almost equal (with a good precision)
Let us now define an operator that will work almost as the normalized Fourier transform, except that it will be applied to the extended signals and preserve the $x_E$ part. This can be easily done using the following block matrix:
Using this definition, multiplying an extended signal $\tilde{x}$ by $\tilde{W}$ will result in a vector containing the Fourier transform of $x$ followed by $x_E$, preserving the known initial values.
11. Recalling the low-pass matrix $P$ defined previously, build $\tilde{P}$ s.t. when applied to $\tilde{W}\tilde{x}$ it results in a vector containing the filtered values (still in the Fourier domain) followed by $x_E$.
13. We can now use all defined above to implement a function that performs an iteration, taking as input the extension of $x$ as defined above, that does the following:
- compute the filtered version of the signal using $\tilde{W}^H \tilde{P}\tilde{W}$ (i.e. projecting on the space of "smooth signals")
- restore the known values in the signal by using $D_E = C_E\tilde{B}_E$ (i.e project back on the space of signals having the known initial values we have set)
Make sure to force the result to real values by using `np.real` before returning
Depending on $M$, you will find that the method can reconstruct perfectly the signal or not.
%% Cell type:markdown id:94962762 tags:
%% Cell type:markdown id:ccf83efe tags:
<!-- END QUESTION -->
%% Cell type:markdown id:738294c5 tags:
## Submission
Make sure you have run all cells in your notebook in order before running the cell below, so that all images/graphs appear in the output. The cell below will generate a zip file for you to submit. **Please save before exporting!**
%% Cell type:code id:7d10d7ed tags:
``` python
# Save your notebook first, then run this cell to export your submission.
