Deconvolution/2D+: Difference between revisions
m
syntax highlighting fixup automation
SqrtNegInf (talk | contribs) m (→{{header|Raku}}: sync rref/trim-system with other tasks, general tidying) |
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
||
Line 82:
=={{header|C}}==
Very tedious code: unpacks 2D or 3D matrix into a vector with padding, do 1D FFT, then pack result back into matrix.
<
#include <stdlib.h>
#include <math.h>
Line 299:
printf("\n");
}
}*/</
-6 -8 -5 9
-7 9 -6 -8
Line 316:
8 5 8
-2 -6 -4 </
=={{header|D}}==
<
class M(T) {
Line 509:
writeln(" f = ", f);
writeln("g deconv h = ", g.deconvolute(h));
}</
''todo(may be not :): pretty print & convert to normal D array''
{{out}}
Line 520:
=={{header|Go}}==
{{trans|C}}
<
import (
Line 711:
}
}
}</
{{out}}
Line 740:
Actually it is a matter of setting up the linear equations and then solving them.
'''Implementation''': <
sz =. x >:@-&$ y NB. shape of z
poi =. ,<"1 ($y) ,"0/&(,@i.) sz NB. pair of indexes
Line 748:
sz $ 1e_8 round ({:"1 %. }:"1) T1
)
round=: [ * <.@%~</
'''Data''': <
f1=: 6 _9 _7 _5
g1=: _48 84 _16 95 125 _70 7 29 54 10
Line 806:
_42 _31 _103 _30 _23 _8
6 4 _26 _10 26 12
)</
'''Tests''': <
1
h2 -: g2 deconv3 f2
1
h3 -: g3 deconv3 f3 NB. -: checks for matching structure and data
1</
=={{header|Julia}}==
Julia has a deconv() function that works on Julia's builtin multidimensional arrays, but not on the nested type 2D and 3D arrays used in the task. So, the solution function, deconvn(), sets up repackaging for 1D fft. The actual solving work is done on one line of ifft/fft, and the rest of the code is merely to repackage the nested arrays.
<
const h1 = [-8, 2, -9, -2, 9, -8, -2]
Line 908:
println(deconvn(f3nested, g3nested,
(length(h3nested), length(h3nested[1]), length(h3nested[1][1])))) # 3D
</
<pre>
Array{Array{Int64,1},1}[[[-6, -8, -5, 9], [-7, 9, -6, -8], [2, -7, 9, 8]], [[7, 4, 4, -6], [9, 9, 4, -4], [-3, 7, -2, -3]]]
Line 914:
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<
Round[ListDeconvolve[{{-5, 2, -2, -6, -7}, {9, 7, -6, 5, -7}, {1, -1, 9, 2, -7}, {5, 9, -9, 2, -5}, {-8, 5, -2, 8, 5}},
Line 927:
{-19, 29, 35, -148, -11, 45}}, {{-55, -147, -146, -31, 55, 60}, {-88, -45, -28, 46, -26, -144},
{-12, -107, -34, 150, 249, 66}, {11, -15, -34, 27, -78, -50}}, {{56, 67, 108, 4, 2, -48}, {58, 67, 89, 32, 32, -8},
{-42, -31, -103, -30, -23, -8}, {6, 4, -26, -10, 26, 12}}}, Method -> "Wiener"]]</
{{out}}
Line 945:
=={{header|Nim}}==
{{trans|D}}
<
type Size = uint64
Line 1,155:
echo "f == g deconv h ? ", f == g.deconvolute(h)
echo " f = ", f
echo "g deconv f = ", g.deconvolute(h)</
{{out}}
Line 1,167:
{{libheader|ntheory}}
{{trans|Raku}}
<
use ntheory qw/forsetproduct/;
Line 1,317:
pretty_print(0,@h);
print "\nff =\n";
pretty_print(0,@ff);</
{{out}}
<pre>3D arrays:
Line 1,354:
Quite frankly I'm fairly astonished that it actually works...<br>
(be warned this contains an exciting mix of 0- and 1- based indexes)
<!--<
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Deconvolution.exw</span>
<span style="color: #008080;">with</span> <span style="color: #000000;">javascript_semantics</span>
Line 1,572:
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">f3</span><span style="color: #0000FF;">))</span>
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deconvolve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">h3</span><span style="color: #0000FF;">))</span>
<!--</
{{out}}
<pre>
Line 1,600:
https://math.stackexchange.com/questions/380720/is-deconvolution-simply-division-in-frequency-domain
<
"""
Line 1,731:
pprint.pprint(deconv(g,f))
</syntaxhighlight>
Output:
Line 1,750:
Translation of Tcl.
<syntaxhighlight lang="raku"
sub deconvolve-N ( @g, @f ) {
my @hsize = @g.shape »-« @f.shape »+» 1;
Line 1,888:
my @ff = deconvolve-N( @g, @h-shaped );
say "\nff =";
pretty-print( @ff );</
Output:
Line 1,925:
The trick to doing this (without using a library to do all the legwork for you) is to recast the higher-order solutions into solutions in the 1D case. This is done by regarding an ''n''-dimensional address as a coding of a 1-D address.
<
namespace path {::tcl::mathfunc ::tcl::mathop}
Line 2,089:
return $h
}</
Demonstrating how to use for the 3-D case:
<
proc pretty matrix {
set size [rank $matrix]
Line 2,142:
# Now do the deconvolution and print it out
puts h:\ [pretty [deconvolve $g $f]]</
Output:
<pre>
Line 2,158:
the [http://www.netlib.org/lapack/lug/node27.html <code>dgelsd</code>] function from the Lapack library.
The <code>break</code> function breaks a long list into a sequence of sublists according to a given template, and the <code>band</code> function is taken from the [[Deconvolution/1D]] solution.
<
#import nat
Line 2,171:
lapack..dgelsd^^(
(~&||0.!**+ ~&B^?a\~&Y@a ^lriFhNSS2iDrlYSK7LS2SL2rQ/~&alt band@alh2faltPrDPMX)^|\~&+ gang,
@t =>~&l ~&L+@r))</
The equations tend to become increasingly sparse in higher dimensions,
so the following alternative implementation uses the sparse matrix
Line 2,177:
instead of Lapack, which is also callable in Ursala, adjusted as shown for the different
[http://www.basis.uklinux.net/avram/refman/umf-input-parameters.html calling convention].
<
~&?\math..div! iota; ~&!*; @h|\; -+
Line 2,189:
~&rFS+ num*rSS+ zipt^*D/~&r ^lrrSPT/~<K33tx zipt^/~&r ~&lSNyCK33+ zipp0,
^/~&rx ~&B->NlNSPC ~&bt+-+-+-,
@t =>~&l ~&L+@r)+-</
UMFPACK doesn't solve systems with more equations than unknowns, so
the system is pruned to a square matrix by first selecting an equation
Line 2,200:
However, some improvement may be possible by averaging the results over several runs.
Here is a test program.
<
f = <<<-9.,5.,-8.>,<3.,5.,1.>>,<<-1.,-7.,2.>,<-5.,-6.,6.>>,<<8.,5.,8.>,<-2.,-6.,-4.>>>
Line 2,234:
'h': deconv3(g,f),
'f': deconv3(g,h)>
</syntaxhighlight>
output:
<pre style="height:25ex;overflow:scroll">
Line 2,287:
{{libheader|Wren-complex}}
{{libheader|Wren-fmt}}
<
import "/fmt" for Fmt
Line 2,463:
}
if (i < kx-1) System.print()
}</
{{out}}
|