Deconvolution/1D: Difference between revisions

=={{header|11l}}==

{{trans|D}}

V result = [0]*(g.len - f.len + 1)

L(&e) result

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

{{out}}

<pre>

=={{header|Ada}}==

This is a translation of the '''D''' solution.

procedure Main is

print (deconv (g, h));

end Main;

{{output}}

<pre>

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

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.

</pre>

'''Library gonum/mat:'''

import (

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

{{out}}

<pre>

=={{header|Haskell}}==

deconv1d xs ys = takeWhile (/= 0) $ deconv xs ys

where

main :: IO ()

{{Out}}

<pre>True</pre>

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>

Integer inputs may need to be converted and copied to floating point to use deconv().

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]

fanswer = deconv(float.(g), float.(h))

{{output}}

=={{header|Kotlin}}==

{{trans|Go}}

fun deconv(g: DoubleArray, f: DoubleArray): DoubleArray {

println("${f.map { it.toInt() }}")

println("${deconv(g, h).map { it.toInt() }}")

{{out}}

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

=={{header|Mathematica}} / {{header|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%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.

=={{header|Nim}}==

var h: seq[float] = newSeq[float](len(g) - len(f) + 1)

for n in 0..<len(h):

echo deconv(g, f)

echo f

{{out}}

<pre>

Using <code>rref</code> routine from [[Reduced row echelon form#Perl|Reduced row echelon form]] task.

{{trans|Raku}}

sub deconvolve {

print ' conv(f,h) = g = ' . join(' ', my @g = convolve(\@f, \@h)) . "\n";

print 'deconv(g,f) = h = ' . join(' ', deconvolve(\@g, \@f)) . "\n";

{{out}}

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

=={{header|Phix}}==

{{trans|D}}

<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>

<span style="color: #008080;">function</span> <span style="color: #000000;">deconv</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">g</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">)</span>

<span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"deconv(g,f)"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"h"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">deconv</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f</span><span style="color: #0000FF;">),</span><span style="color: #000000;">h</span><span style="color: #0000FF;">)</span>

<span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"deconv(g,h)"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"f"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">deconv</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">h</span><span style="color: #0000FF;">),</span><span style="color: #000000;">f</span><span style="color: #0000FF;">)</span>

{{out}}

<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

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

Based on the R version.

import numpy

print(deconv(g,h))

Output

* 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|Raku}}==

Translation of Python, using a modified version of the subroutine <code>rref</code> from [[Reduced row echelon form#Raku| Reduced row echelon form]] task.

my \h = 1 + @g - @f;

my @m;

.say for ~@g, ~convolve(@f, @h),'';

.say for ~@h, ~deconvolve(@g, @f),'';

{{out}}

=={{header|REXX}}==

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"

if @.$1.t= @.$2.t then iterate /*create array list. */

say "***error*** arrays" $1 ' and ' $2 "aren't equal."

{{out|output|text=&nbsp; when using the default internal inputs:}}

<pre>

=={{header|Scala}}==

{{Out}}Best seen running in your browser either by [https://scalafiddle.io/sf/ENWyl3Z/0 ScalaFiddle (ES aka JavaScript, non JVM)] or [https://scastie.scala-lang.org/bFag8sS1Qr2Z062LN8dr6A Scastie (remote JVM)].

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)

println(sb.result())

=={{header|Swift}}==

{{trans|Kotlin}}

let fs = f.count

var ret = [Double](repeating: 0, count: g.count - fs + 1)

print("\(f.map({ Int($0) }))")

{{out}}

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

<pre>deconv(g,f) = h = [-8,-9,-3,-1,-6,7]

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>

=={{header|Vlang}}==

{{trans|Go}}

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]

}

return h

{{out}}

=={{header|Wren}}==

{{trans|Go}}

var h = List.filled(g.count - f.count + 1, 0)

for (n in 0...h.count) {

System.print(deconv.call(g, f))

System.print(f)

{{out}}

=={{header|zkl}}==

Using GNU Scientific Library:

fcn dconv1D(f,g){

fsz,hsz:=f.len(), g.len() - fsz +1;

h:=A.AxEQb(g);

h

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

f:=dconv1D(h,g);

{{out}}

<pre>

Or, using lists:

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