Deconvolution/1D: Difference between revisions

g = [24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7]

</code>

 

=={{header|BBC BASIC}}==

{{works with|BBC BASIC for Windows}}

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

DIM h(5), f(15), g(20)

h() = -8,-9,-3,-1,-6,7

a$ += STR$(a(i%)) + ", "

NEXT

{{out}}

<pre>

=={{header|C}}==

Using [[FFT]]:

#include <stdlib.h>

#include <math.h>

for (int i = 0; i < lh; i++) printf(" %g", h2[i]);

printf("\n");

{{out}}<pre>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</pre>

 

=={{header|Common Lisp}}==

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

 

(defun make-conv-matrix (x m n)

(let ((lx (cadr (array-dimensions x)))

(A (make-conv-matrix f lg lh)))

 

 

Example:

 

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

(-2.0000000000000004)

(-7.000000000000001)

 

=={{header|D}}==

int flen = f.length;

int glen = g.length;

writeln(deconv(g, f) == h, " ", deconv(g, f));

writeln(deconv(g, h) == f, " ", deconv(g, h));

{{out}}

<pre>true [-8, -9, -3, -1, -6, 7]

=={{header|Fortran}}==

This solution uses the LAPACK95 library.

! Build

! Windows: ifort /I "%IFORT_COMPILER11%\mkl\include\ia32" deconv1d.f90 "%IFORT_COMPILER11%\mkl\ia32\lib\*.lib"

 

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

=={{header|Go}}==

{{trans|D}}

 

import "fmt"

}

return h

{{out}}

<pre>

 

{{trans|C}}

 

import (

y[k], y[k+n/2] = y[k]+tf, y[k]-tf

}

{{out}}

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.

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

</pre>

 

=={{header|Haskell}}==

 

h, f, g :: [Double]

 

 

 

 

=={{header|J}}==

This solution borrowed from [[Formal_power_series#J|Formal power series]]:

 

 

Sample data:

 

f=: _3 _6 _1 8 _6 3 _1 _9 _9 3 _2 5 2 _2 _7 _1

 

Example use:

_8 _9 _3 _1 _6 7

g divide h

 

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:

 

 

=={{header|Java}}==

{{trans|Go}}

 

public class Deconvolution1D {

System.out.println(sb.toString());

}

{{out}}

<pre>

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]

</pre>

 

=={{header|Lua}}==

Using metatables:

local h = setmetatable({}, {__index = function(self, n)

if n == 1 then self[1] = g[1] / f[1]

local _ = h[#g - #f + 1]

return setmetatable(h, nil)

 

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

 

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%2B]]

deconv[f_List, g_List] :=

Module[{A =

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:

<pre>

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:

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

ans =

 

 

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

 

 

 

 

 

}

}

 

}

}

 

 

 

 

 

 

 

{{out}}

<pre>

</pre>

 

=={{header|PicoLisp}}==

 

(de deconv (G F)

(dec 'H

(*/ M (get F (- N I)) 1.0) ) )

Test:

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

 

(test H (deconv G F))

 

=={{header|Python}}==

 

Inspired by the TCL solution, and using the <code>ToReducedRowEchelonForm</code> function to reduce to row echelon form from [[Reduced row echelon form#Python|here]]

if not M: return

lead = 0

M[i] = [ iv - lv*rv for rv,iv in zip(M[r],M[i])]

lead += 1

def pmtx(mtx):

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

 

=={{header|R}}==

* solution is ifft(fft(a)*fft(b)), truncated.

 

p <- length(a)

q <- length(b)

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)

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

 

=={{header|Racket}}==

#lang racket

(require math/matrix)

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

(deconvolve g f)

(deconvolve g h)

{{out}}

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

 

=={{header|Tcl}}==

{{works with|Tcl|8.5}}

This builds the a command, <code>1D</code>, with two subcommands (<code>convolve</code> and <code>deconvolve</code>) for performing convolution and deconvolution of these kinds of arrays. The deconvolution code is based on a reduction to [[Reduced row echelon form#Tcl|reduced row echelon form]].

namespace eval 1D {

namespace ensemble create; # Will be same name as namespace

return $result

}

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

proc pp {name nlist} {

set sep ""

pp "deconv(g,f) = h" [1D deconvolve $g $f]

pp "deconv(g,h) = f" [1D deconvolve $g $h]

{{out}}

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

 

=={{header|Ursala}}==

the same length by appending zeros to the short ones).

 

#import nat

 

 

deconv = lapack..dgelsd^\~&l ~&||0.!**+ band

test program:

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': deconv(g,f),

'f': deconv(g,h)>

{{out}}

<pre>

-7.000000e+00,

-1.000000e+00>>

</pre>

 

=={{header|zkl}}==

{{trans|D}}

flen, glen, delta:=f.len(), g.len(), glen - flen + 1;

result:=List.createLong(delta); // allocate list with space for items

}

result;

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

{{out}}

<pre>