function	VBfreq = prepare_online_vb_filter(proj_root,VBfilt,curr_parm)
%
%  VBfreq = prepare_online_vb_filter(VBfilt,Fs,T,freq_range)
%
% VBfilt.KW = VB inverse filter  (Njact,Nch,Nwindow,Ntask)
% --- Input
% curr_parm.filterfile
% curr_parm.megfile
% curr_parm.twin_meg
% curr_parm.Tperiod
% curr_parm.Tnext
% curr_parm.freq_range : band frequency for feature extraction [Nband x 2]
%           freq_range(n ,:)  = [freq_low freq_high] Hz (for n-th band)
% --- output variable
% 
% VBfreq.KW     = VB inverse filter  (Njact,Nch)
% VBfreq.fxlist = index list for band frequency
% VBfreq.Wfreq  = Sparse matrix to get sum of frequency component in each band
% VBfreq.NFFT   =  number of time points to calculate FFT
% VBfreq.Twindow = Time window index;
% VBfreq.freq_range = freq_range;
% 
% 
% 
% Masa-aki Sato 2009-2-6

%%%%%% load VBMEG estimated result
%filterfile = [proj_root '/' curr_parm.filterfile];  % not necessary? mod yu 090330
%load(filterfile, 'VBfilt');       % mod yu 090330

megfile = [proj_root, curr_parm.megfile{1}];
%MEGinfo = BMI_load_eeginfo(megfile);
load(megfile,'MEGinfo');

% Fs : sampling frequency [Hz]
Fs = MEGinfo.SampleFreq;%MEGinfo.SampleFrequency;

if ~isfield(curr_parm,'freq_range'), 
	freq_range = [1 30];
else
	freq_range = curr_parm.freq_range
end;
if ~isfield(curr_parm,'fft_window'), 
	fft_window = 0;
else
	fft_window = curr_parm.fft_window;
end

[Njact,Nch,Nw,Ntask] = size(VBfilt.KW);

if Nw==1 && Ntask==1, 
	KW = VBfilt.KW(:,:,1,1) ;
else
	fprintf('(Nwindow=%d ,Ntask=%d)\n',Nw,Ntask)
	fprintf('Take average of VB filters\n')
	KW = squeeze(mean( mean( VBfilt.KW,3) ,4));
end

Tsample = curr_parm.twin_meg(2) - curr_parm.twin_meg(1) + 1;
Twindow = calc_timewindow(Tsample, curr_parm.Tperiod, curr_parm.Tnext);
Tperiod = Twindow(:,2) - Twindow(:,1) + 1;

T = max(Tperiod);
NFFT = 2^nextpow2(T);

% frequency for FFT
df=Fs/NFFT;
freq=0:df:df*(NFFT/2);

if max(freq_range(:)) > Fs/2, 
	error('max freq is larger than Nyquist freq')
end

% band frequency
[Nband ,D] = size(freq_range)
if D ~=2, error('freq_range should be [Nband x 2]'); end;

% index list for band frequency
flist  = zeros(Nband ,2);

for n=1:Nband*2
	[tmp, ix] = min(abs(freq - freq_range(n)));
	flist(n) = ix;
end

fxlist = [];
ixlist = [];
nxlist = [];
Nlast  = 0;
%Ilist  = zeros(Nband,2);

for n=1:Nband
	Nlist  = flist(n,2) - flist(n,1) + 1;
	Nstart = Nlast + 1;
	Nlast  = Nlast + Nlist;
	
	fxlist = [fxlist, flist(n,1):flist(n,2)];
	ixlist = [ixlist, Nstart:Nlast];
	nxlist = [nxlist, repmat(n, [1, Nlist])];
	
%	Ilist(n,:) = [Nstart Nlast];
end

% Sparse matrix to get sum of frequency component in each band
Nall = length(ixlist);
Wfreq = sparse( ixlist , nxlist , 1 , Nall, Nband) ;

% Set parameters for VB-filter to calculate frequency feature
VBfreq.KW     = KW;
VBfreq.fxlist = fxlist;
VBfreq.Wfreq  = Wfreq;
VBfreq.NFFT   = NFFT;
VBfreq.Twindow = Twindow;
VBfreq.freq_range = freq_range;

% Bf = fft( B, NFFT, 2);
% Zf = Bf(:,fxlist) ).^2 * Wfreq 
% Zf(:,n) : n-th band power 
%           summed over [flist(n,1):flist(n,2)] in freq_range(n)

% Zf = vb_feature_calc(Bt,KW,fxlist,Wfreq,NFFT);

