File: pramode.html

package info (click to toggle)
lg-issue98 1-1
links: PTS
area: main
in suites: sarge
size: 1,180 kB
ctags: 91
sloc: sh: 78; makefile: 34; perl: 23
file content (681 lines) | stat: -rw-r--r-- 15,134 bytes

<html>
<head>
<link href="../lg.css" rel="stylesheet" type="text/css" media="screen, projection"  />
<title>Mathematical Explorations with Scilab/Linux LG #98</title>

<style type="text/css" media="screen, projection">
<!--


.articlecontent {
	position:absolute;
	top:143px;
}


-->
</style>


</head>

<body>


<img src="../gx/2003/newlogo-blank-200-gold2.jpg" id="logo" alt="Linux Gazette"/>
<p id="fun">...making Linux just a little more fun!</p>


<div class="content articlecontent">

<div id="previousnexttop">
<A HREF="orr.html" >&lt;-- prev</A> | <A HREF="pranevich.html" >next --&gt;</A>
</div>



<h1>Mathematical Explorations with Scilab/Linux</h1>
<p id="by"><b>By <A HREF="../authors/pramode.html">Pramode C.E</A></b></p>

<p>
Little would
<a href="http://www-gap.dcs.st-and.ac.uk/~history/Mathematicians/Fourier.html">
Jean Baptiste Joseph Fourier</a>,
the 18th century French mathematician and revolutionary, have imagined that
the analytical techniques he had invented to study 
the behaviour of mathematical functions would someday become one
of the most powerful tools in the hands of scientists and
engineers working in disciplines as diverse as neurophysiology and
digital communication. 

<P> As I was fast sliding into the depths
of mathematical ignorance, I thought maybe I would refresh some
high school memories by trying to understand a bit of Fourier's 
math. Much of what I read flew far above my head - my only 
consolation was that I discovered Linux to be an ideal platform 
not only for Operating System hacking but also for
mathematical recreation and research.

<P> I came upon a great tool called Scilab and also a nice little tutorial on
Fourier
Math by <a href="http://www.ibiblio.org/obp/py4fun/wave/wave.html">Chris Meyers</a> which demonstrated
some interesting sine-wave combination/analysis stuff using
Python code. This article demonstrates a few simple Scilab 
tricks and reimplements Chris's code in Scilab's native
scripting language. Readers looking for mathematical wisdom are
warned not to rely too much on what I say here!


<h2>What is Scilab?</h2>
<p>
<a href="http://www.scilab.org">Scilab</a> is a powerful, free
environment for mathematical computation. It provides an extensible
framework for general matrix manipulation and `toolboxes' for doing
stuff like control system design, digital signal processing etc.
The C/Fortran source code is available for download from the project
home page - I had absolutely no difficulty in building the system -
the standard `configure; make; make install' magic worked perfectly.
<p>
Here is a screen shot of Scilab running on my Linux box:
<p>
<img src="misc/pramode/sci1.png">

<h2>Simple math</h2>
<p>
Let's get started by doing a few simple matrix manipulations. A 3-by-3
matrix is created by simply typing, at the Scilab prompt:

<pre>
--&gt;a = [1,10,20; 5,6,7; 12,11,45]
 a  =
 
!   1.     10.    20. !
!   5.     6.     7.  !
!   12.    11.    45. !
 
--&gt;
</pre>

It's easy to get the transposed matrix:
<pre>
---&gt;a'
ans  =
 
!   1.     5.    12. !
!   10.    6.    11. !
!   20.    7.    45. !
 
--&gt;

</pre>
A few other functions:

<pre>

--&gt;sum(a, 'c')
 ans  =
 
!   31. !
!   18. !
!   68. !
 
--&gt;sum(a, 'r')
 ans  =
 
!   18.    27.    72. !
 
--&gt;diag(a)
 ans  =
 
!   1.  !
!   6.  !
!   45. !
 
--&gt;

</pre>

<p>
Elements can be extracted from matrices in many different
ways - the simplest is the standard indexing procedure. 
Writing a(1,2) would yield the element at row 1 and column 2 (note
that the index starts at 1). Indexing a matrix beyond its bound
will result in an error. Writing to a non-existent index will
result in the matrix growing dynamically.

<pre>

--&gt;a(3,4) = 3
 a  =
 
!   1.     10.    20.    0. !
!   5.     6.     7.     0. !
!   12.    11.    45.    3. !
 
--&gt;

</pre>
</p>

<h2>The 'colon' operator</h2>
<p>

The 'colon' is a cute little operator. We can
create a vector of numbers 1,2,3 ... 10 by just
writing:

<pre>

--&gt;a = 1:10
 a  =
 
!   1.    2.    3.    4.    5.    6.    7.    8.    9.    10. !
 
--&gt;

</pre>

Many other tricks are possible:

<pre>

--&gt;b
 b  =
 
!   1.    2.    3. !
!   4.    5.    6. !
!   7.    8.    9. !
 

--&gt;b(1:3,2:3)
ans  =
 
!   2.    3. !
!   5.    6. !
!   8.    9. !
 
--&gt;1:2:10
 ans  =
 
!   1.    3.    5.    7.    9. !
 
--&gt;

</pre>
<p>
Note that 1:2:10 means create a vector starting from
1, each successive element being computed by adding 
2, until the value becomes greater than 10.

<h2>Simple plotting</h2>
<p>
Let's look at an example of a simple sine wave
plot. We want one full cycle of the sine curve
(from 0 to 2*PI) - let's take 240 points in
between, so each division would be 2*PI/240. First
we create a vector containing all the angle values
in this range and then we plot it (%pi is a constant 
standing for the value of PI):

<pre>
--&gt; = 0:(2*%pi)/240:2*%pi
 x  =
         column 1 to 5
 
!   0.    0.0261799    0.0523599    0.0785398    0.1047198 !
 
         column 6 to 9
 
!   0.1308997    0.1570796    0.1832596    0.2094395 !
 
         column 10 to 13
 
!   0.2356194    0.2617994    0.2879793    0.3141593 !
 
         column 14 to 17
 
!   0.3403392    0.3665191    0.3926991    0.4188790 !
[More (y or n ) ?] 

</pre>

<p>
Now, we use a simple plot function:
<pre>
--&gt;plot(x, sin(x))
</pre>

<p>
<img src="misc/pramode/sci2.png">

<h2>Writing Scilab scripts</h2>
<p>
Writing Scilab scripts is simple. Here is an example of
a 'for' loop which can be entered at the Scilab prompt itself:

<pre>
--&gt;s = 0         
 s  =
 
    0.  
 
--&gt;for i=1:3:10
--&gt;  s = s + i
--&gt;end
 s  =
 
    1.  
 s  =
 
    5.  
 s  =
 
    12.  
 s  =
 
    22.  
[More (y or n ) ?] 
 
</pre>

<h2>Defining functions</h2>
<p>
The function definition syntax is a wee bit tricky. Here is
a simple example:

<pre>

--&gt;function [r] = my_sqr(x)
--&gt;  r = x * x
--&gt;endfunction
 
--&gt;my_sqr(3)
 ans  =
 
    9.  
--&gt;

</pre>

<p>

After the keyword 'function', we supply a list of `output values'.
Any value written to an `output' value will be `returned' by the
function. The argument 'x' is of course the input argument to the
function. The function returns the value 'r' which is the square of
'x'.

<p>
The question obviously is what if we want to return two values. We
try the following at the Scilab prompt:

<pre>

--&gt;function [r1, r2] = foo (x, y)
--&gt;  r1 = x + y
--&gt;  r2 = x - y
--&gt;endfunction
 
--&gt;[p, q] = foo(10, 20)
 q  =
 
  - 10.  
 p  =
 
    30.  
 
--&gt;

</pre>

Note the special way we call the function. The value of r1 will get
transferred to 'p' and value of r2 to 'q'.

<p>
The following invocations of 'foo' demonstrates the fact that
the language is dynamically typed. 

<pre>
--&gt;[p, q] = foo([1,2], 1)
 q  =
 
!   0.    1. !
 p  =
 
!   2.    3. !
 
--&gt;[p, q] = foo([1,2], [3,4,5])
 !--error     8 
inconsistent addition
at line       2 of function foo   called by :  
[p, q] = foo([1,2], [3,4,5])
 
--&gt;
</pre>

<p>
It is possible to store function definitions in a file and
load them at a later time. Suppose the above function definition
is stored in a file called 'fun.sci'. We need to simply
invoke, at the Scilab prompt:

<pre>
--&gt;exec('fun.sci')
</pre>


<h2>Enter Fourier!</h2>
<p>
<img src="misc/pramode/fourier.jpg">

<p>
We encounter 'signals' everywhere. The PC speaker generates sound by
converting electrical signals to vibrations. We see objects around
us because these objects bounce back light signals to our eyes. Our
TV and radio receive electromagnetic signals. We are immersed in a 
'sea of signals' ! Analysis of signals is therefore of central
importance in most branches of science and engineering.

<p>
The basic Unix philosophy is `Keep it Simple, stupid'. Physicists
(and most other scientists and engineers) often can't stick to this
dictum when they start analysing stuff, simply because the phenomena they are studying have
awesome complexity. But it seems that most complex things in this
world can be explained on the basis of simpler things. Joseph Fourier's
insight was that complex time varying signals can be expressed as
a combination of simple sin/cos curves of varying frequency and
amplitude. We will verify this assumption by plotting a few simple
equations with the help of Scilab.

<p>
Let's start with a simple sum of two 'sin' signals.

<pre>

--&gt;delta = (2*%pi)/240
 delta  =
 
    0.0261799  

--&gt;x = 0:delta:2*%pi

--&gt;a = sin(x) - (1/2)*sin(2*x)
--&gt;plot(x, a)
 
</pre>
Here is the plot:
<p>
<img src="misc/pramode/sci3.png">
<p>

There is very little indication here that something interesting is
going to happen. Next, we try plotting.

<pre>
b = sin(x) - (1/2)*sin(2*x) + (1/3)*sin(3*x)
</pre>

We keep on adding terms to the series, the next term would be 
-(1/4)*sin(4*x), the next one +(1/5)*sin(5*x) and so on. Here is what I got when I plotted this
series with 200 terms in it (you will have to write
a function to do this for you):

<p>
<img src="misc/pramode/sci4.png">

<p>
Seems like magic! The sin curve has vanished completely and we 
have a brand new signal! How exactly Mr.Fourier 'knew' such a
series would ultimately give us something totally different
from the sum of its parts would be more appropriately dealt
with in a mathematics class(Do I hear you yawn? Do we have a
case for a more `practical' math education with students being
given access to Linux boxes running Scilab, Python(Numeric),
and a whole lot of other free, educational tools?)


<h2>Determining the components of a signal</h2>
<p>
We have seen that adding together sines of different
frequency and amplitude gives us signals which look
totally different. Now the question is, given some
numbers which represent a complex waveform, will we
be able to say what combination of sine's (frequency
and amplitude) gave rise to that particular signal? Let's
try.

<p>
Let's first write a function which performs simple numerical
'integration' over the range 0 to 2*PI. We divide the area
under our curve into tiny strips, each of width say 2*PI/240.
The area of a strip at point 'x' (0 &lt; x &lt; 2*PI) will
be its height multiplied by the width, which will be
sin(x) * (2*PI/240). This is the idea behind the integration
function, which can be typed at the Scilab prompt. The argument
to integrate is a vector of sin values in the range 0 to 2*PI-delta
where delta is (2*PI)/240. The difference between two successive
values in the vector is 'delta'.

<pre>
--&gt;function [r] = integrate(a)
--&gt;  r = sum(a)*(2*%pi)/240
--&gt;endfunction
</pre>

<p>

Let's try integrating the simple sin function, sin(x).

<pre>
 
--&gt;x = 0:delta:(2*%pi-delta)

--&gt;integrate(sin(x))
 ans  =
 
    3.837E-16  
 
</pre>

We see that the integral is zero. The sin curve has equal area above and below 
the zero-point.

<p>
Let's try plotting sin(x).*(-sin(x)) (Note that the .* operator performs
memberwise multiplication of two vectors):
<p>
<img src="misc/pramode/sci5.png">

<p>
We see that the function has been shifted completely below the zero-point.
It should now definitely have a non-zero area. 

<pre>

--&gt;integrate(sin(x).*(-sin(x))) 
 ans  =
 
  - 3.1415927  
 
--&gt;

</pre>

Scilab tells us it is -PI. Let's now try plotting sin(2*x).*(-sin(x))

<p>
<img src="misc/pramode/sci6.png">

<p>
The graph tells us that the integral should be zero. We verify this:

<pre>
--&gt;integrate(sin(2*x).*(-sin(x)))  
 ans  =
 
    3.977E-16  
</pre>

We are now beginning to get a 'feel' of the idea we would employ to
separate out the components of our complex signal. Multiplying a sine
with negative of a sine of a different frequency gives us zero - 
only when the frequencies match do we get non zero results. Say our
complex signal is:

<pre>
sin(x) - (1/2)*sin(2*x) + (1/3)*sin(3*x) - (1/5)*sin(5*x)
</pre>

If we multiply this with -sin(x), what we get is:

<pre>
sin(x).*(-sin(x)) - (1/2)*sin(2*x).*(-sin(x)) + 
(1/3)*sin(3*x).*(-sin(x)) - (1/5)*sin(5*x).*(-sin(x))
</pre>

The first term gives us -PI, all other terms become zero. The fact
that we are getting a non zero value tells us that sin(x) is 
present in the signal. Now we multiply the signal with -sin(2*x).
If we get a non-zero result, that means that sin(2*x) is present
in the signal. We repeat this process as many times as we wish.

<p>
How do we get the amplitude of each component? Let's try out
another experiment:

<pre>
 
--&gt;b = sin(x) - (1/2)*sin(2*x) + (1/3)*sin(3*x) - (1/4)*sin(4*x)
 
--&gt;integrate(b.*(-sin(x)))                                      
 ans  =
 
  - 3.1415927  
 
--&gt;integrate(b.*(-sin(2*x)))
 ans  =
 
    1.5707963  
 
--&gt;integrate(b.*(-sin(3*x)))
 ans  =
 
  - 1.0471976  
 
--&gt;integrate(b.*(-sin(4*x)))
 ans  =
 
    0.7853982  
 
</pre>

We see that dividing each result by -PI gives us the amplitude of
each component of the signal.

<h2>Conclusion</h2>
<p>
Very high quality proprietary tools exist for doing
numeric/symbolic math - but they are sometimes priced
beyond the reach of the student or the hobbyist. I hope
this article has convinced you that Free Software alternatives
do exist. Kindly let me know about any inaccuracies you find
in this document. I can be contacted via my home page at
<a href="http://pramode.net">pramode.net</a>. 


<h2>Acknowledgements</h2>
<p>
Thanks to the Scilab team for creating such a wonderful 
tool and also documenting it thoroughly. This article
would not have been written without the help of <a href="http://www.ibiblio.org/obp/py4fun">
Chris Meyers</a> document explaining Fourier's math -
Chris has also written some other very interesting Python
programs which you are sure to enjoy. A big Thank You to
him!

</p>


<!-- *** BEGIN author bio *** -->
<P>&nbsp;
<P>
<!-- *** BEGIN bio *** -->
<P>
<img ALIGN="LEFT" ALT="[BIO]" SRC="../gx/2002/note.png">
<em>
I am an instructor working for IC Software in Kerala, India. I would have loved
becoming an organic chemist, but I do the second best thing possible, which is
play with Linux and teach programming!
</em>
<br CLEAR="all">
<!-- *** END bio *** -->

<!-- *** END author bio *** -->

<div id="articlefooter">

<p>
Copyright &copy; 2004, Pramode C.E. Copying license 
<a href="http://linuxgazette.net/copying.html">http://linuxgazette.net/copying.html</a>
</p>

<p>
Published in Issue 98 of Linux Gazette, January 2004
</p>

</div>


<div id="previousnextbottom">
<A HREF="orr.html" >&lt;-- prev</A> | <A HREF="pranevich.html" >next --&gt;</A>
</div>


</div>






<div id="navigation">

<a href="../index.html">Home</a>
<a href="../faq/index.html">FAQ</a>
<a href="../lg_index.html">Site Map</a>
<a href="../mirrors.html">Mirrors</a>
<a href="../mirrors.html">Translations</a>
<a href="../search.html">Search</a>
<a href="../archives.html">Archives</a>
<a href="../authors/index.html">Authors</a>
<a href="../contact.html">Contact Us</a>

</div>



<div id="breadcrumbs">

<a href="../index.html">Home</a> &gt; 
<a href="index.html">January 2004 (#98)</a> &gt; 
Article

</div>





<img src="../gx/2003/sit3-shine.7-2.gif" id="tux" alt="Tux"/>




</body>
</html>