Note that Tcllib provides math::fourier.
See also Fast Fourier Transform.
Arjen Markus 2004-03-22: 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.
Note: I am not satisfied with the current representation of the Fourier transform that needs to distinguish between even and odd numbered data series via an empty last list element. IMO it is not very elegant. I want to revise it and use a tagged list instead.
# 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]TV 2004-03-22: I just tried
puts [::Fourier::dft {1.0 1.0 1.0 1.0}]as a first check which gives:
2.0 0.0 -9.18454765367e-017 1.11022302463e-016 0.0 {}shouldn't a constant signal just give a transform result with only on component
1.0 * cos(0*x)
which represents the dc component so to speak,being the first component of the result vector? I seem to remember (without looking anything up to check though) there was something with the first component.
To test fft a bit, lets generate a single harmonic, of wavelength equal to the equivalent fft interval:
set l {}
for {set i 0.0} {$i < [expr 2.0*$::Fourier::PI]} {set i [expr $i+2*$::Fourier::PI/5]} {
lappend l [expr sin($i)]
}
puts [::Fourier::dft $l]-4.08277855968e-012 0.0 6.1242344529e-012 0.999999999999 -4.08268974184e-012 -2.14650519581e-012
Changing phase to the other basis functions (sine above becomes cos):
puts [::Fourier::dft $l]
5.61934943022e-012 0.0 1.0 -1.02068131724e-011 -3.47304407455e-012 1.7763568394e-016
No dc component, 1 principle of double equiv wavelength 0, second is the tone put in at amplitude 1, rest zero.
AM: I do not take any rounding errors into account - hence the very small but non-zero amplitudes in the above results. And there is something non-intuitive about discrete Fourier transforms: for even numbers of data you have an extra cosine ...
TV: To me that doesn't explain a constant signal transform of double amplitude in the FFT domain.
Lars H: That's because the code above doesn't normalize the constant term coefficient in the DFT conversion. (The constant function in the trigonometric ON basis for functions on an interval is exceptional in that its square norm is twice that of the other functions. This often shows up in formulae as an extra factor 1/2 for that term.)
Note that the above isn't FFT (Fast Fourier Transform, which is an algorithm), but merely the straightforward implementation of DFT, having O(n^2) execution time rather than the O(n \log n) of FFT. Rather than speaking of the FFT domain, one should speak of the frequency domain.
And concerning AM's trouble with odd/even number of data points: You may want to consider using the discrete cosine transform instead. That always gives the same number of coefficients as there are data points.
TV: Well, normalizing I guess would be reasonable as normal part of the transform, or the backtransform or subsequent analysis may well be of.
Apart from whether the perfect shuffle with reused elements is being used (which I didn't check), calling the transform of any of these type of limited sum intervals the Frequency Domain is dangerous, because it is only a fixed number of components which one computes in this discrete frequency transform domain, not like with the real Fourier Transform, where theoretically one does compute actual frequency components.
That's why often averaging is used of various repeated fourier transforms, and why 'frequency analysis' of many programs based on the (real or not) FFT is not the same as actual frequency analysis, and based on the observation of the limitations of the Niquist rate even possibly ambiguous.
I'd not compare the DCT as used in Jpeg compression directly. There the aforementioned (in the case of jpeg 2 dimenstional) is applied in a special way, and when you would include the phase in the transform (either by using an orthogonal basis of sine and cosine like in regular fourier analyis, or by transforming to norm/angle complex numbers), the results would become comparable.
This might fit nicely into tcllib. (KBK: Please, not in its current form. People will expect a Fast Fourier Transform - which this isn't.)
TV: As a background note: try : [L1 ] and especially for instance the Common Waveforms --> Sawtooth Perfect last 5 submenu entries for examples of what a (fast or not) fourier transform of a limited interval can do to perfectly cyclic functions and their transform, when the interval doesn't match the wavelength (as is often the case). Don't forget to first do (from the left javascript menu) Waveform Lab --> Theory --> View Applet to see the transform applet. (Requires Java 1.1 IIRC)
AM For an analoguous technique using polynomials, see Polynomial fitting
AK: See also http://www.dspguide.com for an online book (pdf) on digital signal processing.