function hpiloc_largeHpiArray( array_config, meg_array, hpi_moment, noise, time, f_smpl , path ) %hpiloc_varyNcoil Summary of this function goes here % Detailed explanation goes here %% -- Definitions -- load([path, 'grads.mat']); load([path, 'electrode_positions.mat']); load([path, 'headSurface.mat']); hpigrad.tra = diag(ones(1,10)); if strcmp(meg_array,'elekta') meggrad = meggrad_elekta; elseif strcmp(meg_array,'onscalp') meggrad = meggrad_onscalp; else disp('ERROR: wrong MEG array') meggrad = []; end headmodel = []; headmodel.unit = 'm'; headmodel.type = 'infinite'; %hpicoilfreq = [ 370 420 470 520 570 620 670 720 770 820 ]; %hpicoilfreq = [ 218 225 232 239 246 253 260 267 274 281 ]; hpi_freq = 218+(1:32)*7; locError_mean = zeros(length(array_config),1); locError_median = zeros(length(array_config),1); %% -- Calculations -- meg_error = zeros(102,1); n_HPI = zeros(length(array_config),1); for m = 1:length(array_config) n_HPI(m) = length(array_config{m}); %% HPI array hpi_array_pos = zeros(n_HPI(m),3); hpi_array_ori = zeros(n_HPI(m),3); hpigrad_array = []; hpigrad_array.unit = 'm'; hpigrad_array.tra = zeros(n_HPI(m)); for i = 1:n_HPI(m) i_tmp = find(strncmp(elec.label,array_config{m}{i},5)); tmp = elec.elecpos(i_tmp,:); [hpi_array_pos(i,:), hpi_array_ori(i,:)] = snap2surf(tmp, bnd.pos, bnd.tri, bnd.nn); hpigrad_array.label{i} = meggrad.label{i}; hpigrad_array.chanpos(i,:) = hpi_array_pos(i,:); hpigrad_array.coilpos(i,:) = hpi_array_pos(i,:); hpigrad_array.coilori(i,:) = hpi_array_ori(i,:); hpigrad_array.tra(i,i) = hpigrad.tra(1,1); end %% - Generate data - amp = hpi_moment; t = (1:time)/f_smpl; data = []; data.time{1} = t; data.trial{1} = zeros(102,time); data.label = meggrad.label; for i = 1:n_HPI(m) signal = amp*cos(2*pi*hpi_freq(i)*t); lf = ft_compute_leadfield(hpigrad_array.coilpos(i,:),meggrad,headmodel); lf = lf * hpigrad_array.coilori(i,:)'; data.trial{1} = data.trial{1}+lf*signal; end data.trial{1} = data.trial{1}+noise.*randn(102,time); % Add noise %% - FFT - cfg = []; cfg.method = 'mtmfft'; cfg.output = 'fourier'; cfg.taper = 'boxcar'; cfg.foi = hpi_freq(1:n_HPI(m))'; freq = ft_freqanalysis(cfg, data); %% - Freq reverse - freqrev = []; freqrev.label = hpigrad_array.label; freqrev.dimord = freq.dimord; freqrev.fourierspctrm = nan(1, n_HPI(m), 102); freqrev.freq = (1:102)'; % all frequencies combined into one reverse HPI topography per channel freqrev.fourierspctrm = permute(freq.fourierspctrm(:,:,:), [1 3 2]); %% - Fit mags - tmp_coilpos = meggrad.coilpos; parfor i=1:102 cfg = []; cfg.channel = 'all'; cfg.headmodel = headmodel; cfg.gridsearch = 'no'; cfg.dip.pos = tmp_coilpos(i,:); cfg.dip.unit = 'm'; cfg.model = 'regional'; cfg.frequency = i; cfg.grad = hpigrad_array; megfit_individual = ft_dipolefitting(cfg, freqrev); meg_error(i)=norm(tmp_coilpos(i,:)-megfit_individual.dip.pos); clc; m end locError_mean(m) = mean(meg_error); locError_max(m) = max(meg_error); locError_min(m) = min(meg_error); locError_std(m) = std(meg_error); locError_median(m) = median(meg_error); end %% - Temporarily save results - res = []; res.array_config = array_config; res.n_hpi = n_HPI; res.c_ring = meggrad.coilpos; res.noise = noise; res.meg_array = meg_array; res.hpi_moment = hpi_moment; res.locError_mean = locError_mean; res.locError_max = locError_max; res.locError_min = locError_min; res.locError_std = locError_std; res.locError_median = locError_median; %% -- Save results -- save([path 'largeHpiArray_' meg_array '.mat'],'res'); %% -- Plot -- f = figure(); hold on plot(n_HPI,locError_mean.*1e3,'*-') hold off xlabel('Number of HPI coils') ylabel('Localization error [mm]') title('Localization error vs. number of HPI coils in full-head array') saveas(f, [path 'Loc_largeHpiArray/largeHpiArray_' meg_array '_mean'], 'fig'); f = figure(); hold on plot(n_HPI,locError_median.*1e3,'*-') hold off xlabel('Number of HPI coils') ylabel('Localization error [mm]') title('Localization error vs. number of HPI coils in full-head array') saveas(f, [path 'Loc_largeHpiArray/largeHpiArray_' meg_array '_median'], 'fig'); end