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]


 
! Build
! Windows: ifort /I "%IFORT_COMPILER11%\mkl\include\ia32" deconv1d.f90 "%IFORT_COMPILER11%\mkl\ia32\lib\*.lib"
! Linux:
 
program deconv
  ! Use gelsd from LAPACK95.
  use mkl95_lapack, only : gelsd
 
  implicit none
  real(8), allocatable :: g(:), href(:), A(:,:), f(:)
  real(8), pointer     :: h(:), r(:)
  integer              :: N
  character(len=16)    :: cbuff
  integer              :: i
  intrinsic            :: nint
 
  ! Allocate data arrays
  allocate(g(21),f(16))
  g = [24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7]
  f = [-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1]
 
  ! Calculate deconvolution
  h => deco(f, g)
 
  ! Check result against reference
  N = size(h)
  allocate(href(N))
  href = [-8,-9,-3,-1,-6,7]
  cbuff = ' '
  write(cbuff,'(a,i0,a)') '(a,',N,'(i0,a),i0)'
  if (any(abs(h-href) > 1.0d-4)) then
     write(*,'(a)') 'deconv(f, g) - FAILED'
  else
     write(*,cbuff) 'deconv(f, g) = ',(nint(h(i)),', ',i=1,N-1),nint(h(N))
  end if
 
  ! Calculate deconvolution
  r => deco(h, g)
 
  cbuff = ' '
  N = size(r)
  write(cbuff,'(a,i0,a)') '(a,',N,'(i0,a),i0)'
  if (any(abs(r-f) > 1.0d-4)) then
     write(*,'(a)') 'deconv(h, g) - FAILED'
  else
     write(*,cbuff) 'deconv(h, g) = ',(nint(r(i)),', ',i=1,N-1),nint(r(N))
  end if
 
contains
  function deco(p, q)
    real(8), pointer    :: deco(:)
    real(8), intent(in) :: p(:), q(:)
 
    real(8), allocatable, target :: r(:)
    real(8), allocatable         :: A(:,:)
    integer :: N
 
    ! Construct derived arrays
    N = size(q) - size(p) + 1
    allocate(A(size(q),N),r(size(q)))
    A = 0.0d0
    do i=1,N
       A(i:i+size(p)-1,i) = p
    end do
 
    ! Invoke the LAPACK routine to do the work
    r = q
    call gelsd(A, r)
 
    deco => r(1:N)
  end function deco
 
end program deconv
 