MEG activity of Jean Louis was recorded while he was trying to fall asleep. For this purpose he was sleep deprived the previous night.
I then took a ~24 minute interval (starting at about 25 minutes after start) of recording and tried to sonify his oscillatory brain patterns. In other words, find a way to listen to his brainwaves in ‘real’ time. With ‘real’ time I mean not time-dilated or speeded up. Brain activity that is measured with EEG or MEG oscillates at frequencies below the hearing range, say between about 0.5 to 80 Hz. For this purpose of transposing brain activity to the human audible spectrum I did the following:
- Calculate the hilbert envelope or three frequencies: theta, alpha and gamma, around 7, 10 and 40Hz respectively.
- Average channels for five different regions of the MEG helmet: front, back, left, right and center.
- Multiply the averaged hilvert envelope with the amplitude of three different tones for each frequency (2000, 3000 and 4000Hz).
- Upsampled to 10000Hz and write to WAV files.
You can listen to some examples here:
Alpha at the back of helmet:
Theta at the front of the helmet
Gamma at the top of the helmet
Note the dynamics! After so many years of analysing neuronal oscillations I never listened to them. It certainly is inspiring thoughts about the brain.
You can see the code here:
cfg = [];
cfg.dataset = 'sleep-1_tsss_mc.fif';
cfg.continuous = 'yes';
cfg.hpfilter = 'no';
cfg.detrend = 'no';
cfg.continuous = 'yes';
cfg.demean = 'yes';
cfg.dftfilter = 'yes';
cfg.dftfreq = [50 100 150];
cfg.lpfilter = 'yes';
cfg.lpfreq = 150;
cfg.hpfilter = 'yes';
cfg.hpfreq = 2;
data = ft_preprocessing(cfg);
cfg = [];
cfg.viewmode = 'vertical';
cfg.channel = 30:31;
cfg.layout = 'neuromag306all';
ft_databrowser(cfg,data);
cfg = [];
date_timelocked = ft_timelockanalysis(cfg,data);
clear data;
chanlist_front = {'MEG0311', 'MEG0322', 'MEG0323', 'MEG0321', 'MEG0333', 'MEG0332', 'MEG0331', 'MEG0513', 'MEG0512', 'MEG0511', 'MEG0523', 'MEG0521', 'MEG0532', 'MEG0533', 'MEG0531', 'MEG0542', 'MEG0543', 'MEG0541', 'MEG0613', 'MEG0612', 'MEG0611', 'MEG0622', 'MEG0623', 'MEG0621', 'MEG0642', 'MEG0643', 'MEG0641', 'MEG0813', 'MEG0822', 'MEG0823', 'MEG0821', 'MEG0913', 'MEG0911', 'MEG0923', 'MEG0922', 'MEG0921', 'MEG0932', 'MEG0933', 'MEG0931', 'MEG0942', 'MEG0943', 'MEG0941', 'MEG1013', 'MEG1012', 'MEG1011', 'MEG1023', 'MEG1022', 'MEG1021', 'MEG1032', 'MEG1033', 'MEG1031', 'MEG1213', 'MEG1212', 'MEG1232', 'MEG1233', 'MEG1243', 'MEG1242', 'MEG1241'};
chanlist_top = {'MEG0422', 'MEG0421', 'MEG0432', 'MEG0433', 'MEG0431', 'MEG0633', 'MEG0631', 'MEG0713', 'MEG0712', 'MEG0711', 'MEG0723', 'MEG0722', 'MEG0721', 'MEG0733', 'MEG0732', 'MEG0731', 'MEG0743', 'MEG0742', 'MEG0741', 'MEG1043', 'MEG1041', 'MEG1112', 'MEG1111', 'MEG1142', 'MEG1143', 'MEG1141', 'MEG1822', 'MEG1823', 'MEG1821', 'MEG1832', 'MEG1833', 'MEG1831', 'MEG1843', 'MEG1842', 'MEG1841', 'MEG2013', 'MEG2012', 'MEG2011', 'MEG2023', 'MEG2022', 'MEG2021', 'MEG2212', 'MEG2213', 'MEG2211', 'MEG2233', 'MEG2232', 'MEG2231', 'MEG2242', 'MEG2243', 'MEG2241', 'MEG2313'};
chanlist_left = {'MEG0132', 'MEG0131', 'MEG0143', 'MEG0142', 'MEG0141', 'MEG0213', 'MEG0212', 'MEG0211', 'MEG0222', 'MEG0223', 'MEG0221', 'MEG0232', 'MEG0233', 'MEG0231', 'MEG0243', 'MEG0242', 'MEG0241', 'MEG0413', 'MEG0412', 'MEG0443', 'MEG0442', 'MEG1512', 'MEG1513', 'MEG1511', 'MEG1522', 'MEG1523', 'MEG1521', 'MEG1533', 'MEG1532', 'MEG1531', 'MEG1543', 'MEG1542', 'MEG1541', 'MEG1613', 'MEG1612', 'MEG1611', 'MEG1622', 'MEG1623', 'MEG1621', 'MEG1632', 'MEG1633', 'MEG1631', 'MEG1643', 'MEG1642', 'MEG1641', 'MEG1722', 'MEG1723', 'MEG1721', 'MEG1813', 'MEG1812', 'MEG1943', 'MEG1942', 'MEG1941'};
chanlist_right = {'MEG1312', 'MEG1323', 'MEG1321', 'MEG1333', 'MEG1332', 'MEG1331', 'MEG1342', 'MEG1343', 'MEG1341', 'MEG1433', 'MEG1432', 'MEG1431', 'MEG1442', 'MEG2321', 'MEG2412', 'MEG2413', 'MEG2411', 'MEG2423', 'MEG2422', 'MEG2421', 'MEG2433', 'MEG2432', 'MEG2431', 'MEG2441', 'MEG2522', 'MEG2523', 'MEG2521', 'MEG2612', 'MEG2613', 'MEG2611', 'MEG2623', 'MEG2622', 'MEG2621', 'MEG2633', 'MEG2632', 'MEG2631', 'MEG2642', 'MEG2643', 'MEG2641'};
chanlist_back = {'MEG1632', 'MEG1641', 'MEG1732', 'MEG1733', 'MEG1731', 'MEG1743', 'MEG1742', 'MEG1741', 'MEG1912', 'MEG1913', 'MEG1911', 'MEG1923', 'MEG1922', 'MEG1921', 'MEG1932', 'MEG1933', 'MEG1931', 'MEG1943', 'MEG1942', 'MEG1941', 'MEG2013', 'MEG2012', 'MEG2011', 'MEG2023', 'MEG2022', 'MEG2021', 'MEG2032', 'MEG2033', 'MEG2031', 'MEG2042', 'MEG2043', 'MEG2041', 'MEG2113', 'MEG2112', 'MEG2111', 'MEG2122', 'MEG2123', 'MEG2121', 'MEG2133', 'MEG2132', 'MEG2131', 'MEG2143', 'MEG2142', 'MEG2141', 'MEG2312', 'MEG2313', 'MEG2311', 'MEG2323', 'MEG2322', 'MEG2321', 'MEG2332', 'MEG2333', 'MEG2331', 'MEG2343', 'MEG2342', 'MEG2341', 'MEG2433', 'MEG2432', 'MEG2431', 'MEG2442', 'MEG2512', 'MEG2513', 'MEG2511', 'MEG2543', 'MEG2542', 'MEG2541'};
cfg = [];
cfg.avgoverchan = 'yes';
cfg.channel = chanlist_front;
data_front = ft_selectdata(cfg,data_timelocked);
cfg.channel = chanlist_top;
data_top = ft_selectdata(cfg,data_timelocked);
cfg.channel = chanlist_left;
data_left = ft_selectdata(cfg,data_timelocked);
cfg.channel = chanlist_right;
data_right = ft_selectdata(cfg,data_timelocked);
cfg.channel = chanlist_back;
data_back = ft_selectdata(cfg,data_timelocked);
data_virtual_channels = ft_appenddata([],data_front, data_top, data_left, data_right, data_back);
data_virtual_channels.label = {'front','top','left','right','back'};
clear data_left data_right data_top data_back data_front data_timelocked
carrierfreq{1} = 1000;
carrierfreq{2} = 2000;
carrierfreq{3} = 3000;
sourcebp{1} = [1 7];
sourcebp{2} = [10 11];
sourcebp{3} = [40 60];
freqname = {'theta', 'alpha', 'gamma'};
wavFs = 10000; % samplerate for WAV file
for ifreq = 1 : 3
disp('Bandpass filtering');
cfg = [];
if sourcebp{ifreq}(1) == 0
cfg.lpfilter = 'yes';
cfg.lpfreq = sourcebp{ifreq}(2);
else
cfg.bpfilter = 'yes';
cfg.bpfreq = sourcebp{ifreq};
cfg.bpfiltord = 3;
end
cfg.demean = 'yes';
data_bp{ifreq} = ft_preprocessing(cfg,data_virtual_channels);
figure; plot(data_bp{ifreq}.trial{1}(1,1:10000));
disp('Hilbert transform');
cfg = [];
cfg.hilbert = 'abs';
data_hilbert{ifreq} = ft_preprocessing(cfg,data_bp{ifreq});
disp('Upsampling');
data_hilbert_upsampled{ifreq} = removefields(data_hilbert{ifreq},'trial');
for ichan = 5 : 5
data_hilbert_upsampled{ifreq}.trial{1}(ichan,:) = interp(data_hilbert{ifreq}.trial{1}(ichan,:),wavFs / data_hilbert{ifreq}.fsample,1);
end
data_AM{ifreq} = data_hilbert_upsampled{ifreq};
for ichan = 5 : 5
disp('Multiplying with carrier wave');
isin = 0;
for x = 1 : size(data_hilbert_upsampled{ifreq}.trial{1},2)
isin = isin + 2*pi/wavFs*carrierfreq{ifreq};
carrierwave{ifreq}(ichan,x) = sin(isin);
data_AM{ifreq}.avg(ichan,x) = sin(isin) .* data_hilbert_upsampled{ifreq}.trial{1}(ichan,x);
end;
data_AM{ifreq}.avg(ichan,:) = data_AM{ifreq}.avg(ichan,:) * (1 / max(data_AM{ifreq}.avg(ichan,:))) * 0.7; % change volume
wavwrite(data_AM{ifreq}.avg(ichan,:),wavFs, [ freqname{ifreq} '_' data_AM{ifreq}.label{ichan} '_' int2str(carrierfreq{ifreq}) 'Hz']);
end;
end;
The sound is horrible to listen to though. It works better to use more pleasant frequencies. This was not the initial purpose of the exercise, but Robert was quick to play along and write some MATLAB code to transpose the WAV files (He did not have access to the initial data, and it takes quite a while to run the initial script). Next time we’ll do it like that from the start.
%% preprocessing and visualisation
cfg = [];
cfg.dataset = pwd;
cfg.dataformat = 'combined_ds';
cfg.headerformat = 'combined_ds';
data = ft_preprocessing(cfg);
cfg = [];
cfg.viewmode = 'vertical';
ft_databrowser(cfg, data);
%% transposition
% get the frequency for all notes
scientific_pitch_notation
% alpha at 2000 Hz
for chan=1:5
disp(chan)
% this takes about 10 seconds per channel for 15430000 samples
s = data.trial{1}(chan,:);
t = data.time{1};
h = hilbert(s);
p = exp((C(4)-2000)*2*pi*t*1i); % phase rotation, linear over time
s = real(h.*p); % reconstruct amplitide from phase rotated signal
data.trial{1}(chan,:) = s;
end
% gamma at 4000 Hz
for chan=6:9
disp(chan)
% this takes about 10 seconds per channel for 15430000 samples
s = data.trial{1}(chan,:);
t = data.time{1};
h = hilbert(s);
p = exp((E(4)-3000)*2*pi*t*1i); % phase rotation, linear over time
s = real(h.*p); % reconstruct amplitide from phase rotated signal
data.trial{1}(chan,:) = s;
end
% theta at 1000 Hz
for chan=10:14
disp(chan)
% this takes about 10 seconds per channel for 15430000 samples
s = data.trial{1}(chan,:);
t = data.time{1};
h = hilbert(s);
p = exp((G(3)-1000)*2*pi*t*1i); % phase rotation, linear over time
s = real(h.*p); % reconstruct amplitide from phase rotated signal
data.trial{1}(chan,:) = s;
end
%% quick and dirty analysis, play one second of sound
psd(spectrum.welch, s(1,1:100000), 'Fs', data.fsample)
soundsc(s(1,1:100000), data.fsample)
%% export each channel to a separate file to disk
hdr = [];
hdr.Fs = data.fsample;
hdr.nChans = 1;
ft_write_data('transposed/alpha_back_C4', data.trial{1}(1,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/alpha_front_C4', data.trial{1}(2,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/alpha_left_C4', data.trial{1}(3,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/alpha_right_C4', data.trial{1}(4,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/alpha_top_C4', data.trial{1}(5,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/gamma_front_E4', data.trial{1}(6,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/gamma_left_E4', data.trial{1}(7,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/gamma_right_E4', data.trial{1}(8,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/gamma_top_E4', data.trial{1}(9,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/theta_back_G3', data.trial{1}(10,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/theta_front_G3', data.trial{1}(11,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/theta_left_G3', data.trial{1}(12,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/theta_right_G3', data.trial{1}(13,:), 'header', hdr, 'dataformat', 'riff_wave')
ft_write_data('transposed/theta_top_G3', data.trial{1}(14,:), 'header', hdr, 'dataformat', 'riff_wave')
Robert also looked up this nice overview table:
% I removed the first column to make the indexing consistent with the notation % http://en.wikipedia.org/wiki/Scientific_pitch_notation C = [32.703 65.406 130.81 261.63 523.25 1046.5 2093.0 4186.0 8372.0 16744.0]; Db = [34.648 69.296 138.59 277.18 554.37 1108.7 2217.5 4434.9 8869.8 17739.7]; D = [36.708 73.416 146.83 293.66 587.33 1174.7 2349.3 4698.6 9397.3 18794.5]; Eb = [38.891 77.782 155.56 311.13 622.25 1244.5 2489.0 4978.0 9956.1 19912.1]; E = [41.203 82.407 164.81 329.63 659.26 1318.5 2637.0 5274.0 10548.1 21096.2]; F = [43.654 87.307 174.61 349.23 698.46 1396.9 2793.8 5587.7 11175.3 22350.6]; Gb = [46.249 92.499 185.00 369.99 739.99 1480.0 2960.0 5919.9 11839.8 23679.6]; G = [48.999 97.999 196.00 392.00 783.99 1568.0 3136.0 6271.9 12543.9 25087.7]; Ab = [51.913 103.83 207.65 415.30 830.61 1661.2 3322.4 6644.9 13289.8 26579.5]; A = [55.000 110.00 220.00 440.00 880.00 1760.0 3520.0 7040.0 14080.0 28160.0]; Bb = [58.270 116.54 233.08 466.16 932.33 1864.7 3729.3 7458.6 14917.2 29834.5]; B = [61.735 123.47 246.94 493.88 987.77 1975.5 3951.1 7902.1 15804.3 31608.5];