1. Introduction
This document describes how to control an external model using the orientation of your Muse device using the Bluetooth Low Energy (BLE) functionalities provided by MATLAB®.
This example is based on the standard Muse communication protocol [A] and MATLAB® official documentation, as described in our quick start guide [B]:
2. Collect Data from your Muse
This example shows how to collect and plot orientation quaternion from your Muse device using the BLE communication support provided by MATLAB® functions.
First, check that the BLE device is powered on and support connections by finding it in MATLAB® using the blelist function. After the device has been found in MATLAB®, connect to it by calling the ble function. Specify the name of the device if the name is unique or specify the device address.
%% Discover and connect to Muse device
devlist = blelist;
muse = ble("muse_roberto");
Access the Muse v3 Custom Service and relative command and data characteristics, to get the corresponding UUIDs.
%% Access custom service and characteristics of interest cmd = characteristic(muse, "C8C0A708-E361-4B5E-A365-98FA6B0A836F", ... "D5913036-2D8A-41EE-85B9-4E361AA5C8A7"); dat = characteristic(muse, "C8C0A708-E361-4B5E-A365-98FA6B0A836F", ... "09BF2C52-D1D9-C0B7-4145-475964544307");
Next, enable subscription to both command and data characteristics.
%% Subscribe to characteristic notify subscribe(cmd); subscribe(dat);
Check device status and get its configuration, before starting acquisition.
%% Get device status (i.e., just to check setup consistency) % GET_STATE command code (hex): 0x82 write(cmd, [130 0]); cmd_response = read(cmd); state = displayDeviceState(cmd_response);
Then start orientation data acquisition as follows.
%% Start data acquisition in STREAMING mode
% To stop acquisition at runtime, write the following command on Command
% Window: write(cmd, [2 1 2]). It set the device in SYS_IDLE state.
if (state == 2)
% SET_STATE command code (hex): 0x02
% streaming type: 0x08 (continuous streaming)
% acquisition mode: 0x02 (accelerometer only)
% acquisition freq: 0x01 (25 Hz)
write(cmd, [2 5 8 16 0 0 1]);
else
disp('Acquisition already running.');
end
The data acquisition is managed through a dedicated callback function assigned to the DataAvailableFcn property of the data characteristic.
%% Manage data acquisition, decoding and visualization
% Setup figure for real time plot
fig=figure;
% Fix the axes scaling, and set a nice view angle
axis('image'); xlabel('X'), ylabel('Y'), zlabel('Z')
view([40 20]);
% grid on, grid minor,
ax = gca;
% axis equal
ax.Color = 'none';
ax.XColor = 'none';
ax.YColor = 'none';
ax.ZColor = 'none';
arrX = mArrow3([0 0 0],[20 0 0], 'facealpha', 0.5, 'color', 'b', 'stemWidth', 0.4);
arrY = mArrow3([0 0 0],[0 20 0], 'facealpha', 0.5, 'color', 'r', 'stemWidth', 0.4);
arrZ = mArrow3([0 0 0],[0 0 20], 'facealpha', 0.5, 'color', 'y', 'stemWidth', 0.4);
ax.XLim(1) = -25;
ax.XLim(2) = 25;
ax.YLim(1) = -25;
ax.YLim(2) = 25;
ax.ZLim(1) = -25;
ax.ZLim(2) = 25;
hold all
% Assign data callback
dat.DataAvailableFcn = @displayCharacteristicData;
The result will be a reference system represented as a set of three arrows. Such a representation of the reference frame will move based on device orientation.
3. Insights
As regard the message decoding, the use of dedicated standard function based on specific commands as well as callback functions are the most straight forward approach that can be adopted. As an example, we provide in the following some code snippets that show how to decode orientation quaternion received from the device.
Refers to [C] for any further detail about raw data decoding and management.
%% Callback function on data characteristic to manage quaternion decoding
function displayCharacteristicData(src,evt)
global q_curr q_prev q_curr_list i
global arrX arrY arrZ
% Read data
data = read(src,'oldest');
% Decode imaginary components (in this example we ignore the first
% 8 bytes that represents the timestamp)
qi = [data(9) data(10)];
qi_val = typecast(uint8(qi),'int16');
qi_val = cast(qi_val,'single') / 32767;
qj = [data(11) data(12)];
qj_val = typecast(uint8(qj),'int16');
qj_val = cast(qj_val,'single') / 32767;
qk = [data(13) data(14)];
qk_val = typecast(uint8(qk),'int16');
qk_val = cast(qk_val,'single') / 32767;
% Compute (unit) quaternion real component
qw_val = 1.0;
if ((1 - ((qi_val*qi_val) + (qj_val*qj_val) + (qk_val*qk_val))) > 0)
qw_val = sqrt(cast(1 - ((qi_val*qi_val) + (qj_val*qj_val) + (qk_val*qk_val)),'single'));
end
% Update orientation quaternion variables
q_prev = q_curr;
q_curr = [qw_val, qi_val, qj_val, qk_val];
q_curr_list(i,:) = q_curr;
i = i + 1;
% compute angular difference from quaternions
qdelta = quaternionProduct(q_curr,quaternionConjugate(q_prev));
Rquat = quat2rotationmatrix(qdelta);
% apply transformation to the model
tmp = (Rquat*((arrX.Vertices - arrX.Vertices(1,:)).')).';
set(arrX,'Vertices',tmp);
tmp = (Rquat*((arrY.Vertices - arrY.Vertices(1,:)).')).';
set(arrY,'Vertices',tmp);
tmp = (Rquat*((arrZ.Vertices - arrZ.Vertices(1,:)).')).';
set(arrZ,'Vertices',tmp);
pause(0.005);
% Display and/or plot data
% Uncomment to display the numeric values
% fprintf('%f\t%f\t%f\t%f\n',qw_val,qi_val,qj_val,qk_val);
end