101 lines
2.4 KiB
Tcl
101 lines
2.4 KiB
Tcl
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
|
|
}
|
|
}
|