vb_observation_noise_specification

 Specify observation noise covariance structure and its magnitude

 --- Syntax
 function [COV,sx0] = vb_observation_noise_specification(noise_parm)

 --- Input
 (Fields of noise_parm)
 .noise_model     : Integer number for specifying noise model
 .megfile_baseline: MEG file for noise calculation
 .twin_noise      : Time window for noise estimation
 .noise_reg       : Regularization parameter for noise covariance
 .temporal_filter (optional): 
 .trial_average   (optional): 

  1.  trial average (optional)
  2.  temporal filter (optional)
  3.  calculate BB' for each session
  4.  calculate sx0  

 --- History
 2005-08-22 
   ver.30b
 2005-12-22 
   bug fix (megfile --> megnoisefile)
 2008-06-27 Taku Yoshioka 
   (megnoisefile -> megfile_baseline)
 2008-11-28 Taku Yoshioka
   Use vb_disp() for displaying message

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [COV,sx0] = vb_observation_noise_specification(noise_parm)
0002 % Specify observation noise covariance structure and its magnitude
0003 %
0004 % --- Syntax
0005 % function [COV,sx0] = vb_observation_noise_specification(noise_parm)
0006 %
0007 % --- Input
0008 % (Fields of noise_parm)
0009 % .noise_model     : Integer number for specifying noise model
0010 % .megfile_baseline: MEG file for noise calculation
0011 % .twin_noise      : Time window for noise estimation
0012 % .noise_reg       : Regularization parameter for noise covariance
0013 % .temporal_filter (optional):
0014 % .trial_average   (optional):
0015 %
0016 %  1.  trial average (optional)
0017 %  2.  temporal filter (optional)
0018 %  3.  calculate BB' for each session
0019 %  4.  calculate sx0
0020 %
0021 % --- History
0022 % 2005-08-22
0023 %   ver.30b
0024 % 2005-12-22
0025 %   bug fix (megfile --> megnoisefile)
0026 % 2008-06-27 Taku Yoshioka
0027 %   (megnoisefile -> megfile_baseline)
0028 % 2008-11-28 Taku Yoshioka
0029 %   Use vb_disp() for displaying message
0030 %
0031 % Copyright (C) 2011, ATR All Rights Reserved.
0032 % License : New BSD License(see VBMEG_LICENSE.txt)
0033 
0034 const = vb_define_verbose;
0035 VERBOSE_LEVEL_NOTICE = const.VERBOSE_LEVEL_NOTICE;
0036 
0037 % field names used in this function
0038 vb_struct2vars(noise_parm,{'noise_model','megfile_baseline','twin_noise',...
0039                     'noise_reg','temporal_filter','trial_average'}); 
0040 
0041 if isempty(temporal_filter), temporal_filter = false; end
0042 if isempty(trial_average), trial_average = false; end
0043 
0044 megfile = megfile_baseline;
0045 
0046 % MEG data file names given as a cell array
0047 if iscell(megfile),
0048   Nfile = length(megfile); 
0049 else
0050   tmp = megfile;
0051   megfile = {tmp};
0052   Nfile = 1;
0053 end;
0054 
0055 switch noise_model
0056  case 1, 
0057   vb_disp('Noise model: spherical',VERBOSE_LEVEL_NOTICE);
0058  case 2, 
0059   vb_disp('Noise model: diagonal',VERBOSE_LEVEL_NOTICE);
0060  case 3, 
0061   vb_disp('Noise model: full covariance',VERBOSE_LEVEL_NOTICE);
0062 end
0063 
0064 COV = cell(Nfile,1); % Noise covariance
0065 sx0 = 0;
0066 Ndata = 0;
0067 
0068 for n = 1:Nfile
0069   % Load MEG data
0070   vb_disp(['Load MEG data for noise estimate [' megfile{n} ']'], ...
0071           VERBOSE_LEVEL_NOTICE);
0072   bexp = vb_load_meg_data(megfile{n});
0073   MEGinfo = vb_load_meg_info(megfile{n});
0074   
0075   % MEG data averaging
0076   if trial_average, bexp = mean(bexp,3); end; 
0077   
0078   % Variable preparation
0079   Nch       = size(bexp,1);
0080   Ntrial  = size(bexp,3);
0081   Cov = zeros(Nch,Nch);       % Noise covariance matrix
0082   
0083   % Time window for noise evaluation
0084   Tnoise = twin_noise(1):twin_noise(2); 
0085   Nt = length(Tnoise);
0086   
0087   % Noise model
0088   Tstart = vb_index_to_time(Tnoise(1), MEGinfo);
0089   Tend = vb_index_to_time(Tnoise(end), MEGinfo);
0090   vb_disp(['Session ' num2str(n) ': ' num2str(Tstart,'%4.1f') ...
0091            ' - ' num2str(Tend,'%4.1f') '[ms]'],VERBOSE_LEVEL_NOTICE);     
0092   b = bexp(:,Tnoise,:);
0093   
0094   % temporal filter
0095   if temporal_filter
0096     for s = 1:Ntrial
0097       b(:,:,s) = vb_temporal_smooth(b(:,:,s));
0098     end
0099   end
0100   
0101   B = reshape(b, [Nch, Nt*Ntrial]);
0102   BB = B*B';
0103   
0104   switch noise_model
0105    case 1,
0106     Cov = eye(Nch,Nch); 
0107     invCov = eye(Nch,Nch);
0108     
0109    case 2,
0110     % covariance calculation
0111     Cov = diag(diag(BB));
0112     % normalization
0113     Cov = Cov/(trace(Cov)/Nch);
0114     invCov = pinv(Cov);
0115     
0116    case 3,
0117     % covariance calculation (+regularization)
0118     Cov = BB+(noise_reg*trace(BB)/Nch)*eye(Nch,Nch);        
0119     % normalization
0120     Cov = Cov/(trace(Cov)/Nch);
0121     invCov = pinv(Cov);
0122   end
0123   COV{n} = Cov;
0124   
0125   Ndata = Ndata + Nch*Ntrial*Nt;
0126   sx0 = sx0 + sum(sum(invCov.*(BB)));;
0127 end
0128 sx0 = sx0/Ndata;