%% 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: %% 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 % % % % April 2016