Sampling Signals with Finite Python-Mini-Project Sampling Si

Python-Mini-Project Sampling Signals with Finite Rate of Innovation with an Application to Image Super-Resolution Communications and Signal Processing Group Department of Electrical and Electronic Engineering Imperial College London The coursework consists of a Python Miniproject. It counts 25% of the ffnal mark. Students should produce a written ffnal report along with code which is to be submitted separately. We expect the code to be neatly written and well documented. 1 Submission details and Guidelines The deliverable for this course consists of two parts: • The ffnal written report must be submitted through blackboard by the deadline indicated on the module page • The code submission, consisting of a single python script called submission.py, must be submitted to a separate submission page on blackboard. A blank submission.py already exists in the project folder. • The project ffles including skeleton code and data can be found on blackboard. 1.1 Report submission The report should give full answers to each question including plots and derivations, however, details of code implementation should be conffned to the code submission and not the report. 1.2 Code submission Code submissions for this coursework must be written in Python. Python has been chosen over MATLAB since it has more widespread use outside of academic circles. Furthermore, it is free to use without a license. Code should be clear and well documented. In cases where code produces the wrong result and is also convoluted or hard to understand, no marks will be given. 1.2.1 Allowed Python modules Python scripts must utilise only the following modules:

PyWavelets for generating wavelets
numpy for general matrix and vector manipulation
matplotlib for plotting 1The following libraries will also be used by the image registration skeleton code: • Pillow • scipy • scikit-image 1.2.2 Submission format Each question that requires a code submission will specify a prototype of a function (i.e. the name, inputs, and outputs) that will be required to answer the question. For each input and output the type will also be given. For example: name dummy_exercise arg1 1D Array of signal samples : numpy.ndarray arg2 Sample rate : float output 1D array of innovation times : numpy.ndarray Functions that do not adhere to these prototypes will lose a portion of the marks allocated for code correctness. In total, the code submission for this coursework should consist of a single python script called ‘submission.py’ containing function deffnitions for each of the requested functions. 1.2.3 Project ffles contents The supplied project folder contains a number of ffles listed here for reference:
submission.py - The empty python script to be fflled with requested functions
data.py - Python module that contains some data required for the exercises
image fusion.py - Python module containing functions and variables needed for the ffnal exercise on super resolution
Pipffle and Pipffle.lock - Pipenv ffles. These can be used for setting a python environment with all the 代写Sampling Signals with Finite appropriate modules installed. See Pipenv for more details.
project ffles & data - Folder containing data for the project. It is not necessary to access this folder manually. A Python “pipenv” is provided in the project repo should students wish to use them. Use of any of this helper tool is not required but provided for the convenience of students. 1.2.4 Supplementary code Any code not explicitly asked for in the code submission may be placed in an appendix of the main report if it is a useful aid to understanding your report. You should avoid placing all your code in the appendices. 22 Preliminaries The classic sampling theorem states that any bandlimited function x(t) such that X(ω) = 0 ∀ |ω| > ωmax can be exactly recovered from its samples given that the rate 2π/T is greater or equal to twice the highest frequency component. The continuous time signal is recovered with x(t) = P n∈Z y[n]sinc(t/T − n) where sinc(t) = sin(πt)/πt and y[n] = x(nT). A fun extension of this theorem has recently been developed, where it is made possible to sample signals that are not bandlimited but completely determined by a ffnite number of parameters. Such signals are called Finite Rate of Innovation (FRI) signals. In particular, it is made possible to sample and exactly recover streams of Diracs, differentiated Diracs and piecewise polynomials using a smoothing kernel φ(t) that satisffes the Strang-Fix conditions. In this case, the samples are given by y[n] = ⟨x(t), φ(t − n)⟩. The main aim of this project is to make the student familiar with these recently developed sampling theorems. Notice that this project refers mainly to Chapter 5 of the lecture notes. 33 Strang-Fix conditions One of the main conditions on the sampling kernel φ(t) is that it should be able to reproduce polynomials. Therefore, there exist coefffcients cm,n such that X n∈Z cm,nφ(t − n) = t m m = 0, 1, . . . , N. (1) In our context, we assume that {φ(t − n)}n∈Z spans a Riesz space. Therefore there exists a dual-basis {φ˜(t − n)}n∈Z that satisffes ⟨φ(t), φ˜(t − n)⟩ = δ0. The coefffcients cm,n can thus be found by applying the inner product with the dual-basis φ˜(t − n) on both sides of equation (1). It results that cm,n = ⟨t m, φ˜(t−n)⟩. Notice that a scaling function that generates wavelets with N + 1 vanishing moments is able to reproduce polynomials of degree N. It can therefore be used in our case as a sampling kernel. Exercise 1. The Daubechies scaling function a. Which Daubechies scaling function should be used in order to reproduce polynomials of maximum degree 3? b. Write a function that calculates coefffcients to reproduce polynomials of degree 0 to 3: name get_cmn_db arg1 shift number : Int arg2 polynomial order : Int output coefficient : float c. Compute the coefffcients to reproduce polynomials of degree 0 to 3 for n = 0, . . . , 32−L where L is the support of the kernel. Plot the results with the shifted kernels and the reproduced polynomials as in Figure 5.2 of the lecture notes. Guideline: • Choose the right Daubechies scaling function. For example, assume it is ‘dB2’. • To generate the φ(t) with a resolution 1 2 J = 1 64 , do as follows:

import pywt db2 = pywt.Wavelet(’db2’) phi, psi, x = db2.wavefun(level=6) phi = np.pad(phi,(0,1984)) • The array phi is in fact a 6-iteration approximation of φ(t) padded with zeros to a length of 2048. The shifted version φ(t − nT) where T = 1 is obtained by shifting phi by 64n points. To do this you can either reassign with index slices or by using np.roll, though be aware that this rotates the array like circshift in MATLAB, not simply shifting it. • Note that a period of T = 1 is equivalent to 64 indices due to the chosen resolution. 4Exercise 2. The B-Spline scaling function a. Which B-Spline scaling function should be used in order to reproduce polynomials of maximum degree 3? b. Find the dual-basis of the scaling function and write a function that calculates coefffcients to reproduce polynomials of degree 0 to 3: name get_cmn_bspline arg1 shift number : Int arg2 polynomial order : Int output coefficient : float c. Compute the coefffcients to reproduce polynomials of degree 0 to 3 for n = 0, . . . , 32−L where L is the support of the kernel. Plot the results with the shifted kernels and the reproduced polynomials as in the previous exercise. Hint: • Use the Daubechies formula for spectral factorization in order to ffnd the dual-basis. 4 The annihilating fflter method Assume we observe a discrete signal τ [m] of the form τ [m] = PK−1 k=0 akt m k and assume we want to retrieve the locations tk and the weights ak. One efffcient way of determining these parameters is through the annihilating fflter. Consider the discrete fflter h[m] given by z-transform H(z) = K Y−1 k=0 (1 − tkz −1 ). (2) Clearly, the convolution between h and τ will result in a zero output as the tks are the roots of the fflter (for more details, we refer to Section 5.2.1 of the lecture notes). Therefore h is called annihilating fflter of τ . In practice, the fflter can be found by posing h[0] = 1 and solving the following system of equations The knowledge of the fflter is sufffcient to determine the tks as they are the roots of the polynomial in (2). The weights ak can be determined using K samples of τ [m]. These form a classic Vandermonde system which has unique solution given that the tks are distinct. 5Exercise 3. The annihilating filter a. Write a function that finds the annihilating filter h[m] for a given set of signal moments τ [m]: name annihilating_filter arg1 1D array of signal moments : np.ndarray arg2 filter order, K : Int output annihilating filter, h[m] : np.ndarray You may assume that there are always at least 2K signal moments b. Write a function that determines the innovation parameters ak and tk for k = 1, . . . , K from the annihilating filter generated by the function in part (a): name innovation_params arg1 annihilating filter, h[m] : np.ndarray arg2 filter order, K : Int output 2D array, a_k along first row, t_k along second row : np.ndarray c. Apply your functions to the signal moments, τ [m], provided in the project folder. Retrieve the locations tk and the amplitudes ak. Remark: we have set K = 2. You can load moments as a numpy array from the data module using data.tau. 5 Sampling Diracs All the necessary pieces for implementing an end-to-end algorithm capable of recovering sampled streams of Diracs have been written. Indeed, consider a signal of type x(t) = PK−1 k=0 akδ(t−tk). Then the locations tk and the amplitudes ak can be recovered from the samples by applying the annihilating filter method to where the last equality derives from the fact that x(t) is a stream of Diracs and R ∞ −∞ f(t)δ(t − tk)dt = f(tk). Exercise 4. Sampling Diracs a. Create a stream of Diracs where K = 2. Filter the Diracs with an appropriate Daubechies scaling function, then sample the resulting signal. For the scaling function, using a sampling period of T = 1. Plot the resulting samples. b. Reconstruct the continuous time Diracs from the samples. Compare this result with the ground truth. 6Exercise 5. Reconstruction We have sampled a stream of Diracs (K = 2, T = 1) with a ‘db4’ Daubechies scaling function. The observed samples are available from the data module as data.samples. a. Use the reconstruction algorithm from the previous question to exactly recover the locations and the amplitudes of all the Diracs. 6 Reconstruction in the presence of noise The methods developed so far are not resilient to noise, so we want to use total-least squares and the Cadzow routine to make them more resilient to noise. Exercise 6. Sampling in noise a. Using the routine you developed in Exercise 4, create a stream of Diracs where K = 2. Sample the signal with an appropriate Daubechies scaling function that can reproduce polynomial up to degree N − 1 and N > 2K using a sampling period of T = 1. Compute the moments s[m] = P n cm,ny[n] with m = 0, 1, .., N − 1 and N > 2K. Add zero-mean Gaussian noise to the moments s[m] leading to sˆ[m] = s[m] + ϵ[m]. Use the command np.random.randn to generate Gaussian noise. Choose different values for the variance σ 2 . Guideline: • To generate zero-mean Gaussian noise you can use the numpy randn function: import numpy as np noise = np.random.randn(N) where N is the number of samples of noise Exercise 7. Reconstruction in noise

Write a function that performs reconstruction on a set of signal moments using TLS denoising: name tls arg1 1D array of signal moments : np.ndarray arg2 Number of Diracs, K : Int output Dirac times and magnitudes : np.ndarray * np.ndarray
Write a function that performs reconstruction on a set of signal moments using Cadzow denoising: name cadzow arg1 1D array of signal moments : np.ndarray arg2 Number of Diracs, K : Int arg3 Number of iterations : Int output Dirac times and magnitudes : np.ndarray * np.ndarray
Given the noisy moments generated in the previous exercise, use the TLS and Cadzow methods to retrieve the original two Diracs. Compare and contrast your results. Try different values of N. 77 Application: image super-resolution You have just implemented a new sampling algorithm capable of perfectly recovering certain classes of signals that have been sampled with a rate below the Nyquist rate. This new scheme has a wide variety of applications ranging from wideband communications to image superresolution. It is this last application that we are now asking you to implement. Suppose that we have access to N low-resolution cameras that are all pointing to the same fixed scene. The positions of each camera is unknown. The goal of image super-resolution is to generate a single high-resolution image, called the super-resolved image, by fusing the set of low-resolution images acquired by these cameras. The super-resolved image (SR) has therefore more pixels than any of the low-resolution image (LR), but more importantly, it has also more details (i.e. it is not a simple interpolation). This is made possible because each camera observes the scene from a different viewpoint. Image super-resolution can usually be decomposed into two main steps (see Figure 1):
Image registration: the geometric transformations between the LR images are estimated so that they can be overlaid accurately.
Image fusion and restoration: the data from each LR image is fused on a single highresolution grid and missing data is estimated. Restoration finally enhances the resulting image by removing any blur and noise. In the next exercise, you are asked to write in Matlab the image registration function of an image super-resolution algorithm. The approach used to register the LR images is based on the notion of continuous moments and is explained next. The fusion/restoration step is already written for you. Suppose that the i-th camera observes the scene: fi (x, y) = f1 (x + dxi , y + dyi) i = 1, . . . , N, where f1 (x, y) is the scene observed by the first camera taken as the reference. We therefore have (dx1, dy1) T = (0, 0)T . We have implicitly assumed here that the geometric transformations between the views are 2-D translations. At each camera, the point spread function φ (x, y) of the lens is modeled by a 2-D cubic B-spline and blurs the observed scene fi (x, y) (see Figure 2). B-splines are functions with compact support and satisfy Strang-Fix conditions. Therefore, similarly to Equation (1), there exists a set of coefficients n The blurred scene fi (x, y) ∗φ (−x, −y) is then uniformly sampled by the CCD array which acts as an analog to discrete converter. The output discrete image of the i-th camera has then the samples S (i) m,n: S (i) m,n = ⟨fi (x, y), φ (x − m, y − n)⟩ By using a similar development as written in Section 5, it can be shown that: mp,q = Z Z f (x, y) x m,n Sm,n (4) Equation (3) is the definition of the geometric moments mp,q of order (p + q), p, q ∈ N, of a 2-D continuous function f (x, y). Equation (4) shows that it is possible to retrieve the exact 8moments of the observed continuous scene from its blurred and sampled version by computing a linear combination with the proper coefficients. Knowing those moments, we can compute the barycenter of f (x, y): can be retrieved from the difference between the barycenter of f1 and the barycenter of fi Exercise 8. Image super-resolution In this exercise, you have access to N = 40 low-resolution color images of 64x64 pixels. The goal is to generate a super-resolved color image of 512x512 pixels. The images observe the same scene and have been taken from different viewpoints. We assume that the geometric transformations between the views are 2-D translations. a. Write a function that registers each image with respect to the first image. This function should return the estimated shift of each image: name get_barycenters output barycenters_x and barycenters_y relative : np.ndarray * np.ndarray to first low resolution image Note: the function has no arguments. barycenterx[k, l] should be the barycenter in the x axis for the lth layer of the kth low resolution image relative to the 1st low resolution image. Thus barycenterx[0, l] = 0. As stated in the guidelines, you should threshold each low resolution image with an appropriate threshold to increase SNR (use pywt.threshold with mode=’hard’). You can see the results of your registration by running image fusion.superres(get barycenters) b. What is the correct threshold for pixel values and how did you determine this? Use a principled way of determining the correct threshold. What is the PSNR achieved with your implementation? c. Can you pre/postprocess the data to achieve a higher PSNR? You must submit the original implementation of get barycenters but you can discuss any improvement you make in the written report Guidelines and tips:
The main python module is image fusion.py. This module handles loading relevant data, image fusion, and plotting. You can edit this file if you want but all submissions will be tested against the original.
Functions written in your submission script can use variables and methods from image - fusion.py by accessing them with the following convention image fusion.x where x is a variable in image fusion or image fusion.f() where f is a function in image fusion
The variables and functions from image fusion of interest to you are: 9name Purpose n sensors Number of low resolution images n layers Number of layers per image (RGB) lr ims Low resolution images in a list, each image has the layer axis last hr ref im High resolution reference image, layer axis last lr size Width/height of low resolution images hr size Width/height of high resolution reference image and target resolution coeff 0 0 Three matrices of size 64x64: You will use these to compute the continuous moments from the samples of your LR images. Coef 0 0 is used to compute the continuous moment m0,0, Coef 1 0 is for m1,0 and Coef 0 1 is for m0,1. coeff 0 1 coeff 1 0 superres(get barycenters) A function that performs image fusion. To test your barycenter finding function, pass it to superres and superres will generate the corresponding reconstruction and plot it. The PSNR and the super resolved image are returned from this function. These can be accessed using dot notation. For example, if you want to use the superres function you can call it like so: image fusion.superres(get barycenters)
Consider first the Red layer of each LR images, compute the continuous moments and find the 2-D translations. Then only, consider the Green layer, and then the Blue one. For each LR image, you should then have retrieved three similar 2-D translations.
In order to reduce the noise created by the background of the image and get accurate moments, it is required that you select a threshold for each layer and set to 0 all the samples that are below it.
All images are scaled such that sample values are between 0.0 and 1.0 10Restoration Super-resolved image Set of low-resolution images Image Registration HR grid estimation LR image 0 ... LR image k Figure 1: Image super-resolution principles S (i) fi(x, y) T m,n φ(x, y) (a) Acquisition model. (b) 2-D B-spline φ (x, y). Figure 2: Model for each of the N cameras.

WX：codinghelp