Arjen Markus (22 march 2004) Fourier series and transforms are a powerful mathematical technique. The script below implements the so-called discrete Fourier transform. It might grow into a complete package for dealing with Fourier transforms one day.
# dft.tcl --
# Discrete Fourier transformation
#
#
namespace eval ::Fourier {
variable PI [expr {4.0*atan(1.0)}]
}
# dft --
# Perform a discrete Fourier transformation and return the
# results as a list
# Arguments:
# data The data to transform (a simple list)
# Result:
# A list of cosine and sine amplitudes
# Note:
# The base period is implicitly taken as 1, so the period
# of component j is 1/j (j>=1).
# To preserve information about the length of the original
# list of data, the last element (which has no meaning any
# other way) is set to {} for an even number of data.
#
proc ::Fourier::dft { data } {
variable PI
set result {}
set length [llength $data]
set nocomps [expr {$length/2}]
set even [expr {$length%2==0}]
for { set j 0 } { $j <= $nocomps } { incr j } {
set omega [expr {2.0*$PI*$j/double($length)}]
set sumcos 0.0
set sumsin 0.0
set phi 0.0
foreach d $data {
set sumcos [expr {$sumcos+$d*cos($omega*$phi)}]
set sumsin [expr {$sumsin+$d*sin($omega*$phi)}]
set phi [expr {$phi+1.0}]
}
if { $j < $nocomps || ! $even } {
lappend result [expr {2.0*$sumcos/$length}] [expr {2.0*$sumsin/$length}]
} else {
lappend result [expr {2.0*$sumcos/$length}] {}
}
}
return $result
}
# estimate --
# Estimate the value at a given coordinate
# Arguments:
# comps The amplitude of the components (a simple list)
# xval The value of the independent coordinate
# Result:
# Estimated value
# Note:
# The method is numerically unstable, but it should be okay with
# small series
#
proc ::Fourier::estimate { comps xval } {
variable PI
set value 0.0
set length [llength $comps]
set nodata [expr {$length-1}]
set even [expr {[lindex $comps end] == {}}]
if { $even } {incr nodata -1}
set omega [expr {2.0*$PI*$xval}]
set sumcos 0.0
set sumsin 0.0
set phi 1.0
set sumcos [expr {[lindex $comps 0]/2.0}]
set last end
if { $even } { set last end-2 }
foreach {c s} [lrange $comps 2 $last] {
set sumcos [expr {$sumcos+$c*cos($omega*$phi)}]
set sumsin [expr {$sumsin+$s*sin($omega*$phi)}]
set phi [expr {$phi+1.0}]
}
if { $even } {
foreach {c s} [lrange $comps end-1 end] {
set sumcos [expr {$sumcos+$c*cos($omega*$phi)/2.0}]
}
}
set value [expr {$sumcos+$sumsin}]
return $value
}
# energy --
# Determine the "energy" density or the cumulative energy of a
# Fourier series
# Arguments:
# option Type of output: -density (default) or -cumulative
# comps Fourier series
# Result:
# Energy density or the cumulative energy
#
proc ::Fourier::energy { args } {
variable PI
if { [llength $args] == 0 || [llength $args] > 2} {
return -code error "Should be: energy ?option? comps"
}
if { [llength $args] == 1 } {
set option "-density"
set comps [lindex $args 0]
} else {
set option [lindex $args 0]
set comps [lindex $args 1]
}
set length [llength $comps]
set nodata [expr {$length-1}]
set even [expr {[lindex $comps end] == {}}]
if { $even } {incr nodata -1}
set result [list [expr {$nodata*pow([lindex $comps 0]/2.0,2)}]]
set sum 0.0
if { $option == "-cumulative" } { set sum [lindex $result 0] }
foreach {c s} [lrange $comps 2 end] {
if { $s == {} } {
set s 0.0
set c [expr {$c/sqrt(2.0)}]
}
set sum [expr {$sum+0.5*$nodata*($c*$c+$s*$s)}]
lappend result $sum
if { $option == "-density" } {
set sum 0.0
}
}
return $result
}
# truncate --
# Truncate the Fourier series (all components after a given one are
# set to zero)
# Arguments:
# last Index of the last component
# comps Fourier series
# Result:
# List with trailing zeroes
#
proc ::Fourier::truncate { last comps } {
set length [llength $comps]
set nodata [expr {$length-1}]
set even [expr {[lindex $comps end] == {}}]
if { $even } {incr nodata -1}
set result {}
set idx 0
foreach {c s} $comps {
if { $idx > $last } {
set c 0.0
if { $s != {} } {
set s 0.0
}
}
lappend result $c $s
incr idx
}
return $result
}
# main --
# Small test of the procedures
#
puts [::Fourier::dft {1.0 0.0 0.0}]
puts [::Fourier::dft {1.0 0.0 0.0 0.0}]
puts [::Fourier::dft {1.0 1.0 1.0 1.0}]
puts [::Fourier::dft {1.0 1.0 -1.0 -1.0}]
puts [::Fourier::invdft [::Fourier::dft {1.0 0.0 0.0}]]
puts "Energy:"
puts [set series1 [::Fourier::dft {1.0 0.0 0.0 0.0 0.0 0.0 0.0}]]
puts [set series2 [::Fourier::dft {1.0 0.0 1.0 0.0 1.0 0.0}]]
puts "Density: [::Fourier::energy $series1]"
puts "Cumulative: [::Fourier::energy -cumulative $series1]"
puts "Density: [::Fourier::energy $series2]"
puts "Cumulative: [::Fourier::energy -cumulative $series2]"
puts "Truncate:"
puts [::Fourier::truncate 2 $series1]
puts [set series3 [::Fourier::truncate 2 $series2]]
puts "Energy of truncated series: [::Fourier::energy -cumulative $series3]"
puts "Estimates:"
puts [::Fourier::estimate $series2 0.0]
puts [::Fourier::estimate $series2 0.1]
puts [::Fourier::estimate $series2 0.16667]
puts [::Fourier::estimate $series2 0.2]
puts [::Fourier::estimate $series2 0.3]
puts [::Fourier::estimate $series2 0.33333]
puts [::Fourier::estimate $series2 0.4]]