[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]

----
[[
[Category Mathematics]
| Category Numerical Analysis]
]]