[FieldTrip] help with inverse computing

[Ludwing Torres]

Hello, I have the next script that computes the forward model of the brain activity:
 
clear all 
close all
clc
% load electrodes positions
elec = [];
load electrodes1020
% array of electrode positions
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];
% load concentric spheres or bem 
load standard_vol.mat
% vol.cond = [1 1/80 1];
% load spheric_vol_162.mat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
cfg = [];
cfg.vol = vol; 
% cfg.elec = elec;
cfg.elec.pnt = elec.pnt(i_elec,1:3);
cfg.elec.label = elec.label(i_elec);
fs=250;
ns=250;
cfg.fsample = fs; 
time = (1:ns)/fs; 
signal1 = 1*sin(2*time*6*pi);
signal2 = 1*cos(2*time*12*pi);
signal3 = 1*cos(2*time*8*pi);

% head model
[vols, sens, cfg] = prepare_headmodel(cfg, []);

dipin=[];
 
[x,y,z]=meshgrid(linspace(-40,40,4), linspace(-40,40,4), linspace(0,30,3) );
dipin.pos=[x(:) y(:) z(:)];
mom_t=zeros(3*length(dipin.pos),ns);
 
x_orig=mom_t;
n_orig=mom_t;
n_orig(20,:)=signal1;
n_orig(80,:)=signal2;
n_orig(40,:)=signal3;
a11=1.2;
b11=0.05;
M11=3*length(dipin.pos);
A = a11*diag(ones(M11,1)) + b11*(diag(ones(M11-1,1),1) + diag(ones(M11-1,1),-1));
b12=0;
b13=0;
c11=-0.9;
for b=1:ns,
if b==1,
x_orig(:,b)=n_orig(:,b);
end
if b==2,
x_orig(:,b)=A*x_orig(:,b-1)+b12*x_orig(:,b-1).^2-b13*x_orig(:,b-1).^3+n_orig(:,b);
end
if b>2,
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);
end
end
C=compute_leadfield(dipin.pos, sens, vols);
dat=C*x_orig;
 
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.
For that, I made the next script (After runing the previous) In which I use different functions of FieldTrip to recover the x data:
 
 
 
[vol1,sens1]=prepare_vol_sens(vol, sens);
inside=inside_vol(dipin.pos,vol1); 
 
%After that I get 48 dipoles inside and 0 outside
dipin.outside=[];
dipin.inside=1:48;

 
%Inverse computing using MNE
dipout=minimumnormestimate(dipin,sens1,vol1,dat); 
 

% I store the dipole moments obtained as a matrix
Xres=[];
for ii=1:48
Xres=[Xres;dipout.mom{ii}];
end
 
% I compute a relative error between the original x and the obtained x to compare

error1=norm(Xres-x_orig,'fro')/norm(x_orig,'fro')
 

%Using Residual variance
[dipout2] = residualvariance(dipin, sens1, vol1, dat);

Xres2=[];
for ii=1:48
Xres2=[Xres2;dipout2.mom{ii}];
end

error2=norm(Xres2-x_orig,'fro')/norm(x_orig,'fro')
 


%Using dipole fitting
[dipout4] = dipole_fit(dipin, sens1, vol1, dat)

error4=norm(dipout4.mom-x_orig,'fro')/norm(x_orig,'fro')
 
The error obtained using residual variance is very different to the error of the two other methods.
Besides, I get a very big error (92%) with the other two methods.
I don't know if I'm doing something wrong, please If anyone could help me, I could pay you.
Thanks
 
C=compute_leadfield(dipin.pos, sens, vols);
dat=C*x_orig; 		 	   		  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20110127/4e39e4ff/attachment-0001.html>