From cb30899b1a80ce45c0545b04b7394f9ad0d25e7f Mon Sep 17 00:00:00 2001 From: Robert Seymour Date: Tue, 23 Feb 2021 14:28:28 +0000 Subject: [PATCH] VERY HACKY - needs work #36 --- data_functions/ft_opm_create_chunks.m | 256 ++++++++++++++++++++++++++ preprocessing/plot_powspctrm.m | 6 +- trialfun/OPM_trialfun_usemat.m | 64 +++++++ 3 files changed, 323 insertions(+), 3 deletions(-) create mode 100644 data_functions/ft_opm_create_chunks.m create mode 100644 trialfun/OPM_trialfun_usemat.m diff --git a/data_functions/ft_opm_create_chunks.m b/data_functions/ft_opm_create_chunks.m new file mode 100644 index 0000000..d7643d0 --- /dev/null +++ b/data_functions/ft_opm_create_chunks.m @@ -0,0 +1,256 @@ +function [rawData] = ft_opm_create_chunks(cfg) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% WARNING: USE THIS AT YOUR OWN RISK - NOT FULLY TESTED +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Function to read optically-pumped magnetencephalography (OPMEG) data +% acquired from the UCL Wellcome Centre for Neuroimaging in chunks of Ns +% and save to disk +% +% Hopefully this will be useful for very long OPM datasets. +% +% EXAMPLE USEAGE: data = ft_opm_create_chunks(cfg) +% ...where, cfg is the input structure +% +%%%%%%%%%%% +% Inputs: +%%%%%%%%%%% +% cfg.chunksize = number in seconds for the length of the chunk +% (default = 300s) +% cfg.downsample = 'yes' or 'no' to 1000Hz (default = 'no') +% +%%%%%%%%%%% +% Option 1 +%%%%%%%%%%% +% cfg.data = path to raw .bin file +% (requires .json and _channels.tsv in same folder +% with same name as the .bin file) +%%%%%%%%%%% +% Option 2 +%%%%%%%%%%% +% cfg.folder = path to folder containing data +% organised in BIDS format +% cfg.bids.task = 'AEF' +% cfg.bids.sub = '001' +% cfg.bids.ses = '003' +% cfg.bids.run = '002' +%__________________________________________________________________________ +% Copyright (C) 2020 Wellcome Trust Centre for Neuroimaging +% Adapted from spm_opm_create (Tim Tierney) + +% Authors: Nicholas Alexander (n.alexander@ucl.ac.uk) +% Robert Seymour (rob.seymour@ucl.ac.uk) +%__________________________________________________________________________ + +%% Set default values +if ~isfield(cfg, 'fs') + cfg.fs = 6000; +end +if ~isfield(cfg, 'nSamples') + cfg.nSamples = 1000; +end + +if ~isfield(cfg, 'precision') + cfg.precision = 'double'; +end + +if ~isfield(cfg, 'chunksize') + cfg.chunksize = 300; +end + +if ~isfield(cfg, 'downsample') + cfg.downsample = 'no'; +end + +%% Determine whether to use cfg.data or cfg.folder +if ~isfield(cfg, 'data') + use_bids = 1; +else + path_to_bin_file = cfg.data; + use_bids = 0; +end + +% If the user has specified both cfg.data and cfg.folder default to +% cfg.data, but warn the user +if isfield(cfg, 'data') && isfield(cfg, 'folder') + ft_warning(['Both cfg.data and cfg.folder have been supplied.'... + ' Defaulting to cfg.data']); +end + +%% If using option 2 = BIDS! +if use_bids + try + file_name_bids = ['sub-' cfg.bids.sub '_ses-' cfg.bids.ses ... + '_task-' cfg.bids.task '_run-' cfg.bids.run '_meg.bin']; + catch + error('Did you specify all the required cfg.bids information?') + end + path_to_bin_file = fullfile(cfg.folder,['sub-' cfg.bids.sub],... + ['ses-' cfg.bids.ses],'meg',file_name_bids); +end + + +%% Read Binary File +try % to read data + [direc, dataFile] = fileparts(path_to_bin_file); + dat = fopen(path_to_bin_file); + % data_raw = fread(dat,Inf,cfg.precision,0,'b'); + %fclose(dat); + binData = 1; +catch % if not readable check if it is numeric + if ~isa(path_to_bin_file,'numeric') % if not numeric throw error + error(['Cannot read the file: ' path_to_bin_file]) + end + binData = 0; + direc = pwd(); + dataFile = cfg.fname; +end + +%% +% Identify potential BIDS Files +base = strsplit(dataFile,'meg'); +chanFile = fullfile(direc,[base{1},'channels.tsv']); +megFile = fullfile(direc,[base{1},'meg.json']); +posFile = fullfile(direc,[base{1},'positions.tsv']); +coordFile = fullfile(direc,[base{1},'coordsystem.json']); + +%% Load a channels file. +try + channels = ft_read_chan_tsv(chanFile); +catch + error(['Cannot read the file: ' chanFile]); +end + +% Reformat data according to channel info. +numChan = size(channels.name,1); + +%% Get header information from json files +% Check for MEG Info +try + meg = ft_read_json(megFile); +catch + error(['ERROR. Attempted to read: ' megFile]); +end + +try + coordsystem = ft_read_json(coordFile); +catch + ft_warning('Not loaded any coordinate system'); + coordsystem = []; +end + +% Position File check +try % to load a channels file + posOri = ft_read_pos_tsv(posFile); + positions = 1; +catch + try % to load a BIDS channel file + posOri = ft_read_pos_tsv(posFile); + positions = 1; + catch + positions = 0; + end +end + +%% +time_for_loop = 0; +sampleinfo_for_loop = 0; +count = 1; + +%% Chunk it up +% compute the number of samples and chunks + +% Read file in chunks to save memory +Chunk = numChan*meg.SamplingFrequency*cfg.chunksize; +Count = Chunk; % Dummy value to satisfy WHILE condition + +while Count == Chunk + [data_raw, Count] = fread(dat,Chunk,cfg.precision,0,'b'); + + + %% Load in channels and reorganise the raw data + + if binData + data = reshape(data_raw,numChan,numel(data_raw)/numChan); + elseif numChan ~= size(data_raw,1) + error('number of channels in cfg.data different to cfg.channels') + else + data = data_raw; + end + + clear data_raw + + %% Make the main fieldtrip structure + % Main structure + rawData = []; + rawData.label = channels.name; + rawData.fsample = meg.SamplingFrequency; + rawData.time{1} = ((0:1:(size(data,2)-1))./rawData.fsample)+time_for_loop; + rawData.trial{1} = data; + rawData.sampleinfo = [sampleinfo_for_loop sampleinfo_for_loop+length(data(1,:))]; + + % Save for next loop + time_for_loop = rawData.time{1}(end)+(1/meg.SamplingFrequency); + sampleinfo_for_loop = rawData.sampleinfo(2)+1; + + % Construct a cfg file + %%%%%%%%% + % TO DO + %%%%%%%%%% + %rawData.cfg = cfg; + + % Now do the same for the header information + rawData.hdr = []; + rawData.hdr.label = channels.name; + rawData.hdr.chantype = channels.type; + rawData.hdr.chanunit = channels.units; + rawData.hdr.dimord = 'chan_time'; + rawData.hdr.Fs = meg.SamplingFrequency; + rawData.hdr.nSamples = length(data(1,:)); + rawData.hdr.nTrials = 1; + rawData.hdr.nSamplesPre = 0; + rawData.hdr.fieldori = channels.fieldori; + + if positions + rawData.grad.chanpos = posOri.chanpos; + rawData.grad.chanori = posOri.chanori; + rawData.grad.label = posOri.label; + rawData.grad.coilpos = posOri.chanpos; + rawData.grad.coilori = posOri.chanori; + + % Find the position of the channels in overall channel list + find_chan = contains(channels.name,posOri.label); + rawData.grad.chanunit = channels.units(find_chan); + rawData.grad.tra = diag(ones(1,length(rawData.grad.label))); + + end + +% % Resample +% sampleinfo = rawData.sampleinfo; + + if strcmp(cfg.downsample,'yes') + %% Resample to 1000Hz + disp(['Resampling Chunk... ' num2str(count)]); + cfg2 = []; + cfg2.resamplefs = 1000; + [rawData] = ft_resampledata(cfg2, rawData); + %rawData.sampleinfo = sampleinfo./10; + end + + %% Save + disp(['Saving Chunk... ' num2str(count)]); + save(['rawData_' num2str(count)],'rawData'); + count = count+1; + +end + +end + + + + + + + diff --git a/preprocessing/plot_powspctrm.m b/preprocessing/plot_powspctrm.m index c13d85b..e72b69f 100644 --- a/preprocessing/plot_powspctrm.m +++ b/preprocessing/plot_powspctrm.m @@ -34,7 +34,7 @@ function plot_powspctrm(po,cfg,label,freq) hold on; % Plot the mean in black -plot(freq,mean_po,'-k','LineWidth',3); +plot(freq,mean_po,'-k','LineWidth',2); hold on xp2 =0:round(freq(end)); yp2=ones(1,round(freq(end))+1)*15; @@ -59,9 +59,9 @@ function plot_powspctrm(po,cfg,label,freq) set(gca, 'YScale', 'log'); hold on; -xlabel('Frequency (Hz)','FontSize',25) +xlabel('Frequency (Hz)','FontSize',20) labY = ['$$PSD (' 'fT' ' \sqrt[-1]{Hz}$$)']; -ylabel(labY,'interpreter','latex','FontSize',25) +ylabel(labY,'interpreter','latex','FontSize',20) % Adjust limits based on cfg.foi xlim([cfg.foi(1), cfg.foi(end)]); diff --git a/trialfun/OPM_trialfun_usemat.m b/trialfun/OPM_trialfun_usemat.m new file mode 100644 index 0000000..f1563ed --- /dev/null +++ b/trialfun/OPM_trialfun_usemat.m @@ -0,0 +1,64 @@ +function trl = OPM_trialfun_usemat(cfg) + + +pos_of_trig = contains(cfg.rawData.hdr.label,cfg.trialdef.trigchan); + +% % Read header +% hdr = ft_read_header(strcat(cfg.dataset(1:end-3),'mat')); +% +% % Read continuous event channel. +% cfgEvent = []; +% cfgEvent.continuous = 'yes'; +% cfgEvent.datafile = cfg.dataset; +% cfgEvent.headerfile = strcat(strcat(cfg.dataset(1:end-3),'mat')); +% cfgEvent.channel = cfg.trialdef.trigchan; +% cfgEvent.dataformat = 'spmeeg_mat'; +% eventContStruct = ft_preprocessing(cfgEvent); + +% % Downsample from 6000Hz +% cfgDS = []; +% cfgDS.resamplefs = cfg.trialdef.downsample; +% cfgDS.detrend = 'no'; +% eventContStruct = ft_resampledata(cfgDS,eventContStruct); + +eventCont = cfg.rawData.trial{1}(pos_of_trig,:); + +time = cfg.rawData.time{1}; +sample = 1:length(time); +% convert to binary +eventDisc = (eventCont > 3); + +% find first sample. first find differences +tmp = eventDisc - [0,eventDisc(1:end-1)]; +tmp2 = (tmp == 1); +events = sample(tmp2)'; +events = round(events); + +if ~isempty(cfg.correct_time) + disp(['Correcting timing by ' num2str(cfg.correct_time) 's']); + time_to_correct = cfg.rawData.fsample*cfg.correct_time; + events = events+time_to_correct; +end + + +if isfield(cfg.trialdef,'downsample') + events = round((events/(cfg.rawData.fsample/cfg.trialdef.downsample))); +end + +% create trl structure based upon the events +trl = []; + +for j = 1:length(events) + if isfield(cfg.trialdef,'downsample') + trlbegin = (events(j) - (cfg.trialdef.prestim*cfg.trialdef.downsample)); + trlend = (events(j) + (cfg.trialdef.poststim*cfg.trialdef.downsample)); + offset = (-cfg.trialdef.prestim*cfg.trialdef.downsample); + trl(end+1, :) = ([trlbegin trlend offset j]); + else + trlbegin = events(j) - cfg.trialdef.prestim*cfg.rawData.hdr.Fs; + trlend = events(j) + cfg.trialdef.poststim*cfg.rawData.hdr.Fs; + offset = -cfg.trialdef.prestim*cfg.rawData.hdr.Fs; + trl(end+1, :) = ([trlbegin trlend offset j]); + end +end +end \ No newline at end of file