Towards Generalized FRI Sampling with an Application to Source Resolution in Radioastronomy
BC 322 (Bâtiment BC)
It is a classic problem to estimate continuous time sparse signals, like point sources in a direction-of-arrival problem, or pulses in a time-of-flight measurement. The earliest occurrence is the estimation of sinusoids in time series using Prony’s method. This is at the root of a substantial line of work on high resolution spectral estimation. The estimation of continuous-time sparse signals from discretetime samples is the goal of the sampling theory for finite rate of innovation (FRI) signals. Both spectral estimation and FRI sampling usually assume uniform sampling.
But not all measurements are obtained uniformly, as exemplified by a concrete radioastronomy problem we set out to solve. Thus, we develop the theory and algorithm to reconstruct sparse signals, typically sum of sinusoids, from non-uniform samples. We achieve this by identifying a linear transformation that relates the unknown uniform samples of sinusoids to the given measurements. These uniform samples are known to satisfy the annihilation equations. A valid solution is then obtained by solving a constrained minimization such that the reconstructed signal is consistent with the given measurements and satisfies the annihilation constraint.
Thanks to this new approach, we unify a variety of FRI based methods. We demonstrate the versatility and robustness of the proposed approach with five FRI reconstruction problems, namely Dirac reconstructions with irregular time or Fourier domain samples, FRI curve reconstructions, Dirac reconstructions on the sphere and point source reconstructions in radioastronomy. The proposed algorithm improves substantially over state of the art methods and is able to reconstruct point sources accurately from irregularly sampled Fourier measurements under severe noise conditions.
This repository contains all the code to reproduce the results of the paper Towards Generalized FRI Sampling with an Application to Source Resolution in Radioastronomy. It contains a python implementation of the proposed algorithm.
A number of scripts were written to exemplify the application of the proposed finite rate of innovation (FRI) reconstruction algorithm on solving various FRI reconstruction problems, including
- Dirac reconstruction from non-uniformly sampled time domain samples (
- Dirac reconstruction from non-uniformly sampled Fourier domain samples (
- FRI curve reconstruction (
- Point source reconstruction in radioastronomy (
We are available for questions or requests relating to either the code or the theory behind it.
Recreate the figures
Start an ipython/jupyter qt-console
and then type in the following commands in an ipython shell.
# Start inline plotting mode %matplotlib inline # Plot the radioastronomy example in the introduction %run fig2_radio_ast_eg1.py # Dirac reconstruction from non-uniformly sampled # time domain measurements (noiseless and noisy cases) %run fig6a_irreg_time_samp.py %run fig6b_irreg_time_samp.py # Dirac reconstruction from non-uniformly sampled # Fourier domain measurements # 1) the spectrum satisfies the periodic assumption %run fig7_irreg_fourier_samp.py # 2) the spectrum DOES NOT satisfy the periodic assumption # (noiseless and noisy cases) %run fig8a_irreg_fourier_samp.py %run fig8b_irreg_fourier_samp.py # FRI curve reconstruction # 1) Monte Carlo simulation for different noise levels %run fig9_fri_curve.py # 2) A visual example for FRI curve reconstruction %run fig10_fri_curve.py # Point source reconstruction for radioastronomy %run fig11abc_radio_ast_eg2.py %run fig11d_radio_ast_eg2.py
Data used in the paper
The randomly generated signal parameters are saved under the folder
# Simulation with non-uniform sampled time domain measurements data/Dirac_Data_20-12_02_22.npz # Simulation with non-uniform sampled Fourier measurements # periodic spectrum data/freq_Dirac_Data_p20-12_11_48.npz # aperiodic spectrum data/freq_Dirac_Data_ap20-12_11_55.npz # Experiment with FRI curve data/coef.npz # Experiments for radioastronomy data/Dirac_Data_8-2_19_47.npz data/Dirac_Data_4-2_02_28.npz
Saved simulation results
The Monte Carlo simulation results are saved in folder
avg_perf_K_3_L_3_all.npz: FRI curve reconstruction with different levels of noise
radio_ast_batch_res1000.npz: Point source reconstruction (3 sources)
radio_ast_eg3_batch_res1000.npz: Point source reconstruction (5 sources)
Remove the saved results to trigger new Monte Carlo simulation.
Overview of results
We apply the proposed algorithmic framework for several FRI reconstruction problems.
Dirac reconstruction from time domain samples (Fig. 6)
We demonstrate how Diracs are reconstructed from non-uniformly sampled time domain measurements.
Dirac reconstruction from Fourier domain samples (Fig. 7 and Fig. 8)
We demonstrate how Diracs are reconstructed from non-uniformly sampled time domain measurements. We consider both cases where the spectrum is periodic and aperiodic. The difference is whether we have an exact mapping from a set of unknown uniform samples of sinusoids to the given non-uniformly sampled measurements.
FRI curve reconstruction (Fig. 9 and Fig. 10)
Visiual comparison for (SNR = 5dB)
|(a) Noisy Samples||(b) Cadzow's Method|
|(c) Structured Low-rank Approximation||(d) Proposed Approach|
Point source reconstruction in radioastronomy (Fig. 2 and Fig. 11)
We demonstrate one application of the proposed FRI-based reconstruction algorithm on reconstructing point sources in radioastronomy.
Resolve closely located sources (SNR = 5dB)
(a) Inverse Fourier Transform of the Raw Measurements (b) Statistics of Retrieved Sources using FRI (c) Zoom-in of (a) (d) Zoom-in of (b)
Comparison with classic approaches
(a) The Given Measurements (b) ℓ1 Minimization Result (c) Continuous Domain Reconstruction (FRI) (d) Statistics of the Retrieved Sources Locations (FRI)
- A working distribution of Python 3.
- Numpy, Scipy
- We use the distribution anaconda to simplify the setup of the environment.
- We use the MKL extension of Anaconda to speed things up. There is a free license for academics.
- We use joblib for parallel computations.
- matplotlib for plotting the results.
List of standard packages needed
numpy, scipy, matplotlib, mkl, joblib
|Machine||MacBook Pro Retina 15-inch, Mid 2014|
|System||macOS Sierra 10.12|
|CPU||Intel Core i7|
System Info: ------------ Darwin Kernel Version 16.0.0: Mon Aug 29 17:56:20 PDT 2016; root:xnu-3789.1.32~3/RELEASE_X86_64 x86_64 Python Info: ------------ Python 3.5.2 :: Anaconda custom (x86_64) Python Packages Info (conda) ---------------------------- # packages in environment at /Users/pan/anaconda: accelerate 2.3.0 np111py35_3 accelerate_cudalib 2.0 0 anaconda custom py35_0 anaconda-client 1.4.0 py35_0 anaconda-navigator 1.3.1 py35_0 joblib 0.9.4 py35_0 jupyter 1.0.0 py35_3 jupyter_client 4.3.0 py35_0 jupyter_console 4.1.1 py35_0 jupyter_core 4.1.0 py35_0 matplotlib 1.5.3 np111py35_0 mkl 11.3.3 0 mkl-service 1.1.2 py35_2 nb_anacondacloud 1.1.0 py35_0 numpy 1.11.1 py35_0 pyqt 5.6.0 py35_0 qt 5.6.0 0 qtawesome 0.3.3 py35_0 qtconsole 4.2.1 py35_1 qtpy 1.1.2 py35_0 scipy 0.18.1 np111py35_0
Copyright (c) 2016, Hanjie Pan, Thierry Blu and Martin Vetterli
The source code is released under the MIT license.