A script for solving a Lorenz-like initial value problem and export the time series of each coordinate as a sound file

hearChaos.m

%% Hearing chaos
% How chaos sounds? In this script we will create a chaotic time series and
% translate it to a sound signal
%
% More information at: <http://mappingignorance.org/2016/05/09/the-sound-of-chaos/>

%% Harvest chaos from the Lorenz equation
% The Lorenz equation is the classical example of deterministic dynamical
% system giving rise to chaos. The equation reads:
% 
% $$\dot x = d \left( y - x \right)$$
%
% $$\dot y = r x - y - e x z$$
%
% $$\dot z = -bz + e x y $$
%
% With the following set of parameters
%
% $$d = 10, r = 28, b = \frac{8}{3}, e = 1$$
%
% we will find chaos for almost all initial conditions, being the only
% exception the coordinate origin (an unstable steady state).
%
% So now we can pose the problem in Matlab code

d = 10;
r = 28;
b = 8/3;
e = 1;

fun = @(t,x) [d*(x(2) - x(1));
              r*x(1) - x(2) - e.*x(1).*x(3);
              -b*x(3) + e.*x(1).*x(2)];

%% Solve numerically
% The Lorenz equation is not analytically solvable, so we will solve it
% numerically for a given initial condition
x0 = [0; 1; 0]; % Set initial condition
tSpan = [0 250]; % Set time span
sol = ode45(fun, tSpan, x0); % Create solver object

% Resample with the desired time step
tStep = .5e-2;
t = tSpan(1):tStep:tSpan(end);
x = deval(sol, t);

%% Plot the Lorenz attractor
% We can use our numerical solution to plot the world-famous Lorenz
% attractor
plot3(x(1,:), x(2,:), x(3,:)); xlabel('x'); ylabel('y'); zlabel('z'); axis equal;
view(0, 0); title('Frontal view'); snapnow;
view(90,0); title('Side view'); snapnow;
view(0,90); title('Upper view'); snapnow;

%% Generate an animated version
h = plot3(x(1,:), x(2,:), x(3,:));
xlabel('x'); ylabel('y'); zlabel('z');
axis vis3d;
set(gca,'xticklabel',[]); set(gca,'yticklabel',[]); set(gca,'zticklabel',[]);
m = struct('cdata', [], 'colormap', []);
for i = 1:1:360
    view(i-1, 30);

    im = frame2im(getframe);
    [A, map] = rgb2ind(im, 256); 
	if i == 1;
		imwrite(A,map, 'lorenz.gif', 'gif', 'LoopCount', Inf, 'DelayTime', .1);
	else
		imwrite(A, map, 'lorenz.gif', 'gif', 'WriteMode', 'append', 'DelayTime', .1);
	end
end

%% Extract a one-dimensional time series
% An easy way to extract a one-dimensional time series out of the Lorenz
% attractor is focusing only on one coordinate:
timeSeriesX = x(1,:);
timeSeriesY = x(2,:);
timeSeriesZ = x(3,:);

%% View the time series
plot(timeSeriesX); xlabel('t'); ylabel('x'); snapnow;
plot(timeSeriesY); xlabel('t'); ylabel('y'); snapnow;
plot(timeSeriesZ); xlabel('t'); ylabel('z'); snapnow;

%% Export to audio files
%
% The _sonification_ process follows this steps:
% 
% * Normalize the signal to the range [-1,1]
% * Pick a sampling frequency
% * Export as _.wav_
%

% Normalization
timeSeriesXNormalized = timeSeriesX./max(abs(timeSeriesX));
timeSeriesYNormalized = timeSeriesY./max(abs(timeSeriesY));
timeSeriesZNormalized = timeSeriesZ./max(abs(timeSeriesZ));

% Sampling
sampling_freq = 30e3;

% Exporting
audiowrite('x.wav', timeSeriesXNormalized, sampling_freq);
audiowrite('y.wav', timeSeriesYNormalized, sampling_freq);
audiowrite('z.wav', timeSeriesZNormalized, sampling_freq);

%% For direct playback use
% 
%   soundsc(timeSeriesX, sampling_freq);
% 

%% Credits
% Pablo Rodríguez-Sánchez 
% 
% <https://sites.google.com/site/pablorodriguezsanchez/>
%
% April 2016