vb_feature_calc

 Calculate bandpass power of estimated current
 --- Usage
   T = size(Bt,2);
   VBfreq = prepare_online_vb_filter(VBfilt,Fs,T,freq_range)

   Zpw = vb_feature_calc(Bt, VBfreq);
 --- Input
 Bt = MEG/EEG data (Nch,Tsample)
 KW = VB inverse filter  (Njact,Nch)

 freq_range : band frequency for feature extraction [Nband x 2]
 freq_range(n ,:)  = [freq_low freq_high] Hz (for n-th band)
 Fs     : sampling frequency [Hz]
 Wfreq  : Sparse matrix to get sum of frequency component in each band
 fxlist : index list for bandpass components in FFT-signal
 NFFT   : number of time points to calculate FFT
 --- Output
 Zf(:,n) : n-th band power summed over [freq_range(n,1) freq_range(n,2)]
 Bt(:,t) : signal at time index t (t=1:Nt)
 Bf = fft( vb_repmultiply(Bt, hanning(Nt)') ,NFFT, 2);
 Zf = Bf(:,fxlist) ).^2 * Wfreq 
 Zf = current of Njact vertices (Njact,)

 Masa-aki Sato 2008-12-2

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

0001 function    Zpw = vb_feature_calc(B, VBfreq)
0002 % Calculate bandpass power of estimated current
0003 % --- Usage
0004 %   T = size(Bt,2);
0005 %   VBfreq = prepare_online_vb_filter(VBfilt,Fs,T,freq_range)
0006 %
0007 %   Zpw = vb_feature_calc(Bt, VBfreq);
0008 % --- Input
0009 % Bt = MEG/EEG data (Nch,Tsample)
0010 % KW = VB inverse filter  (Njact,Nch)
0011 %
0012 % freq_range : band frequency for feature extraction [Nband x 2]
0013 % freq_range(n ,:)  = [freq_low freq_high] Hz (for n-th band)
0014 % Fs     : sampling frequency [Hz]
0015 % Wfreq  : Sparse matrix to get sum of frequency component in each band
0016 % fxlist : index list for bandpass components in FFT-signal
0017 % NFFT   : number of time points to calculate FFT
0018 % --- Output
0019 % Zf(:,n) : n-th band power summed over [freq_range(n,1) freq_range(n,2)]
0020 % Bt(:,t) : signal at time index t (t=1:Nt)
0021 % Bf = fft( vb_repmultiply(Bt, hanning(Nt)') ,NFFT, 2);
0022 % Zf = Bf(:,fxlist) ).^2 * Wfreq
0023 % Zf = current of Njact vertices (Njact,)
0024 %
0025 % Masa-aki Sato 2008-12-2
0026 %
0027 % Copyright (C) 2011, ATR All Rights Reserved.
0028 % License : New BSD License(see VBMEG_LICENSE.txt)
0029 
0030 % Current
0031 
0032 Twindow = VBfreq.Twindow;
0033 Nwindow = size(Twindow,1);
0034 Nact    = size(VBfreq.KW,1);
0035 Nband   = size(VBfreq.Wfreq,2);
0036 
0037 Zpw = zeros(Nact,Nband,Nwindow);
0038 
0039 if isfield(VBfreq,'fft_window') 
0040     fft_window = VBfreq.fft_window;
0041 else
0042     fft_window = 0;
0043 end
0044 
0045 for n=1:Nwindow
0046     % FFT of MEG/EEG data with Hanning window
0047     Bt = B(:,Twindow(n,1):Twindow(n,2));
0048     
0049     Nt = size(Bt,2);
0050     
0051     if fft_window==1,
0052         Bt = vb_repmultiply(Bt, hanning(Nt)');
0053     end
0054     
0055     % FFT
0056     Bf  = fft( Bt , VBfreq.NFFT, 2);
0057     
0058     % select frequency components
0059     Bfx = Bf(:, VBfreq.fxlist);
0060     
0061     % Current reconstruction in frequency space
0062     Zf = abs( VBfreq.KW * Bfx ).^2 ;
0063     
0064     % sum over frequency bands
0065     Zpw(:,:,n) = Zf * VBfreq.Wfreq;
0066 end
0067 
0068