Deconvolution/1D: Difference between revisions
(→{{header|Ursala}}: Added zkl) |
|||
(42 intermediate revisions by 25 users not shown) | |||
Line 79: | Line 79: | ||
g = [24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7] |
g = [24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7] |
||
</code> |
</code> |
||
=={{header|11l}}== |
|||
{{trans|D}} |
|||
<syntaxhighlight lang="11l">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))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
[-8, -9, -3, -1, -6, 7] |
|||
[-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1] |
|||
</pre> |
|||
=={{header|Ada}}== |
|||
This is a translation of the '''D''' solution. |
|||
<syntaxhighlight lang="ada">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; |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
(-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) |
|||
</pre> |
|||
=={{header|BBC BASIC}}== |
=={{header|BBC BASIC}}== |
||
{{works with|BBC BASIC for Windows}} |
{{works with|BBC BASIC for Windows}} |
||
As several others, this is a translation of the '''D''' solution. |
As several others, this is a translation of the '''D''' solution. |
||
< |
<syntaxhighlight lang="bbcbasic"> *FLOAT 64 |
||
DIM h(5), f(15), g(20) |
DIM h(5), f(15), g(20) |
||
h() = -8,-9,-3,-1,-6,7 |
h() = -8,-9,-3,-1,-6,7 |
||
Line 120: | Line 208: | ||
a$ += STR$(a(i%)) + ", " |
a$ += STR$(a(i%)) + ", " |
||
NEXT |
NEXT |
||
= LEFT$(LEFT$(a$))</ |
= LEFT$(LEFT$(a$))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 129: | Line 217: | ||
=={{header|C}}== |
=={{header|C}}== |
||
Using [[FFT]]: |
Using [[FFT]]: |
||
< |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include <stdlib.h> |
#include <stdlib.h> |
||
#include <math.h> |
#include <math.h> |
||
Line 220: | Line 308: | ||
for (int i = 0; i < lh; i++) printf(" %g", h2[i]); |
for (int i = 0; i < lh; i++) printf(" %g", h2[i]); |
||
printf("\n"); |
printf("\n"); |
||
}</ |
}</syntaxhighlight> |
||
{{out}}<pre>f[] data is : -3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1 |
{{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 |
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 |
h[] data is : -8 -9 -3 -1 -6 7 |
||
deconv(g, f): -8 -9 -3 -1 -6 7</pre> |
deconv(g, f): -8 -9 -3 -1 -6 7</pre> |
||
=={{header|C++}}== |
|||
<syntaxhighlight lang="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)); |
|||
} |
|||
</syntaxhighlight> |
|||
{{ out }} |
|||
<pre> |
|||
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] |
|||
</pre> |
|||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |
||
Uses the routine (lsqr A b) from [[Multiple regression]] and (mtp A) from [[Matrix transposition]]. |
Uses the routine (lsqr A b) from [[Multiple regression]] and (mtp A) from [[Matrix transposition]]. |
||
< |
<syntaxhighlight lang="lisp">;; Assemble the mxn matrix A from the 2D row vector x. |
||
(defun make-conv-matrix (x m n) |
(defun make-conv-matrix (x m n) |
||
(let ((lx (cadr (array-dimensions x))) |
(let ((lx (cadr (array-dimensions x))) |
||
Line 250: | Line 387: | ||
(A (make-conv-matrix f lg lh))) |
(A (make-conv-matrix f lg lh))) |
||
(lsqr A (mtp g))))</ |
(lsqr A (mtp g))))</syntaxhighlight> |
||
Example: |
Example: |
||
< |
<syntaxhighlight lang="lisp">(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 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))) |
(setf g #2A((24 75 71 -34 3 22 -45 23 245 25 52 25 -67 -96 96 31 55 36 29 -43 -7))) |
||
Line 282: | Line 419: | ||
(-2.0000000000000004) |
(-2.0000000000000004) |
||
(-7.000000000000001) |
(-7.000000000000001) |
||
(-0.9999999999999994))</ |
(-0.9999999999999994))</syntaxhighlight> |
||
=={{header|D}}== |
=={{header|D}}== |
||
< |
<syntaxhighlight lang="d">T[] deconv(T)(in T[] g, in T[] f) pure nothrow { |
||
int flen = f.length; |
int flen = f.length; |
||
int glen = g.length; |
int glen = g.length; |
||
Line 307: | Line 444: | ||
writeln(deconv(g, f) == h, " ", deconv(g, f)); |
writeln(deconv(g, f) == h, " ", deconv(g, f)); |
||
writeln(deconv(g, h) == f, " ", deconv(g, h)); |
writeln(deconv(g, h) == f, " ", deconv(g, h)); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>true [-8, -9, -3, -1, -6, 7] |
<pre>true [-8, -9, -3, -1, -6, 7] |
||
Line 314: | Line 451: | ||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
This solution uses the LAPACK95 library. |
This solution uses the LAPACK95 library. |
||
< |
<syntaxhighlight lang="fortran"> |
||
! Build |
! Build |
||
! Windows: ifort /I "%IFORT_COMPILER11%\mkl\include\ia32" deconv1d.f90 "%IFORT_COMPILER11%\mkl\ia32\lib\*.lib" |
! Windows: ifort /I "%IFORT_COMPILER11%\mkl\include\ia32" deconv1d.f90 "%IFORT_COMPILER11%\mkl\ia32\lib\*.lib" |
||
Line 388: | Line 525: | ||
end program deconv |
end program deconv |
||
</syntaxhighlight> |
|||
</lang> |
|||
Results: |
Results: |
||
< |
<syntaxhighlight lang="fortran"> |
||
deconv(f, g) = -8, -9, -3, -1, -6, 7 |
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 |
deconv(h, g) = -3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1 |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|FreeBASIC}}== |
|||
<syntaxhighlight lang="vbnet">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</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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]</pre> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
{{trans|D}} |
{{trans|D}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import "fmt" |
import "fmt" |
||
Line 425: | Line 625: | ||
} |
} |
||
return h |
return h |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 435: | Line 635: | ||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 499: | Line 699: | ||
y[k], y[k+n/2] = y[k]+tf, y[k]-tf |
y[k], y[k+n/2] = y[k]+tf, y[k]-tf |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{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. |
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. |
||
Line 507: | Line 707: | ||
[-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] |
||
[-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> |
|||
'''Library gonum/mat:''' |
|||
<syntaxhighlight lang="go">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))) |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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⎦ |
|||
</pre> |
</pre> |
||
=={{header|Haskell}}== |
=={{header|Haskell}}== |
||
<syntaxhighlight lang="haskell">deconv1d :: [Double] -> [Double] -> [Double] |
|||
<lang haskell>import Data.List |
|||
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, f, g :: [Double] |
||
h = [-8,-9,-3,-1,-6,7] |
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] |
|||
f = [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1] |
|||
scale x ys = map (x*) ys |
|||
deconv1d :: (Fractional a) => [a] -> [a] -> [a] |
|||
deconv1d xs ys = takeWhile (/=0) $ deconv xs ys |
|||
where [] `deconv` _ = [] |
|||
(0:xs) `deconv` (0:ys) = xs `deconv` ys |
|||
(x:xs) `deconv` (y:ys) = |
|||
q : zipWith (-) xs (scale q ys ++ repeat 0) `deconv` (y:ys) |
|||
where q = x / y</lang> |
|||
g = |
|||
Check: |
|||
[ 24 |
|||
<lang haskell>*Main> h == deconv1d g f |
|||
, 75 |
|||
True |
|||
, 71 |
|||
, -34 |
|||
, 3 |
|||
, 22 |
|||
, -45 |
|||
, 23 |
|||
, 245 |
|||
, 25 |
|||
, 52 |
|||
, 25 |
|||
, -67 |
|||
, -96 |
|||
, 96 |
|||
, 31 |
|||
, 55 |
|||
, 36 |
|||
, 29 |
|||
, -43 |
|||
, -7 |
|||
] |
|||
main :: IO () |
|||
*Main> f == deconv1d g h |
|||
main = print $ (h == deconv1d g f) && (f == deconv1d g h)</syntaxhighlight> |
|||
True</lang> |
|||
{{Out}} |
|||
<pre>True</pre> |
|||
=={{header|J}}== |
=={{header|J}}== |
||
Line 538: | Line 825: | ||
This solution borrowed from [[Formal_power_series#J|Formal power series]]: |
This solution borrowed from [[Formal_power_series#J|Formal power series]]: |
||
< |
<syntaxhighlight lang="j">Ai=: (i.@] =/ i.@[ -/ i.@>:@-)&# |
||
divide=: [ +/ .*~ [:%.&.x: ] +/ .* Ai</ |
divide=: [ +/ .*~ [:%.&.x: ] +/ .* Ai</syntaxhighlight> |
||
Sample data: |
Sample data: |
||
< |
<syntaxhighlight lang="j">h=: _8 _9 _3 _1 _6 7 |
||
f=: _3 _6 _1 8 _6 3 _1 _9 _9 3 _2 5 2 _2 _7 _1 |
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</ |
g=: 24 75 71 _34 3 22 _45 23 245 25 52 25 _67 _96 96 31 55 36 29</syntaxhighlight> |
||
Example use: |
Example use: |
||
< |
<syntaxhighlight lang="j"> g divide f |
||
_8 _9 _3 _1 _6 7 |
_8 _9 _3 _1 _6 7 |
||
g divide h |
g divide h |
||
_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</syntaxhighlight> |
||
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: |
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: |
||
< |
<syntaxhighlight lang="j">divide=: [ +/ .*~ [:%. ] +/ .* Ai</syntaxhighlight> |
||
=={{header|Java}}== |
=={{header|Java}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
< |
<syntaxhighlight lang="java">import java.util.Arrays; |
||
public class Deconvolution1D { |
public class Deconvolution1D { |
||
Line 587: | Line 874: | ||
System.out.println(sb.toString()); |
System.out.println(sb.toString()); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 594: | Line 881: | ||
f = [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1] |
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] |
deconv(g, h) = [-3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1] |
||
</pre> |
|||
=={{header|jq}}== |
|||
{{trans|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''' |
|||
<syntaxhighlight lang="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) |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
[-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] |
|||
</pre> |
|||
=={{header|Julia}}== |
|||
The deconv function for floating point data is built into Julia, though <code>using DSP</code> is required with version 1.0. |
|||
Integer inputs may need to be converted and copied to floating point to use deconv(). |
|||
<syntaxhighlight lang="julia">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")</syntaxhighlight> |
|||
{{output}} |
|||
<pre>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]</pre> |
|||
=={{header|Kotlin}}== |
|||
{{trans|Go}} |
|||
<syntaxhighlight lang="scala">// 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() }}") |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
[-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] |
|||
</pre> |
</pre> |
||
=={{header|Lua}}== |
=={{header|Lua}}== |
||
Using metatables: |
Using metatables: |
||
< |
<syntaxhighlight lang="lua">function deconvolve(f, g) |
||
local h = setmetatable({}, {__index = function(self, n) |
local h = setmetatable({}, {__index = function(self, n) |
||
if n == 1 then self[1] = g[1] / f[1] |
if n == 1 then self[1] = g[1] / f[1] |
||
Line 612: | Line 999: | ||
local _ = h[#g - #f + 1] |
local _ = h[#g - #f + 1] |
||
return setmetatable(h, nil) |
return setmetatable(h, nil) |
||
end</ |
end</syntaxhighlight> |
||
Tests: |
Tests: |
||
< |
<syntaxhighlight lang="lua"> |
||
local f = {-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1} |
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 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} |
local h = {-8,-9,-3,-1,-6,7} |
||
print(unpack(deconvolve(f, g))) --> -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</ |
print(unpack(deconvolve(h, g))) --> -3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1</syntaxhighlight> |
||
=={{header|Mathematica}}== |
=={{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]] |
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]] |
||
<syntaxhighlight lang="mathematica"> |
|||
<lang Mathematica> |
|||
deconv[f_List, g_List] := |
deconv[f_List, g_List] := |
||
Module[{A = |
Module[{A = |
||
Line 630: | Line 1,017: | ||
Table[Band[{n, 1}] -> f[[n]], {n, 1, Length[f]}], {Length[g], Length[f] - 1}]}, |
Table[Band[{n, 1}] -> f[[n]], {n, 1, Length[f]}], {Length[g], Length[f] - 1}]}, |
||
Take[LinearSolve[A, g], Length[g] - Length[f] + 1]] |
Take[LinearSolve[A, g], Length[g] - Length[f] + 1]] |
||
</syntaxhighlight> |
|||
</lang> |
|||
Usage: |
Usage: |
||
<pre> |
<pre> |
||
Line 643: | Line 1,030: | ||
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". |
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: |
To test that this operates according to the task spec we can test the criteria above: |
||
< |
<syntaxhighlight lang="matlab">>> 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]; |
>> 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]; |
>> f = [-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1]; |
||
Line 656: | Line 1,043: | ||
ans = |
ans = |
||
-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</syntaxhighlight> |
||
Therefore, "deconv(a,b)" behaves as expected. |
Therefore, "deconv(a,b)" behaves as expected. |
||
=={{header| |
=={{header|Nim}}== |
||
<syntaxhighlight lang="nim">proc deconv(g, f: openArray[float]): seq[float] = |
|||
{{works with|Rakudo 2010.12}} |
|||
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] |
|||
Translation of Python, using a modified version of the Reduced Row Echelon Form subroutine <code>rref()</code> from [[Reduced row echelon form#Perl 6|here]]. |
|||
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)</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
[-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] |
|||
</pre> |
|||
=={{header|Perl}}== |
|||
<lang perl6>sub deconvolve (@g, @f) { |
|||
Using <code>rref</code> routine from [[Reduced row echelon form#Perl|Reduced row echelon form]] task. |
|||
my $h = 1 + @g - @f; |
|||
{{trans|Raku}} |
|||
my @m; |
|||
<syntaxhighlight lang="perl">use v5.36; |
|||
@m[^@g]>>.[^$h] >>+=>> 0; |
|||
use Math::Cartesian::Product; |
|||
@m[^@g]>>.[$h] >>=<< @g; |
|||
for ^$h -> $j { for @f.kv -> $k, $v { @m[$j + $k][$j] = $v } } |
|||
return rref( @m )[^$h]>>.[$h]; |
|||
} |
|||
sub |
sub deconvolve($g,$f) { |
||
my @g = |
my @g = @{$g}; |
||
my @f = @{$f}; |
|||
my(@m,@d); |
|||
} |
|||
my $h = 1 + @g - @f; |
|||
# Reduced Row Echelon Form simultaneous equation solver. |
|||
push @m, [(0) x $h, $g[$_]] for 0..$#g; |
|||
# Can handle over-specified systems of equations. |
|||
for my $j (0..$h-1) { |
|||
for my $k (0..$#f) { |
|||
sub rref ($m is rw) { |
|||
$m[$j + $k][$j] = $f[$k] |
|||
return unless $m; |
|||
my ($lead, $rows, $cols) = 0, +$m, +$m[0]; |
|||
# Trim off over specified rows if they exist. |
|||
# Not strictly necessary, but can save a lot of |
|||
# redundant calculations. |
|||
if $rows >= $cols { |
|||
$m = trim_system($m); |
|||
$rows = +$m; |
|||
} |
|||
for ^$rows -> $r { |
|||
$lead < $cols or return $m; |
|||
my $i = $r; |
|||
until $m[$i][$lead] { |
|||
++$i == $rows or next; |
|||
$i = $r; |
|||
++$lead == $cols and return $m; |
|||
} |
} |
||
$m[$i, $r] = $m[$r, $i] if $r != $i; |
|||
my $lv = $m[$r][$lead]; |
|||
$m[$r] >>/=>> $lv; |
|||
for ^$rows -> $n { |
|||
next if $n == $r; |
|||
$m[$n] >>-=>> $m[$r] >>*>> $m[$n][$lead]; |
|||
} |
|||
++$lead; |
|||
} |
} |
||
rref(\@m); |
|||
push @d, @{ $m[$_] }[$h] for 0..$h-1; |
|||
@d; |
|||
} |
|||
sub convolve($f,$h) { |
|||
# Reduce a system of equations to n equations with n unknowns. |
|||
my @f = @{$f}; |
|||
# Looks for an equation with a true value for each position. |
|||
my @h = @{$h}; |
|||
# If it can't find one, assumes that it has already taken one |
|||
my @i; |
|||
# and pushes in the first equation it sees. This assumtion |
|||
for my $x (cartesian {@_} [0..$#f], [0..$#h]) { |
|||
# will alway be successful except in some cases where an |
|||
push @i, @$x[0]+@$x[1]; |
|||
# under-specified system has been supplied, in which case, |
|||
# it would not have been able to reduce the system anyway. |
|||
sub trim_system ($m is rw) { |
|||
my ($vars, @t) = +$m[0]-1, (); |
|||
for ^$vars -> $lead { |
|||
for ^$m -> $row { |
|||
@t.push( $m.splice( $row, 1 ) ) and last if $m[$row][$lead]; |
|||
} |
|||
} |
|||
while (+@t < $vars) and +$m { @t.push( $m.splice( 0, 1 ) ) }; |
|||
return @t; |
|||
} |
} |
||
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) { |
|||
my @h = (-8,-9,-3,-1,-6,7); |
|||
$lead < $cols or return; |
|||
my @f = (-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1); |
|||
my $i = $r; |
|||
my @g = (24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7); |
|||
until ($m[$i][$lead]) { |
|||
++$i == $rows or next; |
|||
$i = $r; |
|||
++$lead == $cols and return; |
|||
} |
|||
@m[$i, $r] = @m[$r, $i]; |
|||
.say for ~@g, ~convolve(@f, @h),''; |
|||
my $lv = $m[$r][$lead]; |
|||
$_ /= $lv foreach @{ $m[$r] }; |
|||
my @mr = @{ $m[$r] }; |
|||
.say for ~@h, ~deconvolve(@g, @f),''; |
|||
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>; |
|||
.say for ~@f, ~deconvolve(@g, @h),'';</lang> |
|||
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";</syntaxhighlight> |
|||
{{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 |
|||
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</pre> |
|||
=={{header|Phix}}== |
|||
{{trans|D}} |
|||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<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: #004080;">integer</span> <span style="color: #000000;">lf</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">lg</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">lh</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lg</span><span style="color: #0000FF;">-</span><span style="color: #000000;">lf</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lh</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">lh</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">e</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">g</span><span style="color: #0000FF;">[</span><span style="color: #000000;">n</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">lf</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">2</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">e</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">h</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">[</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #000000;">h</span><span style="color: #0000FF;">[</span><span style="color: #000000;">n</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">e</span><span style="color: #0000FF;">/</span><span style="color: #000000;">f</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">h</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">conv</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</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: #004080;">integer</span> <span style="color: #000000;">lf</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">lh</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">h</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">lg</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lf</span><span style="color: #0000FF;">+</span><span style="color: #000000;">lh</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">g</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lg</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">lh</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">lf</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">j</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> |
|||
<span style="color: #000000;">g</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">h</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">g</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span><span style="color: #000000;">7</span><span style="color: #0000FF;">},</span> |
|||
<span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span> |
|||
<span style="color: #000000;">g</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">24</span><span style="color: #0000FF;">,</span><span style="color: #000000;">75</span><span style="color: #0000FF;">,</span><span style="color: #000000;">71</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">34</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">22</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">45</span><span style="color: #0000FF;">,</span><span style="color: #000000;">23</span><span style="color: #0000FF;">,</span><span style="color: #000000;">245</span><span style="color: #0000FF;">,</span><span style="color: #000000;">25</span><span style="color: #0000FF;">,</span><span style="color: #000000;">52</span><span style="color: #0000FF;">,</span><span style="color: #000000;">25</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">67</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">96</span><span style="color: #0000FF;">,</span><span style="color: #000000;">96</span><span style="color: #0000FF;">,</span><span style="color: #000000;">31</span><span style="color: #0000FF;">,</span><span style="color: #000000;">55</span><span style="color: #0000FF;">,</span><span style="color: #000000;">36</span><span style="color: #0000FF;">,</span><span style="color: #000000;">29</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">43</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">7</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #004080;">string</span> <span style="color: #000000;">desc</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">eq</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">e</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s (%ssame as %s): %V\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">desc</span><span style="color: #0000FF;">,</span><span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">==</span><span style="color: #000000;">e</span><span style="color: #0000FF;">?</span><span style="color: #008000;">""</span><span style="color: #0000FF;">:</span><span style="color: #008000;">"**NOT** "</span><span style="color: #0000FF;">),</span><span style="color: #000000;">eq</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
<span style="color: #000000;">test</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" conv(h,f)"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"g"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">conv</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><span style="color: #000000;">g</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> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
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} |
|||
24 75 71 -34 3 22 -45 23 245 25 52 25 -67 -96 96 31 55 36 29 -43 -7 |
|||
deconv(g,h) (same as f): {-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1} |
|||
-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 |
|||
</pre> |
</pre> |
||
=={{header|PicoLisp}}== |
=={{header|PicoLisp}}== |
||
< |
<syntaxhighlight lang="picolisp">(load "@lib/math.l") |
||
(de deconv (G F) |
(de deconv (G F) |
||
Line 767: | Line 1,211: | ||
(dec 'H |
(dec 'H |
||
(*/ M (get F (- N I)) 1.0) ) ) |
(*/ M (get F (- N I)) 1.0) ) ) |
||
(link (*/ H 1.0 A)) ) ) ) )</ |
(link (*/ H 1.0 A)) ) ) ) )</syntaxhighlight> |
||
Test: |
Test: |
||
< |
<syntaxhighlight lang="picolisp">(setq |
||
F (-3. -6. -1. 8. -6. 3. -1. -9. -9. 3. -2. 5. 2. -2. -7. -1.) |
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.) |
G (24. 75. 71. -34. 3. 22. -45. 23. 245. 25. 52. 25. -67. -96. 96. 31. 55. 36. 29. -43. -7.) |
||
Line 775: | Line 1,219: | ||
(test H (deconv G F)) |
(test H (deconv G F)) |
||
(test F (deconv G H))</ |
(test F (deconv G H))</syntaxhighlight> |
||
=={{header|Python}}== |
=={{header|Python}}== |
||
Line 781: | Line 1,225: | ||
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]] |
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]] |
||
< |
<syntaxhighlight lang="python">def ToReducedRowEchelonForm( M ): |
||
if not M: return |
if not M: return |
||
lead = 0 |
lead = 0 |
||
Line 805: | Line 1,249: | ||
M[i] = [ iv - lv*rv for rv,iv in zip(M[r],M[i])] |
M[i] = [ iv - lv*rv for rv,iv in zip(M[r],M[i])] |
||
lead += 1 |
lead += 1 |
||
return M |
|||
def pmtx(mtx): |
def pmtx(mtx): |
||
Line 833: | Line 1,278: | ||
g = [24,75,71,-34,3,22,-45,23,245,25,52,25,-67,-96,96,31,55,36,29,-43,-7] |
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 convolve(f,h) == g |
||
assert deconvolve(g, f) == h</ |
assert deconvolve(g, f) == h</syntaxhighlight> |
||
Based on the R version. |
|||
<syntaxhighlight lang="python"> |
|||
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)) |
|||
</syntaxhighlight> |
|||
Output |
|||
<pre> |
|||
[ 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.] |
|||
</pre> |
|||
=={{header|R}}== |
=={{header|R}}== |
||
Line 842: | Line 1,341: | ||
* solution is ifft(fft(a)*fft(b)), truncated. |
* solution is ifft(fft(a)*fft(b)), truncated. |
||
< |
<syntaxhighlight lang="r">conv <- function(a, b) { |
||
p <- length(a) |
p <- length(a) |
||
q <- length(b) |
q <- length(b) |
||
Line 859: | Line 1,358: | ||
return(y[1:n]) |
return(y[1:n]) |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
To check : |
To check : |
||
<syntaxhighlight lang="r"> |
|||
<lang R> |
|||
h <- c(-8,-9,-3,-1,-6,7) |
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) |
f <- c(-3,-6,-1,8,-6,3,-1,-9,-9,3,-2,5,2,-2,-7,-1) |
||
Line 871: | Line 1,370: | ||
max(abs(deconv(g,f) - h)) |
max(abs(deconv(g,f) - h)) |
||
max(abs(deconv(g,h) - f)) |
max(abs(deconv(g,h) - f)) |
||
</syntaxhighlight> |
|||
</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. |
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. |
||
Line 877: | Line 1,376: | ||
R has also a function convolve, |
R has also a function convolve, |
||
<syntaxhighlight lang="r"> |
|||
<lang R> |
|||
conv(a, b) == convolve(a, rev(b), type="open") |
conv(a, b) == convolve(a, rev(b), type="open") |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Racket}}== |
=={{header|Racket}}== |
||
< |
<syntaxhighlight lang="racket"> |
||
#lang racket |
#lang racket |
||
(require math/matrix) |
(require math/matrix) |
||
Line 901: | Line 1,400: | ||
(define lh (+ (- lg lf) 1)) |
(define lh (+ (- lg lf) 1)) |
||
(least-square (convolution-matrix f lg lh) g)) |
(least-square (convolution-matrix f lg lh) g)) |
||
</syntaxhighlight> |
|||
</lang> |
|||
Test: |
Test: |
||
< |
<syntaxhighlight lang="racket"> |
||
(define f (col-matrix [-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1])) |
(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 h (col-matrix [-8 -9 -3 -1 -6 7])) |
||
Line 910: | Line 1,409: | ||
(deconvolve g f) |
(deconvolve g f) |
||
(deconvolve g h) |
(deconvolve g h) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
< |
<syntaxhighlight lang="racket"> |
||
#<array '#(6 1) #[-8 -9 -3 -1 -6 7]> |
#<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]> |
#<array '#(16 1) #[-3 -6 -1 8 -6 3 -1 -9 -9 3 -2 5 2 -2 -7 -1]> |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Raku}}== |
|||
(formerly Perl 6) |
|||
Translation of Python, using a modified version of the subroutine <code>rref</code> from [[Reduced row echelon form#Raku| Reduced row echelon form]] task. |
|||
<syntaxhighlight lang="raku" line>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),'';</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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</pre> |
|||
=={{header|REXX}}== |
|||
<syntaxhighlight lang="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. */</syntaxhighlight> |
|||
{{out|output|text= when using the default internal inputs:}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|RPL}}== |
|||
{{trans|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''' |
|||
≫ ≫ '<span style="color:blue">DECONV</span>' 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 <span style="color:blue">DECONV</span> h == |
|||
g h <span style="color:blue">DECONV</span> f == AND |
|||
≫ ≫ ‘<span style="color:blue">TASK</span>’ STO |
|||
{{out}} |
|||
<pre> |
|||
1: 1 |
|||
</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)]. |
|||
<syntaxhighlight lang="scala">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()) |
|||
}</syntaxhighlight> |
|||
=={{header|Swift}}== |
|||
{{trans|Kotlin}} |
|||
<syntaxhighlight lang="swift">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) }))")</syntaxhighlight> |
|||
{{out}} |
|||
<pre>[-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]</pre> |
|||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
||
{{works with|Tcl|8.5}} |
{{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]]. |
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]]. |
||
< |
<syntaxhighlight lang="tcl">package require Tcl 8.5 |
||
namespace eval 1D { |
namespace eval 1D { |
||
namespace ensemble create; # Will be same name as namespace |
namespace ensemble create; # Will be same name as namespace |
||
Line 1,020: | Line 1,740: | ||
return $result |
return $result |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
To use the above code, a simple demonstration driver (which solves the specific task): |
To use the above code, a simple demonstration driver (which solves the specific task): |
||
< |
<syntaxhighlight lang="tcl"># Simple pretty-printer |
||
proc pp {name nlist} { |
proc pp {name nlist} { |
||
set sep "" |
set sep "" |
||
Line 1,039: | Line 1,759: | ||
pp "deconv(g,f) = h" [1D deconvolve $g $f] |
pp "deconv(g,f) = h" [1D deconvolve $g $f] |
||
pp "deconv(g,h) = f" [1D deconvolve $g $h] |
pp "deconv(g,h) = f" [1D deconvolve $g $h] |
||
pp " conv(f,h) = g" [1D convolve $f $h]</ |
pp " conv(f,h) = g" [1D convolve $f $h]</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>deconv(g,f) = h = [-8,-9,-3,-1,-6,7] |
|||
<pre> |
|||
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] |
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] |
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]</pre> |
||
</pre> |
|||
=={{header|Ursala}}== |
=={{header|Ursala}}== |
||
Line 1,055: | Line 1,773: | ||
the same length by appending zeros to the short ones). |
the same length by appending zeros to the short ones). |
||
< |
<syntaxhighlight lang="ursala">#import std |
||
#import nat |
#import nat |
||
Line 1,061: | Line 1,779: | ||
deconv = lapack..dgelsd^\~&l ~&||0.!**+ band |
deconv = lapack..dgelsd^\~&l ~&||0.!**+ band |
||
</syntaxhighlight> |
|||
</lang> |
|||
test program: |
test program: |
||
< |
<syntaxhighlight 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.> |
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.> |
g = <24.,75.,71.,-34.,3.,22.,-45.,23.,245.,25.,52.,25.,-67.,-96.,96.,31.,55.,36.,29.,-43.,-7.> |
||
Line 1,074: | Line 1,792: | ||
'h': deconv(g,f), |
'h': deconv(g,f), |
||
'f': deconv(g,h)> |
'f': deconv(g,h)> |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,102: | Line 1,820: | ||
-7.000000e+00, |
-7.000000e+00, |
||
-1.000000e+00>> |
-1.000000e+00>> |
||
</pre> |
|||
=={{header|V (Vlang)}}== |
|||
{{trans|Go}} |
|||
<syntaxhighlight lang="v (vlang)">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 |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
[-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] |
|||
</pre> |
|||
=={{header|Wren}}== |
|||
{{trans|Go}} |
|||
<syntaxhighlight lang="wren">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))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
[-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] |
|||
</pre> |
</pre> |
||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
Using GNU Scientific Library: |
|||
<syntaxhighlight lang="zkl">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 |
|||
}</syntaxhighlight> |
|||
<syntaxhighlight lang="zkl">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();</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
-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 |
|||
</pre> |
|||
Or, using lists: |
|||
{{trans|D}} |
{{trans|D}} |
||
< |
<syntaxhighlight lang="zkl">fcn deconv(g,f){ |
||
flen, glen, delta:=f.len(), g.len(), glen - flen + 1; |
flen, glen, delta:=f.len(), g.len(), glen - flen + 1; |
||
result:=List.createLong(delta); // allocate list with space for items |
result:=List.createLong(delta); // allocate list with space for items |
||
Line 1,116: | Line 1,926: | ||
} |
} |
||
result; |
result; |
||
}</ |
}</syntaxhighlight> |
||
< |
<syntaxhighlight lang="zkl">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); |
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, |
g:=T(24,75,71,-34,3,22,-45,23,245,25,52,25,-67, |
||
-96,96,31,55,36,29,-43,-7); |
-96,96,31,55,36,29,-43,-7); |
||
println(deconv(g, f) == h, " ", deconv(g, f)); |
println(deconv(g, f) == h, " ", deconv(g, f)); |
||
println(deconv(g, h) == f, " ", deconv(g, h));</ |
println(deconv(g, h) == f, " ", deconv(g, h));</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Latest revision as of 08:12, 22 April 2024
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
anddeconv
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
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
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 -1deconv(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]
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
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]
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
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
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
// 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]
Perl
Using rref
routine from Reduced row echelon form task.
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
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
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
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
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
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/~<K33tx 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)
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
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:
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)