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:
https://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