Deconvolution/1D: Difference between revisions

Content added Content deleted

Inline

Revision as of 14:28, 18 March 2010

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

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

for all integers ${\mathit {n}}$ . Assume $F(n)$ can be non-zero only for $0$ ≤ ${\mathit {n}}$ ≤ $|{\mathit {F}}|$ , where $|{\mathit {F}}|$ is the "length" of ${\mathit {F}}$ , and similarly for ${\mathit {G}}$ and ${\mathit {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 $|{\mathit {F}}|=6$ and $|{\mathit {H}}|=5$ would determine the following value of ${\mathit {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 ${\mathit {h}}$ given ${\mathit {f}}$ and ${\mathit {g}}$ .

The function should work for ${\mathit {G}}$ of arbitrary length (i.e., not hard coded or constant) and ${\mathit {F}}$ of any length up to that of ${\mathit {G}}$ . Note that $|{\mathit {H}}|$ will be given by $|{\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$ 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]

Fortran

This solution uses the LAPACK95 library. <lang fortran> ! 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 </lang> Results: <lang fortran> deconv(f, g) = -8, -9, -3, -1, -6, 7 deconv(h, g) = -3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1 </lang>

J

This solution borrowed from Formal power series:

<lang J>Ai=: (i.@] =/ i.@[ -/ i.@>:@-)&# divide=: [ +/ .*~ [:%.&.x: ] +/ .* Ai</lang>

Sample data:

Example use: <lang J> g divide f _8 _9 _3 _1 _6 7

  g divide h

_3 _6 _1 8 _6 3 _1 _9 _9 3 _2 5 2 _2 _7 _1</lang>

That said, note that this particular implementation is slow since it uses extended precision intermediate results. It will run quite a bit faster for this example with no notable loss of precision if floating point is used. In other words:

<lang J>divide=: [ +/ .*~ [:%. ] +/ .* Ai</lang>

Lua

Using metatables: <lang lua>function deconvolve(f, g)

  local h = setmetatable({}, {__index = function(self, n)
     if n == 1 then self[1] = g[1] / f[1]
     else
        self[n] = g[n]
        for i = 1, n - 1 do
           self[n] = self[n] - self[i] * (f[n - i + 1] or 0)
        end
        self[n] = self[n] / f[1]
     end
     return self[n]
  end})
  local _ = h[#g - #f + 1]
  return setmetatable(h, nil)

end</lang>

Tests: <lang lua> local f = {-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1} local g = {24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7} local h = {-8,-9,-3,-1,-6,7} print(unpack(deconvolve(f, g))) --> -8 -9 -3 -1 -6 7 print(unpack(deconvolve(h, g))) --> -3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1</lang>

Python

Works with: Python version 3.x

Inspired by the TCL solution, and using the ToReducedRowEchelonForm function to reduce to row echelon form from here <lang python>def ToReducedRowEchelonForm( M ):

   if not M: return
   lead = 0
   rowCount = len(M)
   columnCount = len(M[0])
   for r in range(rowCount):
       if lead >= columnCount:
           return
       i = r
       while M[i][lead] == 0:
           i += 1
           if i == rowCount:
               i = r
               lead += 1
               if columnCount == lead:
                   return
       M[i],M[r] = M[r],M[i]
       lv = M[r][lead]
       M[r] = [ mrx / lv for mrx in M[r]]
       for i in range(rowCount):
           if i != r:
               lv = M[i][lead]
               M[i] = [ iv - lv*rv for rv,iv in zip(M[r],M[i])]
       lead += 1

def pmtx(mtx):

   print ('\n'.join(.join(' %4s' % col for col in row) for row in mtx))

def convolve(f, h):

   g = [0] * (len(f) + len(h) - 1)
   for hindex, hval in enumerate(h):
       for findex, fval in enumerate(f):
           g[hindex + findex] += fval * hval
   return g

def deconvolve(g, f):

   lenh = len(g) - len(f) + 1
   mtx = [[0 for x in range(lenh+1)] for y in g]
   for hindex in range(lenh):
       for findex, fval in enumerate(f):
           gindex = hindex + findex
           mtx[gindex][hindex] = fval
   for gindex, gval in enumerate(g):        
       mtx[gindex][lenh] = gval
   ToReducedRowEchelonForm( mtx )
   return [mtx[i][lenh] for i in range(lenh)]  # h

if __name__ == '__main__':

   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]
   assert convolve(f,h) == g
   assert deconvolve(g, f) == h</lang>

R

Here we won't solve the system but use the FFT instead. The method :

extend vector arguments so that they are the same length, a power of 2 larger than the length of the solution,
solution is ifft(fft(a)*fft(b)), truncated.

<lang R> nextpow <- function(n, s=2) { r <- 1 while(r<n) r <- r*s r }

conv <- function(a, b) { p <- length(a) q <- length(b) n <- p + q - 1 r <- nextpow(n) y <- fft(fft(c(a, rep(0, r-p))) * fft(c(b, rep(0, r-q))), inverse=TRUE)/r y[1:n] }

deconv <- function(a, b) { p <- length(a) q <- length(b) n <- p - q + 1 r <- nextpow(max(p, q)) y <- fft(fft(c(a, rep(0, r-p))) / fft(c(b, rep(0, r-q))), inverse=TRUE)/r return(y[1:n]) } </lang>

To check :

<lang R> h <- c(-8,-9,-3,-1,-6,7) f <- c(-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1) g <- c(24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7)

max(abs(conv(f,h) - g)) max(abs(deconv(g,f) - h)) max(abs(deconv(g,h) - f)) </lang>

This solution often introduces complex numbers, with null or tiny imaginary part. If it hurts in applications, type Re(conv(f,h)) and Re(deconv(g,h)) instead, to return only the real part. It's not hard-coded in the functions, since they may be used for complex arguments as well.

R has also a function convolve, <lang R> conv(a, b) == convolve(a, rev(b), type="open") </lang>

Tcl

Works with: Tcl version 8.5

This builds the a command, 1D, with two subcommands (convolve and deconvolve) for performing convolution and deconvolution of these kinds of arrays. The deconvolution code is based on a reduction to reduced row echelon form. <lang tcl>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

}</lang> To use the above code, a simple demonstration driver (which solves the specific task): <lang tcl># Simple pretty-printer proc pp {name nlist} {

   set sep ""
   puts -nonewline "$name = \["
   foreach n $nlist {

puts -nonewline [format %s%g $sep $n] set sep ,

   }
   puts "\]"

}

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

pp "deconv(g,f) = h" [1D deconvolve $g $f] pp "deconv(g,h) = f" [1D deconvolve $g $h] pp " conv(f,h) = g" [1D convolve $f $h]</lang> Output:

deconv(g,f) = h = [-8,-9,-3,-1,-6,7]
deconv(g,h) = f = [-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1]
  conv(f,h) = 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

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.>

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>>