Version 1 of Fiddling with Fourier series

Updated 2003-04-16 07:23:39

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 rightaway, 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
 }

 #
 # Some simple tests
 #
 if {[file tail $::argv0] == [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 vv [toTimeseries $f 21]
    puts $vv
 }

---

[ Category Mathematics ]