225 lines
6.8 KiB
Plaintext
225 lines
6.8 KiB
Plaintext
-- demo\rosetta\Deconvolution.exw
|
|
|
|
function m_size(sequence m)
|
|
--
|
|
-- returns the size of a matrix as a list of lengths
|
|
--
|
|
sequence res = {}
|
|
object me = m
|
|
while sequence(me) do
|
|
res &= length(me)
|
|
me = me[1]
|
|
end while
|
|
return res
|
|
end function
|
|
|
|
function product(sequence s)
|
|
--
|
|
-- multiply all elements of s together
|
|
--
|
|
integer res = s[1]
|
|
for i=2 to length(s) do
|
|
res *= s[i]
|
|
end for
|
|
return res
|
|
end function
|
|
|
|
function make_coordset(sequence size)
|
|
--
|
|
-- returns all points in the matrix, in zero-based indexes,
|
|
-- eg {{0,0,0}..{3,3,5}} for a 4x4x6 matrix [96 in total]
|
|
--
|
|
sequence res = {}
|
|
integer count = product(size)
|
|
for i=0 to count-1 do
|
|
sequence coords = {}
|
|
integer j = i
|
|
for s=length(size) to 1 by -1 do
|
|
integer dimension = size[s]
|
|
coords &= mod(j,dimension)
|
|
j = floor(j/dimension)
|
|
end for
|
|
coords = reverse(coords)
|
|
res = append(res,coords)
|
|
end for
|
|
return res
|
|
end function
|
|
|
|
function row(sequence g, f, gs, gc, fs, hs)
|
|
--
|
|
--# Assembles a row, which is one of the simultaneous equations that needs
|
|
--# to be solved by reducing the whole set to reduced row echelon form. Note
|
|
--# that each row describes the equation for a single cell of the 'g' function.
|
|
--#
|
|
--# Arguments:
|
|
--# g The "result" matrix of the convolution being undone.
|
|
--# h The known "input" matrix of the convolution being undone.
|
|
--# gs The size descriptor of 'g', passed as argument for efficiency.
|
|
--# gc The coordinate in 'g' that we are generating the equation for.
|
|
--# fs The size descriptor of 'f', passed as argument for efficiency.
|
|
--# hs The size descriptor of 'h' (the unknown "input" matrix), passed
|
|
--# as argument for efficiency.
|
|
--
|
|
sequence row = {},
|
|
coords = make_coordset(hs)
|
|
for i=1 to length(coords) do
|
|
sequence hc = coords[i]
|
|
object fn = f
|
|
for k=1 to length(gc) do
|
|
integer d = gc[k]-hc[k]
|
|
if d<0 or d>=fs[k] then
|
|
fn = 0
|
|
exit
|
|
end if
|
|
fn = fn[d+1]
|
|
end for
|
|
row = append(row,fn)
|
|
end for
|
|
object gn = g
|
|
for i=1 to length(gc) do
|
|
gn = gn[gc[i]+1]
|
|
end for
|
|
row = append(row,gn)
|
|
return row
|
|
end function
|
|
|
|
function toRREF(sequence m)
|
|
--
|
|
-- [renamed] copy of Reduced_row_echelon_form.htm#Phix
|
|
-- plus one small tweak, as noted below, exit->return,
|
|
-- not that said seems to make any actual difference.
|
|
--
|
|
integer lead = 1,
|
|
rows = length(m),
|
|
cols = length(m[1])
|
|
for r=1 to rows do
|
|
if lead>=cols then exit end if
|
|
integer i = r
|
|
while m[i][lead]=0 do
|
|
i += 1
|
|
if i=rows then
|
|
i = r
|
|
lead += 1
|
|
-- if lead=cols then exit end if
|
|
if lead=cols then return m end if
|
|
end if
|
|
end while
|
|
-- nb m[i] is assigned before m[r], which matters when i=r:
|
|
{m[r],m[i]} = {sq_div(m[i],m[i][lead]),m[r]}
|
|
for j=1 to rows do
|
|
if j!=r then
|
|
m[j] = sq_sub(m[j],sq_mul(m[j][lead],m[r]))
|
|
end if
|
|
end for
|
|
lead += 1
|
|
end for
|
|
return m
|
|
end function
|
|
|
|
function lset(sequence h, sequence idx, object v)
|
|
-- helper routine: store v somewhere deep inside h
|
|
integer i1 = idx[1]+1
|
|
if length(idx)=1 then
|
|
h[i1] = v
|
|
else
|
|
h[i1] = lset(h[i1],idx[2..$],v)
|
|
end if
|
|
return h
|
|
end function
|
|
|
|
function deconvolve(sequence g, f)
|
|
--
|
|
--# Deconvolve a pair of matrixes. Solves for 'h' such that 'g = f convolve h'.
|
|
--#
|
|
--# Arguments:
|
|
--# g The matrix of data to be deconvolved.
|
|
--# f The matrix describing the convolution to be removed.
|
|
--
|
|
-- Compute the sizes of the various matrixes involved.
|
|
sequence gsize = m_size(g),
|
|
fsize = m_size(f),
|
|
hsize = sq_add(sq_sub(gsize,fsize),1)
|
|
|
|
-- Prepare the set of simultaneous equations to solve
|
|
sequence toSolve = {},
|
|
coords = make_coordset(gsize)
|
|
for i=1 to length(coords) do
|
|
toSolve = append(toSolve,row(g,f,gsize,coords[i],fsize,hsize))
|
|
end for
|
|
|
|
-- Solve the equations
|
|
sequence solved = toRREF(toSolve)
|
|
|
|
-- Create a result matrix of the right size
|
|
object h = 0
|
|
for i=length(hsize) to 1 by -1 do
|
|
h = repeat(h,hsize[i])
|
|
end for
|
|
|
|
-- Fill the results from the equations into the result matrix
|
|
coords = make_coordset(hsize)
|
|
for i=1 to length(coords) do
|
|
h = lset(h,coords[i],solved[i][$])
|
|
end for
|
|
return h
|
|
end function
|
|
|
|
constant f1 = { 6, -9, -7, -5},
|
|
g1 = {-48, 84, -16, 95, 125, -70, 7, 29, 54, 10},
|
|
h1 = {-8, 2, -9, -2, 9, -8, -2}
|
|
|
|
if deconvolve(g1, f1)!=h1 then ?9/0 end if
|
|
if deconvolve(g1, h1)!=f1 then ?9/0 end if
|
|
|
|
constant f2 = {{-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}},
|
|
g2 = {{ 40, -21, 53, 42, 105, 1, 87, 60, 39, -28},
|
|
{-92, -64, 19,-167, -71, -47, 128,-109, 40, -21},
|
|
{ 58, 85,-93, 37, 101, -14, 5, 37, -76, -56},
|
|
{-90,-135, 60,-125, 68, 53, 223, 4, -36, -48},
|
|
{ 78, 16, 7,-199, 156,-162, 29, 28,-103, -10},
|
|
{-62, -89, 69, -61, 66, 193, -61, 71, -8, -30},
|
|
{ 48, -6, 21, -9,-150, -22, -56, 32, 85, 25}},
|
|
h2 = {{-8, 1,-7,-2,-9, 4},
|
|
{ 4, 5,-5, 2, 7,-1},
|
|
{-6,-3,-3,-6, 9, 5}}
|
|
|
|
if deconvolve(g2, f2)!=h2 then ?9/0 end if
|
|
if deconvolve(g2, h2)!=f2 then ?9/0 end if
|
|
|
|
constant f3 = {{{-9, 5, -8}, { 3, 5, 1}},
|
|
{{-1, -7, 2}, {-5, -6, 6}},
|
|
{{ 8, 5, 8}, {-2, -6, -4}}},
|
|
g3 = {{{ 54, 42, 53, -42, 85, -72},
|
|
{ 45,-170, 94, -36, 48, 73},
|
|
{-39, 65,-112, -16, -78, -72},
|
|
{ 6, -11, -6, 62, 49, 8}},
|
|
{{-57, 49, -23, 52,-135, 66},
|
|
{-23, 127, -58, -5,-118, 64},
|
|
{ 87, -16, 121, 23, -41, -12},
|
|
{-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}}},
|
|
h3 = {{{ -6, -8, -5, 9},
|
|
{ -7, 9, -6, -8},
|
|
{ 2, -7, 9, 8}},
|
|
{{ 7, 4, 4, -6},
|
|
{ 9, 9, 4, -4},
|
|
{ -3, 7, -2, -3}}}
|
|
|
|
if deconvolve(g3, f3)!=h3 then ?9/0 end if
|
|
if deconvolve(g3, h3)!=f3 then ?9/0 end if
|
|
|
|
ppOpt({pp_Nest,2,pp_IntFmt,"%3d"})
|
|
pp(deconvolve(g3, f3))
|
|
pp(deconvolve(g3, h3))
|