ImaGIN_artefact_fit.m
3.25 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
function [sig_corr,rc] = ImaGIN_artefact_fit(Signal, mmArt, flagRC, IntCor)
% -=============================================================================
% This function is part of the ImaGIN software:
% https://f-tract.eu/
%
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE AUTHORS
% DO NOT ASSUME ANY LIABILITY OR RESPONSIBILITY FOR ITS USE IN ANY CONTEXT.
%
% Copyright (c) 2000-2018 Inserm U1216
% =============================================================================-
%
% Authors: Olivier David
%correct assuming the same RC on the series of pulses and select the RC
%value with the largest artefact power
signal=Signal(:,1:IntCor);
Uns=ones(1,size(mmArt,3));
% Power=zeros(size(signal,1),1);
if strcmp(flagRC,'same')
R2=zeros(size(signal,1),size(mmArt,2));
Power=inf*ones(size(signal,1),size(mmArt,2));
% Power=zeros(size(signal,1),size(mmArt,2));
Artefact=zeros(size(Signal,1),size(Signal,2),size(mmArt,2));
for i1=1:size(mmArt,2) %loop on RC
Art=squeeze(mmArt(:,i1,:));
m=max(abs(Art),[],2);
index=find(m>0);
Art=Art(index,:)./(m(index)*Uns);
Art=unique(Art,'rows');
for i2=1:size(Art,1)
X=[Uns(1:IntCor);Art(i2,1:IntCor)];
XX=X*X';
if rcond(XX)>10^(-15)
P=XX\X;
beta=P*signal';
artefactCut=beta(2,:)'*Art(i2,1:IntCor);
Y=repmat(mean(signal), size(signal,1),1);
% r2=1-(sum((signal-artefactCut).^2,2))./(sum((signal-Y).^2,2));
pow=sum(abs(artefactCut-signal),2);
index=find(pow<Power(:,i1));
artefact=beta(2,:)'*Art(i2,:);
if ~isempty(index)
Power(index,i1)=pow(index);
Artefact(index,:,i1)=artefact(index,:);
end
end
end
end
%select the RC with the best fit on all pulses
Power=mean(Power);
rc=min(find(Power==min(Power)));
ArtefactFit=squeeze(Artefact(:,:,rc));
% [~,rc]=max(median(R2,1));
% ArtefactFit(:,:)=squeeze(Artefact(:,:,rc));
else
%no constraint of same RC
Power=inf*ones(size(signal,1),1);
Artefact=zeros(size(signal),size(Signal,2));
Art=reshape(mmArt,size(mmArt,1)*size(mmArt,2),size(mmArt,3));
m=max(abs(Art),[],2);
index=find(m>0);
Art=Art(index,:)./(m(index)*Uns);
Art=unique(Art,'rows');
for i2=1:size(Art,1)
X=[Uns(1:IntCor);Art(i2,1:IntCor)];
P=inv(X*X')*X;
beta=P*signal';
artefactCut=beta(2,:)'*Art(i2,1:IntCor);
pow=sum(abs(artefactCut-signal),2);
index=find(pow<Power);
artefact=beta(2,:)'*Art(i2,:);
% pow=sum(abs(artefact),2);
% index=find(pow>Power);
if ~isempty(index)
Power(index)=pow(index);
Artefact(index,:)=artefact(index,:);
end
end
ArtefactFit=Artefact;
rc=0;
end
sig_corr=Signal-ArtefactFit;
end