Skip to content

Instantly share code, notes, and snippets.

@PabRod PabRod/hearChaos.m
Last active May 12, 2016

Embed
What would you like to do?
A script for solving a Lorenz-like initial value problem and export the time series of each coordinate as a sound file
%% 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.