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 UpdateIn 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