<html>
<head>
<style><!--
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Tahoma
}
--></style>
</head>
<body class='hmmessage'>
Hello, I have the next script that computes the forward model of the brain activity:<BR>
<BR><FONT face="Courier New" size=2>
clear </FONT><FONT face="Courier New" color=#a020f0 size=2>all</FONT><FONT face="Courier New" size=2> <BR>
close </FONT><FONT face="Courier New" color=#a020f0 size=2>all<BR></FONT><FONT face="Courier New" size=2>
clc</FONT><FONT face="Courier New" color=#228b22 size=2><BR>
% load electrodes positions<BR></FONT><FONT face="Courier New" size=2>
elec = [];<BR>
load </FONT><FONT face="Courier New" color=#a020f0 size=2>electrodes1020<BR></FONT><FONT face="Courier New" color=#228b22 size=2>
% array of electrode positions<BR></FONT><FONT face="Courier New" size=2>
i_elec=[4 6 18 19 21 23 25 27 28 31 33 35 37 40 41 43 45 47 49 50 53 55 57 59 62 63 65 67 69 71 72 84 86];<BR></FONT><FONT face="Courier New" color=#228b22 size=2>
% load concentric spheres or bem <BR></FONT><FONT face="Courier New" size=2>
load </FONT><FONT face="Courier New" color=#a020f0 size=2>standard_vol.mat</FONT><FONT face="Courier New" color=#228b22 size=2><BR>
% vol.cond = [1 1/80 1];<BR>
% load spheric_vol_162.mat<BR>
<BR>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% <BR></FONT><FONT face="Courier New" size=2>
cfg = [];<BR>
cfg.vol = vol; <BR></FONT><FONT face="Courier New" color=#228b22 size=2>
% cfg.elec = elec;<BR></FONT><FONT face="Courier New" size=2>
cfg.elec.pnt = elec.pnt(i_elec,1:3);<BR>
cfg.elec.label = elec.label(i_elec);<BR>
fs=250;<BR>
ns=250;<BR>
cfg.fsample = fs; <BR>
time = (1:ns)/fs; <BR>
signal1 = 1*sin(2*time*6*pi);<BR>
signal2 = 1*cos(2*time*12*pi);<BR>
signal3 = 1*cos(2*time*8*pi);<BR>
<BR></FONT><FONT face="Courier New" color=#228b22 size=2>
% head model<BR></FONT><FONT face="Courier New" size=2>
[vols, sens, cfg] = prepare_headmodel(cfg, []);<BR>
<BR>
dipin=[];<BR></FONT><FONT face="Courier New" color=#228b22 size=2>
<BR></FONT><FONT face="Courier New" size=2>
[x,y,z]=meshgrid(linspace(-40,40,4), linspace(-40,40,4), linspace(0,30,3) );<BR>
dipin.pos=[x(:) y(:) z(:)];<BR>
mom_t=zeros(3*length(dipin.pos),ns);<BR></FONT><FONT face="Courier New" color=#228b22 size=2>
<BR></FONT><FONT face="Courier New" size=2>
x_orig=mom_t;<BR>
n_orig=mom_t;<BR>
n_orig(20,:)=signal1;<BR>
n_orig(80,:)=signal2;<BR>
n_orig(40,:)=signal3;<BR>
a11=1.2;</FONT><FONT face="Courier New" color=#228b22 size=2><BR></FONT><FONT face="Courier New" size=2>
b11=0.05;</FONT><FONT face="Courier New" color=#228b22 size=2><BR></FONT><FONT face="Courier New" size=2>
M11=3*length(dipin.pos);<BR>
A = a11*diag(ones(M11,1)) + b11*(diag(ones(M11-1,1),1) + diag(ones(M11-1,1),-1));<BR>
b12=0;</FONT><FONT face="Courier New" color=#228b22 size=2><BR></FONT><FONT face="Courier New" size=2>
b13=0;</FONT><FONT face="Courier New" color=#228b22 size=2><BR></FONT><FONT face="Courier New" size=2>
c11=-0.9;</FONT><FONT face="Courier New" color=#228b22 size=2><BR></FONT><FONT face="Courier New" color=#0000ff size=2>
for</FONT><FONT face="Courier New" size=2> b=1:ns,<BR>
</FONT><FONT face="Courier New" color=#0000ff size=2>if</FONT><FONT face="Courier New" size=2> b==1,<BR>
x_orig(:,b)=n_orig(:,b);<BR>
</FONT><FONT face="Courier New" color=#0000ff size=2>end<BR></FONT><FONT face="Courier New" size=2>
</FONT><FONT face="Courier New" color=#0000ff size=2>if</FONT><FONT face="Courier New" size=2> b==2,<BR>
x_orig(:,b)=A*x_orig(:,b-1)+b12*x_orig(:,b-1).^2-b13*x_orig(:,b-1).^3+n_orig(:,b);<BR>
</FONT><FONT face="Courier New" color=#0000ff size=2>end<BR></FONT><FONT face="Courier New" size=2>
</FONT><FONT face="Courier New" color=#0000ff size=2>if</FONT><FONT face="Courier New" size=2> b>2,<BR>
x_orig(:,b)=A*x_orig(:,b-1)+b12*x_orig(:,b-1).^2+b13*x_orig(:,b-1).^3+c11*x_orig(:,b-2)+n_orig(:,b);<BR>
</FONT><FONT face="Courier New" color=#0000ff size=2>end<BR>
end<BR><FONT face="Courier New" size=2>
C=compute_leadfield(dipin.pos, sens, vols);<BR>
dat=C*x_orig;<BR>
<BR>
now I have to compute the inverse model, that is, I have to obtain the x_orig data, knowing only the data stored in dat and C.<BR>
For that, I made the next script (After runing the previous) In which I use different functions of FieldTrip to recover the x data:<BR>
<BR>
<BR><FONT face="Courier New" size=2>
<BR>
[vol1,sens1]=prepare_vol_sens(vol, sens);<BR>
inside=inside_vol(dipin.pos,vol1); <BR>
<BR>
%After that I get 48 dipoles inside and 0 outside<BR>
</FONT><FONT face="Courier New" size=2>dipin.outside=[];<BR>
dipin.inside=1:48;<BR>
<BR></FONT><FONT face="Courier New" size=2>
<BR>
%Inverse computing using MNE<BR>
dipout=minimumnormestimate(dipin,sens1,vol1,dat); <BR>
<BR>
<BR></FONT><FONT face="Courier New" color=#228b22 size=2>
% I store the dipole moments obtained as a matrix<BR></FONT><FONT face="Courier New" size=2>
Xres=[];<BR></FONT><FONT face="Courier New" color=#0000ff size=2>
for</FONT><FONT face="Courier New" size=2><FONT color=#000000> ii=1:48</FONT><BR>
Xres=[Xres;dipout.mom{ii}];<BR></FONT><FONT face="Courier New" color=#0000ff size=2>
end<BR>
<BR>
% I compute a relative error between the original x and the obtained x to compare<BR>
<BR></FONT><FONT face="Courier New" size=2>
error1=norm(Xres-x_orig,</FONT><FONT face="Courier New" color=#a020f0 size=2>'fro'</FONT><FONT face="Courier New" size=2>)/norm(x_orig,</FONT><FONT face="Courier New" color=#a020f0 size=2>'fro'</FONT><FONT face="Courier New" size=2>)</FONT><BR>
<FONT face="Courier New" size=2> <BR>
<BR></FONT><FONT face="Courier New" color=#228b22 size=2>
%Using Residual variance<BR></FONT><FONT face="Courier New" size=2>
[dipout2] = residualvariance(dipin, sens1, vol1, dat);<BR>
<BR>
Xres2=[];<BR></FONT><FONT face="Courier New" color=#0000ff size=2>
for</FONT><FONT face="Courier New" size=2><FONT color=#000000> ii=1:48</FONT><BR>
Xres2=[Xres2;dipout2.mom{ii}];<BR></FONT><FONT face="Courier New" color=#0000ff size=2>
end<BR>
<BR></FONT><FONT face="Courier New" size=2>
error2=norm(Xres2-x_orig,</FONT><FONT face="Courier New" color=#a020f0 size=2>'fro'</FONT><FONT face="Courier New" size=2>)/norm(x_orig,</FONT><FONT face="Courier New" color=#a020f0 size=2>'fro'</FONT><FONT face="Courier New" size=2>)</FONT><BR>
<FONT face="Courier New" size=2> <BR>
<BR>
<BR></FONT><FONT face="Courier New" color=#228b22 size=2>
%Using dipole fitting<BR></FONT><FONT face="Courier New" size=2>
[dipout4] = dipole_fit(dipin, sens1, vol1, dat)<BR>
<BR>
error4=norm(dipout4.mom-x_orig,</FONT><FONT face="Courier New" color=#a020f0 size=2>'fro'</FONT><FONT face="Courier New" size=2>)/norm(x_orig,</FONT><FONT face="Courier New" color=#a020f0 size=2>'fro'</FONT><FONT face="Courier New" size=2>)</FONT><BR>
<FONT face="Courier New" size=2></FONT> <BR>
<FONT face="Courier New" size=2>The error obtained using residual variance is very different to the error of the two other methods.</FONT><BR>
Besides, I get a very big error (92%) with the other two methods.<BR>
<FONT face="Courier New" size=2>I don't know if I'm doing something wrong, please If anyone could help me, I could pay you.</FONT><BR>
<FONT face="Courier New" size=2>Thanks<BR></FONT></FONT>
<BR></FONT><FONT face="Courier New" size=2>
C=compute_leadfield(dipin.pos, sens, vols);<BR>
dat=C*x_orig;<BR></FONT> </body>
</html>