Jump to content

Deconvolution/1D

From Rosetta Code
Task
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 and of an integer variable is defined as the function satisfying

for all integers . Assume can be non-zero only for ≤ ≤ , where is the "length" of , and similarly for and , so that the functions can be modeled as finite sequences by identifying with , etc. Then for example, values of and would determine the following value of by definition.

We can write this in matrix form as:

or

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 for given and .

  • The function should work for of arbitrary length (i.e., not hard coded or constant) and of any length up to that of . Note that will be given by .
  • 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 and deconv 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]

11l

Translation of: D
F deconv(g, f)
   V result = [0]*(g.len - f.len + 1)
   L(&e) result
      V n = L.index
      e = g[n]
      V lower_bound = I n >= f.len {n - f.len + 1} E 0
      L(i) lower_bound .< n
         e -= result[i] * f[n - i]
      e /= f[0]
   R result

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

Ada

This is a translation of the D solution. 

with Ada.Text_IO; use Ada.Text_IO;

procedure Main is
   package real_io is new Float_IO (Long_Float);
   use real_io;

   type Vector is array (Natural range <>) of Long_Float;

   function deconv (g, f : Vector) return Vector is
      len : Positive :=
        Integer'Max ((g'Length - f'length), (f'length - g'length));
      h     : Vector (0 .. len);
      Lower : Natural := 0;
   begin
      for n in h'range loop
         h (n) := g (n);
         if n >= f'length then
            Lower := n - f'length + 1;
         end if;
         for i in Lower .. n - 1 loop
            h (n) := h (n) - (h (i) * f (n - i));
         end loop;
         h (n) := h (n) / f (0);
      end loop;
      return h;
   end deconv;

   procedure print (v : Vector) is
   begin
      Put ("(");
      for I in v'range loop
         Put (Item => v (I), Fore => 1, Aft => 1, Exp => 0);
         if I < v'Last then
            Put (" ");
         else
            Put_Line (")");
         end if;
      end loop;
   end print;

   h : Vector := (-8.0, -9.0, -3.0, -1.0, -6.0, 7.0);
   f : Vector :=
     (-3.0, -6.0, -1.0, 8.0, -6.0, 3.0, -1.0, -9.0, -9.0, 3.0, -2.0, 5.0, 2.0,
      -2.0, -7.0, -1.0);
   g : Vector :=
     (24.0, 75.0, 71.0, -34.0, 3.0, 22.0, -45.0, 23.0, 245.0, 25.0, 52.0, 25.0,
      -67.0, -96.0, 96.0, 31.0, 55.0, 36.0, 29.0, -43.0, -7.0);
begin
   print (h);
   print (deconv (g, f));
   print (f);
   print (deconv (g, h));
end Main;
Output:
(-8.0 -9.0 -3.0 -1.0 -6.0 7.0)
(-8.0 -9.0 -3.0 -1.0 -6.0 7.0)
(-3.0 -6.0 -1.0 8.0 -6.0 3.0 -1.0 -9.0 -9.0 3.0 -2.0 5.0 2.0 -2.0 -7.0 -1.0)
(-3.0 -6.0 -1.0 8.0 -6.0 3.0 -1.0 -9.0 -9.0 3.0 -2.0 5.0 2.0 -2.0 -7.0 -1.0)

BBC BASIC

Works with: BBC BASIC for Windows

As several others, this is a translation of the D solution. 

      *FLOAT 64
      DIM h(5), f(15), g(20)
      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
      
      PROCdeconv(g(), f(), x())
      PRINT "deconv(g,f) = " FNprintarray(x())
      x() -= h() : IF SUM(x()) <> 0 PRINT "Error!"
      
      PROCdeconv(g(), h(), y())
      PRINT "deconv(g,h) = " FNprintarray(y())
      y() -= f() : IF SUM(y()) <> 0 PRINT "Error!"
      END
      
      DEF PROCdeconv(g(), f(), RETURN h())
      LOCAL f%, g%, i%, l%, n%
      f% = DIM(f(),1) + 1
      g% = DIM(g(),1) + 1
      DIM h(g% - f%)
      FOR n% = 0 TO g% - f%
        h(n%) = g(n%)
        IF n% < f% THEN l% = 0 ELSE l% = n% - f% + 1
        IF n% THEN
          FOR i% = l% TO n% - 1
            h(n%) -= h(i%) * f(n% - i%)
          NEXT
        ENDIF
        h(n%) /= f(0)
      NEXT n%
      ENDPROC
      
      DEF FNprintarray(a())
      LOCAL i%, a$
      FOR i% = 0 TO DIM(a(),1)
        a$ += STR$(a(i%)) + ", "
      NEXT
      = LEFT$(LEFT$(a$))
Output:
deconv(g,f) = -8, -9, -3, -1, -6, 7
deconv(g,h) = -3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1

C

Using FFT: 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
 
double PI;
typedef double complex cplx;
 
void _fft(cplx buf[], cplx out[], int n, int step)
{
	if (step < n) {
		_fft(out, buf, n, step * 2);
		_fft(out + step, buf + step, n, step * 2);
 
		for (int i = 0; i < n; i += 2 * step) {
			cplx t = cexp(-I * PI * i / n) * out[i + step];
			buf[i / 2]     = out[i] + t;
			buf[(i + n)/2] = out[i] - t;
		}
	}
}
 
void fft(cplx buf[], int n)
{
	cplx out[n];
	for (int i = 0; i < n; i++) out[i] = buf[i];
	_fft(buf, out, n, 1);
}

/* pad array length to power of two */
cplx *pad_two(double g[], int len, int *ns)
{
	int n = 1;
	if (*ns) n = *ns;
	else while (n < len) n *= 2;

	cplx *buf = calloc(sizeof(cplx), n);
	for (int i = 0; i < len; i++) buf[i] = g[i];
	*ns = n;
	return buf;
}

void deconv(double g[], int lg, double f[], int lf, double out[]) {
	int ns = 0;
	cplx *g2 = pad_two(g, lg, &ns);
	cplx *f2 = pad_two(f, lf, &ns);

	fft(g2, ns);
	fft(f2, ns);

	cplx h[ns];
	for (int i = 0; i < ns; i++) h[i] = g2[i] / f2[i];
	fft(h, ns);

	for (int i = 0; i >= lf - lg; i--)
		out[-i] = h[(i + ns) % ns]/32;
	free(g2);
	free(f2);
}

int main()
{
	PI = atan2(1,1) * 4;
	double g[] = {24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7};
	double f[] = { -3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1 };
	double h[] = { -8,-9,-3,-1,-6,7 };

	int lg = sizeof(g)/sizeof(double);
	int lf = sizeof(f)/sizeof(double);
	int lh = sizeof(h)/sizeof(double);

	double h2[lh];
	double f2[lf];

	printf("f[] data is : ");
	for (int i = 0; i < lf; i++) printf(" %g", f[i]);
	printf("\n");

	printf("deconv(g, h): ");
	deconv(g, lg, h, lh, f2);
	for (int i = 0; i < lf; i++) printf(" %g", f2[i]);
	printf("\n");

	printf("h[] data is : ");
	for (int i = 0; i < lh; i++) printf(" %g", h[i]);
	printf("\n");

	printf("deconv(g, f): ");
	deconv(g, lg, f, lf, h2);
	for (int i = 0; i < lh; i++) printf(" %g", h2[i]);
	printf("\n");
}
Output:
f[] data is :  -3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1

deconv(g, h):  -3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1
h[] data is :  -8 -9 -3 -1 -6 7


deconv(g, f):  -8 -9 -3 -1 -6 7

C++ 

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <vector>

void print_vector(const std::vector<int32_t>& list) {
	std::cout << "[";
	for ( uint64_t i = 0; i < list.size() - 1; ++i ) {
		std::cout << list[i] << ", ";
	}
	std::cout << list.back() << "]" << std::endl;
}

std::vector<int32_t> deconvolution(const std::vector<int32_t>& a, const std::vector<int32_t>& b) {
	std::vector<int32_t> result(a.size() - b.size() + 1, 0);
	for ( uint64_t n = 0; n < result.size(); n++ ) {
		result[n] = a[n];
		uint64_t start = std::max((int) (n - b.size() + 1), 0);
		for ( uint64_t i = start; i < n; i++ ) {
			result[n] -= result[i] * b[n - i];
		}
		result[n] /= b[0];
	}
	return result;
}

int main() {
	const std::vector<int32_t> h = { -8, -9, -3, -1, -6, 7 };
	const std::vector<int32_t> f = { -3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1 };
	const std::vector<int32_t> g = { 24, 75, 71, -34, 3, 22, -45, 23, 245, 25, 52,
                                     25, -67, -96, 96, 31, 55, 36, 29, -43, -7 };

    std::cout << "h =                   "; print_vector(h);
	std::cout << "deconvolution(g, f) = "; print_vector(deconvolution(g, f));
	std::cout << "f =                   "; print_vector(f);
	std::cout << "deconvolution(g, h) = "; print_vector(deconvolution(g, h));
}
Output:
h =                   [-8, -9, -3, -1, -6, 7]
deconvolution(g, f) = [-8, -9, -3, -1, -6, 7]
f =                   [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]
deconvolution(g, h) = [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]

Common Lisp

Uses the routine (lsqr A b) from Multiple regression and (mtp A) from Matrix transposition. 

;; Assemble the mxn matrix A from the 2D row vector x.
(defun make-conv-matrix (x m n)
  (let ((lx (cadr (array-dimensions x)))
        (A  (make-array `(,m ,n) :initial-element 0)))

    (loop for j from 0 to (- n 1) do
         (loop for i from 0 to (- m 1) do
              (setf (aref A i j)
                    (cond ((or (< i j) (>= i (+ j lx)))
                           0)
                          ((and (>= i j) (< i (+ j lx)))
                           (aref x 0 (- i j)))))))
    A))

;; Solve the overdetermined system A(f)*h=g by linear least squares.
(defun deconv (g f)
  (let* ((lg (cadr (array-dimensions g)))
         (lf (cadr (array-dimensions f)))
         (lh (+ (- lg lf) 1))
         (A  (make-conv-matrix f lg lh)))

    (lsqr A (mtp g))))

Example: 

(setf f #2A((-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1)))
(setf h #2A((-8 -9 -3 -1 -6 7)))
(setf g #2A((24 75 71 -34 3 22 -45 23 245 25 52 25 -67 -96 96 31 55 36 29 -43 -7)))

(deconv g f)
#2A((-8.0)
    (-9.000000000000002)
    (-2.999999999999999)
    (-0.9999999999999997)
    (-6.0)
    (7.000000000000002))

(deconv g h)
#2A((-2.999999999999999)
    (-6.000000000000001)
    (-1.0000000000000002)
    (8.0)
    (-5.999999999999999)
    (3.0000000000000004)
    (-1.0000000000000004)
    (-9.000000000000002)
    (-9.0)
    (2.9999999999999996)
    (-1.9999999999999991)
    (5.0)
    (1.9999999999999996)
    (-2.0000000000000004)
    (-7.000000000000001)
    (-0.9999999999999994))

D 

T[] deconv(T)(in T[] g, in T[] f) pure nothrow {
    int flen = f.length;
    int glen = g.length;
    auto result = new T[glen - flen + 1];
    foreach (int n, ref e; result) {
        e = g[n];
        immutable lowerBound = (n >= flen) ? n - flen + 1 : 0;
        foreach (i; lowerBound .. n)
                e -= result[i] * f[n - i];
        e /= f[0];
    }
    return result;
}

void main() {
    import std.stdio;
    immutable h = [-8,-9,-3,-1,-6,7];
    immutable f = [-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1];
    immutable g = [24,75,71,-34,3,22,-45,23,245,25,52,25,-67,
                   -96,96,31,55,36,29,-43,-7];
    writeln(deconv(g, f) == h, " ", deconv(g, f));
    writeln(deconv(g, h) == f, " ", deconv(g, h));
}
Output:
true [-8, -9, -3, -1, -6, 7]
true [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]

EasyLang

Translation of: Java
func[] deconv g[] f[] .
   len h[] len g[] - len f[] + 1
   for n = 1 to len h[]
      h[n] = g[n]
      low = higher 1 (n - len f[] + 1)
      for i = low to n - 1
         h[n] -= h[i] * f[n - i + 1]
      .
      h[n] /= f[1]
   .
   return h[]
.
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 ]
print h[]
print deconv g[] f[]
print f[]
print deconv g[] h[]
Output:
[ -8 -9 -3 -1 -6 7 ]
[ -8 -9 -3 -1 -6 7 ]
[ -3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1 ]
[ -3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1 ]

Fortran

This solution uses the LAPACK95 library. 

! 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

Results: 

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

FreeBASIC 

Sub Deconv(g() As Double, f() As Double, h() As Double)
    Dim As Integer n, i, lower
    Dim As Integer hCount = Ubound(g) - Ubound(f) + 2
    Redim h(hCount - 1)
    
    For n = 0 To hCount - 1
        h(n) = g(n)
        lower = Iif(n >= Ubound(f) + 1, n - Ubound(f), 0)
        i = lower
        While i < n
            h(n) -= h(i) * f(n - i)
            i += 1
        Wend
        h(n) /= f(0)
    Next n
End Sub

Dim As Integer i
Dim As Double h(5) = {-8, -9, -3, -1, -6, 7}
Dim As Double f(15) = {-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1}
Dim As Double g(20) = {24, 75, 71, -34, 3, 22, -45, 23, 245, 25, 52, 25, -67, -96, 96, 31, 55, 36, 29, -43, -7}
Dim As Double result()

Print !"h:\n[";
For i = Lbound(h) To Ubound(h)
    Print h(i); ",";
Next i
Print Chr(8) & !"]\n";

Deconv(g(), f(), result())
Print !"\deconv(g, f):\n[";
For i = Lbound(result) To Ubound(result)-1
    Print result(i); ",";
Next i
Print Chr(8) & !"]\n";

Print
Print !"f:\n[";
For i = Lbound(f) To Ubound(f)
    Print f(i); ",";
Next i
Print Chr(8) & !"]\n";

Deconv(g(), h(), result())
Print !"\deconv(g, h):\n[";
For i = Lbound(result) To Ubound(result)-1
    Print Using "##_,"; result(i);
Next i
Print Chr(8) & !"]\n";
Sleep
Output:
h:
[-8,-9,-3,-1,-6, 7]
deconv(g, f):
[-8,-9,-3,-1,-6, 7]

f:
[-3,-6,-1, 8,-6, 3,-1,-9,-9, 3,-2, 5, 2,-2,-7,-1]
deconv(g, h):
[-3,-6,-1, 8,-6, 3,-1,-9,-9, 3,-2, 5, 2,-2,-7,-1]

Go

Translation of: D
package main

import "fmt"

func main() {
    h := []float64{-8, -9, -3, -1, -6, 7}
    f := []float64{-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1}
    g := []float64{24, 75, 71, -34, 3, 22, -45, 23, 245, 25, 52, 25, -67, -96,
        96, 31, 55, 36, 29, -43, -7}
    fmt.Println(h)
    fmt.Println(deconv(g, f))
    fmt.Println(f)
    fmt.Println(deconv(g, h))
}

func deconv(g, f []float64) []float64 {
    h := make([]float64, len(g)-len(f)+1)
    for n := range h {
        h[n] = g[n]
        var lower int
        if n >= len(f) {
            lower = n - len(f) + 1
        }
        for i := lower; i < n; i++ {
            h[n] -= h[i] * f[n-i]
        }
        h[n] /= f[0]
    }
    return h
}
Output:
[-8 -9 -3 -1 -6 7]
[-8 -9 -3 -1 -6 7]
[-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1]
[-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1]
Translation of: C
package main

import (
    "fmt"
    "math"
    "math/cmplx"
)

func main() {
    h := []float64{-8, -9, -3, -1, -6, 7}
    f := []float64{-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1}
    g := []float64{24, 75, 71, -34, 3, 22, -45, 23, 245, 25, 52, 25, -67, -96,
        96, 31, 55, 36, 29, -43, -7}
    fmt.Printf("%.1f\n", h)
    fmt.Printf("%.1f\n", deconv(g, f))
    fmt.Printf("%.1f\n", f)
    fmt.Printf("%.1f\n", deconv(g, h))
}

func deconv(g, f []float64) []float64 {
    n := 1
    for n < len(g) {
        n *= 2
    }
    g2 := make([]complex128, n)
    for i, x := range g {
        g2[i] = complex(x, 0)
    }
    f2 := make([]complex128, n)
    for i, x := range f {
        f2[i] = complex(x, 0)
    }
    gt := fft(g2)
    ft := fft(f2)
    for i := range gt {
        gt[i] /= ft[i]
    }
    ht := fft(gt)
    it := 1 / float64(n)
    out := make([]float64, len(g)-len(f)+1)
    out[0] = real(ht[0]) * it
    for i := 1; i < len(out); i++ {
        out[i] = real(ht[n-i]) * it
    }
    return out
}

func fft(in []complex128) []complex128 {
    out := make([]complex128, len(in))
    ditfft2(in, out, len(in), 1)
    return out
}

func ditfft2(x, y []complex128, n, s int) {
    if n == 1 {
        y[0] = x[0]
        return
    }
    ditfft2(x, y, n/2, 2*s)
    ditfft2(x[s:], y[n/2:], n/2, 2*s)
    for k := 0; k < n/2; k++ {
        tf := cmplx.Rect(1, -2*math.Pi*float64(k)/float64(n)) * y[k+n/2]
        y[k], y[k+n/2] = y[k]+tf, y[k]-tf
    }
}
Output:

Some results have errors out in the last decimal place or so. Only one decimal place shown here to let results fit in 80 columns. 

[-8.0 -9.0 -3.0 -1.0 -6.0 7.0]
[-8.0 -9.0 -3.0 -1.0 -6.0 7.0]
[-3.0 -6.0 -1.0 8.0 -6.0 3.0 -1.0 -9.0 -9.0 3.0 -2.0 5.0 2.0 -2.0 -7.0 -1.0]
[-3.0 -6.0 -1.0 8.0 -6.0 3.0 -1.0 -9.0 -9.0 3.0 -2.0 5.0 2.0 -2.0 -7.0 -1.0]

Library gonum/mat: 

package main

import (
    "fmt"

    "gonum.org/v1/gonum/mat"
)

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

func band(g, f []float64) *mat.Dense {
    nh := len(g) - len(f) + 1
    b := mat.NewDense(len(g), nh, nil)
    for j := 0; j < nh; j++ {
        for i, fi := range f {
            b.Set(i+j, j, fi)
        }
    }
    return b
}

func deconv(g, f []float64) mat.Matrix {
    z := mat.NewDense(len(g)-len(f)+1, 1, nil)
    z.Solve(band(g, f), mat.NewVecDense(len(g), g))
    return z
}

func main() {
    fmt.Printf("deconv(g, f) =\n%.1f\n\n", mat.Formatted(deconv(g, f)))
    fmt.Printf("deconv(g, h) =\n%.1f\n", mat.Formatted(deconv(g, h)))
}
Output:
deconv(g, f) =
⎡-8.0⎤
⎢-9.0⎥
⎢-3.0⎥
⎢-1.0⎥
⎢-6.0⎥
⎣ 7.0⎦

deconv(g, h) =
⎡-3.0⎤
⎢-6.0⎥
⎢-1.0⎥
⎢ 8.0⎥
⎢-6.0⎥
⎢ 3.0⎥
⎢-1.0⎥
⎢-9.0⎥
⎢-9.0⎥
⎢ 3.0⎥
⎢-2.0⎥
⎢ 5.0⎥
⎢ 2.0⎥
⎢-2.0⎥
⎢-7.0⎥
⎣-1.0⎦

Haskell 

deconv1d :: [Double] -> [Double] -> [Double]
deconv1d xs ys = takeWhile (/= 0) $ deconv xs ys
  where
    [] `deconv` _ = []
    (0:xs) `deconv` (0:ys) = xs `deconv` ys
    (x:xs) `deconv` (y:ys) =
      let q = x / y
      in q : zipWith (-) xs (scale q ys ++ repeat 0) `deconv` (y : ys)

scale :: Double -> [Double] -> [Double]
scale = map . (*)

h, f, g :: [Double]
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
  ]

main :: IO ()
main = print $ (h == deconv1d g f) && (f == deconv1d g h)
Output:
True

J

This solution borrowed from Formal power series: 

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

Sample data: 

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

Example use: 

   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

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: 

divide=: [ +/ .*~ [:%. ] +/ .* Ai

Java

Translation of: Go
import java.util.Arrays;

public class Deconvolution1D {
    public static int[] deconv(int[] g, int[] f) {
        int[] h = new int[g.length - f.length + 1];
        for (int n = 0; n < h.length; n++) {
            h[n] = g[n];
            int lower = Math.max(n - f.length + 1, 0);
            for (int i = lower; i < n; i++)
                h[n] -= h[i] * f[n - i];
            h[n] /= f[0];
        }
        return h;
    }

    public static void main(String[] args) {
        int[] h = { -8, -9, -3, -1, -6, 7 };
        int[] f = { -3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1 };
        int[] g = { 24, 75, 71, -34, 3, 22, -45, 23, 245, 25, 52, 25, -67, -96,
                96, 31, 55, 36, 29, -43, -7 };

        StringBuilder sb = new StringBuilder();
        sb.append("h = " + Arrays.toString(h) + "\n");
        sb.append("deconv(g, f) = " + Arrays.toString(deconv(g, f)) + "\n");
        sb.append("f = " + Arrays.toString(f) + "\n");
        sb.append("deconv(g, h) = " + Arrays.toString(deconv(g, h)) + "\n");
        System.out.println(sb.toString());
    }
}
Output:
h = [-8, -9, -3, -1, -6, 7]
deconv(g, f) = [-8, -9, -3, -1, -6, 7]
f = [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]
deconv(g, h) = [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]

jq

Translation of: Wren

Works with jq, the C implementation of jq

Works with gojq, the Go implementation of jq

Works with jaq, the Rust implementation of jq 

def deconv($g; $f):
  { h: [range(0; ($g|length) - ($f|length) + 1) | 0] } 
  | reduce range  ( 0;.h|length) as $n (.;
      .h[$n] = $g[$n]
      | (if $n >= ($f|length) then $n - ($f|length) + 1 else 0 end) as $lower
      | .i = $lower
      | until(.i >= $n;
          .h[$n] -= .h[.i] * $f[$n - .i]
          | .i += 1 )
      | .h[$n] /= $f[0] )
  | .h ;

### The tasks

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

h,
deconv(g; f),
f,
deconv(g; h)
Output:
[-8,-9,-3,-1,-6,7]
[-8,-9,-3,-1,-6,7]
[-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1]
[-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1]

Julia

The deconv function for floating point data is built into Julia, though using DSP is required with version 1.0. Integer inputs may need to be converted and copied to floating point to use deconv(). 

h = [-8, -9, -3, -1, -6, 7]
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]

hanswer = deconv(float.(g), float.(f))
println("The deconvolution deconv(g, f) is $hanswer, which is the same as h = $h\n")

fanswer = deconv(float.(g), float.(h))
println("The deconvolution deconv(g, h) is $fanswer, which is the same as f = $f\n")
Output:
The deconvolution deconv(g, f) is [-8.0, -9.0, -3.0, -1.0, -6.0, 7.0], 
which is the same as h = [-8, -9, -3, -1, -6, 7]

The deconvolution deconv(g, h) is [-3.0, -6.0, -1.0, 8.0, -6.0, 3.0, -1.0, -9.0, -9.0, 3.0, -2.0, 5.0, 2.0, -2.0, -7.0, -1.0], 
which is the same as f = [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]

Kotlin

Translation of: Go
// version 1.1.3

fun deconv(g: DoubleArray, f: DoubleArray): DoubleArray {
    val fs = f.size   
    val h = DoubleArray(g.size - fs + 1)
    for (n in h.indices) {
        h[n] = g[n]
        val lower = if (n >= fs) n - fs + 1 else 0
        for (i in lower until n) h[n] -= h[i] * f[n -i]
        h[n] /= f[0] 
    }
    return h
}

fun main(args: Array<String>) {
    val h = doubleArrayOf(-8.0, -9.0, -3.0, -1.0, -6.0, 7.0)
    val f = doubleArrayOf(-3.0, -6.0, -1.0,  8.0, -6.0,  3.0, -1.0, -9.0, 
                          -9.0,  3.0, -2.0,  5.0,  2.0, -2.0, -7.0, -1.0)
    val g = doubleArrayOf(24.0,  75.0, 71.0, -34.0,  3.0,  22.0, -45.0, 
                          23.0, 245.0, 25.0,  52.0, 25.0, -67.0, -96.0,
                          96.0,  31.0, 55.0,  36.0, 29.0, -43.0,  -7.0)
    println("${h.map { it.toInt() }}")
    println("${deconv(g, f).map { it.toInt() }}")
    println()
    println("${f.map { it.toInt() }}")
    println("${deconv(g, h).map { it.toInt() }}")
}
Output:
[-8, -9, -3, -1, -6, 7]
[-8, -9, -3, -1, -6, 7]

[-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]
[-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]

Lua

Using metatables: 

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

Tests: 

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

Mathematica / Wolfram Language

This function creates a sparse array for the A matrix and then solves it with a built-in function. It may fail for overdetermined systems, though. Fast approximate methods for deconvolution are also built into Mathematica. See Deconvolution/2D+ 

deconv[f_List, g_List] := 
 Module[{A = 
    SparseArray[
     Table[Band[{n, 1}] -> f[[n]], {n, 1, Length[f]}], {Length[g], Length[f] - 1}]}, 
  Take[LinearSolve[A, g], Length[g] - Length[f] + 1]]

Usage: 

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};
deconv[f,g]
Output:
{-8, -9, -3, -1, -6, 7}

MATLAB

The deconvolution function is built-in to MATLAB as the "deconv(a,b)" function, where "a" and "b" are vectors storing the convolved function values and the values of one of the deconvoluted vectors of "a". To test that this operates according to the task spec we can test the criteria above: 

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

ans =

   -8.0000   -9.0000   -3.0000   -1.0000   -6.0000    7.0000

>> deconv(g,h)

ans =

    -3    -6    -1     8    -6     3    -1    -9    -9     3    -2     5     2    -2    -7    -1

Therefore, "deconv(a,b)" behaves as expected.

Nim 

proc deconv(g, f: openArray[float]): seq[float] =
  var h: seq[float] = newSeq[float](len(g) - len(f) + 1)
  for n in 0..<len(h):
    h[n] = g[n]
    var lower: int
    if n >= len(f):
      lower = n - len(f) + 1
    for i in lower..<n:
      h[n] -= h[i] * f[n - i]
    h[n] /= f[0]
  h

let h = [-8'f64, -9, -3, -1, -6, 7]
let f = [-3'f64, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]
let g = [24'f64, 75, 71, -34, 3, 22, -45, 23, 245, 25, 52, 25, -67, -96, 
         96, 31, 55, 36, 29, -43, -7]
echo h
echo deconv(g, f)
echo f
echo deconv(g, h)
Output:
[-8.0, -9.0, -3.0, -1.0, -6.0, 7.0]
@[-8.0, -9.0, -3.0, -1.0, -6.0, 7.0]
[-3.0, -6.0, -1.0, 8.0, -6.0, 3.0, -1.0, -9.0, -9.0, 3.0, -2.0, 5.0, 2.0, -2.0, -7.0, -1.0]
@[-3.0, -6.0, -1.0, 8.0, -6.0, 3.0, -1.0, -9.0, -9.0, 3.0, -2.0, 5.0, 2.0, -2.0, -7.0, -1.0]

PascalABC.NET

Translation of: Nim
##
function deconv(g, f: array of real): array of real;
begin
  var h: array of real; 
  setlength(h, g.length - f.length + 1);
  for var n := 0 to h.length - 1 do
  begin
    h[n] := g[n];
    var lower: integer;
    if n >= f.length then 
      lower := n - f.length + 1;
    for var i := lower to n - 1 do
      h[n] -= h[i] * f[n - i];
    h[n] /= f[0];
  end;
  result := h;
end;

var h := |-8.0, -9, -3, -1, -6, 7|;
var f := |-3.0, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1|;
var g := |24.0, 75, 71, -34, 3, 22, -45, 23, 245, 25, 52, 25, -67, -96, 
          96, 31, 55, 36, 29, -43, -7|;
h.Println;
deconv(g, f).Println;
f.Println;
deconv(g, h).Println;
Output:
-8 -9 -3 -1 -6 7
-8 -9 -3 -1 -6 7
-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1
-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1

Perl

Using rref routine from Reduced row echelon form task.

Translation of: Raku
use v5.36;
use Math::Cartesian::Product;

sub deconvolve($g,$f) {
    my @g = @{$g};
    my @f = @{$f};
    my(@m,@d);

    my $h = 1 + @g - @f;
    push @m, [(0) x $h, $g[$_]] for 0..$#g;
    for my $j (0..$h-1) {
        for my $k (0..$#f) {
            $m[$j + $k][$j] = $f[$k]
        }
    }
    rref(\@m);
    push @d, @{ $m[$_] }[$h] for 0..$h-1;
    @d;
}

sub convolve($f,$h) {
    my @f = @{$f};
    my @h = @{$h};
    my @i;
    for my $x (cartesian {@_} [0..$#f], [0..$#h]) {
        push @i, @$x[0]+@$x[1];
    }
    my $cnt = 0;
    my @g = (0) x (@f + @h - 1);
    for my $x (cartesian {@_} [@f], [@h]) {
        $g[$i[$cnt++]] += @$x[0]*@$x[1];
    }
    @g;
}

sub rref($m) {
  my @m = @{$m};
  @m or return;
  my ($lead, $rows, $cols) = (0, scalar(@m), scalar(@{$m[0]}));

  for my $r (0 .. $rows - 1) {
      $lead < $cols or return;
      my $i = $r;

      until ($m[$i][$lead]) {
          ++$i == $rows or next;
          $i = $r;
          ++$lead == $cols and return;
      }

      @m[$i, $r] = @m[$r, $i];
      my $lv = $m[$r][$lead];
      $_ /= $lv foreach @{ $m[$r] };

      my @mr = @{ $m[$r] };
      for my $i (0 .. $rows - 1) {
          $i == $r and next;
          ($lv, my $n) = ($m[$i][$lead], -1);
          $_ -= $lv * $mr[++$n] foreach @{ $m[$i] };
      }
      ++$lead;
   }
}

my @h = qw<-8 -9 -3 -1 -6 7>;
my @f = qw<-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1>;
print '  conv(f,h) = g = ' . join(' ', my @g = convolve(\@f, \@h)) . "\n";
print 'deconv(g,f) = h = ' . join(' ', deconvolve(\@g, \@f)) . "\n";
print 'deconv(g,h) = f = ' . join(' ', deconvolve(\@g, \@h)) . "\n";
Output:
  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
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

Phix

Translation of: D
with javascript_semantics
function deconv(sequence g, f)
    integer lf = length(f), lg = length(g), lh = lg-lf+1
    sequence h = repeat(0,lh)
    for n=1 to lh do
        atom e = g[n]
        for i=max(n-lf,0) to n-2 do
            e -= h[i+1] * f[n-i]
        end for
        h[n] = e/f[1]
    end for
    return h
end function

function conv(sequence f, h)
    integer lf = length(f), lh = length(h), lg = lf+lh-1
    sequence g = repeat(0,lg)
    for i=1 to lh do
        for j=1 to lf do
            integer k = i+j-1
            g[k] += f[j] * h[i]
        end for
    end for
    return g
end function

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

procedure test(string desc, eq, sequence r, e)
    printf(1,"%s (%ssame as %s): %V\n",{desc,iff(r==e?"":"**NOT** "),eq,r})
end procedure

test("  conv(h,f)","g",  conv(h,f),g)
test("deconv(g,f)","h",deconv(g,f),h)
test("deconv(g,h)","f",deconv(g,h),f)
Output:
  conv(h,f) (same as g): {24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7}
deconv(g,f) (same as h): {-8,-9,-3,-1,-6,7}
deconv(g,h) (same as f): {-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1}

PicoLisp 

(load "@lib/math.l")

(de deconv (G F)
   (let A (pop 'F)
      (make
         (for (N . H) (head (- (length F)) G)
            (for (I . M) (made)
               (dec 'H
                  (*/ M (get F (- N I)) 1.0) ) )
            (link (*/ H 1.0 A)) ) ) ) )

Test: 

(setq
   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.)
   H (-8. -9. -3. -1. -6. 7.) )

(test H (deconv G F))
(test F (deconv G H))

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 

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

Based on the R version. 

import numpy

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]

# https://stackoverflow.com/questions/14267555/find-the-smallest-power-of-2-greater-than-n-in-python

def shift_bit_length(x):
    return 1<<(x-1).bit_length()

def conv(a, b):
    p = len(a)
    q = len(b)
    n = p + q - 1
    r = shift_bit_length(n)
    y = numpy.fft.ifft(numpy.fft.fft(a,r) * numpy.fft.fft(b,r),r)
    return numpy.trim_zeros(numpy.around(numpy.real(y),decimals=6))

def deconv(a, b):
    p = len(a)
    q = len(b)
    n = p - q + 1
    r = shift_bit_length(max(p, q))
    y = numpy.fft.ifft(numpy.fft.fft(a,r) / numpy.fft.fft(b,r), r)
    return numpy.trim_zeros(numpy.around(numpy.real(y),decimals=6))
    
# should return g
   
print(conv(h,f))

# should return h

print(deconv(g,f))

# should return f

print(deconv(g,h))

Output 

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

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.
conv <- function(a, b) {
	p <- length(a)
	q <- length(b)
	n <- p + q - 1
	r <- nextn(n, f=2)
	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 <- nextn(max(p, q), f=2)
	y <- fft(fft(c(a, rep(0, r-p))) / fft(c(b, rep(0, r-q))), inverse=TRUE)/r
	return(y[1:n])
}

To check : 

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

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, 

conv(a, b) == convolve(a, rev(b), type="open")

Racket 

#lang racket
(require math/matrix)
(define T matrix-transpose)

(define (convolution-matrix f m n)
  (define l (matrix-num-rows f))
  (for*/matrix m n ([i (in-range 0 m)] [j (in-range 0 n)])
      (cond [(or  (< i j) (>= i (+ j l)))  0]
            [(matrix-ref f (- i j) 0)])))
 
(define (least-square X y)
  (matrix-solve (matrix* (T X) X) (matrix* (T X) y)))

(define (deconvolve g f)
  (define lg (matrix-num-rows g))
  (define lf (matrix-num-rows f))
  (define lh (+ (- lg lf) 1))
  (least-square (convolution-matrix f lg lh) g))

Test: 

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

(deconvolve g f)
(deconvolve g h)
Output:
#<array '#(6 1) #[-8 -9 -3 -1 -6 7]>
#<array '#(16 1) #[-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1]>

Raku

(formerly Perl 6)

Translation of Python, using a modified version of the subroutine rref from Reduced row echelon form task. 

sub deconvolve (@g, @f) {
    my \h = 1 + @g - @f;
    my @m;
    @m[^@g;^h] »+=» 0;
    @m[^@g; h] »=«  @g;
    for ^h -> \j { for @f.kv -> \k, \v { @m[j+k;j] = v } }
    (rref @m)[^h;h]
}

sub convolve (@f, @h) {
    my @g = 0 xx + @f + @h - 1;
    @g[^@f X+ ^@h] »+=« (@f X× @h);
    @g
}
 
# Reduced Row Echelon Form simultaneous equation solver
# Can handle over-specified systems of equations (N unknowns in N + M equations)
sub rref (@m) {
    @m = trim-system @m;
    my ($lead, $rows, $cols) = 0, @m, @m[0];
    for ^$rows -> $r {
        return @m unless $lead < $cols;
        my $i = $r;
        until @m[$i;$lead] {
            next unless ++$i == $rows;
            $i = $r;
            return @m if ++$lead == $cols;
        }
        @m[$i, $r] = @m[$r, $i] if $r != $i;
        @m[$r] »/=» $ = @m[$r;$lead];
        for ^$rows -> $n {
            next if $n == $r;
            @m[$n] »-=» @m[$r] »×» (@m[$n;$lead] // 0);
        }
        ++$lead;
    }
    @m
 }

# Reduce to N equations in N unknowns; a no-op unless rows > cols
sub trim-system (@m) {
    return @m unless @m ≥ @m[0]; 
    my (\vars, @t) = @m[0] - 1;
    for ^vars -> \lead {
        for ^@m -> \row {
            @t.append: @m.splice(row, 1) and last if @m[row;lead];
        }
    }
    while @t < vars and @m { @t.push: shift @m }
    @t
}
  
my @h = (-8,-9,-3,-1,-6,7);
my @f = (-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1);
my @g = (24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7);
 
.say for ~@g,   ~convolve(@f, @h),'';
.say for ~@h, ~deconvolve(@g, @f),'';
.say for ~@f, ~deconvolve(@g, @h),'';
Output:
24 75 71 -34 3 22 -45 23 245 25 52 25 -67 -96 96 31 55 36 29 -43 -7
24 75 71 -34 3 22 -45 23 245 25 52 25 -67 -96 96 31 55 36 29 -43 -7

-8 -9 -3 -1 -6 7
-8 -9 -3 -1 -6 7

-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1
-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1

REXX 

/*REXX pgm performs deconvolution of two arrays:    deconv(g,f)=h   and   deconv(g,h)=f */
call make 'H',  "-8 -9 -3 -1 -6 7"
call make 'F',  "-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1"
call make 'G',  "24 75 71 -34 3 22 -45 23 245 25 52 25 -67 -96 96 31 55 36 29 -43 -7"
call show 'H'                                    /*display the elements of array  H.    */
call show 'F'                                    /*   "     "     "      "   "    F.    */
call show 'G'                                    /*   "     "     "      "   "    G.    */
call deco 'G',  "F", 'X'                         /*deconvolution of  G  and  F  ───►  X */
call test 'X',  "H"                              /*test: is array  H  equal to array  X?*/
call deco 'G',  "H", 'Y'                         /*deconvolution of  G  and  H  ───►  Y */
call test 'F',  "Y"                              /*test: is array  F  equal to array  Y?*/
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
deco: parse arg $1,$2,$r;   b= @.$2.# + 1;   a= @.$1.# + 1      /*get sizes of array 1&2*/
      @.$r.#= a - b                                             /*size of return array. */
            do n=0  to a-b                                      /*define  return array. */
            @.$r.n= @.$1.n                                      /*define RETURN element.*/
            if n<b  then L= 0                                   /*define the variable L.*/
                    else L= n - b + 1                           /*   "    "     "     " */
            if n>0  then do j=L  to n-1;                _= n-j  /*define elements > 0.  */
                         @.$r.n= @.$r.n - @.$r.j * @.$2._       /*compute   "     " "   */
                         end   /*j*/                            /* [↑] subtract product.*/
            @.$r.n= @.$r.n / @.$2.0                             /*divide array element. */
            end   /*n*/;                     return
/*──────────────────────────────────────────────────────────────────────────────────────*/
make: parse arg $,z;                     @.$.#= words(z) - 1    /*obtain args; set size.*/
            do k=0  to @.$.#;            @.$.k= word(z, k + 1)  /*define array element. */
            end   /*k*/;                     return             /*array starts at unity.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
show: parse arg $,z,_;    do s=0  to @.$.#;  _= strip(_ @.$.s)  /*obtain the arguments. */
                          end   /*s*/                           /* [↑]  build the list. */
      say 'array' $": " _;                   return             /*show the list;  return*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
test: parse arg $1,$2;    do t=0  to max(@.$1.#, @.$2.#)        /*obtain the arguments. */
                          if @.$1.t= @.$2.t  then iterate       /*create array list.    */
                          say "***error*** arrays"   $1    ' and '    $2   "aren't equal."
                          end   /*t*/;       return             /* [↑]  build the list. */
output   when using the default internal inputs:
array H:  -8 -9 -3 -1 -6 7
array F:  -3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1
array G:  24 75 71 -34 3 22 -45 23 245 25 52 25 -67 -96 96 31 55 36 29 -43 -7

RPL

Translation of: D

When translating to RPL, it is mandatory to take into account that:

  • array indexes start at 1
  • For loop variables, j shall be preferred to i, which is the name of the internal constant that equals √-1
  • FOR..NEXT loops are executed at least once
≪ → g f
   ≪ g SIZE f SIZE - 1 + 1 →LIST 0 CON
      1 g 1 GET f 1 GET / PUT
      2 OVER SIZE FOR n
         g n GET
         1 n f SIZE - 0 MAX +
         n 1 - FOR j
            OVER j GET
            f n j - 1 + GET * -
         NEXT 
         f 1 GET / n SWAP PUT
      NEXT
≫ ≫ 'DECONV' STO

≪ [-8 -9 -3 -1 -6 7]
   [-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1]
   [24 75 71 -34 3 22 -45 23 245 25 52 25 -67 -96 96 31 55 36 29 -43 -7]
  → h f g
  ≪ g f DECONV h ==
    g h DECONV f == AND
≫ ≫ ‘TASK’ STO
Output:
1: 1

Scala

Output:

Best seen running in your browser either by ScalaFiddle (ES aka JavaScript, non JVM) or Scastie (remote JVM). 

object Deconvolution1D extends App {
  val (h, f) = (Array(-8, -9, -3, -1, -6, 7), Array(-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1))
  val g = Array(24, 75, 71, -34, 3, 22, -45, 23, 245, 25, 52, 25, -67, -96, 96, 31, 55, 36, 29, -43, -7)
  val sb = new StringBuilder

  private def deconv(g: Array[Int], f: Array[Int]) = {
    val h = Array.ofDim[Int](g.length - f.length + 1)

    for (n <- h.indices) {
      h(n) = g(n)
      for (i <- math.max(n - f.length + 1, 0) until n) h(n) -= h(i) * f(n - i)
      h(n) /= f(0)
    }
    h
  }

  sb.append(s"h = ${h.mkString("[", ", ", "]")}\n")
    .append(s"deconv(g, f) = ${deconv(g, f).mkString("[", ", ", "]")}\n")
    .append(s"f = ${f.mkString("[", ", ", "]")}\n")
    .append(s"deconv(g, h) = ${deconv(g, h).mkString("[", ", ", "]")}")
  println(sb.result())

}

Swift

Translation of: Kotlin
func deconv(g: [Double], f: [Double]) -> [Double] {
  let fs = f.count
  var ret = [Double](repeating: 0, count: g.count - fs + 1)

  for n in 0..<ret.count {
    ret[n] = g[n]
    let lower = n >= fs ? n - fs + 1 : 0

    for i in lower..<n {
      ret[n] -= ret[i] * f[n - i]
    }

    ret[n] /= f[0]
  }

  return ret
}

let h = [-8.0, -9.0, -3.0, -1.0, -6.0, 7.0]
let f = [-3.0, -6.0, -1.0,  8.0, -6.0,  3.0, -1.0, -9.0,
         -9.0,  3.0, -2.0,  5.0,  2.0, -2.0, -7.0, -1.0]
let g = [24.0,  75.0, 71.0, -34.0,  3.0,  22.0, -45.0,
         23.0, 245.0, 25.0,  52.0, 25.0, -67.0, -96.0,
         96.0,  31.0, 55.0,  36.0, 29.0, -43.0,  -7.0]

print("\(h.map({ Int($0) }))")
print("\(deconv(g: g, f: f).map({ Int($0) }))\n")


print("\(f.map({ Int($0) }))")
print("\(deconv(g: g, f: h).map({ Int($0) }))")
Output:
[-8, -9, -3, -1, -6, 7]
[-8, -9, -3, -1, -6, 7]

[-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]
[-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]

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. 

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

To use the above code, a simple demonstration driver (which solves the specific task): 

# 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]
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). 

#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

test program: 

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

V (Vlang)

Translation of: Go
fn main() {
    h := [f64(-8), -9, -3, -1, -6, 7]
    f := [f64(-3), -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]
    g := [f64(24), 75, 71, -34, 3, 22, -45, 23, 245, 25, 52, 25, -67, -96,
        96, 31, 55, 36, 29, -43, -7]
    println(h)
    println(deconv(g, f))
    println(f)
    println(deconv(g, h))
}
 
fn deconv(g []f64, f []f64) []f64 {
    mut h := []f64{len: g.len-f.len+1}
    for n in 0..h.len {
        h[n] = g[n]
        mut lower := 0
        if n >= f.len {
            lower = n - f.len + 1
        }
        for i in lower..n {
            h[n] -= h[i] * f[n-i]
        }
        h[n] /= f[0]
    }
    return h
}
Output:
[-8, -9, -3, -1, -6, 7]
[-8, -9, -3, -1, -6, 7]
[-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]
[-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]

Wren

Translation of: Go
var deconv = Fn.new { |g, f|
    var h = List.filled(g.count - f.count + 1, 0)
    for (n in 0...h.count) {
        h[n] = g[n]
        var lower = (n >= f.count) ? n - f.count + 1 : 0
        var i = lower
        while (i < n) {
            h[n] = h[n] - h[i]*f[n-i]
            i = i + 1
        }
        h[n] = h[n] / f[0]
    }
    return h
}

var h = [-8, -9, -3, -1, -6, 7]
var f = [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]
var g = [24, 75, 71, -34, 3, 22, -45, 23, 245, 25, 52, 25, -67, -96, 96, 31, 55, 36, 29, -43, -7]
System.print(h)
System.print(deconv.call(g, f))
System.print(f)
System.print(deconv.call(g, h))
Output:
[-8, -9, -3, -1, -6, 7]
[-8, -9, -3, -1, -6, 7]
[-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]
[-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1]

zkl

Using GNU Scientific Library: 

var [const] GSL=Import("zklGSL");	// libGSL (GNU Scientific Library)
fcn dconv1D(f,g){
   fsz,hsz:=f.len(), g.len() - fsz +1;
   A:=GSL.Matrix(g.len(),hsz);
   foreach n,fn in ([0..].zip(f)){ foreach rc in (hsz){ A[rc+n,rc]=fn } }
   h:=A.AxEQb(g);
   h
}
f:=GSL.VectorFromData(-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1);
g:=GSL.VectorFromData(24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7);
h:=dconv1D(f,g);
h.format().println();

f:=dconv1D(h,g);
f.format().println();
Output:
-8.00,-9.00,-3.00,-1.00,-6.00,7.00
-3.00,-6.00,-1.00,8.00,-6.00,3.00,-1.00,-9.00,-9.00,3.00,-2.00,5.00,2.00,-2.00,-7.00,-1.00

Or, using lists:

Translation of: D
fcn deconv(g,f){
   flen, glen, delta:=f.len(), g.len(), glen - flen + 1;
   result:=List.createLong(delta); // allocate list with space for items
   foreach n in (delta){
      e:=g[n];
      lowerBound:=(if (n>=flen) n - flen + 1 else 0);
      foreach i in ([lowerBound .. n-1]){ e-=result[i]*f[n - i]; }
      result.append(e/f[0]);
    }
    result;
}
h:=T(-8,-9,-3,-1,-6,7);
f:=T(-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1);
g:=T(24,75,71,-34,3,22,-45,23,245,25,52,25,-67,
                   -96,96,31,55,36,29,-43,-7);
println(deconv(g, f) == h, " ", deconv(g, f));
println(deconv(g, h) == f, " ", deconv(g, h));
Output:
True L(-8,-9,-3,-1,-6,7)
True L(-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1)
Retrieved from "https://rosettacode.org/wiki/Deconvolution/1D?oldid=370736"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information