The convolution of two functions 





F




{\displaystyle {\mathit {F}}}

 and 





H




{\displaystyle {\mathit {H}}}

 of
an integer variable is defined as the function 





G




{\displaystyle {\mathit {G}}}


satisfying





G
(
n
)
=

∑

m
=
−
∞


∞


F
(
m
)
H
(
n
−
m
)


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


for all integers 





n




{\displaystyle {\mathit {n}}}

. Assume 



F
(
n
)


{\displaystyle F(n)}

 can be non-zero only for 



0


{\displaystyle 0}

 ≤ 





n




{\displaystyle {\mathit {n}}}

 ≤ 




|



F



|



{\displaystyle |{\mathit {F}}|}

, where 




|



F



|



{\displaystyle |{\mathit {F}}|}

 is the "length" of 





F




{\displaystyle {\mathit {F}}}

, and similarly for 





G




{\displaystyle {\mathit {G}}}

 and 





H




{\displaystyle {\mathit {H}}}

, so that the functions can be modeled as finite sequences by identifying 




f

0


,

f

1


,

f

2


,
…


{\displaystyle f_{0},f_{1},f_{2},\dots }

 with 



F
(
0
)
,
F
(
1
)
,
F
(
2
)
,
…


{\displaystyle F(0),F(1),F(2),\dots }

, etc.
Then for example, values of 




|



F



|

=
6


{\displaystyle |{\mathit {F}}|=6}

 and 




|



H



|

=
5


{\displaystyle |{\mathit {H}}|=5}

 would determine the following value of 





g




{\displaystyle {\mathit {g}}}

 by definition.










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








{\displaystyle {\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}}}


We can write this in matrix form as:






(





g

0







g

1







g

2







g

3







g

4







g

5







g

6







g

7







g

8







g

9






)

=

(





f

0







f

1





f

0







f

2





f

1





f

0







f

3





f

2





f

1





f

0







f

4





f

3





f

2





f

1





f

0







f

5





f

4





f

3





f

2





f

1








f

5





f

4





f

3





f

2









f

5





f

4





f

3










f

5





f

4











f

5






)



(





h

0







h

1







h

2







h

3







h

4






)



{\displaystyle \left({\begin{array}{l}g_{0}\\g_{1}\\g_{2}\\g_{3}\\g_{4}\\g_{5}\\g_{6}\\g_{7}\\g_{8}\\g_{9}\\\end{array}}\right)=\left({\begin{array}{lllll}f_{0}\\f_{1}&f_{0}\\f_{2}&f_{1}&f_{0}\\f_{3}&f_{2}&f_{1}&f_{0}\\f_{4}&f_{3}&f_{2}&f_{1}&f_{0}\\f_{5}&f_{4}&f_{3}&f_{2}&f_{1}\\&f_{5}&f_{4}&f_{3}&f_{2}\\&&f_{5}&f_{4}&f_{3}\\&&&f_{5}&f_{4}\\&&&&f_{5}\end{array}}\right)\;\left({\begin{array}{l}h_{0}\\h_{1}\\h_{2}\\h_{3}\\h_{4}\\\end{array}}\right)}


or





g
=
A

h


{\displaystyle g=A\;h}


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 represented by the above matrix 



A


{\displaystyle A}

 for 





h




{\displaystyle {\mathit {h}}}

 given 





f




{\displaystyle {\mathit {f}}}

 and 





g




{\displaystyle {\mathit {g}}}

.

 The function should work for 





G




{\displaystyle {\mathit {G}}}

 of arbitrary length (i.e., not hard coded or constant) and 





F




{\displaystyle {\mathit {F}}}

 of any length up to that of 





G




{\displaystyle {\mathit {G}}}

. Note that 




|



H



|



{\displaystyle |{\mathit {H}}|}

 will be given by 




|



G



|

−

|



F



|

+
1


{\displaystyle |{\mathit {G}}|-|{\mathit {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


{\displaystyle (g,f)=h}

 and deconv



(
g
,
h
)
=
f


{\displaystyle (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]


package require Tcl 8.5
namespace eval 1D {
    namespace ensemble create;   # Will be same name as namespace
    namespace export convolve deconvolve
    # Access core language math utility commands
    namespace path {::tcl::mathfunc ::tcl::mathop}
 
    # Utility for converting a matrix to Reduced Row Echelon Form
    # From http://rosettacode.org/wiki/Reduced_row_echelon_form#Tcl
    proc toRREF {m} {
	set lead 0
	set rows [llength $m]
	set cols [llength [lindex $m 0]]
	for {set r 0} {$r < $rows} {incr r} {
	    if {$cols <= $lead} {
		break
	    }
	    set i $r
	    while {[lindex $m $i $lead] == 0} {
		incr i
		if {$rows == $i} {
		    set i $r
		    incr lead
		    if {$cols == $lead} {
			# Tcl can't break out of nested loops
			return $m
		    }
		}
	    }
	    # swap rows i and r
	    foreach j [list $i $r] row [list [lindex $m $r] [lindex $m $i]] {
		lset m $j $row
	    }
	    # divide row r by m(r,lead)
	    set val [lindex $m $r $lead]
	    for {set j 0} {$j < $cols} {incr j} {
		lset m $r $j [/ [double [lindex $m $r $j]] $val]
	    }
 
	    for {set i 0} {$i < $rows} {incr i} {
		if {$i != $r} {
		    # subtract m(i,lead) multiplied by row r from row i
		    set val [lindex $m $i $lead]
		    for {set j 0} {$j < $cols} {incr j} {
			lset m $i $j \
			    [- [lindex $m $i $j] [* $val [lindex $m $r $j]]]
		    }
		}
	    }
	    incr lead
	}
	return $m
    }
 
    # How to apply a 1D convolution of two "functions"
    proc convolve {f h} {
	set g [lrepeat [+ [llength $f] [llength $h] -1] 0]
	set fi -1
	foreach fv $f {
	    incr fi
	    set hi -1
	    foreach hv $h {
		set gi [+ $fi [incr hi]]
		lset g $gi [+ [lindex $g $gi] [* $fv $hv]]
	    }
	}
	return $g
    }
 
    # How to apply a 1D deconvolution of two "functions"
    proc deconvolve {g f} {
	# Compute the length of the result vector
	set hlen [- [llength $g] [llength $f] -1]
 
	# Build a matrix of equations to solve
	set matrix {}
	set i -1
	foreach gv $g {
	    lappend matrix [list {*}[lrepeat $hlen 0] $gv]
	    set j [incr i]
	    foreach fv $f {
		if {$j < 0} {
		    break
		} elseif {$j < $hlen} {
		    lset matrix $i $j $fv
		}
		incr j -1
	    }
	}
 
	# Convert to RREF, solving the system of simultaneous equations
	set reduced [toRREF $matrix]
 
	# Extract the deconvolution from the last column of the reduced matrix
	for {set i 0} {$i<$hlen} {incr i} {
	    lappend result [lindex $reduced $i end]
	}
	return $result
    }
}