Deconvolution/1D

Deconvolution/1D
You are encouraged to solve this task according to the task description, using any language you may know.

The convolution of two functions $F$ and $H$ of an integer variable is defined as the function $G$ satisfying

$G(n)=\sum _{m=-\infty }^{\infty }F(m)H(m-n)$

for all integers $n$. Assume $F(n)$ can be non-zero only for $0$ ≤ $n$ < $|F|$, where $|F|$ is the "length" of $F$, and similarly for $G$ and $H$, so that the functions can be modeled as finite sequences by identifying $f_{0},f_{1},f_{2}\dots$ with $F(0),F(1),F(2)\dots$, etc.. Then for example, values of $|F|=6$ and $|H|=5$ would determine the following value of $g$ by definition.

${\begin{array}{lllllllllll}g_{0}&=&f_{0}h_{0}\\g_{1}&=&f_{1}h_{0}&+&f_{0}h_{1}\\g_{2}&=&f_{2}h_{0}&+&f_{1}h_{1}&+&f_{0}h_{2}\\g_{3}&=&f_{3}h_{0}&+&f_{2}h_{1}&+&f_{1}h_{2}&+&f_{0}h_{3}\\g_{4}&=&f_{4}h_{0}&+&f_{3}h_{1}&+&f_{2}h_{2}&+&f_{1}h_{3}&+&f_{0}h_{4}\\g_{5}&=&f_{5}h_{0}&+&f_{4}h_{1}&+&f_{3}h_{2}&+&f_{2}h_{3}&+&f_{1}h_{4}\\g_{6}&=&&&f_{5}h_{1}&+&f_{4}h_{2}&+&f_{3}h_{3}&+&f_{2}h_{4}\\g_{7}&=&&&&&f_{5}h_{2}&+&f_{4}h_{3}&+&f_{3}h_{4}\\g_{8}&=&&&&&&&f_{5}h_{3}&+&f_{4}h_{4}\\g_{9}&=&&&&&&&&&f_{5}h_{4}\end{array}}$

For this task, implement a function (or method, procedure, subroutine, etc.) deconv to perform deconvolution (i.e., the inverse of convolution) by constructing and solving such a system of equations for $h$ given $f$ and $g$.

• The function should work for $G$ of arbitrary length (i.e., not hard coded or constant) and $F$ of any length up to that of $G$. Note that $|H|$ will be given by $|G|-|F|+1$.
• There may be more equations than unknowns. If convenient, use a function from a library that finds the best fitting solution to an overdetermined system of linear equations (as in the Multiple regression task). Otherwise, prune the set of equations as needed and solve as in the Reduced row echelon form task.
• Test your solution on the following data. Be sure to verify both that deconv$(g,f)=h$ and deconv$(g,h)=f$ and display the results in a human readable form.

h = [-8,-9,-3,-1,-6,7]
f = [-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1]
g = [24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7]

Ursala

The user defined function band constructs the required matrix as a list of lists given the pair of sequences to be deconvolved, and the lapack..dgelsd function solves the system. Some other library functions used are zipt (zipping two unequal length lists by truncating the longer one) zipp0 (zipping unequal length lists by padding the shorter with zeros) and pad0 (making a list of lists all the same length by appending zeros to the short ones).

<lang Ursala>#import std

1. import nat

band = pad0+ ~&rSS+ zipt^*D(~&r,^lrrSPT/~&ltK33tx zipt^/~&r ~&lSNyCK33+ zipp0)^/~&rx ~&B->NlNSPC ~&bt

deconv = lapack..dgelsd^\~&l ~&||0.!**+ band </lang> test program: <lang Ursala>h = <-8.,-9.,-3.,-1.,-6.,7.> f = <-3.,-6.,-1.,8.,-6.,3.,-1.,-9.,-9.,3.,-2.,5.,2.,-2.,-7.,-1.> g = <24.,75.,71.,-34.,3.,22.,-45.,23.,245.,25.,52.,25.,-67.,-96.,96.,31.,55.,36.,29.,-43.,-7.>

1. cast %eLm

test =

<

  'h': deconv(g,f),
  'f': deconv(g,h)>

</lang> output:

<
   'h': <
      -8.000000e+00,
      -9.000000e+00,
      -3.000000e+00,
      -1.000000e+00,
      -6.000000e+00,
      7.000000e+00>,
   'f': <
      -3.000000e+00,
      -6.000000e+00,
      -1.000000e+00,
      8.000000e+00,
      -6.000000e+00,
      3.000000e+00,
      -1.000000e+00,
      -9.000000e+00,
      -9.000000e+00,
      3.000000e+00,
      -2.000000e+00,
      5.000000e+00,
      2.000000e+00,
      -2.000000e+00,
      -7.000000e+00,
      -1.000000e+00>>