Fiddling with Fourier series

Arjen Markus (16 April 2003) This is just a small experiment in dealing with Fourier series via Tcl. I thought I put it on the Wiki right away, so that someone who is interested can have a look at it.

To make it more useful, a lot of things need to added :) Feel free to add comments and suggestions - you might come up with things I have not yet thought about ...

# fourier.tcl --
#    Package to manipulate Fourier series
#
# Author: Arjen Markus ([email protected])
#
# 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

 Category Mathematics