function read_epochs_noise_reduction(trigger_channels, epoch_start, epoch_end)
% 
% read_epochs_noise_reduction(trigger_channels, epoch_start, epoch_end)
%
% read epochs from SQD file, given trigger, and epoch start and end time
% and perform multi-channel-noise-reduction
% It will generate two directories: 
%  1. 'matdata' contains a set of matlab .mat file for each trigger.
%  These are original data, without noise reduction
%  2. 'matdata_nr' contains one .mat file for each trigger, after noise reduction.  
%
%
% Usage:
%
% trigger_channels, example: [165], [165:166], or [165:166, 171] 
% epoch_start: trial starting in miliseconds, according to the trigger, example:  -3600 or 0;
% epoch_end: trial end in miliseconds, according to the trigger, example:  0 or 400;
%
% exmaple: read_epochs_noise_reduction([160:166],0, 1000);
% 

%   clear;
% 	trigger_channels = [160:166];
% 	epoch_start = 0;
% 	epoch_end = 1000;
% 

    output_dir_or = 'matdata';   % original data, without noise reduction
    output_dir_nr = 'matdata_nr';  % after noise reduction
 % downsampling = 1;
[fn, pn] = uigetfile('*.sqd','SQD File Selector');

filename = [pn fn];

% make output directory for orignal data
[s,m,mid] = mkdir(pn, output_dir_or);
output_path_or = [pn output_dir_or];


[s,m,mid] = mkdir(pn, output_dir_nr);
[s,m,mid] = mkdir([pn output_dir_nr '\'], 'temp');
    
output_path_nr = [pn output_dir_nr];
output_path_temp = [pn output_dir_nr '\temp'];

info = sqdread(filename,'Info');

Fs = info.SampleRate;


trigger = find_trigger(filename, trigger_channels);

nEpoch = size(trigger,1);   %if nEpoch ~= nEpochExpected disp('warning, number of original epochs is not expected'); end;

nTrigger = size(trigger,2);

epoch_start_sampling = epoch_start * Fs/1e3;
epoch_end_sampling = epoch_end * Fs/1e3 ;
lenEpoch = epoch_end_sampling - epoch_start_sampling;

data_ref = sqdread(filename,'Channels',[157:159]);


disp('                   ');
disp('Begin noise reduction for channel 0 to channel 156...');
disp('                   ');

% read data, filtering and save it 
for  channel = 0:156
    channel
    data = sqdread(filename,'Channels',channel);
    
    data2 = multi_reference_adaptive_filter(data, data_ref);

    savefile = [output_path_temp '\temp_nr_' num2str(channel)];
 	save(savefile,'data2');
end;

disp('                   ');
disp('Read data for each trigger...');
disp('                   ');


for iTrigger = 1:nTrigger
    sprintf('Now processing trigger %g \n',  iTrigger);
    epoch = zeros(nEpoch*lenEpoch,157);

	% put triggers from same together
    con_trigger_index = zeros(lenEpoch*nEpoch,1);
	for iEpoch = 1: nEpoch
            con_trigger_index(lenEpoch*(iEpoch-1)+(1:lenEpoch)) = trigger(iEpoch, iTrigger)+ [epoch_start_sampling: epoch_end_sampling-1];
	end;
    
    
    for iChannel = 0:156
     	loadfile = [output_path_temp '\temp_nr_' num2str(iChannel)];
        load(loadfile);
        epoch(:,iChannel+1) = data2(con_trigger_index);
    end;
	savefile = [output_path_nr '\trigger' num2str(trigger_channels(iTrigger))];
	save(savefile,'epoch');
    
end;
    
[s,m] = rmdir(output_path_temp,'s');

disp('                   ');
disp('Read data for each trigger without noise reduction...');
disp('                   ');




for TriggerNo = 1:length(trigger_channels);

	
	
    trigger = find_trigger(filename, trigger_channels(TriggerNo));
    
	nEpoch = size(trigger,1);   %if nEpoch ~= nEpochExpected disp('warning, number of original epochs is not expected'); end;
	nTrigger = size(trigger,2);
	
    epoch_start_sampling = epoch_start * Fs/1e3;
	epoch_end_sampling = epoch_end * Fs/1e3 ;
	lenEpoch = epoch_end_sampling - epoch_start_sampling;
	
	epoch = zeros(nEpoch*lenEpoch,157);
	
	for iEpoch = 1: nEpoch
        epoch(lenEpoch*(iEpoch-1)+(1:lenEpoch),:) = sqdread(filename,'Channels',[0:156],'SAMPLES',[trigger(iEpoch)+1 trigger(iEpoch) + lenEpoch]);
	end;
	
	savefile = [output_path_or '\trigger' num2str(trigger_channels(TriggerNo))];
	save(savefile,'epoch');
end;
