From 1eedf033c5a5dc3f75e6b15ba91d064e7cb66278 Mon Sep 17 00:00:00 2001
From: Claude HeilandAllen
Date: Sat, 23 Mar 2013 19:51:14 +0000
Subject: [PATCH] paper updates

doc/lyapunovfm.tex  153 +++++++++++++++++++++++++++
1 file changed, 79 insertions(+), 74 deletions()
diff git a/doc/lyapunovfm.tex b/doc/lyapunovfm.tex
index f450813..29f0d26 100644
 a/doc/lyapunovfm.tex
+++ b/doc/lyapunovfm.tex
@@ 25,16 +25,15 @@
\begin{abstract}
\begin{contentsmall}
I consider two coupled oscillators, each modulating the other's
+Consider two coupled oscillators, each modulating the other's
frequency. This system is governed by four parameters: the base
frequency and modulation index for each oscillator. For some
parameter values the system becomes unstable. I use the Lyapunov
exponent to measure the instability. I generate images of the
parameter space, implementing the number crunching on graphics
hardware using OpenGL. I link the mouse position over the displayed
image to realtime audio output, creating an audiovisual browser
+parameter values the system becomes unstable. The Lyapunov
+exponent is used to measure the instability. Images of the
+parameter space are generated, with the number crunching implemented on graphics
+hardware using OpenGL. The mouse position over the displayed
+image is linked to realtime audio output, creating an audiovisual browser
for the 4D parameter space.

\end{contentsmall}
\end{abstract}
@@ 46,25 +45,23 @@ chaos, DSP, GPU
\begin{figure}[H]
\centering
\includegraphics[width=6.25in]{example.png}
+\includegraphics[width=\textwidth]{example.png}
\caption{Example output.}
\label{example}
\end{figure}
\section{Introduction}

In Soft Rock EP~\cite{SoftRockEP} and Soft Rock DVD~\cite{SoftRockDVD} I
+Soft Rock EP~\cite{SoftRockEP} and Soft Rock DVD~\cite{SoftRockDVD}
explored the transitions between order and chaos in coupled FM
oscillators, implemented in Puredata using GridFlow for
visualization of the output waveforms over time. More recently
I developed the idea to map the parameter space on a perceptually
relevant level and in performance choose parameters on the basis
+oscillators. A more recent continuation of
+this project is to make a map of the parameter space of coupled FM oscillators on a perceptually
+relevant level and use it in live performance, choosing parameters on the basis
of desired sound character.
Seeing a bifurcation diagram produced by an analogue Moog
synthesizer~\cite{Slater98} and images of Lyapunov fractals~\cite{Dewdney91},
I decided to apply the latter technique to coupled FM oscillators
in the digital realm.
+A bifurcation diagram produced by an analogue Moog
+synthesizer~\cite{Slater98} and images of Lyapunov fractals~\cite{Dewdney91}
+were inspiration to apply the latter technique to the parameter space
+of coupled FM oscillators in the digital realm.
\newpage
@@ 78,7 +75,11 @@ in the digital realm.
\label{fig:pd}
\end{figure}
Consider the two coupled oscillators in Figure~\ref{fig:pd}.
I formulate this into a mutual recurrence relation:
+Puredata's model of interconnected components each with their
+own internal state maps poorly to GPU architecture. Considering
+the whole system as one and flattening the internal state into
+a single phase space vector leads to the following formulation
+as a mutual recurrence relation:
\begin{equation}\begin{aligned}
x_{n+1} &= \wrap(x_{n} + \inc(f_x + m_x \cos (2 \pi y_{nd}))) \\
y_{n+1} &= \wrap(y_{n} + \inc(f_y + m_y \cos (2 \pi x_{nd})))
@@ 92,15 +93,17 @@ modulation index of each oscillator as a MIDI note number. $\wrap(t)$
performs wrapping into $[0,1)$, with $\lfloor t \rfloor$ being the
flooring operation (the largest integer not greater than $t$).
I'll write the fourdimensional parameter space vector
+The fourdimensional parameter space vector will be written
\begin{equation*}a=(f_x,f_y,m_x,m_y)\end{equation*}
and the $(2d+2)$dimensional phase space vector
\begin{equation*}z=(x_n,y_n,x_{n1},y_{n1},\ldots,x_{nd},y_{nd})\end{equation*}
with sample rate $\SR = 48000$. I'll fix $d=1$ for reasons
explained in Section~\ref{sec:glissues}.
+with sample rate $\SR = 48000$. For reasons
+explained in Section~\ref{sec:glissues}, $d=1$ will be fixed.
\subsection{Lyapunov Exponents}
+Lyapunov exponents can be used to measure the stability (or otherwise)
+of a dynamical system.
A good introduction is found in Chapter~4.3 \emph{Lyapunov Exponent}~\cite{Elert07}.
The definition is covered in Chapter~13.7 \emph{Liapounov exponents and entropies}~\cite{Falconer03}
which also relates it to measures of fractal dimension.
@@ 112,13 +115,16 @@ z_1(0) \to z_0(0)\end{array}} \frac{1}{t} \log{\frac{\leftz_1(t)  z_0(t)\right
An attracting orbit has $\lambda < 0$ and a divergent (chaotic) orbit
has $\lambda > 0$.
The wrapping of phase into $[0,1)$ requires a modified norm:
+A modified norm is required to take into account
+the wrapping of phase into $[0,1)$:
\begin{equation*}\leftz\right_{\wrap} = \sqrt{\sum_{i} {\left(\min(\wrap(z_i), 1  \wrap(z_i))\right)^2}}\end{equation*}
+For example the distance between $0.1$ and $0.9$ is properly $0.2$ (not $0.8$)
+because $0.1$ can be phaseunwrapped to $1.1$.
\subsection{Viewing Planes}
An image is 2D, which requires choosing a subset of the 4D parameter
space to visualize. I chose two particular planes:
+space to visualize. Two particular planes were chosen:
\begin{equation}\begin{aligned}
A_+(a_0, r_0) &= a_0 + r_0 \left(\begin{array}{rr}
1 & 0 \\
@@ 140,15 +146,15 @@ flexible enough to explore the whole 4D space.
The $A_+$ plane varies both oscillators in the same direction, while
the $A_$ plane varies each oscillator in opposite directions. To
center on a particular target point $(f_x, f_y, m_x, m_y)$ one might
use the $A_+$ plane to center on
+use the $A_+$ plane to center on the midpoint
\begin{equation*}\left(\frac{f_x + f_y}{2}, \frac{f_x + f_y}{2}, \frac{m_x + m_y}{2}, \frac{m_x + m_y}{2}\right)\end{equation*}
and then switch to the $A_$ plane to break the $(x,y)$ symmetry.
\section{Implementation}
I used OpenGL~\cite{OpenGL} and OpenGL Shading Language~\cite{GLSL}, with
GLUT~\cite{GLUT} for windowing and input event handling, and JACK~\cite{JACK} for
audio output.
+The implementation uses OpenGL~\cite{OpenGL} and OpenGL Shading Language~\cite{GLSL}
+for computation and graphical rendering, GLUT~\cite{GLUT} for windowing
+and input event handling, and JACK~\cite{JACK} for audio output.
\subsection{Introduction to Modern OpenGL}
@@ 167,34 +173,35 @@ a framebuffer.
\subsection{Computation Overview}
To render an image I first fill a texture with $(u,v)$ coordinates
using a framebuffer object and a fragment shader. I copy this texture
to a vertex buffer object, interleaving an initial phase space vector
$z=(0,0,0,0)$ and Lyapunov exponent statistics vector $l=(0,0,0,0)$.
+To render an image a texture is first filled with $(u,v)$ coordinates
+using a framebuffer object and a fragment shader. This texture is copied
+to a vertex buffer object, interleaved with an initial phase space vector
+$z=(0,0,0,0)$ and Lyapunov exponent statistics vector $l=(0,0,0,0)$ for
+each point.
Using a vertex shader, I
calculate $a$ from $(u,v)$ using Equation~\ref{eq:view}, and then
compute a rough estimate of the Lyapunov exponent using
+Using a vertex shader, $a$ is calculated from $(u,v)$ using Equation~\ref{eq:view}, and then
+a rough estimate of the Lyapunov exponent is computed using
Equation~\ref{eq:lexp} by perturbing $z_1(0) = z_0(0) + \delta$ with
$\delta$ small and performing $t = 256$ iterations of Equation~\ref{eq:osc}.
I discard the first few repetitions and those resulting in $\infty$,
and accumulate the rough $\lambda$ estimates in $l$.
+The first few repetitions are discarded, along with those resulting in $\infty$,
+and the rough $\lambda$ estimates are accumulated in $l$.
Between each repetition I compact the working set using a geometry
shader, plotting points whose mean Lyapunov exponent estimate changed
very little during the previous step. I keep the others to refine
+Between each repetition the working set is compacted using a geometry
+shader. Points whose mean Lyapunov exponent estimate changed
+very little during the previous step are plotted and removed from the
+working set. The other points are kept to be refined
further, directing the computational effort on the points that
need it most:~those slow to converge.
To ensure user interface responsiveness, the computation is amortized
over several frames. I divide the target frame period by the measured
+over several frames. The target frame period is divided by the measured
time for one repetition to compute how many repetitions to perform that
frame. The repetitionsperframe increases as the working set becomes
smaller.
\subsection{Noise Increases Stability}
At the end of each repetition I keep $z_1$ instead of $z_0$. This
+At the end of each repetition $z_1$ is kept instead of $z_0$. This
effectively adds a small amount of noise, counterintuitively
increasing stability. Noise allows more of the phase space to be
explored, and makes it more likely for the perturbed orbit to reach an
@@ 202,13 +209,13 @@ attracting part of the phase space.
\subsection{Dither Increases Quality}
To reduce grid sampling artifacts, I perturb $(u,v)$ within the bounds
+To reduce grid sampling artifacts, $(u,v)$ is perturbed within the bounds
of its corresponding pixel before calculating the $a$ parameter vector
for each repetition.
\section{Results}
\begin{figure*}[ht]
+\begin{figure*}[ht!]
\centering
\begin{tabular}{cc}
\subfigure[$A_+((120, 120, 0, 0), 72)$]{\includegraphics[width=3in]{example01.png}\label{ex:01}} &
@@ 216,11 +223,22 @@ for each repetition.
\subfigure[$A_+((95.2, 95.2, 32.6, 32.6), 4.5)$]{\includegraphics[width=3in]{example03.png}\label{ex:03}} &
\subfigure[$A_((95.2, 95.2, 32.6, 32.6), 4.5)$]{\includegraphics[width=3in]{example04.png}\label{ex:04}} \\
\subfigure[$A_+((151.57, 151.57, 1.64, 1.64), 0.07)$]{\includegraphics[width=3in]{example05.png}\label{ex:05}} &
\subfigure[$A_((151.57, 151.57, 1.64, 1.64), 0.07)$]{\includegraphics[width=3in]{example06.png}\label{ex:06}}
+\subfigure[$A_((151.57, 151.57, 1.64, 1.64), 0.07)$]{\includegraphics[width=3in]{example06.png}\label{ex:06}} \\
+\subfigure[$A_((117.0, 148.4, 20.4, 2.7), 1.8)$]{\includegraphics[width=3in]{example07.png}\label{ex:07}} &
+\subfigure[$A_+((103.65, 108.41, 33.42, 10.93), 0.14)$]{\includegraphics[width=3in]{example09.png}\label{ex:09}}
\end{tabular}
\caption{Example images. Darker shades are stable, lighter shades chaotic.}
\label{fig:examples1}
\end{figure*}
+%\begin{figure}[ht]
+%\centering
+%\subfigure[$A_((117.0, 148.4, 20.4, 2.7), 1.8)$]{\includegraphics[width=3in]{example07.png}\label{ex:07}}
+%\subfigure[$A_((141.46, 146.22, 22.76, 0.27), 0.14)$]{\includegraphics[width=3in]{example08.png}\label{ex:08}}
+%\subfigure[$A_+((103.65, 108.41, 33.42, 10.93), 0.14)$]{\includegraphics[width=3in]{example09.png}\label{ex:09}}
+%\subfigure[$A_((89.8, 137.5, 17.5, 7.1), 3.7)$]{\includegraphics[width=3in]{example10.png}\label{ex:10}}
+%\caption{More examples.}
+%\label{fig:examples2}
+%\end{figure}
\subsection{Examples}
@@ 241,22 +259,12 @@ When $f_x=f_y$ and $m_x=m_y$ the $A_+$ plane has mirror symmetry about its
horizontal axis, and the $A_$ plane has twofold rotational symmetry
about its centre.
Breaking the symmetry and setting $f_x \ne f_y$ or $m_x \ne m_y$ leads
to diverse forms, shown in Figure~\ref{fig:examples2}. In particular Figure~\ref{ex:09} has shapes that
+to diverse forms. In particular Figure~\ref{ex:09} has shapes that
resemble those of Lyapunov space images of the logistic map.
\begin{figure}[ht]
\centering
\subfigure[$A_((117.0, 148.4, 20.4, 2.7), 1.8)$]{\includegraphics[width=3in]{example07.png}\label{ex:07}}
\subfigure[$A_((141.46, 146.22, 22.76, 0.27), 0.14)$]{\includegraphics[width=3in]{example08.png}\label{ex:08}}
\subfigure[$A_+((103.65, 108.41, 33.42, 10.93), 0.14)$]{\includegraphics[width=3in]{example09.png}\label{ex:09}}
\subfigure[$A_((89.8, 137.5, 17.5, 7.1), 3.7)$]{\includegraphics[width=3in]{example10.png}\label{ex:10}}
\caption{More examples.}
\label{fig:examples2}
\end{figure}

\subsection{Interactive Explorer}
I implemented an interactive audiovisual explorer for the parameter
+The implementation is an interactive audiovisual explorer for the parameter
space of coupled FM oscillators.
Clicking with the mouse zooms the view about the clicked point.
The left button (or scroll up) zooms in, the right button (or scroll
@@ 273,7 +281,7 @@ a reference frame for chosing parameters to audition by moving the mouse.
\subsection{Original Intent}
I had originally experimented with one Puredata batch mode instance per CPU
+Earlier experiments used one Puredata batch mode instance per CPU
core each sending analysis data to a realtime Puredata instance.
The analysis used various methods (including FFT for spectral
statistics and the sigmund external for pitch tracking) to classify
@@ 292,14 +300,14 @@ attribute to four components with the maximum number of attributes
typically limited to sixteen. This totals 64 floats per vertex, 6 of
which are needed for the pixel coordinates and Lyapunov exponent
statistics accumulation. Therefore using OpenGL imposes a limit
$d < 28$. For comparison my original experiments in
+$d < 28$. For comparison the original experiments in
\emph{Soft Rock EP} used Puredata's default block size of 64, with $d = 32$.
Moreover, increasing $d$ increases video memory consumption. With
the maximum $d = 27$, browsing at $1920 \times 1080$ resolution would
require over 1GB.
My future work on this project will look into using OpenCL, which
provides a heterogenous CPU and GPU computation framework. I hope
+Future work on this project will look into using OpenCL, which
+provides a heterogenous CPU and GPU computation framework, in the hope that
it will avoid the inherent awkwardness of abusing OpenGL shaders to
perform calculations.
@@ 308,13 +316,13 @@ perform calculations.
While the implementation works as intended, with $d=1$
the sound is nowhere near as rich and varied as with $d=32$. With
small $d$ there is much more very high frequency content in
interestinglooking regions. I've not found any regions of the
+interestinglooking regions. There seem to be few if any regions of the
parameter space with both interesting appearance and palatable audio
frequencies at $d=1$, while with high $d$ there are
parameters that generate sounds that fluctuate intermittently between
smooth tones and noise. I haven't been able to visualize
with high $d$ to know whether their neighbourhoods
look as interesting as they sound.
+smooth tones and noise. Visualization with high $d$ has not been
+possible so far, so whether their neighbourhoods
+look as interesting as they sound remains an openquestion .
Unfortunately, heavy use of the GPU in the interactive browser can
block the operating system for too long and cause audible glitches
@@ 325,24 +333,21 @@ improve, allowing use of the browser in a live situation.
Despite these shortcomings, I think the images look good.
I plan to render a selection at high resolution and print
postcards and posters. For huge images I divide the image plane into tiles and
+postcards and posters. For huge images it is possible to divide the image plane into tiles and
compute each tile in succession, finally combining the pieces
into one large picture.
There is also scope for video work, moving and rotating the viewing
plane through the 4D parameter space, with different shapes forming
and collapsing over time. My rough benchmarks take 510 seconds
per frame at $1920 \times 1080$, so for now I'll wait until faster cheaper
+and collapsing over time. Rough benchmarks take 510 seconds
+per frame at $1920 \times 1080$, so it seems sensible to wait until faster cheaper
graphics cards become available.
\section{Obtaining the Implementation}
I wrote the implementation on GNU/Linux Debian Wheezy using
proprietary drivers running on a quadcore AMD64 processor
with NVIDIA GTX~550Ti graphics card.

The source code for the
implementation is available at:
+The implementation was written on GNU/Linux Debian Wheezy running on a quadcore AMD64 processor
+with NVIDIA GTX~550Ti graphics card using
+proprietary drivers. The source code is available at:
\url{https://gitorious.org/maximus/lyapunovfm}

2.1.4