# fourier.tcl --
# Package to manipulate Fourier series
#
# Author: Arjen Markus (arjen.markus@wldelft.nl)
#
# Notes on the implementation:
# 1. The Fourier series is represented as a tagged list, where
# the tag is either FOURIER or FOURIER-HALF. Aftre the tag
# follow two values, the amplitudes for the cosine and sine
# components.
# 2. The interval of interest is always [-1,1]
# 3. The function is supposed to be periodic on an interval [-1,1]
# or on the interval [-2,2], indicated as FOURIER-HALF
#
# math::fourier --
# Namespace for the commands
#
namespace eval ::math::fourier {
namespace export sineFourier toTimeseries
# Possible extension: plenty!
}
# sineFourier --
# Construct a Fourier series based on sines only from an expression
# for the amplitudes (so an antisymmetric function results)
#
# Arguments:
# ncomps Number of components
# var Name of the variable in the expression (the index of
# the component)
# expression Expression to determine the amplitudes
# option Whether the interval [-1,1] covers the whole
# period or only half of it (-half). Optional
# Return value:
# Tagged list
#
proc ::math::fourier::sineFourier { ncomps var expression {option {}} } {
if { $option == "-half" } {
set series "FOURIER-HALF"
} else {
set series "FOURIER"
}
#
# TODO: properly evaluate other variables in caller's scope
#
lappend series 0.0 0.0
for { set $var 0 } { [set $var] < $ncomps } { incr $var } {
set ampl [expr $expression]
lappend series 0.0 $ampl
}
return $series
}
# toTimeseries --
# Return the values for equally spaced points
#
# Arguments:
# series Fourier series to use
# npoints Number of points
# Return value:
# List of values
#
proc ::math::fourier::toTimeseries { series npoints } {
set tag [lindex $series 0]
if { $tag == "FOURIER-HALF" } {
set period -0.5
} elseif { $tag == "FOURIER" } {
set period 0.0
} else {
error "Series is not a Fourier series"
}
if { $npoints < 2 } {
error "Number of points must be at least 2"
}
set ncomps [expr {([llength $series]-1)/2}]
set delx [expr {2.0/($npoints-1)}]
set a0 [lindex $series 1]
for { set j 0 } { $j < $npoints } { incr j } {
lappend result $a0
}
set elem 3
for { set i 1 } { $i < $ncomps } { incr i } {
set acos [lindex $series $elem]
incr elem
set asin [lindex $series $elem]
incr elem
set period [expr {$period+1.0}]
set newresult {}
for { set j 0 } { $j < $npoints } { incr j } {
set x [expr {-1.0+$j*$delx}]
set resx [lindex $result $j]
set acosx [expr {$acos * cos(3.1415926*$period*$x)}]
set asinx [expr {$asin * sin(3.1415926*$period*$x)}]
set resx [expr {$resx + $acosx + $asinx}]
lappend newresult $resx
}
set result $newresult
}
return $result
}
# plotTimeseries --
# Plot a timeseries
#
# Arguments:
# series Fourier series to use
# npoints Number of points
# colour Which colour to use
# Return value:
# List of values
#
proc ::math::fourier::plotTimeseries { series npoints colour } {
set yvalues [toTimeseries $series $npoints]
set x0 -1.0
set dx [expr {2.0/($npoints-1)}]
set xvalues {}
set x1 $x0
set y1 [lindex $yvalues 0]
for { set i 1 } { $i < $npoints } { incr i } {
set x2 [expr {$x0+$dx*$i}]
set y2 [lindex $yvalues $i]
foreach {sx1 sy1} [scaleToCanvas $x1 $y1] {break}
foreach {sx2 sy2} [scaleToCanvas $x2 $y2] {break}
.c create line $sx1 $sy1 $sx2 $sy2 -fill $colour
set x1 $x2
set y1 $y2
}
}
# scaleToCanvas --
# Scale points for drawing on the canvas
#
# Arguments:
# x X coordinate
# y Y coordinate
# Return value:
# List of scaled x and y
#
proc scaleToCanvas { x y } {
global scale
set sx [expr {$scale(pxmin)+($x-$scale(xmin))*$scale(xfactor)}]
set sy [expr {$scale(pymax)-($scale(ymax)-$y)*$scale(yfactor)}]
return [list $sx $sy]
}
# setCanvasScales --
# Set the drawing area and the extreme coordinates for drawing
#
# Arguments:
# area List of pixel coordinates: xmin, ymin, xmax, ymax
# coords List of world coordinates: xmin, ymin, xmax, ymax
# Return value:
# None
# Side effect:
# Calculation of coordinate transformation
#
proc setCanvasScale { area coords } {
global scale
set scale(pxmin) [lindex $area 0]
set scale(pymax) [lindex $area 1] ;# Because of inversion of y-coordinate
set scale(xmin) [lindex $coords 0]
set scale(ymax) [lindex $coords 3]
set xmin [lindex $coords 0]
set xmax [lindex $coords 2]
set pxmax [lindex $area 2]
set xfactor [expr {($pxmax-$scale(pxmin))/($xmax-$xmin)}]
set ymin [lindex $coords 1]
set ymax [lindex $coords 3]
set pymin [lindex $area 3]
set yfactor [expr {($scale(pymax)-$pymin)/($ymax-$ymin)}]
set scale(xfactor) $xfactor
set scale(yfactor) $yfactor
return
}
#
# Some simple tests
#
catch { console show }
if {[file tail $::argv0] == [file tail [info script]]} {
namespace import ::math::fourier::*
# This Fourier series converges to the function:
# f(x) = -1/2 if x < 0
# f(x) = 1/2 if x >= 0
set f [sineFourier 20 k {1.0/($k+0.5)/3.1415926} -half]
puts $f
set values [toTimeseries $f 41]
puts $values
if { [info exists tk_version] } {
canvas .c -width 400 -height 300
pack .c -fill both
setCanvasScale {50 50 380 280} {-1.0 -1.0 1.0 1.0}
foreach {xaxis1 yaxis1} [scaleToCanvas -1.0 -1.0] {break}
foreach {xaxis2 yaxis2} [scaleToCanvas 1.0 1.0] {break}
foreach {xaxis3 yaxis3} [scaleToCanvas 0.0 0.0] {break}
.c create line $xaxis1 $yaxis3 $xaxis2 $yaxis3 -fill black
.c create line $xaxis3 $yaxis1 $xaxis3 $yaxis2 -fill black
foreach {ncomps colour} {
1 blue 2 green 3 yellow 4 red 5 magenta 20 black} {
set series [sineFourier $ncomps k {1.0/($k+0.5)/3.1415926} -half]
plotTimeseries $series 100 $colour
}
}
}TV Time for some graphics?AM That would surely improve the script :) But then I was thinking also of procs to manipulate the Fourier series (simulate the effect of low- and high-pass filters for instance, do basic arithmetic and so on). Well, it is but a first exercise ...Okay, I enhanced the script with a colourful graph.
AM I have continued the experiment with a Discrete Fourier Transform
