Skip to content
Snippets Groups Projects
Commit f8ae4042 authored by Nicolas Aspert's avatar Nicolas Aspert
Browse files

exercises week3

parent 8164b87b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
exercises/week03/Projections and signal restoration.ipynb 0 → 100644
+ 1089
− 0
View file @ f8ae4042
%% Cell type:code id:9692b307 tags:
``` python
# Initialize Otter
import otter
grader = otter.Notebook("Projections and signal restoration.ipynb")
```
%% Cell type:markdown id:4133723e-952b-4b1c-bb3d-a7f78a957024 tags:
# Matrix Analysis 2025 - EE312
## Week 3 - Signal restoration using projections
[LTS2](https://lts2.epfl.ch)
%% Cell type:markdown id:fb42154f-8832-4960-9fee-3372b473ef69 tags:
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.
%% Cell type:code id:d238cbc9-055f-4d55-a6ae-29d83881e2f5 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
N = 100
k = np.arange(0, N)
w1 = 3
w2 = 7
w3 = 12
# generate simple signals
x1 = np.cos(2*w1*np.pi*k/N) + np.sin(2*w2*np.pi*k/N)
x2 = np.sin(2*w1*np.pi*k/N) + np.cos(2*w3*np.pi*k/N)
```
%% Cell type:code id:59a703e8-f351-4a02-8278-12560a67b023 tags:
``` python
# visualize the base signals
plt.plot(x1)
plt.plot(x2, 'r')
```
%% Cell type:markdown id:09249ca9-e441-4537-9bc5-018be67bc099 tags:
1. Implement a function that returns the normalized Fourier matrix of size N (use your implementation from week 2)
%% Cell type:code id:f886d8a6-4cb7-4187-ae4b-4615a06caa2d tags:otter_answer_cell
``` python
def fourier_matrix(N):
...
```
%% Cell type:code id:e101a935 tags:
``` python
grader.check("q1")
```
%% Cell type:code id:62fae5a8-f3c4-477b-a53f-a0106730b842 tags:
``` python
W_hat = fourier_matrix(N)
```
%% Cell type:markdown id:1e0a2337-446d-49d3-8521-55d07a1539bc tags:
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.
%% Cell type:code id:aeed4067-1130-4e6c-9fbd-ab8d0a3a41c0 tags:otter_answer_cell
``` python
def lowpass_matrix(N, w_c):
"""
Computes the P matrix that will remove high-frequency coefficients in a DFT transformed signal
Parameters
----------
N : length of the input signal
w_c : cutoff frequency
Returns
-------
The P matrix
"""
...
```
%% Cell type:code id:d1f7f40c tags:
``` python
grader.check("q2")
```
%% Cell type:markdown id:fef68750-d97f-4625-bd5e-5f13e7894f7d tags:
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$.
%% Cell type:code id:bc1370d5-c7b2-4598-b573-d57bec56aa1b tags:otter_answer_cell
``` python
w_c = ...
```
%% Cell type:code id:85787e70-1427-4860-8a7b-94d20432bd76 tags:otter_answer_cell
``` python
# Compute the DFT of x1 and x2
X1 = W_hat@x1
X2 = W_hat@x2
# Get the lowpass filter matrix
P = lowpass_matrix(N, w_c)
# Filter X1 and X2 and apply inverse DFT
x1_f = ...
x2_f = ...
```
%% Cell type:markdown id:700eb2e2-3137-4916-a829-9e8e2ba8c6e5 tags:
NB: Make sure the filtered output is **real** (or its imaginary part should be smaller than $10^{-12}$).
%% Cell type:code id:e960728d-c2c0-4b9d-ab73-3db971f559a1 tags:
``` python
# Plot the filtered signals
# x1_f should still contain both frequencies, x2_f only one
plt.plot(x1_f)
plt.plot(x2_f, 'r')
```
%% Cell type:code id:dd892e03 tags:
``` python
grader.check("q3")
```
%% Cell type:markdown id:c34b8fc7-8a51-4fe5-a6ab-1cc27da2cd1e tags:
<!-- BEGIN QUESTION -->
4. Prove that $P$ is a projection.
%% Cell type:markdown id:610729bb tags:otter_answer_cell
_Type your answer here, replacing this text._
%% Cell type:markdown id:7f17f3ae-6a39-4e86-a29c-879c96b4e5fb tags:
<!-- END QUESTION -->
### 2. Signal extension
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.
%% Cell type:markdown id:f9434f5f-ceb9-42bf-9b5e-bd7d7a94dae8 tags:
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.
%% Cell type:code id:dae1c4aa-c6ce-4f35-b309-aa6dd09e9abf tags:otter_answer_cell
``` python
def Bt_E(N, E):
"""
Computes the $\tilde{B}_E$ matrix
Parameters
----------
N : length of the input signal
E : list of suitable indices. e.g. for N=5 a valid E could be [1, 3]
Returns
-------
The $\tilde{B}_E$ matrix
"""
...
```
%% Cell type:code id:2dddfff8 tags:
``` python
grader.check("q5")
```
%% Cell type:markdown id:853c1a3b-9ed9-461c-86f9-12992e07337d tags:
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.$
6. Write a function that generates $C_E$.
%% Cell type:code id:932ae153-f336-4b0f-b7c0-774c02ccaad1 tags:otter_answer_cell
``` python
def C_E(N, E):
"""
Computes the $C_E$ matrix
Parameters
----------
N : length of the input signal
E : list of suitable indices. e.g. for N=5 a valid E could be [1, 3]
Returns
-------
The $C_E$ matrix
"""
...
```
%% Cell type:code id:577b22ec tags:
``` python
grader.check("q6")
```
%% Cell type:markdown id:0e57be3f-1c44-445a-86ad-07e20fbcd1d2 tags:
<!-- BEGIN QUESTION -->
7. What is the effect of $D_E = C_E \tilde{B}_E$ on an extended signal $\tilde{x}$ ?
%% Cell type:markdown id:273c0bc5 tags:otter_answer_cell
_Type your answer here, replacing this text._
%% Cell type:markdown id:3b2caf03-41cb-41f9-93eb-778047b9b944 tags:
<!-- END QUESTION -->
<!-- BEGIN QUESTION -->
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)
%% Cell type:code id:94dbc7f8-c2fc-4337-9828-dd2e3311830d tags:otter_answer_cell
``` python
# Set some parameters
N=5
E=[1,3]
# Compute D_E
D_E = ...
# Now check that D_E is a projector
...
```
%% Cell type:markdown id:7fb30993-5350-4e9b-aa45-b657b450a3a1 tags:
<!-- END QUESTION -->
### 3. Extended signals in the Fourier domain
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:
$\tilde{W} = \begin{pmatrix}\hat{W} & 0 \\0 & I_M\end{pmatrix}$.
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.
%% Cell type:markdown id:3157e22e-f7d7-44b8-a4fd-9ddec1d7aa8a tags:
<!-- BEGIN QUESTION -->
9. Prove that $\tilde{W}$ is orthonormal.
%% Cell type:markdown id:0fd7c9a7 tags:otter_answer_cell
_Type your answer here, replacing this text._
%% Cell type:markdown id:f140c720-eaed-4fd2-8b9f-6b99a1a5b0ff tags:
<!-- END QUESTION -->
10. Write a function that returns $\tilde{W}$.
%% Cell type:code id:b744a698-40d5-41a6-9623-76a19f50bcd2 tags:otter_answer_cell
``` python
def Wt_E(N, E):
"""
Computes the $\tilde{W}_E$ matrix
Parameters
----------
N : length of the input signal
E : list of suitable indices. e.g. for N=5 a valid E could be [1, 3]
Returns
-------
The $\tilde{W}_E$ matrix
"""
...
```
%% Cell type:code id:94532bbb tags:
``` python
grader.check("q10")
```
%% Cell type:markdown id:1afee45d-78e6-498b-905b-cbc97eeaf2d5 tags:
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$.
%% Cell type:code id:4576bfbc-a249-4d89-919d-47763c8a60f9 tags:otter_answer_cell
``` python
def Pt_E(N, M, w_c):
...
```
%% Cell type:markdown id:781d7752-74ae-4ec2-8ea2-4279b54ee5e7 tags:
<!-- BEGIN QUESTION -->
12. A signal $\tilde{x}$ will be filtered by doing $\tilde{W}^H \tilde{P}\tilde{W}\tilde{x}$.
Prove that this is a projection.
%% Cell type:markdown id:53c720a8 tags:otter_answer_cell
_Type your answer here, replacing this text._
%% Cell type:markdown id:1c8c4469-7264-4968-b700-897be635a193 tags:
<!-- END QUESTION -->
### 4. Signal restoration
%% Cell type:markdown id:ac24806d-105d-4fbf-bc98-788dbcc7c0ba tags:
<!-- BEGIN QUESTION -->
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
%% Cell type:code id:155e3277-bcef-44ed-b53f-a1b2e5ad7363 tags:otter_answer_cell
``` python
def restoration_iter(Wt, Pt, Dt, xt):
"""
Performs a full restoration iteration by
- projecting the signal into Fourier, performing low-pass filtering followed by inverse DFT
- restoring the known initial values into the filtered signal
Parameters
----------
Wt : \tilde{W} matrix
Pt : \tilde{P} matrix
Dt : \tilde{D} matrix
xt : \tilde{x} vector
Returns
-------
The new $\tilde{x} vector after the iteration
"""
...
```
%% Cell type:markdown id:758aa9ce-1ebc-49e1-b2f6-ca8e1035b31c tags:
<!-- END QUESTION -->
<!-- BEGIN QUESTION -->
15. Finally we are ready to apply all this to an example.
%% Cell type:code id:b467649e-2a46-40f2-bc98-73c354fb0484 tags:
``` python
# Setup an example
N = 100
w_c = 10 # cut-off
w1 = 3
w2 = 7
E = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95])
# E = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90]) # try with less known points
M = len(E)
x = np.cos(2*w1*np.pi*np.arange(0, N)/N) + np.sin(2*w2*np.pi*np.arange(0,N)/N) # original signal
# Create a signal that is only noise, except at known initial values
y = np.random.rand(N) #
y[E] = x[E] # restore known values
xe = x[E]
```
%% Cell type:code id:510b6f8f-b1c2-48b7-af00-eee8e55b5868 tags:
``` python
plt.plot(y) # plot the "acquired" signal
plt.plot(x, 'r+-', alpha=0.5) # plot the ground truth signal
```
%% Cell type:markdown id:b71f1ff5-73d6-4bc1-ba06-1452c8bf8adb tags:
Generate $\tilde{y}$ (this only need to be done once)
%% Cell type:code id:4e0943bf-0f35-4837-b6db-1615ba466ae5 tags:
``` python
yt = np.append(y, xe) # SOLUTION
```
%% Cell type:markdown id:c25239d9-3015-4a33-b0c0-0c906ee4129e tags:
Run a number of iterations of the method and plot the result:
%% Cell type:code id:6aaf1a71-6c2a-4f43-a51e-82b40cc6d6d0 tags:
``` python
def signal_restore(Wt, Pt, Dt, yt, niter=20):
yr0 = yt # initialize
for k in range(niter):
yr1 = restoration_iter(Wt, Pt, Dt, yr0)
yr0 = yr1
return yr1
```
%% Cell type:code id:951418c2-ee12-4e8d-8f58-cba29d26dd7e tags:
``` python
Wt = Wt_E(N, E)
Dt = C_E(N, E)@Bt_E(N, E)
Pt = Pt_E(N, M, w_c)
yr = signal_restore(Wt, Pt, Dt, yt)
```
%% Cell type:code id:43b3d939-eeb2-4171-9293-0076f1c06c83 tags:
``` python
plt.plot(yr[:N]) # plot reconstructed signal
plt.plot(x, 'r+', alpha=0.7) # plot original for comparison
```
%% Cell type:markdown id:bf876405-b60c-4f7c-9162-56d58cafd19f tags:
Depending on $M$, you will find that the method can reconstruct perfectly the signal or not.
%% Cell type:markdown id:94962762 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.
grader.export(pdf=False, run_tests=True)
```
%% Cell type:markdown id:159a45fb tags:
exercises/week03/images/pocs.png 0 → 100644
+ 0
− 0
View file @ f8ae4042
exercises/week03/images/pocs.png

76.9 KiB

0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment