Table Of Content

Lyapunov Space of Coupled FM Oscillators Claude Heiland-Allen [email protected] Abstract 1 Introduction Soft Rock EP [ClaudiusMaximus, 2005] and Consider two coupled oscillators, each modulating Soft Rock DVD [ClaudiusMaximus, 2006] ex- the other’s frequency. This system is governed by plored the transitions between order and chaos fourparameters: thebasefrequencyandmodulation indexforeachoscillator. Forsomeparametervalues in coupled FM oscillators. A more recent con- the system becomes unstable. The Lyapunov ex- tinuationofthisprojectistomakeamapofthe ponent is used to measure the instability. Images of parameter space of coupled FM oscillators on a theparameterspacearegenerated,withthenumber perceptuallyrelevantlevelanduseitinliveper- crunching implemented on graphics hardware using formance, choosing parameters on the basis of OpenGL.Themousepositionoverthedisplayedim- desired sound character. age is linked to realtime audio output, creating an A bifurcation diagram produced by an ana- audio-visual browser for the 4D parameter space. logue Moog synthesizer [Slater, 1998] and images of Lyapunov fractals [Dewdney, 1991] were Keywords inspiration to apply the latter technique to the parameter space of coupled FM oscillators in chaos, DSP, GPU the digital realm. Figure 1: Example output. 2 Formulation 2.2 Lyapunov Exponents 2.1 Coupled FM Oscillators Lyapunovexponentscanbeusedtomeasurethe stability (or otherwise) of a dynamical system. A good introduction is found in Chapter 4.3 Lyapunov Exponent [Elert, 2007]. The defini- tion is covered in Chapter 13.7 Liapounov exponents and entropies [Falconer, 2003] which also relates it to measures of fractal dimension. The Lyapunov exponent λ measures diver- gence in phase space: |z1(t)−z0(t)| ≈ eλt|z1(0)−z0(0)| 1 |z1(t)−z0(t)| λ = lim log (2) t → ∞ t |z1(0)−z0(0)| z1(0) → z0(0) An attracting orbit has λ < 0 and a divergent (chaotic) orbit has λ > 0. Figure 2: Coupled FM oscillators in Pure-data. A modified norm is required to take into ac- count the wrapping of phase into [0,1): Consider the two coupled oscillators in Fig- ure 2. Pure-data’s model of interconnected components each with their own internal state |z| = (min(%(z ),1−%(z )))2 % i i maps poorly to GPU architecture. Consider- s i X ing the whole system as one and flattening the internal state into a single phase space vector For example the distance between 0.1 and 0.9 is leads to the following formulation as a mutual properly 0.2 (not 0.8) because 0.1 can be phase- recurrence relation: unwrapped to 1.1. xn+1 = %(xn+I(fx+mxcos(2πyn−d))) 2.3 Viewing Planes (1) yn+1 = %(yn+I(fy +mycos(2πxn−d))) An image is 2D, which requires choosing a sub- set of the 4D parameter space to visualize. Two where particular planes were chosen: 440 t−69 %(t) = t−⌊t⌋, I(t) = 2 12 SR 1 0 1 0 u Here xn, yn is the phase of each oscillator at A+(a0,r0) = a0+r0 0 1  v time step n, d is a delay measured in sam- (cid:18) (cid:19) 0 1 ples, f , f is the base frequency of each os-   x y   cillator as a MIDI note number, and m , m 1 0 x y is the modulation index of each oscillator as a −1 0 u MIDI note number. %(t) performs wrapping A−(a0,r0) = a0+r0 0 1  v (cid:18) (cid:19) into [0,1), with ⌊t⌋ being the flooring operation 0 −1   (the largest integer not greater than t).   (3) The four-dimensional parameter space vector where (u,v) is the coordinates of the pixel, a0 will be written is the centre of the view, and r0 is the radius of the view. a = (f ,f ,m ,m ) x y x y These planes were chosen because they are simple, while still being flexible enough to ex- and the (2d+2)-dimensional phase space vector plore the whole 4D space. The A+ plane varies z = (xn,yn,xn−1,yn−1,...,xn−d,yn−d) both oscillators in the same direction, while the A− plane varies each oscillator in opposite di- with sample rate SR = 48000. For reasons ex- rections. To center on a particular target point plained in Section 5.2, d = 1 will be fixed. (fx,fy,mx,my) one might use the A+ plane to center on the midpoint plotted and removed from the working set. The other points are kept to be refined further, di- fx+fy fx+fy mx+my mx+my recting the computational effort on the points , , , 2 2 2 2 that need it most: those slow to converge. (cid:18) (cid:19) To ensure user interface responsiveness, the and then switch to the A− plane to break the computation is amortized over several frames. (x,y) symmetry. The target frame period is divided by the measured time for one repetition to compute how 3 Implementation many repetitions to perform that frame. The The implementation uses OpenGL [Segal, 2013] repetitions-per-frame increases as the working and OpenGL Shading Language [Kessenich, set becomes smaller. 2013] for computation and graphical rendering, 3.3 Noise Increases Stability GLUT [Kilgard, 1996] for windowing and input event handling, and JACK [Davis, 2013] for au- At the end of each repetition z1 is kept instead dio output. of z0. This effectively adds a small amount of noise, counter-intuitively increasing stability. 3.1 Introduction to Modern OpenGL Noise allows more of the phase space to be ex- Modern OpenGL has a programmable shader plored, and makes it more likely for the per- pipeline. Vertex attributes are read from vertex turbed orbit to reach an attracting part of the buffers and processed by vertex shaders. The phase space. outputs of the vertex shader (called varyings) 3.4 Dither Increases Quality are further manipulated by an optional geome- To reduce grid sampling artifacts, (u,v) is per- try shader stage. Geometry shaders can output turbed within the bounds of its corresponding a different vertex count to their input count, pixel before calculating the a parameter vector whereas vertex shaders are one-in one-out. The for each repetition. result of the geometry shader can be captured into another vertex buffer using transform feed- 4 Results back. Following the geometry shader the prim- itives (points or triangles) are rasterized, and 4.1 Examples varyings interpolated across each primitive. Fi- Figure 3(a) shows the initial view on starting nally a fragment shader takes these values and the interactive browser. Low frequencies to the computes the colour at that pixel. The output left are stable even at high modulation index ofafragmentshadercanbecapturedbyattach- away from the central axis. High frequencies to ing a texture to a framebuffer. the right become chaotic at progressively lower modulation index. (b) shows the A− plane at 3.2 Computation Overview the same location. (c) shows bands alternat- To render an image a texture is first filled with ing between stability and chaos. The bands be- (u,v) coordinates using a framebuffer object come distorted and collapse as the modulation andafragmentshader. Thistextureiscopiedto index and frequency increase. (d) shows its A− avertexbufferobject,interleavedwithaninitial plane,bandsbecomerings. Whenthefrequency phasespacevectorz = (0,0,0,0)andLyapunov is greatly increased, the shapes become more exponentstatisticsvectorl = (0,0,0,0)foreach intricate. (e) exhibits spirals of stability, with point. similar spirals in the A− plane in (f). Using a vertex shader, a is calculated from Whenfx = fy andmx = my theA+planehas (u,v) using Equation 3, and then a rough es- mirror symmetry about its horizontal axis, and timate of the Lyapunov exponent is computed the A− plane has two-fold rotational symmetry usingEquation2byperturbingz1(0) = z0(0)+δ about its centre. Breaking the symmetry and with δ small and performing t = 256 iterations setting f 6= f or m 6= m leads to diverse x y x y of Equation 1. The first few repetitions are dis- forms. InparticularFigure3(h)hasshapesthat carded, along with those resulting in −∞, and resemble those of Lyapunov space images of the the rough λ estimates are accumulated in l. logistic map. Between each repetition the working set is 4.2 Interactive Explorer compacted using a geometry shader. Points whose mean Lyapunov exponent estimate The implementation is an interactive audio- changed very little during the previous step are visual explorer for the parameter space of cou- (a) A+((120,120,0,0),72) (b) A−((120,120,0,0),72) (c) A+((95.2,95.2,32.6,32.6),4.5) (d) A−((95.2,95.2,32.6,32.6),4.5) (e) A+((151.57,151.57,1.64,1.64),0.07) (f) A−((151.57,151.57,1.64,1.64),0.07) (g) A−((117.0,148.4,20.4,2.7),1.8) (h) A+((103.65,108.41,33.42,10.93),0.14) Figure 3: Example images. Darker shades are stable, lighter shades chaotic. pled FM oscillators. Clicking with the mouse 5.3 Audio Issues zoomstheviewabouttheclickedpoint. Theleft While the implementation works as intended, button (or scroll up) zooms in, the right button with d = 1 the sound is nowhere near as rich (or scroll down) zooms out, the middle button and varied as with d = 32. With small d there centers the view on the target point. Pressing is much more very high frequency content in the TAB key toggles between the A+ and A− interesting-looking regions. There seem to be planesinEquation3,andF11togglesfullscreen few if any regions of the parameter space with operation. both interesting appearance and palatable au- While the GPU simulates and analyses one diofrequenciesatd = 1, whilewithhighdthere oscillator pair per pixel, the CPU simulates one are parameters that generate sounds that fluc- oscillator pair with a determined from the pixel tuate intermittently between smooth tones and under the mouse pointer. The image acts as a noise. Visualization with high d has not been map, a reference frame for chosing parameters possiblesofar,sowhethertheirneighbourhoods to audition by moving the mouse. look as interesting as they sound remains an openquestion . 5 Conclusions Unfortunately, heavy use of the GPU in the 5.1 Original Intent interactive browser can block the operating sys- Earlier experiments used one Pure-data batch tem for too long and cause audible glitches mode instance per CPU core each sending anal- (JACK xruns). This situation may change as ysis data to a realtime Pure-data instance. The free drivers continue to improve, allowing use of analysis used various methods (including FFT the browser in a live situation. for spectral statistics and the sigmund exter- nal for pitch tracking) to classify points into 5.4 Pretty Pictures pitched (ordered, stable) or unpitched (chaotic, Despite these shortcomings, I think the images unstable) with measures of distortion or noisi- look good. I plan to render a selection at high ness. Sadly this approach proved impractical as resolution and print postcards and posters. For it achieved only tens of pixels per second, even huge images it is possible to divide the image with a fast multi-core CPU, and porting these plane into tiles and compute each tile in succes- signal analysis algorithms to massively-parallel sion, finally combining the pieces into one large programmable graphics hardware seemed to be picture. too difficult. There is also scope for video work, moving 5.2 OpenGL Issues and rotating the viewing plane through the 4D parameter space, with different shapes forming The current implementation is hardcoded with and collapsing over time. Rough benchmarks delay d = 1 and would be very awkward to take 5-10 seconds per frame at 1920 × 1080, generalize. OpenGL architecture limits each so it seems sensible to wait until faster cheaper vertex attribute to four components with the graphics cards become available. maximum number of attributes typically lim- ited to sixteen. This totals 64 floats per ver- 6 Obtaining the Implementation tex, 6 of which are needed for the pixel coordinates and Lyapunov exponent statistics accu- The implementation was written on mulation. Therefore using OpenGL imposes a GNU/Linux Debian Wheezy running on a limit d < 28. For comparison the original ex- quad-core AMD64 processor with NVIDIA perimentsinSoft Rock EP usedPure-data’sde- GTX 550Ti graphics card using propri- fault block size of 64, with d = 32. Moreover, etary drivers. The source code is available at: increasing d increases video memory consump- https://gitorious.org/maximus/lyapunov-fm tion. With the maximum d = 27, browsing at 1920×1080 resolution would require over 1GB. 7 Acknowledgements Future work on this project will look into using OpenCL, which provides a heterogenous Thanks to the anonymous reviewers for their CPU and GPU computation framework, in the constructive criticism on a number of issues, hope that it will avoid the inherent awkward- and to Rob Canning, Adnan Hadzi, and Joanne nessofabusingOpenGLshaderstoperformcal- Seale for their helpful feedback on earlier ver- culations. sions of this paper. References ClaudiusMaximus. 2005. Soft Rock EP. http://archive.org/details/ ClaudiusMaximus - Soft Rock EP. ClaudiusMaximus. 2006. Soft Rock DVD. http://archive.org/details/ ClaudiusMaximus - Soft Rock DVD. Paul Davis. 2013. The JACK Audio Connec- tion Kit. http://jackaudio.org. A. K. Dewdney. 1991. Mathematical Recre- ations: Leaping into Lyapunov Space. Scien- tific American, 265:178–180. Glenn Elert. 2007. The Chaos Hypertextbook. http://hypertextbook.com/chaos/. Kenneth Falconer. 2003. Fractal Geometry: Mathematical Foundations and Applications, Second Edition. Wiley. John Kessenich. 2013. The OpenGL Shading Language. http://www.opengl.org/registry/doc/ GLSLangSpec.4.30.8.pdf. Mark J. Kilgard. 1996. The OpenGL Utility Toolkit (GLUT) Programming Interface. http://www.opengl.org/documentation/ specs/glut/glut-3.spec.pdf. Mark Segal. 2013. The OpenGL Graphics System: A Specification. http://www.opengl.org/registry/doc/ glspec43.core.20130214.pdf. Dan Slater. 1998. Chaotic Sound Synthesis. Computer Music Journal, 22(2):12–19.