Changeset 82


Ignore:
Timestamp:
Jun 16, 2009, 10:54:45 AM (10 years ago)
Author:
tjniskan
Message:
  • Some minor modifications to aedes_readfid.m
  • Added BOLD strength estimation to

plugins/fmri_plugins/basic_fmri_analysis.m

  • Changed the image in about dialog
  • Added support for negative threshold in misclib/fmri_blob_overlay.m

M misclib/fmri_blob_overlay.m
M aedes_helpabout.m
M plugins/fmri_plugins/basic_fmri_analysis.m
M aedes_readfid.m
M aedes_revision.m

Files:
5 edited

Legend:

Unmodified
Added
Removed
  • aedes_helpabout.m

    r80 r82  
    6666imax = axes('parent',fh,...
    6767            'units','pixel',...
    68             'position',[0.05*fig_w 0.12*fig_h 0.9*fig_w 0.9*fig_w],...
     68            'position',[(fig_w-256)/2 0.65*fig_h 256 130],...
    6969            'visible','off',...
    7070            'ydir','reverse',...
     
    7676[fp,fn,fe]=fileparts(fpath);
    7777try
    78   imdata = imread([fp,filesep,'about_aedes.png']);
    79   imdata = imdata(:,:,1);
     78  imdata = imread([fp,filesep,'aedes_logo.png']);
     79  %imdata = imdata(:,:,1);
    8080catch
    8181  delete(fh);
    82   hh=errordlg('Cannot find file "about_aedes.png"!','Error!',...
     82  hh=errordlg('Cannot find file "aedes_logo.png"!','Error!',...
    8383              'modal');
    8484  return
    8585end
     86
    8687sz=size(imdata);
    8788im=image('parent',imax,'cdata',imdata,...
    88          'cdatamapping','scaled');
    89 set(imax,'xlim',[0.5 sz(1)+0.5],...
    90          'ylim',[0.5 sz(2)+0.5],...
    91          'clim',[min(min(imdata)) max(max(imdata))])
    92 set(imax,'clim',[-500 255])
     89  'cdatamapping','scaled');
     90set(imax,'xlim',[0.5 sz(2)+0.5],...
     91  'ylim',[0.5 sz(1)+0.5])
    9392
    94 %% Title text
    95 shadow_tx = text('parent',bgax,...
    96                  'horizontalalign','left',...
    97                  'units','normal',...
    98                  'position',[0.054 0.976],...
    99                  'verticalalign','top',...
    100                  'string','Aedes 1.0',...
    101                  'fontsize',22,...
    102                  'fontweig','bold',...
    103                  'color',[0 0 0]);
    104 title_tx = text('parent',bgax,...
    105                 'horizontalalign','left',...
    106                 'units','normal',...
    107                 'position',[0.05 0.98],...
    108                 'verticalalign','top',...
    109                 'string','Aedes 1.0',...
    110                 'fontsize',22,...
    111                 'fontweig','bold',...
    112                 'color',[0 0 0.85]);
     93% %% Title text
     94% shadow_tx = text('parent',bgax,...
     95%                  'horizontalalign','left',...
     96%                  'units','normal',...
     97%                  'position',[0.054 0.976],...
     98%                  'verticalalign','top',...
     99%                  'string','Aedes 1.0',...
     100%                  'fontsize',22,...
     101%                  'fontweig','bold',...
     102%                  'color',[0 0 0]);
     103% title_tx = text('parent',bgax,...
     104%                 'horizontalalign','left',...
     105%                 'units','normal',...
     106%                 'position',[0.05 0.98],...
     107%                 'verticalalign','top',...
     108%                 'string','Aedes 1.0',...
     109%                 'fontsize',22,...
     110%                 'fontweig','bold',...
     111%                 'color',[0 0 0.85]);
    113112
    114113% Version text
     
    117116                  'horizontalalign','left',...
    118117                  'units','normal',...
    119                   'position',[0.10 0.89],...
     118                  'position',[0.10 0.63],...
    120119                  'verticalalign','top',...
    121120                  'string',...
    122                                   sprintf(['version 1.0 rev %d\n(%s)'],rev,repo),...
     121                                  sprintf(['version 1.0 rev %d\n%s'],rev,'http://aedes.uku.fi'),...
    123122                  'fontsize',10,...
    124123                  'fontweig','bold',...
     
    159158
    160159% Contact info
    161 contact_shadow = text('parent',imax,...
     160contact_shadow = text('parent',bgax,...
    162161                      'horizontalalign','left',...
    163162                      'units','normal',...
    164                       'position',[0.053 0.047],...
     163                      'position',[0.1 0.117],...
    165164                      'verticalalign','bottom',...
    166165                      'string',{'Contact information:','',...
    167166                    'Juha-Pekka Niskanen,',...
    168                     'University of Kuopio,',...
    169                     'Department of Physics,',...
    170                     'PL 1627, 70211 Kuopio, Finland'...
     167                    'Biomedical NMR Group,',...
     168                    'A. I. Virtanen Institute for Molecular Sciences,',...
     169                    'University of Kuopio, Finland'...
    171170                    'Email: Juha-Pekka.Niskanen@uku.fi'},...
    172171                      'fontsize',8,...
    173172                      'fontweig','bold',...
    174173                      'color',[1 1 1]);
    175 contact_tx = text('parent',imax,...
     174contact_tx = text('parent',bgax,...
    176175                  'horizontalalign','left',...
    177176                  'units','normal',...
    178                   'position',[0.05 0.05],...
     177                  'position',[0.1 0.12],...
    179178                  'verticalalign','bottom',...
    180179                  'string',{'Contact information:','',...
    181180                    'Juha-Pekka Niskanen,',...
    182                     'University of Kuopio,',...
    183                     'Department of Physics,',...
    184                     'PL 1627, 70211 Kuopio, Finland'...
     181                    'Biomedical NMR Group,',...
     182                    'A. I. Virtanen Institute for Molecular Sciences,',...
     183                    'University of Kuopio, Finland'...
    185184                    'Email: Juha-Pekka.Niskanen@uku.fi'},...
    186185                  'fontsize',8,...
  • aedes_readfid.m

    r81 r82  
    492492  end
    493493 
     494  % If the sequence is a fast spin echo, try to construct the phasetable
     495  if isempty(phasetable) && isfield(procpar,'petable') && ...
     496      strncmpi(procpar.petable{1},'fse',3)
     497    phasetable = {'t1',l_CreateFsemsPhaseTable(procpar)};
     498  end
     499 
    494500  %% Abort and throw error if phasetable cannot be read and the FID-file
    495501  % has not been sorted
     
    13531359  % patches that hopefully work...
    13541360 
    1355  
    1356   if strcmp(seqcon,'nccnn') && length(size(kspace))==2 && ...
    1357           procpar.ns>=1 && AcqType~=1
     1361  if strcmp(seqcon,'nncnn')
     1362    kspace = permute(kspace,[1 3 2]);
     1363  elseif strcmp(seqcon,'nccnn') && length(size(kspace))==2 && ...
     1364      procpar.ns>=1 && AcqType~=1
    13581365    if ~isempty(Dat.phasetable)
    1359           kssz = size(kspace);
    1360           phsz = size(Dat.phasetable);
    1361           kspace=permute(reshape(kspace,procpar.np/2,...
    1362                 phsz(2),...
    1363                 kssz(2)/phsz(2)),[1 3 2]);
    1364           kspace = permute(reshape(kspace,...
    1365                 procpar.np/2,procpar.ns,kssz(2)/procpar.ns),[1 3 2]);
    1366           if Dat.Sorting
    1367                 kspace(:,Dat.phasetable(:),:)=kspace;
    1368           end
    1369         else
    1370           kspace = reshape(kspace,...
    1371                 [procpar.np/2 procpar.ns ...
    1372                 size(kspace,2)/procpar.ns]);
    1373           kspace=permute(kspace,[1 3 2]);
    1374         end
    1375        
     1366      kssz = size(kspace);
     1367      phsz = size(Dat.phasetable);
     1368      kspace=permute(reshape(kspace,procpar.np/2,...
     1369        phsz(2),...
     1370        kssz(2)/phsz(2)),[1 3 2]);
     1371      kspace = permute(reshape(kspace,...
     1372        procpar.np/2,procpar.ns,kssz(2)/procpar.ns),[1 3 2]);
     1373      if Dat.Sorting
     1374        kspace(:,Dat.phasetable(:),:)=kspace;
     1375      end
     1376    else
     1377      kspace = reshape(kspace,...
     1378        [procpar.np/2 procpar.ns ...
     1379        size(kspace,2)/procpar.ns]);
     1380      kspace=permute(kspace,[1 3 2]);
     1381    end
     1382   
    13761383  elseif strcmp(seqcon,'nscsn') && length(size(kspace))==3 && ...
    13771384          AcqType~=1
     
    15781585tt = tt(:);
    15791586
    1580 phasetable = [reshape(t(1:nv/2),etl*2,[]);...
    1581   flipud(reshape(tt(1:nv/2),etl*2,[]))];
     1587phasetable = [reshape(t(1:nv/2),[],etl);...
     1588  flipud(reshape(tt(1:nv/2),[],etl))];
    15821589phasetable = circshift(fliplr(phasetable),[0 kzero-1]);
    15831590
  • aedes_revision.m

    r81 r82  
    9393% bash-script every time it is called so that this file "aedes_revision.m" is
    9494% always in the list of committed files. DO NOT EDIT THE NEXT LINE!!!
    95 % - SVN Hook -
     95% - Svn Hook -
  • misclib/fmri_blob_overlay.m

    r80 r82  
    2222%       **********              ************
    2323%
     24%       'Deactivation'         'on','off': Select if deactivation blobs are also
     25%                               displayed (default = 'off')
     26%
    2427%       'BlobCMin'              Min color limit for scaling activation
    2528%                               blobs. (default = min(ImBlob))
     
    3033%       'AnatCLim'              A two element vector [CLOW CHIGH] for
    3134%                               scaling anatomical image. (default =
    32 %                               min/max) 
     35%                               min/max)
    3336%
    3437%       'BlobCmap'              String of a valid colormap name
    3538%                               ('jet','hot',etc.) to be used with
    3639%                               activation blobs. (default = 'hot')
     40%
     41%       'BlobNegCmap'           String of a valid colormap name
     42%                               ('cool','jet',etc.) to be used with
     43%                               deactivation blobs. (default = 'cool')
    3744
    3845%       'AnatCmap'               String of a valid colormap name
     
    7279blobcmin = [];
    7380blobcmax = [];
     81blobnegcmin = [];
     82blobnegcmax = [];
    7483anatcmin = [];
    7584anatcmax = [];
     85deactivations = false;
    7686blobcmap = 'hot';
     87blobnegcmap = 'cool';
    7788anatcmap = 'gray';
    7889vnmrfix = false;
     
    95106 
    96107  switch param
     108    case 'deactivation'
     109      if strcmpi(value,'on')
     110        deactivations = true;
     111      end
    97112    case 'blobclim'
    98113      blobcmin = value(1);
     
    100115    case 'blobcmin'
    101116      blobcmin = value;
     117    case 'blobnegcmin'
     118      blobnegcmin = abs(value);
     119    case 'blobnegcmax'
     120      blobnegcmax = abs(value);
    102121    case 'blobcmax'
    103122      blobcmax = value;
     
    107126    case 'blobcmap'
    108127      blobcmap = value;
     128    case 'blobnegcmap'
     129      blobnegcmap = value;
    109130    case 'anatcmap'
    110131      anatcmap = value;
     
    152173end
    153174
    154 
     175if length(thold)==2
     176  if isnan(thold(2))
     177    thold = thold(1);
     178  else
     179    deactivations = true;
     180  end
     181end
     182
     183if deactivations
     184  tholdneg = abs(thold(2));
     185  thold = thold(1);
     186end
    155187
    156188% Check image sizes ------------------------------------
     
    170202end
    171203
     204ImBlobNeg = -ImBlob;
     205
    172206% Threshold blobs
    173207ind_lo = find(ImBlob<=thold);
     
    175209minT = min(ImBlob(:));
    176210maxT = max(ImBlob(:));
     211
     212if deactivations
     213  ind_neg_lo = find(ImBlobNeg<=tholdneg);
     214  ind_neg_hi = find(ImBlobNeg>tholdneg);
     215  minTneg = min(ImBlobNeg(:));
     216  maxTneg = max(ImBlobNeg(:));
     217end
    177218%ImBlob(ImBlob<=thold) = 0;
    178219%ind = find(ImBlob~=0);
     
    184225  ind_hi = find(ImBloIndHi);
    185226end
     227if deactivations
     228  ImBloNegIndHi = ImBlobNeg>tholdneg;
     229  ImBloNegIndHi = bwareaopen(ImBloNegIndHi,voxelsizemin,8);
     230  ind_neg_hi = find(ImBloNegIndHi);
     231end
    186232
    187233% Handle clims for blobs and anatomical image -----------
     
    196242ImBlob=(ImBlob-min(ImBlob(:)));
    197243
     244if deactivations
     245  if ~isempty(blobnegcmin)
     246    ImBlobNeg(ImBlobNeg<blobnegcmin) = blobnegcmin;
     247    minTneg = blobnegcmin;
     248  end
     249  if ~isempty(blobnegcmax)
     250    ImBlobNeg(ImBlobNeg>blobnegcmax) = blobnegcmax;
     251    maxTneg = blobnegcmax;
     252  end
     253  ImBlobNeg=(ImBlobNeg-min(ImBlobNeg(:)));
     254end
     255
    198256if ~isempty(anatcmin)
    199257  ImAnat(ImAnat<anatcmin) = anatcmin;
     
    204262
    205263% Construct indexed image -------------------------------
    206 ImBlob = round(ImBlob./(maxT-minT)*255)+1;
    207 ImAnat = round(ImAnat./max(ImAnat(:))*255)+1;
     264ImBlob = round(ImBlob./(maxT-minT)*255);
     265if deactivations
     266  ImBlobNeg = round(ImBlobNeg./(maxTneg-minTneg)*255);
     267  ImBlobNeg = abs(ImBlobNeg-256);
     268end
     269ImAnat = round(ImAnat./max(ImAnat(:))*255);
    208270ImRes = ImAnat;
    209 ImRes(ind_hi) = ImBlob(ind_hi)+256;
    210 %ImRes = ImRes+ImBlob;
     271if deactivations
     272  ImRes(ind_neg_hi) = ImBlobNeg(ind_neg_hi)+256;
     273  ImRes(ind_hi) = ImBlob(ind_hi)+2*256;
     274else
     275  ImRes(ind_hi) = ImBlob(ind_hi)+256;
     276end
    211277
    212278% Construct joint colormap
    213 ImBlobColormap=eval([blobcmap,'(256)']);
     279if deactivations
     280  tmp = eval([blobcmap,'(256)']);
     281  %tmp2 = jet(512);
     282  %tmp2=tmp2(1:256,:);
     283  tmp2 = eval([blobnegcmap,'(256)']);
     284  ImBlobColormap = [tmp2;tmp];
     285else
     286  ImBlobColormap=eval([blobcmap,'(256)']);
     287end
    214288ImAnatColormap=eval([anatcmap,'(256)']);
    215289StackedColormap = [ImAnatColormap;ImBlobColormap];
     
    219293  cmap = StackedColormap;
    220294  im_rgb = ind2rgb(im,cmap);
    221   minmaxT = [minT maxT];
    222 end
    223 
    224 
     295  if deactivations
     296     minmaxT = [minT maxT -maxTneg -minTneg];
     297  else
     298    minmaxT = [minT maxT];
     299  end
     300end
     301
     302
  • plugins/fmri_plugins/basic_fmri_analysis.m

    r78 r82  
    55
    66if AddInfo.isDataMixed
    7   return
     7  % Make data a 4D-matrix if it is loaded in aedes as a stack of 2D or 3D images
     8  DATA2{1}=DATA{1};
     9  DATA2{1}.FTDATA = zeros([size(DATA{1}.FTDATA,1),...
     10    size(DATA{1}.FTDATA,2),size(DATA{1}.FTDATA,3),length(DATA)],...
     11    class(DATA{1}.FTDATA));
     12  for ii=1:length(DATA)
     13    DATA2{1}.FTDATA(:,:,:,ii)=DATA{ii}.FTDATA;
     14  end
     15  DATA=DATA2;
    816end
    917
     
    5361    'hipass',hipass,'omitvolumes',omitVols,...
    5462    'mask',ROI(roi_ind).voxels{1}(:,:,:,1));
     63  ROI(roi_ind)=[];
    5564else
    5665  [maps_out,fdr_th] = fmri_analysis(DATA{1},'TR',TR,...
     
    93102  DATA{1}.FTDATA(:,:,:,omitVols)=[];
    94103  overlay.FTDATA(:,:,:,omitVols)=[];
     104else
     105  new_onset = onset;
     106  new_durat = durat;
     107end
     108
     109
     110% Estimate BOLD strength for ROIs
     111if ~isempty(ROI)
     112  nScans = size(DATA{1}.FTDATA,4);
     113  fprintf(1,'---------------------------------------\n');
     114  for kk=1:length(ROI)
     115   
     116    % Get time-series indices from ROI
     117    ind = repmat(ROI(kk).voxels{1}(:,:,:,1),[1 1 1 size(DATA{1}.FTDATA,4)]);
     118
     119    % Mean EPI time series
     120    ts_data = reshape(double(DATA{1}.FTDATA(ind)),[],size(DATA{1}.FTDATA,4));
     121    ts_data = ts_data.';
     122    ts_data = detrend(ts_data)+repmat(mean(ts_data),size(ts_data,1),1);
     123   
     124    % Normalize mean to bold-%
     125    mean_ts_data = mean(ts_data,2);
     126   
     127    % estimate BOLD-% for the ROI mean
     128    [bold_prc,z,th,z_hat]=estimate_bold(new_onset,new_durat,TR,nScans,mean_ts_data);
     129   
     130    z_norm = (z./th(2)-1)*100;
     131    z_hat_norm = (z_hat./th(2)-1)*100;
     132   
     133    % Normalize raw timeseries to bold-%
     134    %ts_data_bold = (ts_data./repmat(mean(ts_data(1:20,:)),size(ts_data,1),1)-1)*100;
     135   
     136    % Plot results
     137    figure;
     138    ax=axes;
     139    line(1:length(z),z_norm,'color','k',...
     140      'parent',ax);
     141    line(1:length(z),z_hat_norm,...
     142      'color','r','linewidth',2,'parent',ax);
     143    title(['Time series for ROI: ',ROI(kk).label])
     144    ylabel('BOLD-%');
     145    set(ax,'xlim',[0 length(z)],...
     146      'ylim',[min(z_norm)-min(z_norm)*0.05 ...
     147      max(z_norm)+max(z_norm)*0.05]);
     148   
     149    fprintf(1,'BOLD-%% %s: %.3f\n',ROI(kk).label,bold_prc);
     150  end
     151  fprintf(1,'---------------------------------------\n');
     152else
     153  % Estimate BOLD strength for the time-series corresponding to the
     154  % maximum value in T-map.
     155  nScans = size(DATA{1}.FTDATA,4);
     156  [mx,I_mx]=max(maps_out.tmap(:));
     157  [II_mx,JJ_mx,KK_mx]=ind2sub(size(maps_out.tmap),I_mx);
     158  [mn,I_mn]=min(maps_out.tmap(:));
     159  [II_mn,JJ_mn,KK_mn]=ind2sub(size(maps_out.tmap),I_mn);
     160 
     161  ts_data_max = squeeze(DATA{1}.FTDATA(II_mx,JJ_mx,KK_mx,:));
     162  ts_data_min = squeeze(DATA{1}.FTDATA(II_mn,JJ_mn,KK_mn,:));
     163 
     164  [bold_prc_max,z_max,th_max,z_hat_max]=estimate_bold(new_onset,new_durat,TR,nScans,ts_data_max);
     165  [bold_prc_min,z_min,th_min,z_hat_min]=estimate_bold(new_onset,new_durat,TR,nScans,ts_data_min);
     166 
     167  z_norm_max = (z_max./th_max(2)-1)*100;
     168  z_hat_norm_max = (z_hat_max./th_max(2)-1)*100;
     169 
     170  z_norm_min = (z_min./th_min(2)-1)*100;
     171  z_hat_norm_min= (z_hat_min./th_min(2)-1)*100;
     172 
     173  % Plot results
     174  figure;
     175  ax1=subplot(2,1,1);
     176  line(1:length(z_max),z_norm_max,'color','k',...
     177    'parent',ax1);
     178  line(1:length(z_max),z_hat_norm_max,...
     179    'color','r','linewidth',2,'parent',ax1);
     180  title('Time series for Max T-value');
     181  ylabel('BOLD-%');
     182  set(ax1,'xlim',[0 length(z_max)],...
     183    'ylim',[min(z_norm_max)-min(z_norm_max)*0.05 ...
     184    max(z_norm_max)+max(z_norm_max)*0.05]);
     185 
     186  ax2=subplot(2,1,2);
     187  line(1:length(z_min),z_norm_min,'color','k',...
     188    'parent',ax2);
     189  line(1:length(z_min),z_hat_norm_min,...
     190    'color','r','linewidth',2,'parent',ax2);
     191  title('Time series for Min T-value');
     192  ylabel('BOLD-%');
     193  set(ax2,'xlim',[0 length(z_min)],...
     194    'ylim',[min(z_norm_min)-min(z_norm_min)*0.05 ...
     195    max(z_norm_min)+max(z_norm_min)*0.05]);
     196 
     197  % Print percents to command window
     198  fprintf(1,'---------------------------------------\n');
     199  fprintf(1,'Max T BOLD-%%: %.3f\n',bold_prc_max);
     200  fprintf(1,'Min T BOLD-%%: %.3f\n',bold_prc_min);
     201  fprintf(1,'---------------------------------------\n');
    95202end
    96203
     
    104211
    105212
     213function [bold_prc,z,th,z_hat]=estimate_bold(onset,durat,TR,num_scans,ts_data)
     214%% -------------------------------------------------
     215
     216% Load Rat HRF (8 seconds long) NOTE: DON'T USE FOR HUMAN DATA!!!
     217bf = [
     218  0
     219  0.000008876945146
     220  0.000387621003097
     221  0.003012441669394
     222  0.011548214895067
     223  0.030056522080168
     224  0.061233628208021
     225  0.105349995559045
     226  0.160158500633323
     227  0.221528149631163
     228  0.284405279593501
     229  0.343761753554314
     230  0.395323573307322
     231  0.436002661868442
     232  0.464045998242252
     233  0.478965423658889
     234  0.481327079129854
     235  0.472473786978796
     236  0.454237528221769
     237  0.428680153860941
     238  0.397883140401152
     239  0.363793656845973
     240  0.328124931280142
     241  0.292303457117353
     242  0.257453130446147
     243  0.224406054242129
     244  0.193730677047637
     245  0.165769521647024
     246  0.140680556701329
     247  0.118477987483032
     248  0.099069734765398
     249  0.082290068881472
     250  0.067926763712502
     251  0.055742762013807
     252  0.045492744488986
     253  0.036935220196777
     254  0.029840852217657
     255  0.023997740489123
     256  0.019214335904674
     257  0.015320580887648
     258  0.012167779438633
     259  0.009627606069476
     260  0.007590575511228
     261  0.005964217696074
     262  0.004671136982604
     263  0.003647081075751
     264  0.002839102788446
     265  0.002203865389816
     266  0.001706118264261
     267  0.001317352436400
     268  0.001014633776781
     269  0.000779604140824
     270  0.000597636251764
     271  0.000457125954290
     272  0.000348904854607
     273  0.000265756797153
     274  0.000202022712087
     275  0.000153279812313
     276  0.000116082720651
     277  0.000087755728843
     278  0.000066226941806
     279  0.000049896490408
     280  0.000037532277385
     281  ];
     282
     283
     284%TR=2.039;
     285%onset=[30 75 120]';
     286%durat=[15 15 15]';
     287onset = onset(:);
     288durat = durat(:);
     289
     290%bf = spm_hrf(1/16);
     291
     292% Create stimulus function (32 bin offset)
     293%k = 165; % Number of scans
     294k = num_scans;
     295T     = 16;
     296dt    = TR/T;
     297u     = ones(size(onset));
     298if ~any(durat)
     299  u  = u/dt;
     300end
     301ton = round(onset*TR/dt) + 32;                  % onsets
     302tof = round(durat*TR/dt) + ton + 1;                     % offset
     303sf = zeros((k*T + 128),size(u,2));
     304
     305for j = 1:length(ton)
     306  if numel(sf)>ton(j),
     307    sf(ton(j),:) = sf(ton(j),:) + u(j,:);
     308  end;
     309  if numel(sf)>tof(j),
     310    sf(tof(j),:) = sf(tof(j),:) - u(j,:);
     311  end;
     312end
     313sf        = cumsum(sf);                                 % integrate
     314
     315% Convolve stimulus with the HRF
     316conv_sf = conv(sf,bf);
     317
     318% Resample the convolved stimulus
     319R = conv_sf([0:(k-1)]*16+1+32);
     320R=R(:);
     321
     322% Remove linear trend but keep the mean
     323mean_ts_data = mean(ts_data);
     324ts_data = detrend(ts_data)+mean_ts_data;
     325
     326% Estimate baseline from the paradigm
     327bl_ind = find(sf==0);
     328bl_ind = bl_ind(1:min(length(bl_ind),20));
     329if isempty(bl_ind)
     330  baseline = mean(ts_data);
     331else
     332  baseline = mean(ts_data(bl_ind));
     333end
     334
     335mean_ts_data_bold = (ts_data./baseline-1)*100;
     336
     337
     338% Fit block model to the data
     339%z=mean_ts_data_bold;
     340z=ts_data;
     341H = [R ones(size(R))];
     342th = H\z;
     343z_hat = H*th;
     344
     345% Calculate results
     346bold_prc = ((th(1)*max(R))/th(2))*100;
     347
Note: See TracChangeset for help on using the changeset viewer.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008