function [signal,chan_names,variance, pnts, rate, xmin, xmax]=loadavg_bcl(FILENAME,chanNameList)
	%writed by allen zhang based on EEGLAB's loadavg.m
	%2009-11-25
	% NENU,CHANGCHUN,CHINA 
	% update by Yang Zhang
	%2011-2-28
	% NENU,Changchun, CHINA

	% revised by Yang Zhang
	% 2011-4-6
	%using automatic routine to identifiy the type of the avg file to get the true nsweeps parameter 

	if ~exist('chanNameList','var')
		chanNameList={'all'};
	end
	try
		fid=fopen(FILENAME,'r','ieee-le');
		% read general part of the erp header and set variables
		% -----------------------------------------------------
		fseek(fid, 362, 'cof');% skip the firsts 362 from BOF (368 bytes in orginal comment?)
		% disp(ftell(fid));% 360 bytes
		% hdr.nsweeps = fread(fid, 1, 'ushort');
		% disp(ftell(fid));% 362
		% hdr.compsweeps = fread(fid, 1, 'ushort');% the exact sweep numbers for eeg and single subject avg file| in grand avg file it represents the number of subjects
		% hdr.acceptcnt = fread(fid, 1, 'ushort');% number of accepted sweeps also the exact sweep numbers for grand avg file
		% hdr.rejectcnt = fread(fid, 1, 'ushort');%number of rejected sweeps
		% disp(ftell(fid));% 368

		compsweeps = fread(fid, 1, 'ushort');% the exact sweeps numbers for eeg file| in Grand avg file it represented number of subjects

		acceptcnt = fread(fid, 1, 'ushort');% number of accepted sweeps 
		rejectcnt = fread(fid, 1, 'ushort');%number of rejected sweeps

		% determine the type of avg file and choose the right value for nsweeps
		if (rejectcnt+acceptcnt)~=compsweeps
			disp('It''s a grand average file!!!');
			disp(['Subject number = ',num2str(compsweeps),'; sweeps = ',num2str(acceptcnt)]);
			nsweeps=compsweeps;
		else
			disp('It''s a single subject average file!!!');
			disp(['nsweeps = ',num2str(compsweeps),'; Accepted sweeps = ',num2str(acceptcnt),'; Rejected sweeps = ',num2str(rejectcnt)]);
			nsweeps=acceptcnt;
		end


		pnts=fread(fid, 1, 'ushort');	% number of point per waveform
		chan=fread(fid, 1, 'ushort');  % number of channels
		fseek(fid, 3, 'cof');	% skip 3 bytes
		variance_flag=fread(fid, 1, 'uchar');
		rate=fread(fid, 1, 'ushort');  % sample rate (Hz)
		fseek(fid, 127, 'cof');	% skip 127 bytes
		xmin=fread(fid, 1, 'float32'); % in s
		xmax=fread(fid, 1, 'float32'); % in s
		fseek(fid, 387, 'cof');	% skip 387 bytes

		% read electrode configuration
		% ----------------------------
		for elec = 1:chan
			channel_label_tmp = fread(fid, 10, 'uchar');
			electrodes(elec).tmp=channel_label_tmp;
			chan_names(elec,:) = channel_label_tmp'; %#ok<*AGROW>
			for index = 2:9
				if chan_names(elec,index) == 0
					chan_names(elec,index)=' ';
				end;
			end;
			fseek(fid, 61, 'cof');%skip 61 bytes
			electrodes(elec).calib= fread(fid, 1, 'float32');

			%     	erp = fread(fid, 47-10, 'uchar');
			% 	baseline(elec) = fread(fid, 1, 'ushort');
			% 	erp = fread(fid, 10, 'uchar');
			% 	sensitivity(elec) = fread(fid, 1, 'float32');
			% 	erp = fread(fid, 8, 'uchar');
			% 	    electrodes(elec).calib= fread(fid, 1, 'float32');
			%     
			% 	fprintf('%s: baseline: %d\tsensitivity: %f\tcalibration: %f\n', ...
			% char(chan_names(elec,1:4)), baseline(elec), sensitivity(elec), electrodes(elec).calib);
		end;

		signal   = zeros(pnts, chan);
		variance = zeros(pnts, chan);

		for elec = 1:chan
			% To scale a data point to
			% microvolts, multiply by the channel-specific calibration factor (i.e., for electrode j:
			% channel[j]->calib) and divide by the number of sweeps in the average (i.e.,
			% channel[j]->n);
			% skip sweeps header and read data
			% --------------------------------
			fseek(fid, 5, 'cof');
			signal(:, elec) =fread(fid, pnts, 'float32')*electrodes(elec).calib/nsweeps;

		end;

		if variance_flag
			for elec = 1:chan
				variance(:, elec) = fread(fid, pnts, 'float32')*electrodes(elec).calib/nsweeps;% not sure
			end;
		else
			variance = 'novariance';
		end;
		%%
		if ~strcmpi(chanNameList{1},'all')
			chanIDX=ones(1,chan);
			for ichanList=1:numel(chanNameList)
				for elec=1:chan
					if strcmpi(chanNameList{ichanList},char(chan_names(elec,1:numel(chanNameList{ichanList}))))
						chanIDX(elec)=0;
						break;
					end
				end
			end

			signal=signal(:,~chanIDX);
			if variance_flag
				variance=variance(:,~chanIDX);
			end
			chan_names=chan_names(~chanIDX,:);
		end
		% signal = signal';
		% variance = variance';
		fclose(fid);

	catch errorLOAD
		disp(FILENAME);
		rethrow(errorLOAD);
	end
	return;