Waveguide Synthesis

A small introductory text on the mathematics and methods to deal with waves propagating through (simple) media, and how to represent this in computer form. It is only a working text (not finished, and as yet unchecked).

Working programs to generate sound files and graphs with waveguide based algorithms are already under test, and results can be viewed through (tcl/tk based) animated waveguide viewer (for simultaneously viewing all the data in the wave guide as it changes with time) and through a sample viewing processing package, which includes viewing the data of one element as a sample and through a fft (also 3d, but the package 'wavelab' seems a bit suspect here). Of course samples can also be made audible with a variety of parameters and options (thus far, some interesting percussion-like sounds have been produced already, more later). A three dimensional viewer for the propagation of data through the waveguide is in progress. The 'Mesa' package has been successfully compiled under the cygnus compiler (see other page), which means source code and library for quite powerfull and efficient (3d) rendering are available, including some test programs, which indicate that a serious 3d viewer can be made without too much trouble, and with mainstream rendering quality (goureoud shading, interaction time view changes, etc.) for fee ! Oh, and there are Java interfaces/shared linbs available for the opengl functions as well (which are compatible). Of course the VRML viewers which appear at a rapid pace lately could serve similar purpose, but don't allow one to do one's own programming).

At any rate, the idea is to make some sensible tools for this kind of programming and data processing, to lead up to the non-linear waveguides described elsewhere on my pages, and which supposedly will add new synthesis methods, and general thinking about similar systems (including the brain).



   Waveguide Synthesis

As an experiment to generate sounds from computing wave propagation 
through a medium, a simple waveguide model is derived and implemented 
in C.


   Basic Waveguide

Assuming a 1 dimensional medium, where displacements are linearly 
coupled through the second derivative (i.e. the force on a 
particle is proportional to the displacement), the medium can be 
discretized by the 'spring-weight' analogy, where a piece of the 
mass of the medium is represented by a mass, and the coupling 
with its environment by springs (applying forces linearly 
dependent on their length).

In this picture E1-3 are elements of the medium (lets assumne they 
have equal mass Me), coupled through springs s1,2 with spring 
constants Cs, stretched over distances d1,2:

      -------        -------               -------
      |     |        |     | _   _   _   _ |     |
 . . .| E1  |/\/\/\/\| E2  |/ \_/ \_/ \_/ \| E3  |. . .
      |     |  s1    |     |       s2      |     |
      -------        -------               -------
         |<------------>|<------------------->|
                d1                 d2

Of course only a part of the medium is shown, and it is a 
simplification od most physical media because it is linearly 
coupled has equal mass distribution, has no special boundary 
conditions (yet), is only one dimensional,
and because it is exactly second of order 
per element). The discretization is supposed to be accurate 
enough to introduce no significant sampling error compared to 
a continuous representation of the medium. The drawing assumes 
a longitudinal displacement only, and a rest condition with 
non-zero tension, and coupling only over neighbouring elements,
which can be extended with transversal motion and seperate per element 
zero-pulling springs, and coupling over larger distances.

Losses are not mentioned at all, but can evidently be added at 
all elements or in all springs.



   Mathematical Representation

Assume the following variables and parameters:

e[0]..e[n]      displacement vector
v[0]..v[n]      velocity vector
a[0,1],a[n-1,n] acceleration (force) vector
s               spring constant (coupling)
m               overall element mass

The relation between the elements is given by:

 2
d e   s de
--- = ----
  2   m dx
dt

This is verified by taking f/m=a, with acceleration as the second
derivative of displacement, and f = s de/dx, where the distance 
between the elements is denoted by x. A constant could be added 
to represent the amount of enery residing in the system.

Of course this is a one dimensional version of a special case 
of Schoedingers equation, where e is the wavefunction.



   Computer representation

The above translates almost directly to a set of arrays containing 
displacement, velocity and acceleration data for each element, 
and all that is required now is an evaluation function to 
effectuate the energy transfer by a properly written discrete 
version of the differential equation. Simplifications can be 
made at various points, of which the first is to linearize all 
interactions in small time frames. For each time step dt, one 
computation cycle is performed, which is then repeated over the 
(virtual) computation time interval. Apart from boundary conditions, 
the operations performed at each element e are similar.
The order of computing the various element 
state-variables depends on the the use of temporal variables, 
and how much computations are optimized for memory use, locality
of reference, or repetition of computations.
First the most straightforward evaluation scheme will be used, which is 
not too bad in terms of implementataion efficiency on a von-neumann 
machine (in this case a pc).

Starting from an inital state (by filling the displacement and velocity 
arrays with appropriate data), the forces between the elements can be 
computed by taking the difference in displacement between neighbouring 
elements and appying the spring equation:

F(Ei) = c(E   -E )  + c(E -E   )
           i-1  i        i  i+1

and the acceleration (with non-unity mass):

a(Ei) = F(Ei)/m,

the new (discretized) velocity (simple forward referencing 
integration):

v(Ei,t   ) = v(Ei,t )  +  a dt
      j+i          j



   Computer Evaluation

To 'update' the states of E0-n, that is the displacement and 
velocity of each element, the above equations can be used repeatedly,
after proper initialization:

initialize_states(E)
while simulation not finished
   for each element Ei
      compute force between neigbouring elements
      update velocity v(i)
   for each element Ei
      update displacement e(i)

The acceleration follows from:

a[i] = (c/m) (E[i-1]+e[i+1]+2E[i])

Integrating by simple forward Euler yields the new velocity 
and displacement:

v[i] += a[i] dt
e[i] += v[i] dt,

where dt is appropriately choosen fraction.
Note that the integration step can only be performed after all
accelerations a[i] have been computed, otherwise the displacement 
values need to be stored in an intermediate array.
Using integer math (usually faster), a fixed point format can 
be used to represent these values, by premultiplying them or in  
other words giving them a fixed point interpretation.

Note that the sprig constants and weights are real physical 
entities that can be given any value in the program however, 
to fit the synthesis purposes at hand, and can be combined with 
the integration by simply gathering terms.
Combining this with a simple impulse initialization, and a 
fixed point format with a exponent of 2^10 (1024), a C-program 
to compute the wavepropagation in a length n waveguide 
with circular shape (i.e. the wave is contained in a 'torus') 
looks like:

/************* wave.c: Basic Waveguide Code *****************/

#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <string.h>

#include <sys/types.h>

#define PREMUL 1000            /* Fixed exponent */
#define MAXN   4096            /* Max # nodes    */
#ifndef PI
#define PI 3.1415926535
#endif

int notcl = 0, nosound = 0;    /* to prevent (time expensive) output */

int e[MAXN], v[MAXN];          /* displacement and velocity arrays   */
int u1 = 71, u2 =  71;         /* combined values of (c/m)*PREMUL*dt */
int n;                         /* number of nodes                    */
int steps;                     /* number of (full wave) iter. steps  */
int fd;                        /* file descriptor for sound file     */
int binoutnode;                /* node to output in sound file       */

init(n)
int n;
{
   int i;
   for (i=0; i<n; i++) {
      e[i] = 0; v[i] = 0;      /* init to zero energy / strain       */
   }
   for (i=0; i<n; i++) {
      e[i] = (int) (PREMUL*4 * 
                (sin(
                   8.0 * PI * (double) i/ (double) n) 
                 ) 
                 /exp(
                   pow(
                      (4.0 * (double) (i-n/2)
                         /(double) n
                      ),2)
                 ));
      v[i] = 0; /* (int) (PREMUL/2 * cos(2.0 * PI * (double) i/ (double) n) ); */
   }
   return;

   e[0] = 4 * PREMUL;         /* impulse of 10 units                */
   e[32]= -4 * PREMUL;        /* impulse of -10 units               */
   e[64] = 4 * PREMUL;         /* impulse of 10 units                */
   e[96]= -4 * PREMUL;        /* impulse of -10 units               */
   e[128] = 4 * PREMUL;         /* impulse of 10 units                */
   e[160]= -4 * PREMUL;        /* impulse of -10 units               */
   e[192] = 4 * PREMUL;         /* impulse of 10 units                */
   e[224]= -4 * PREMUL;        /* impulse of -10 units               */
}

inittcl(n)
int n;
{
   int i;
   if (notcl) return;
   printf("if {[winfo exists .c] == 0} {canvas .c; pack .c -expand y -fill both}\n");
   printf("proc graph {n} {global a; for {set i 0} {$i<$n} {set i [expr $i+1]} {.c create line [expr 2*$i] 300 [expr 2*$i] [expr 300-$a($i)/30] -width 1 -fill red}\n.c create line 0 300 [expr 2* $n] 300 -fill green}\n");
   printf("proc graphb {n} {global a; for {set i 1} {$i<$n} {set i [expr $i+1]} {.c create line [expr 2*$i] [expr 300-$a([expr $i-1])/30] [expr 2*$i] [expr 300-$a($i)/30] -width 1 -fill green}\n.c create line 0 300 [expr 2* $n] 300 -fill blue}\n");
   printf("set nnodes %d\n",n);
   printf(".c delete all\n");
   for (i=0; i<n; i++)
      printf("set a(%d) %d\n",i,e[i]);
   printf("graphb $nnodes\n update \n after 2000\nset vt 2\n");
}

initbin(n,s)                     /* could add a .wav header here   */
int n;
char *s;
{
   if (nosound) return;
   fd = open(s,O_WRONLY|O_CREAT|O_TRUNC|O_BINARY,S_IWUSR);
}

outputtcl(n)
int n;
{
   int i;
   if (notcl) return;
   for (i=0; i<n; i++)
      printf("set a(%d) %d\n",i,e[i]);
   printf(".c delete all\n");
   printf("graph $nnodes\n update\n after $vt\n");
}

outputbin(i)
int i;
{
   char *cp;
   if (fd > 0) {
      cp = (char *) &(e[i]);
      write(fd,cp,1);
      write(fd,cp+1,1);
   }
}

propagate(n)                   /* Compute new velocity and           */
int n;                         /* displacements once for all nodes.  */
{                              /* (more optimization possible here)  */
   register int i;
   for (i=0; i<n; i++)
      v[i] += (u1 * (e[(i-1+n)%n]+e[(i+1)%n]-2*e[i]) ) / PREMUL;
   for (i=0; i<n; i++)
      e[i] += (u2 * v[i] ) / PREMUL;
}

simulate(n,tn)
int n, tn;
{
   int t;
   for (t=0; t<tn; t++) {
      propagate(n);
      output(t,n);
      outputbin(binoutnode);
   }
}

output(t,n)
int t,n;
{
   int i;
   outputtcl(n);
   return;


   printf("t=%3d|",t);
   for (i=0; i<n; i++)
      printf(" %3d",e[i]);
   printf("\n");
}

closeup()
{
   if (fd > 0) close(fd);
}


process_args(argc,argv)
int argc;
char *argv[];
{
   int i;
   if ((argc < 3) || ( ((argc-3)%2) != 0 )) {
      printf("Usage %s <#nodes> <#steps> [<option> <value>]...\n",
         argv[0]);
      printf("option = binfile  <file.raw>  | tcl <1|0>        | \n");
      printf("         initstate <file.asc> | outnode <node #>   \n");
      exit(-1);
   }
   sscanf(argv[1],"%d",&n);
   sscanf(argv[2],"%d",&steps);
   for (i=3; i<argc; i++) {
      if (strcmp(argv[i],"binfile") == 0) {
           initbin(n,argv[i+1]); i++;
      }
      if (strcmp(argv[i],"tcl") == 0) {
         if (argv[i+1][0] == '1') notcl=0; else notcl=1; i++;
      }
      if (strcmp(argv[i],"initstate") == 0) {
         printf("Warning: initstate not yet implemented\n"); i++;
      }
      if (strcmp(argv[i],"outnode") == 0) {
         sscanf(argv[i+1],"%d",&binoutnode); i++;
      }
   } 
}

main(argc, argv)
int argc;
char *argv[];
{
   notcl = 1;
   process_args(argc,argv);
   init(n);  inittcl(n);
   simulate(n,steps);
   closeup();
   return 0;
}

/************* end of wave.c *******************************/


Clearly, more sophisticated IO mechanisms and display tools are 
indispensible here, and they are already available and will be made 
public when I find some more webspace.

Latest Update
In this lates update I've added a simple 16 bit raw audio format 
as possible output, and when 'tcl 1' is added on the command line, 
the program will output a tcl/tk file that automatically displays 
an animated version of the whole waveguide during its evaluation.
It suffices to say:

gcc -o sim.exe sim.c -lm
sim 32 2000 tcl 1 > waveani.tcl

under the cygnus compiler/unix environment for windows95 (see elsewhere)
and to start a tcl/tk 'wish' interpreter which then displays 
the waveguide as a bar graph (preceded by the initial values) by:

source waveani.tcl

Wavelab (a demo is on the web) can be used to display an hear 
the waves that appear at a node of the waveguide as time progresses
by using:

  /Theover/Sources $ time sim 128 30000 binfile sound.raw outnode 0
  0.00user 0.00system 0:08.33elapsed 0%CPU (0avgtext+0avgdata 0maxresident)k
  0inputs+0outputs (0major+0minor)pagefaults 0swaps

and load sound.raw as a 16bit, 48kHz, mono sample in wavelab.
On a 100MHz Pentium, computation times for this unoptimized 
version are reasonable, small waveguides could even be used in 
realtime.

Notice that only simple integer arithmetic is used (need not even be 
long words, only multiplies and adds, and even the multiplies 
can be completely alleviated), and that the code is simple 
enough to probably run on DSP's. I did no elaborate testing 
on this version, it was merely intended to illustrate and 
test the principles ( I have other versions), and are provided 
as-is, I'll look more accurately at it as a have time.
In fact this whole text including a compiled version of the 
program was made in about two days.


   Appendix

My homepage contains some references to this and related work at:
http://theover.tripod.com/Dds/index.html

Some html refs to related pages:
http://www.nd.edu/Courses/kantor/cheg258/cheg258-1995/notes/examples/example4_26_95.html
http://www.vision.ee.ethz.ch/~tkoller/java/schroedinger.html