Fractal Images and Music With Perl
The Perl Journal March 2003
By Randy Kobes, Andrea Letkeman, and Olesia Shewchuk
Randy researches and teaches physics at the University of Winnipeg. He can be reached at [email protected]. Andrea currently works for an oil company, and is studying towards her Royal Conservatory of Music Teacher's certification. She can be reached at andrea_letkeman@ hotmail.com. Olesia is a student of science, computers, dance, and music. She can be reached at [email protected].
You may remember taking an IQ test in which you are given some numbers: 1 4 7 10..., and are asked to predict the next number in the sequence. What you end up doing is, at least subconsciously, starting from the first number, looking for a pattern that subsequent numbers fit, and then using that pattern to predict the next (or any other) number in the series. In this article, we'll examine some sequences of numbers that exhibit a fractal nature, and from those, see how Perl can be used to generate images and sounds that represent aspects of this nature.
The Nature of Fractals
A simple example of a fractal sequence can be given by considering the binary form of the first few numbers:
printf("%b ", $_) for (0 .. 10); 0 1 10 11 100 101 110 111 1000 1001 1010
Next, form the sequence by counting the number of times "1" occurs in each number:
for (0 .. 20) { $b = sprintf("%b", $_); $n = $b =~ tr/1//; print "$n "; } 0 1 1 2 1 2 2 3 1 2 2 3 2 3 3 4 1 2 2 3 2
Now take every other number in this sequence:
for (0 .. 20) { next unless $_ % 2 == 0; $b = sprintf("%b", $_); $n = $b =~ tr/1//; print "$n "; } 0 1 1 2 1 2 2 3 1 2 2
We see that, up to the same number of terms considered, the two sequences are the same. This is also true if we had taken every 4th number of the original sequence, and every 8th number, and so on. This self-similarity is a fundamental property of what is known as a "fractal"something that looks similar to itself on different scales. This is why, for example, when we look at a picture filled with clouds or one of a moonscape with many craters, it is hard to judge how far away the picture was takenthe fractal nature of these objects makes it hard for us to associate one particular length scale with them, as they look similar at many different scales. Fractals are also often used in generating patterns for screensavers, some of which are a true source of fascination.
A web search will quickly find the two most popular aspects of fractalsfractal images, which are especially associated with what are called the "Mandelbrot" and "Julia" sets, and fractal music. In this article, we'll describe how easily both of these can be made with Perl. What amazed early workers in the field, and is still a source of wonder today, is just how easily these can be generated, and despite the simplicity of the algorithms used, just how intricate and surprising the results can be.
In the following, we consider particular sequences of numbers to illustrate the results. We imagine the sequences to be stored in some array @x, and denote a particular member by x[n]. The sequence is generated, starting from an initial x[0], through a relation x[n+1]=F(x[n]) for some given function F(x), where n=0, 1, 2,... What this notation means is that we start with some given x[0], and then, with n=0, find x[1]=F(x[0]). Having found x[1], we then set n=1 to find x[2]=F(x[1]), and so on. For example, the sequence 1 4 7 10... can be represented as x[n+1]=x[n]+3, with x[0]=1.
The Fractal Module
The actual work of making the images and sound files in what follows is done through some functions made available in the Fractal module (fractal.pm, available online at http://www.tpj.com/source/). This module contains three functions:
draw()Takes a reference to an array of arrays (for example, [ [x1,y1,i1], [x2,y2,i2] ]), and produces an image using Lincoln Stein's GD module. The (x,y) values represent the horizontal and vertical coordinates of the points to be plotted, while the associated i value is used to help determine the color used to plot this point. This is done by first calculating an index based on the values of (x,y,i), and then using this index to calculate an rgb triplet (r,g,b) of numbers, each between 0 and 255, which specify the amounts of red, green, and blue to use. Default routines to do these are given in the module, or the user can specify his or her own routines by providing code references as draw(rgbindex => \&my_rgbindex, rgb => \&my_rgb).
compose()Takes a reference to an array of points and produces a midi sound file using Sean Burke's MIDI module. Which particular note is associated with each point is determined, as with draw(), by first calculating an index from the point under consideration and then from this index, associating a note. Again, default routines to do this are contained in the module, or the user can specify her or his own as compose(noteindex => \&my_noteindex, note => \&my_note).
min_max()Takes either an array reference or a reference to an array of arrays and calculates the minimum and maximum values. For an array reference, two scalars representing the minimum and maximum values are returned, while for a reference to an array of arrays, two array references are returned containing the corresponding information. This routine is useful in determining the index used in both draw() and compose(), and is also used in draw() in converting the (x,y) points passed to the routine into pixel coordinates based on the size of the image specified.
Fractal Images
With this, we begin by considering fractal images, for which there are two main algorithms. (One can mathematically relate the two algorithms, but it is calculationally convenient to consider them separately.) The first algorithm uses what is called the "escape-time" method. Consider the sequence of numbers (x[n],y[n]):
x[n+1] = x[n]2 - y[n]2 + cx y[n+1] = 2 x[n] y[n] + cy
where cx and cy are constants. (This unlikely looking sequence is more natural in the context of complex numbers, where it can be written as z[n+1]=z[n]2+c, where z[n] and c are complex numbers having a real and an imaginary part (see the documentation on Math:: Complex for details). Starting from some initial points (x[0],y[0]), two distinct possibilities can arisethe numbers in the sequence all remain finite or, at some stage, they start to diverge towards infinity. To avoid dealing with very large numbers, a breakout condition on the size of the numbers is imposed that when met, ends the generation of further numbers in the sequence. This defines the following escape-time algorithm: For each point (i.e., pixel) on a graph of interest, set up the initial conditions and then (up to some maximum number of iterations):
1. Find the next point in the sequence.
2. If the breakout condition on the size of the numbers is met, end the iterations.
3. Plot the point according to some criteria based on when the breakout condition was met.
A typical breakout condition used is to calculate x[n]2+y[n]2 for each point, and then to end the sequence if this is larger than some number.
There are two major types of images typically used in illustrating this algorithm, the difference being in their choices of initial points (x[0],y[0]) and the constants (cx,cy). For the Mandelbrot set, one sets x[0]=0=y[0] and chooses for (cx,cy) the coordinates of the pixel under consideration. For the Julia set, one sets (x[0],y[0]) to the coordinates of the pixel under consideration, and chooses (cx,cy) equal to fixed constants. Listing 1 is a script for generating the Mandelbrot set, while Listing 2 is for generating the Julia set. Running these scripts results in the images in Figures 1(a) and 1(b), respectively.
Apart from producing interesting images, the use of colors gives information about the behavior of these sequences. For the Mandelbrot set, points in the interior (pink, in this case) never satisfy the breakout condition, while points on the edges (shades of green, in this case) quickly do. The interesting region is around the jagged edges, where one can see evidence of the fractal nature through self-similarity. Indeed, if one blows up a region around the edges by changing (xmin,xmax,ymin,ymax) in the mandelbrot.pl script, similar-looking images result. The Julia set shows similar behavior and, in fact, is closely related to the Mandelbrot setinteresting images for the Julia set are found by choosing the constants (cx,cy) near the jagged edges of the Mandelbrot set.
Another type of fractal image is formed from what is called an "iterated function system." For this one, consider sequences of the form:
x[n+1] = a[i] x[n] + b[i] y[n] + e[i] y[n+1] = c[i] x[n] + d[i] y[n] + f[i]
where (a[i],b[i],c[i],d[i],e[i],f[i]), with i=1..N, are sets of constant numbers. (Although theoretically the constants can be chosen arbitrarily, constraints on their magnitude exist if one wishes to generate images of a finite size.) To generate an image from these sequences, one begins with some initial point (x[0],y[0]) and then implements the following algorithm for some number of iterations:
1. Choose, by some criteria, one of the sets of numbers labeled by the index i to generate the next point in the sequence.
2. Find the next point, and plot it.
Listings 3 and 4 produce the images that appear in Figure 2.
As well as the coefficients used in generating the sequences, these differ according to how the index i is chosen. For the snowflake, an i is picked randomly from those available, while for the fern, a weighted set of probabilities is used. For many cases, this latter way of choosing the index produces more complete images in fewer iterations.
Fractal Music
We now come to a second way of viewing fractalsthrough musical scores. Like art, music is a very subjective and personal thingwhat one person passionately likes, another may equally detest. However, there are some common aspects about music that span across genres and time, and one of those is a fractal nature. (The degree of "fractalness" is related to the information content, and entropy, of a system. Lovers of classical music will take pride in knowing that classical music is more "fractal" than most forms of today's popular music; we'll let the reader draw his or her own conclusions.) This again is related to the property of self-similaritypatterns of notes occurring within a few bars are often repeated throughout a score, and also occur (qualitatively) on a larger scale.
Perhaps even more significantly, what makes music interesting and enjoyable is an aspect that was also evident in the images just studiedthat of their relative complexity. Humans like to find patternsa random collection of colored pixels or notes on a scale is unpleasantbut we also like to be surprised. In music, phrases in harmony lead to what is termed "consonance," while the introduction of suspense, or tension, leads to "dissonance." Music needs elements of both to be enjoyable. Fractals, in a sense, live on the border between random noise and dull predictability.
For a simple example of fractal music, consider what is known as the "logistic map":
x[n+1] = A x[n] (1 - x[n])
where A is a constant and each x[n] lies between 0 and 1. (Historically, this sequence was proposed as a simple model of population growth, where x[n] represents the fraction of the maximum that a population attains in year n. The term Ax[n] increases the population, through reproduction, and the -Ax[n]2 term decreases the population, through predators and disease.) The behavior of this sequence depends on the parameter Afor A less than about 3, the sequences of numbers converges to some fixed value; while for A greater than about 3.5 up until 4.0, the numbers appear to be random. (Actually, the behavior is chaotic.) The key word here is "appear"for this range of A, there is actually a fractal nature present.
Sequences such as these are the basic ingredient in making fractal musicthe numbers of the sequence are simply mapped to notes of a musical scale. In Listing 5, we give an example of how this is done for the logistic map; while in Listing 6, we do so for another example called the "Hénon map":
x[n+1] = a - x[n]2 + b y[n] y[n+1] = x[n]
for which we use only the x[n] in generating the notes. For the logistic score, we employ the default methods of the Fractal module to map numbers in the sequence to notes, while for the Hénon score, we specify these within the script.
When played, one can sense these scores as being musical recognizable patterns emerge, and yet there are elements of unpredictability. (If one had to classify fractal scores, a possible candidate might be "new" or "modern" classical music, such as that of Schönberg, in which greater dissonance and a movement away from a fixed key signature plays a prominent role.) One might think it easy to create such scores from an almost arbitrary sequence, but it is really those with a fractal nature that are most pleasing. Try creating a score from just a random sequence of numbers (in the logistic map, this corresponds to setting the parameter A to 4)the difference is quite noticeable. Of course, at this level, fractal music is missing many components of real music such as dynamics, changing tempo, and varied note durations and rests. This, in principle, can be added to the scripts; one learns from these studies just how multifaceted and complicated music really is.
We can compose potentially more interesting scores by obtaining, as for the Mandelbrot set, two related sequences of numbers (x[n],y[n]), which can then be used to make two simultaneous tracks in the score. This is done using the sync() routine of the MIDI module; an example, using a point on the ragged boundary of the Mandelbrot set, appears in Listing 7. Apart from producing richer scores, an interesting philosophical aspect of doing this is the potential for being able to "hear" images and, reversing the process, being able to "see" musical scores.
Conclusion
What do we learn from such studies? The main goal here was to have fun. There's much variety in these algorithmsusing other sequences, changing the ranges of the coordinates used, playing with the indexing routines used to map numbers to pixel colors and notes, and so on. Discovering an unexpectedly rich image or melodic score is a source of amazement, and if one learns some Perl along the way, that can't be a bad thing. This isn't to say, though, that the main use of fractals is in entertainmentthe fern image, for example, indicates there is something "real" about fractals. More generally, science is just discovering that many diverse phenomena, such as weather patterns, fluid flow, heart rates, data- transmission noise, and commodity prices, have a common fractal nature. But even in the context of fractal images and music, something pretty profound is going on here. We are taking a mathematical set of numbers and producing aesthetically pleasing pictures and scores. However, quite a bit of programmer insight is still needed even to create thesesimply using the "wrong" coloring index can turn a beautiful image into an eyesore. Perhaps the lesson to take away from this is that while we are crossing, at least tentatively, into understanding aspects of creativity, we are just beginning to understand how complicated, and wondrous, this creativity is.
Acknowledgments
We'd like to thank the Natural Sciences and Engineering Research Council of Canada for support.
(All listings for this article are also available online at http:// www.tpj.com/source/.)
TPJ
Listing 1
#!/usr/bin/perl # mandelbrot.pl - draw the Mandelbrot set use strict; use warnings; use Fractal qw(draw min_max); # xmin, xmax, ymin, and ymax specify the range of coordinates # steps is the number of points used between (xmin .. xmax), # [which is the same number used between (ymin .. ymax)] # maxit is the maximum number of iterations my ($xmin, $xmax, $ymin, $ymax, $steps, $maxit) = (-1.6, 0.6, -1.2, 1.2, 300, 200); my $points = mandelbrot(); draw(pts => $points, file => 'mandelbrot.png', size => $steps); # Mandelbrot set, as defined by # x[n+1] = x[n]*x[n] - y[n]*y[n] + c_x # y[n+1] = 2*x[n]*y[n] + c_y # where (c_x, c_y) are the pixel coordinates # and the iterations begin at (x, y) = (0, 0) sub mandelbrot { my $pts = []; # determine step size between (xmin .. xmax) and (ymin .. ymax) my $xstep = ($xmax - $xmin) / $steps; my $ystep = ($ymax - $ymin) / $steps; for (my $cx=$xmin; $cx<=$xmax; $cx+=$xstep) { for (my $cy=$ymin; $cy<=$ymax; $cy+=$ystep) { my ($oldx, $oldy, $newx, $newy) = (0,0,0,0); # start at x = y = 0 my $i=1; # iteration number before breakout condition is met for ($i=1; $i<$maxit; $i++) { $newx = $oldx*$oldx - $oldy*$oldy + $cx; $newy = 2*$oldx*$oldy + $cy; last if sqrt($newx*$newx + $newy*$newy) > 2; # breakout condition ($oldx, $oldy) = ($newx, $newy); } push @$pts, [$cx, $cy, $i]; } } return $pts; }
Listing 2
#!/usr/bin/perl # julia.pl - draw the Julia set use strict; use warnings; use Fractal qw(draw min_max); # xmin, xmax, ymin, and ymax specify the range of coordinates # steps is the number of points used between (xmin .. xmax) # the same number is used between (ymin .. ymax) # maxit is the maximum number of iterations my ($xmin, $xmax, $ymin, $ymax, $steps, $maxit) = (-1.3, 0.9, -1.2, 1.2, 300, 200); my ($cx, $cy) = (-0.11, -0.9); # used to define the Julia set my $points = julia(); my ($min, $max) = min_max($points); # used in my_rgbindex and my_rgb draw(pts => $points, file => 'julia.png', size => $steps, rgbindex => \&my_rgbindex, rgb => \&my_rgb); # Julia set, as defined by # x[n+1] = x[n]*x[n] - y[n]*y[n] + c_x # y[n+1] = 2*x[n]*y[n] + c_y # where (c_x, c_y) are fixed constants # and the iterations begin at the pixel coodinates sub julia { my $pts = []; # determine step size between (xmin .. xmax) and (ymin .. ymax) my $xstep = ($xmax - $xmin) / $steps; my $ystep = ($ymax - $ymin) / $steps; for (my $x=$xmin; $x<=$xmax; $x+=$xstep) { for (my $y=$ymin; $y<=$ymax; $y+=$ystep) { my ($oldx, $oldy, $newx, $newy) = ($x, $y, 0, 0); # initialization my $i=1; # iteration number before breakout condition is met for ($i=1; $i<$maxit; $i++) { $newx = $oldx*$oldx - $oldy*$oldy + $cx; $newy = 2*$oldx*$oldy + $cy; last if sqrt($newx*$newx + $newy*$newy) > 2; # breakout ($oldx, $oldy) = ($newx, $newy); } push @$pts, [$x, $y, $i]; } } return $pts; } sub my_rgbindex { my $p = shift; return int($p->[2]/$max->[2]*100); } sub my_rgb { my $index = shift; return map{hex} unpack "a2a2a2", sprintf("%02X", int($index*255**3)); }
Listing 3
#!/usr/bin/perl # snowflake.pl - draw a snowflake use strict; use warnings; use Fractal qw(draw min_max); # specify the maximum number of iterations used, # and the size of the desired image my ($maxit, $size) = (500000, 300); # specify the transformations # x[n+1] = a[i]*x[n] + b[i]*y[n] + e[i] # y[n+1] = c[i]*x[n] + d[i]*y[n] + f[i] # where "i" labels the particular transformation my @a = (.75, .5, .25, .25, .25, .25); my @b = (0, -0.5, 0, 0, 0, 0); my @e = (.125, .5, 0, .75, 0, .75); my @c = (0, .5, 0, 0, 0, 0); my @d = (.75, .5, .25, .25, .25, .25); my @f = (.125, 0, .75, .75, 0, 0); # where to begin the iterations at my ($xstart, $ystart) = (0.6, 0.5); my $points = ifs_random(); my ($min, $max) = min_max($points); draw(pts => $points, file => 'snowflake.png', size => $size, rgbindex => \&my_rgbindex, rgb => \&my_rgb); # routine to find the points of an ifs by # choosing the transformation used randomly sub ifs_random { my $pts = []; push @$pts, [$xstart, $ystart, 0]; my ($oldx, $oldy, $newx, $newy) = ($xstart, $ystart, 0, 0); for (my $count = 1; $count < $maxit; $count++) { # choose a transformation randomly my $index = int(rand @a); $newx = $a[$index]*$oldx + $b[$index]*$oldy + $e[$index]; $newy = $c[$index]*$oldx + $d[$index]*$oldy + $f[$index]; push @$pts, [$newx, $newy, $index]; ($oldx, $oldy) = ($newx, $newy); } return $pts; } sub my_rgbindex { my $p = shift; return int($p->[2]/$max->[2]*100); } sub my_rgb { my $index = shift; return (128, 0, 128); }
Listing 4
#!/usr/bin/perl # fern.pl - draw a fern use strict; use warnings; use Fractal qw(draw min_max); # specify the maximum number of iterations used, # and the size of the desired image my ($maxit, $size) = (100000, 300); # specify the transformations # x[n+1] = a[i]*x[n] + b[i]*y[n] + e[i] # y[n+1] = c[i]*x[n] + d[i]*y[n] + f[i] # where "i" labels the particular transformation my @a = (0.85, 0.20, -0.15, 0); my @b = (0.04, -0.26, 0.28, 0); my @e = (0.075, 0.4, 0.575, 0.5); my @c = (-0.04, 0.23, 0.26, 0); my @d = (0.85, 0.22, 0.24, 0.16); my @f = (0.18, 0.045, -0.086, 0); # specify the probablility for which to choose each transformation # and check that they all add up to 1 my %prob = (0 => 0.77, 1 => 0.12, 2 => 0.1, 3 => 0.01); my $sum = 0; $sum += $_ for values %prob; die 'The probablilities have to add up to 1' unless ( abs($sum - 1) < 0.001); # where to begin the iterations at my ($xstart, $ystart) = (0, 0); my $points = ifs_specified(); my ($min, $max) = min_max($points); draw(pts => $points, file => 'fern.png', size => $size, rgbindex => \&my_rgbindex, rgb => \&my_rgb); # routine to find the points of an ifs by choosing the # transformation as specified by a user-set probability sub ifs_specified { my $pts = []; push @$pts, [$xstart, $ystart, 0]; my ($oldx, $oldy, $newx, $newy) = ($xstart, $ystart, 0, 0); # set up bounds to use in test for deciding which # transformation to use my %bounds; my $last = 0; my @indices = sort {$a <=> $b} keys %prob; for (@indices) { $last += $prob{$_}; $bounds{$_} = $last; } for (my $count=1; $count<$maxit; $count++) { # decide which transformation to use my $index; my $r = rand; for (@indices) { if ($r < $bounds{$_}) { $index = $_; last; } } $newx = $a[$index]*$oldx + $b[$index]*$oldy + $e[$index]; $newy = $c[$index]*$oldx + $d[$index]*$oldy + $f[$index]; push @$pts, [$newx, $newy, $index]; ($oldx, $oldy) = ($newx, $newy); } return $pts; } sub my_rgbindex { my $p = shift; return int($p->[2]/$max->[2]*255); } sub my_rgb { my $index = shift; return ($index*20, 128, (255 - $index*10)); }
Listing 5
#!/usr/bin/perl # logistic.pl - hear the logistic map use strict; use warnings; use Fractal qw(compose); # where to start the map at, the parameter mu, and the number of notes my ($xstart, $A, $maxit) = (0.4, 3.82, 256); my $points = logistic(); compose(pts => $points, file => 'logistic.mid'); sub logistic { my $pts = []; my ($oldx, $newx) = ($xstart, 0); push @$pts, $xstart; for (my $i=1; $i<$maxit; $i++) { $newx = $A*$oldx*(1-$oldx); push @$pts, $newx; $oldx = $newx; } return $pts; }
Listing 6
#!/usr/bin/perl # henon.pl - hear the henon map use strict; use warnings; use Fractal qw(compose min_max); # where to start the map at my ($xstart, $ystart) = (0.4, 0.2); # the parameters a and b, and the number of notes my ($a, $b, $maxit) = (1.4, 0.3, 256); # use a chromatic scale my @notes = (55 .. 89); my $points = henon(); my ($min, $max) = min_max($points); compose(pts => $points, file => 'henon.mid', tempo => 100000, patch => 44, noteindex => \&my_noteindex, note => \&my_note); sub henon { my $pts = []; my ($oldx, $oldy, $newx, $newy) = ($xstart, $ystart, 0, 0); push @$pts, $xstart; for (my $i=1; $i<$maxit; $i++) { $newx = $a - $oldx*$oldx + $b*$oldy; $newy = $oldx; push @$pts, $newx; ($oldx, $oldy) = ($newx, $newy); } return $pts; } sub my_noteindex { my $x = shift; $x = ($x - $min) / ($max - $min); return int($x * $#notes); } sub my_note { my $index = shift; return $notes[$index]; }
Listing 7
#!/usr/bin/perl # two_track.pl - hear a point of the Mandelbrot set use strict; use warnings; use Fractal qw(min_max); use MIDI::Simple; my $maxit = 200; # maximum number of notes my ($cx, $cy) = (-0.1, -0.9); # pick a point on the Mandelbrot set my $file = 'mandelbrot.mid'; my ($tempo, $patch_0, $patch_1) = (90000, 33, 66); my @xnotes = (55, 57, 59, 60, 62, 64, 66, 67); my @ynotes = (67, 69, 71, 72, 74, 76, 78, 79); my $points = mandelbrot(); my ($min, $max) = min_max($points); my ($xmin, $xmax, $ymin, $ymax) = ($min->[0], $max->[0], $min->[1], $max->[1]); new_score(); patch_change 0, $patch_0; patch_change 1, $patch_1; set_tempo($tempo); my @subs = (\&track_0, \&track_1); my ($x, $y); foreach (@$points) { ($x, $y) = ($_->[0], $_->[1]); synch(@subs); } write_score($file); sub mandelbrot { my $pts = []; my ($oldx, $oldy, $newx, $newy) = (0,0,0,0); # start at x = y = 0 push @$pts, [$oldx, $oldy]; my $i=1; # iteration number before breakout condition is met for ($i=1; $i<$maxit; $i++) { $newx = $oldx*$oldx - $oldy*$oldy + $cx; $newy = 2*$oldx*$oldy + $cy; last if sqrt($newx*$newx + $newy*$newy) > 2; # breakout # condition ($oldx, $oldy) = ($newx, $newy); push @$pts, [$newx, $newy]; } return $pts; } sub track_0 { my $it = shift; my $note = xnote(); $it->n('c0', 'hn', $note); } sub track_1 { my $it = shift; my $note = ynote(); $it->n('c1', 'hn', $note); } sub xnote { my $xn = ($x - $xmin) / ($xmax - $xmin); my $index = int($xn * $#xnotes); return $xnotes[$index]; } sub ynote { my $yn = ($ymax - $y) / ($ymax - $ymin); my $index = int($yn * $#ynotes); return $ynotes[$index]; }