Authored by JD

adapt StimDetect method (StimDetectVersion=4)

function [ alignments, alignedStimulations ] = ImaGIN_AlignSegmentsIter( d, stimulations, segmentHWidth, templateHWidth, corrHWidth, do_plot, iterMax )
if size(d,1) ~= 1
error( 'Invalid d' )
end
alignments = nan( length(stimulations), iterMax ) ;
alignedStimulations = stimulations ;
segments = nan( length(alignedStimulations), 2*segmentHWidth+1 ) ;
for i=1:iterMax
fprintf( 1, [ 'alignSegments iter ' num2str(i) '...' '\n' ] ) ;
for s=1:length(alignedStimulations)
segments(s,:) = d( alignedStimulations(s) + (-segmentHWidth:segmentHWidth) ) ;
end
[ alignment, ~ ] = alignSegments( segments, templateHWidth, corrHWidth, do_plot ) ;
alignments(:,i) = alignment ;
if isempty(find(alignment,1))
% fprintf( 1, [ 'alignSegments iter ' num2str(i) ' => BREAK' '\n' ] ) ;
break
else
alignedStimulations = alignedStimulations + alignment' ;
end
% if i == iterMax
% error( [ '>> WARNING: alignment not achieved after iterMax = ' num2str(iterMax) ] ) ;
% % fprintf( 1, [ '>> WARNING: alignment not achieved after iterMax = ' num2str(iterMax) '\n' ] ) ;
% end
end
end
function [ alignment, alignedSegments ] = alignSegments( segments, templateHWidth, corrHWidth, do_plot )
if ~ ( corrHWidth > templateHWidth )
error( 'Invalid corrHWidth' )
end
segmentCount = size(segments,1) ;
segmentWidth = size(segments,2) ;
segmentCenter = round(segmentWidth/2);
templateWidth = 2*templateHWidth+1 ;
corrWidth = 2*corrHWidth+1 ;
corrCount = corrWidth - templateWidth + 1 ;
corr = nan( segmentCount, corrCount ) ;
template = mean( segments(:,segmentCenter+(-templateHWidth:templateHWidth)), 1 ) ;
for s=1:segmentCount
for c=1:corrCount
segStart = segmentCenter - corrHWidth -1 + c ;
segEnd = segStart + templateWidth - 1 ;
r = corrcoef( template, segments(s,segStart:segEnd) ) ;
corr(s,c) = r(2,1) ;
end
end
[ ~, i ] = max( corr,[],2) ;
alignment = segmentCenter - corrHWidth -1 + i + templateHWidth - segmentCenter ;
alignedSegments = applyAlignment( segments, alignment ) ;
if do_plot
plotSegmentsAlignment( segments, alignedSegments, segmentCenter ) ;
end
end
function alignedSegments = applyAlignment( segments, alignment )
alignedSegments = nan( size(segments) ) ;
for s=1:size(segments,1)
if alignment(s) > 0
alignedSegments(s,1:end-alignment(s)) = segments(s,1+alignment(s):end) ;
elseif alignment(s) < 0
alignedSegments(s,1+abs(alignment(s)):end) = segments(s,1:end-abs(alignment(s))) ;
else
alignedSegments(s,:) = segments(s,:);
end
end
end
function plotSegmentsAlignment( segments, alignedSegments, segmentCenter )
figure ;
ax1 = subplot(2,1,1); hold on ;
plot( 1:size(segments,2), segments ) ;
plot( 1:size(segments,2), mean(segments,1), 'color', 'r', 'linewidth', 4 ) ;
line( [ segmentCenter; segmentCenter ], get( gca, 'YLim' )', 'color', 'm', 'linewidth', 2 ) ;
xlim( [ 1 size(segments,2) ] ) ;
title( 'templates' )
ax2 = subplot(2,1,2); hold on ;
plot( 1:size(segments,2), alignedSegments ) ;
plot( 1:size(segments,2), mean(alignedSegments,1), 'color', 'r', 'linewidth', 4 ) ;
line( [ segmentCenter; segmentCenter ], get( gca, 'YLim' )', 'color', 'm', 'linewidth', 2 ) ;
xlim( [ 1 size(segments,2) ] ) ;
title( 'templates aligned' )
linkaxes( [ ax1, ax2 ], 'xy' )
end
... ...
... ... @@ -148,7 +148,9 @@ function ImaGIN_CropManual(D,t,S)
Dnew(:, :, 1) = d;
Dnew = events(Dnew, 1, select_events(Events,[Time(1)/D.fsample+timeonset(D) Time(2)/D.fsample+timeonset(D)]));
Dnew = timeonset(Dnew, Time(1)./D.fsample+D.timeonset);
% Dnew = timeonset(Dnew, Time(1)./D.fsample+D.timeonset);
% first sample has relative time 0
Dnew = timeonset(Dnew, (Time(1)-1)./D.fsample+D.timeonset);
save(Dnew);
end
... ...
... ... @@ -36,8 +36,17 @@ else
thisN = '';
end
% method to be used in ImaGIN_StimDetect
try
StimDetectVersion = S.StimDetectVersion ;
catch
StimDetectVersion = 4 ;
end
D = spm_eeg_load(sFile); % Load the converted file .mat
%% Find existing crop files in DirOut
matExist = dir(fullfile(DirOut,'*.mat')); %Delete all cropped text file
matExist = {matExist.name};
... ... @@ -77,6 +86,11 @@ for c=1:evsize % Navigate all available events
xpr4b = '\w*55.0hz\w*';
xpr5b = '\w*55hz\w*';
xpr6b = '\w*55 hz\w*';
% present in '/gin/data/database/02-raw/0007BUC/SEEG/StimLF_2013-11-13/Electrodes/electrodes_0007BUC_Raw_2.mat' ; 'P10-P11_2mA_1011Hz_3000us'
% while in fact, it is 1Hz
% and using 1011Hz makes a poor detection
xpr6c = '\w*1011Hz\w*';
xpr7 = '\w*alarme\w*';
xpr8 = '\w*SE1Hz\w*';
... ... @@ -95,6 +109,7 @@ for c=1:evsize % Navigate all available events
~isempty(regexpi(Notes{c},xpr12)) || ~isempty(regexpi(Notes{c},xpr4b))|| ...
~isempty(regexpi(Notes{c},xpr12b))|| ~isempty(regexpi(Notes{c},xpr13))|| ...
~isempty(regexpi(Notes{c},xpr5b)) || ~isempty(regexpi(Notes{c},xpr6b))|| ...
~isempty(regexpi(Notes{c},xpr6c)) || ...
~isempty(regexpi(Notes{c},xpr11)) || strcmpi(Notes{c}(1:min([length(Notes{c}) 5])),xpr10)
elseif ~isempty(regexpi(Notes{c},xpr1))
KeepEvent=[KeepEvent c];
... ... @@ -461,8 +476,11 @@ for c = 1:length(KeepEvent) % Navigate all stim events
% S.FindBadChannels=0;
%}
[stimTime,~,~] = ImaGIN_StimDetect(S);
disp(KeepEvent(c)), disp(S.EvtName);
S.StimDetectVersion = StimDetectVersion ;
[stimTime,~,~] = ImaGIN_StimDetect(S);
stimFq = StimFreq;
% if numel(stimTime) >= seuilHz %OD
... ...
... ... @@ -59,6 +59,12 @@ catch
FileOut = [Filename(1:end-3) 'txt'];
end
try
StimDetectVersion=S.StimDetectVersion;
catch
StimDetectVersion=4;
end
% D=timeonset(D,0);
% save(D)
... ... @@ -221,14 +227,226 @@ if 1==1
end
% remove outliers close to other stims
stimulation = removeOutliers( stimulation, Stim ) ;
if ~isempty(stimulation)
%try to correct onset inpired from Lena's method
% alignment of templates
segmentHWidth = ceil(Stim/2) ;
templateHWidth = round( 0.01 * fsample(D) );
corrHWidth = round( 0.02 * fsample(D) ) ;
do_plot = false ;
iterMax = 10 ;
if segmentHWidth > templateHWidth % sometimes it happens to be not the case: 0007BUC/ElectrodeFile#2/S.EvtName = P10-P11_2mA_1011Hz_3000us
[ alignments, alignedStimulations ] = ImaGIN_AlignSegmentsIter( d, stimulation, segmentHWidth, templateHWidth, corrHWidth, do_plot, iterMax ) ;
if 0
[ pathstr, name ] = fileparts(FileOut) ;
fname = fullfile( pathstr, 'templateAlignment', [ name '__templateAlignment.mat' ] ) ;
if ~exist( fileparts(fname), 'dir' ), mkdir( fileparts(fname) ) ; end
fprintf( 1, [ 'Save ' fname '\n' ] ) ;
save( fname, 'stimulation', 'alignments', 'alignedStimulations' ) ;
end
% apply alignment only if convergence and max(abs(alignments)) < 10ms
for i=1:iterMax
% convergence
if isempty(find(alignments(:,i)~=0,1))
% sanity check
if ~isequal( alignedStimulations, stimulation+sum(alignments(:,1:i),2)' ), error( 'Inconsistent alignment' ) ; end
% apply alignment if < 10ms
if max(abs(alignedStimulations-stimulation))/fsample(D) < 10e-3
stimulation = alignedStimulations ;
end
break
end
end
end
% try to correct onset inpired from Lena's method
tmps=sign(mean(Data(:,stimulation),2));
Data2=Data;
for i1=1:size(Data2,1)
Data2(i1,:)=tmps(i1)*Data2(i1,:);
end
if 1==1
TemplateSamples=-segmentHWidth:segmentHWidth ;
TimeTemplate=TemplateSamples./fsample(D);
avg_to_plot = 0 ;
for i1=1:length(stimulation)
avg_to_plot = avg_to_plot + Data(:,stimulation(i1)+TemplateSamples) ;
end
%Same offset for all pulses
% StimDetectVersion = 1 ; % based on template of data + threshold from 10 to 5
% StimDetectVersion = 2 ; % based on template of acceleration d + threshold (perc of the max)
% StimDetectVersion = 3 ; % mixed method: accel used to detect the start of the stim and then data to detect the max (threshold from 10 to 4)
% StimDetectVersion = 4 ; % similar to 3 but offset is computed from data based on a percentage of the max
if StimDetectVersion == 1
TemplateData=0;
for i1=1:length(stimulation)
TemplateData=TemplateData+mean(abs(Data2(:,stimulation(i1)+TemplateSamples)),1);
end
TemplateNorm=ImaGIN_normalisation(TemplateData,2,find(TimeTemplate<-0.02&TimeTemplate>-0.1));
template_to_plot = TemplateNorm ;
for th=10:-1:5
tmpoffset=find( abs(TemplateNorm)>th & abs(TimeTemplate)<0.02);
if ~isempty(tmpoffset)
break
end
end
if isempty(tmpoffset)
offset=nan;
else
i1=min(tmpoffset);
offset=ceil(Stim/2)+1-i1;
end
elseif StimDetectVersion == 2
% use mean template from d instead of data
template = 0 ;
for s=1:length(stimulation)
template = template + d( stimulation(s) + (-segmentHWidth:segmentHWidth) ) ;
end
template_to_plot = template ;
% use a percentage of the max as a threshold
% sometimes there are many peaks (biphasic, high fsample) and the first is not the biggest one
% so we prefer a low threshold on the max
th = 0.3 ;
offsetMaxTime = 0.03 ;
offsetMaxSample = round( offsetMaxTime * fsample(D) ) ;
segmentCenter = segmentHWidth + 1 ;
[ templateMax, iMax ] = max( abs(template(segmentCenter+(-offsetMaxSample:offsetMaxSample))) ) ;
tmpoffset = find( template>th*templateMax & abs(TimeTemplate)<offsetMaxTime, 1 );
if isempty(tmpoffset)
offset=nan;
else
i1=min(tmpoffset);
offset=segmentCenter-i1;
end
elseif StimDetectVersion == 3
% offset is computed from the data
TemplateData=0;
for i1=1:length(stimulation)
TemplateData=TemplateData+mean(abs(Data2(:,stimulation(i1)+TemplateSamples)),1);
end
TemplateNorm=ImaGIN_normalisation(TemplateData,2,find(TimeTemplate<-0.02&TimeTemplate>-0.1));
% try to detect stimulation start to refine detection
% use mean template from d instead of data
template = 0 ;
for s=1:length(stimulation)
template = template + d( stimulation(s) + TemplateSamples ) ;
end
% use a percentage of the max as a threshold
% sometimes there are many peaks (biphasic, high fsample) and the first is not the biggest one
% so we prefer a low threshold on the max
th = 30/100 ;
offsetMaxInterval = [ -0.04 0.02 ];
offsetMaxIntervalSample = round( offsetMaxInterval * fsample(D) ) ;
segmentCenter = segmentHWidth + 1 ;
[ templateMax, ~ ] = max( abs(template(segmentCenter+(offsetMaxIntervalSample(1):offsetMaxIntervalSample(2)))) ) ;
startoffset = find( template>th*templateMax & TimeTemplate>offsetMaxInterval(1) & TimeTemplate<offsetMaxInterval(2), 1 );
if isempty(startoffset)
for th=10:-1:4
tmpoffset=find( abs(TemplateNorm)>th & abs(TimeTemplate)<0.02);
if ~isempty(tmpoffset)
break
end
end
else
startoffset=min(startoffset);
for th=10:-1:4
tmpoffset=find( abs(TemplateNorm)>th & TimeTemplate>=TimeTemplate(startoffset+1) & TimeTemplate<=TimeTemplate(startoffset)+8e-3 );
if ~isempty(tmpoffset)
break
end
end
if isempty(tmpoffset)
% 4ms after startoffset (one sample at 256Hz)
tmpoffset = startoffset + round( fsample(D)*4e-3 ) ;
end
end
if isempty(tmpoffset)
offset=nan;
else
i1=min(tmpoffset);
offset=ceil(Stim/2)+1-i1;
end
elseif StimDetectVersion == 4
% offset is computed from the data
TemplateData=0;
for i1=1:length(stimulation)
TemplateData=TemplateData+mean(abs(Data2(:,stimulation(i1)+TemplateSamples)),1);
end
TemplateNorm=ImaGIN_normalisation(TemplateData,2,find(TimeTemplate<-0.02&TimeTemplate>-0.1));
% try to detect stimulation start to refine detection
% use mean template from d instead of data
template = 0 ;
for s=1:length(stimulation)
template = template + d( stimulation(s) + TemplateSamples ) ;
end
% use a percentage of the max as a threshold
% sometimes there are many peaks (biphasic, high fsample) and the first is not the biggest one
% so we prefer a low threshold on the max
th = 30/100 ;
offsetMaxInterval = [ -0.04 0.02 ];
offsetMaxIntervalSample = round( offsetMaxInterval * fsample(D) ) ;
segmentCenter = segmentHWidth + 1 ;
[ templateMax, ~ ] = max( abs(template(segmentCenter+(offsetMaxIntervalSample(1):offsetMaxIntervalSample(2)))) ) ;
startoffset = find( template>th*templateMax & TimeTemplate>offsetMaxInterval(1) & TimeTemplate<offsetMaxInterval(2), 1 );
if isempty(startoffset)
for th=10:-1:4
tmpoffset=find( abs(TemplateNorm)>th & abs(TimeTemplate)<0.02);
if ~isempty(tmpoffset)
break
end
end
else
startoffset = min(startoffset);
% 6ms allow to move 2 samples away at 256 Hz
tmpMax = max( abs(TemplateNorm(startoffset:startoffset+round(fsample(D)*6e-3))) ) ;
thd = 50/100 ;
tmpoffset = find( abs(TemplateNorm)>thd*tmpMax & TimeTemplate>=TimeTemplate(startoffset) & TimeTemplate<=TimeTemplate(startoffset+round(fsample(D)*6e-3)) );
% we take the sample at the first third of the interval (usefull for bipolar stim at 4096Hz)
tmpoffset = tmpoffset( floor((length(tmpoffset)-1)/3)+1 ) ;
end
if isempty(tmpoffset)
offset=nan;
else
i1=min(tmpoffset);
offset=ceil(Stim/2)+1-i1;
end
elseif 1==1
%Same offset for all pulses
TemplateData=0;
for i1=1:length(stimulation)
... ... @@ -267,7 +485,66 @@ if 1==1
end
end
stimulation=stimulation-offset;
% save and plot for comparison of methods
plot_and_save = false ;
if plot_and_save
if StimDetectVersion == 1 || StimDetectVersion == 2 || StimDetectVersion == 3 || StimDetectVersion == 4
% save th and offset in txt file
[ pathstr, name ] = fileparts(FileOut) ;
fname = fullfile( pathstr, 'offset', [ name '_offset.txt' ] ) ;
if ~exist(fileparts(fname),'dir'), mkdir(fileparts(fname)); end
fprintf( 1, [ 'Save ' fname '\n' ] ) ;
fd = fopen( fname, 'w' ) ;
fprintf( fd, [ num2str(offset) ' ' num2str(offset/fsample(D)) '\n' ] ) ;
fprintf( fd, [ num2str(th) ] ) ;
fclose(fd);
% plot offset
H = figure( 'Position', [189 118 560 700] ) ; %[30 183 560 700] ) ;
subplot(3,1,1); hold on ;
plot( TimeTemplate, avg_to_plot ) ;
line( zeros(2,1), get(gca,'YLim')','Color','k','linewidth',2);
if ~isnan(offset)
line([TimeTemplate(i1);TimeTemplate(i1)],get(gca,'YLim')','color','r','linewidth',2);
end
xlim( [ -.05 .05 ] ) ;
title( strrep( [ name ' - V' num2str(StimDetectVersion) ], '_', ' ' ) )
subplot(3,1,2); hold on ;
plot( TimeTemplate, template,'Color','b','linewidth',3 ) ;
line( zeros(2,1), get(gca,'YLim')','Color','k','linewidth',2);
if ~isempty(startoffset)
line([TimeTemplate(startoffset);TimeTemplate(startoffset)],get(gca,'YLim')','color','b','linewidth',2);
end
xlim( [ -.05 .05 ] ) ;
title( 'template accel' )
subplot(3,1,3); hold on ;
plot( TimeTemplate, TemplateNorm,'Color','g','linewidth',3 ) ;
line( zeros(2,1), get(gca,'YLim')','Color','k','linewidth',2);
if ~isnan(offset)
line([TimeTemplate(i1);TimeTemplate(i1)],get(gca,'YLim')','color','g','linewidth',2);
end
xlim( [ -.05 .05 ] ) ;
title( 'template data' )
fname = fullfile( pathstr, 'offset', [ name '_offset.png' ] ) ;
fprintf( 1, [ 'Save ' fname '...' ] ) ;
saveas(H,fname);
fprintf( 1, [ 'done' '\n' ] ) ;
close(H);
end
end
if ~isnan(offset)
stimulation=stimulation-offset;
end
end
else
... ... @@ -391,20 +668,9 @@ else
end
%remove outliers close to other stims
if length(stimulation)>1
remove=[];
for i1=1:length(stimulation)
difftmp=abs(stimulation-stimulation(i1));
if isempty(find(difftmp>0.95*Stim & difftmp<1.05*Stim))
remove=[remove i1];
end
end
% StimulationFreqU=D.fsample/median(diff(stimulation)); %uncorrected stim frequency
stimulation=stimulation(setdiff(1:length(stimulation),remove));
end
stimulation = removeOutliers( stimulation, Stim ) ;
if isempty(stimulation)
StimulationFreqU=[];
... ... @@ -486,4 +752,22 @@ fclose(fid);
% fprintf(fid,'%f\n',StimulationIndex);
% fclose(fid);
end
%remove outliers close to other stims
function stimulation = removeOutliers( stimulation, Stim )
if length(stimulation)>1
remove=[];
for i1=1:length(stimulation)
difftmp=abs(stimulation-stimulation(i1));
if isempty(find(difftmp>0.95*Stim & difftmp<1.05*Stim,1))
remove=[remove i1];
end
end
% StimulationFreqU=D.fsample/median(diff(stimulation)); %uncorrected stim frequency
stimulation=stimulation(setdiff(1:length(stimulation),remove));
end
end
... ...