%==================================================================
%
%	FilterEgisFile.m
%
%	Markus Junghoefer	[1996]
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      
%
%
%	Function definition
%
%
	function[SaveFilePath]=FilterEgisFile(FilePathMat,FileStatus,...
	ch_FilterFile,TimeStatusVec,ch_ChooseTimeInd,ch_ChooseChannelInd, ...
	PlotStatus,ch_AddSubZeros, ...
	ch_LowPassFilt,LowB,LowA, ...
	ch_HighPassFilt,HighB,HighA, ...
	LowFreqsMat,HighFreqsMat);
%

%==================================================================
if nargin<16; HighFreqsMat=[]; end
if nargin<15; LowFreqsMat=[]; end
if nargin<14; HighA=[]; end
if nargin<13; HighB=[]; end
if nargin<12; ch_HighPassFilt=1; end
if nargin<11; LowA=[]; end
if nargin<10; LowB=[]; end
if nargin<9; ch_LowPassFilt=1; end
if nargin<8; ch_AddSubZeros=[]; end
if nargin<7; PlotStatus=[]; end;
if nargin<6; ch_ChooseChannelInd=[]; end;
if nargin<5; ch_ChooseTimeInd=[]; end;
if nargin<4; TimeStatusVec=[]; end;
if nargin<3; ch_FilterFile=[]; end;
if nargin<2; FileStatus=[]; end;
if nargin<1; FilePathMat=[]; end

if strcmp(FilePathMat,'TestFile')
	TestStatus=1;
	ch_FilterFile=1;
	FileStatus=3;
	File=1;
	NChan=1;
	NPoints=400;
	NCells=1;
	NTrials=1;
	SampRate=500;
	LHeader=2;
	Header=1;
	OrigDataMat=zeros(NChan,NPoints);
	OrigDataMat(1,:)=[1:NPoints]+10.*randn(1,NPoints);
	OrigDataMat(1,100)=[1000];
	TimeStatusVec=[];
	ch_ChooseTimeInd=1;
else
    FileStatus=1
	TestStatus=0;
	DefFile=SetDefPath('r','*.ses*');
	[NFiles,FilePathMat]=ReadFileNames(FilePathMat,DefFile);
	if NFiles==0; return; end
    SetDefPath('w',DefFile);
end
global hPlot2dMenuFig
if ~isempty(hPlot2dMenuFig)
    hPlot2dMenuFig=[];
end

[FileStatus]=IfEmptyInputSpecVal(FileStatus,[1 2],1,...
'Please insert the kind of file:\n',...
'Original               [1]',...
'Averaged               [2]');

for FileIndex=1:NFiles
	[File,Path,FilePath]=GetFileNameOfMat(FilePathMat,FileIndex);
	%==================================================================
	if FileStatus==1
		global NChan
		global NPoints
		global NCells
		global NTrials
		global SampRate
		global LHeader
		[Header,File,Path,FilePath]=ReadEGISHeader(FilePath);
		OriginalFid=fopen(FilePath,'r','b');
		fseek(OriginalFid,LHeader,'bof');
		NTrialsInFile=sum(NTrials);
	elseif FileStatus==2	
		[OrigDataMat,File,Path,FilePath,NTrialAvgVec,StdChanTimeMat,SampRate]=ReadAvgFile(FilePath);
		[NChan,NPoints]=size(OrigDataMat);
		NCells=1;
		NTrials=1;
		NTrialsInFile=1;
	end

	if FileIndex==1
		% ========== Get time filter coeff and time filter channels  ===================
		if ch_LowPassFilt 
			if isempty(LowA) | isempty(LowB) | isempty(LowFreqsMat)
 				[ch_LowPassFilt,LowB,LowA,ch_HighPassFilt,NoUse,NoUse,LowFreqsMat,NoUse] = ...
				GetHighLowFiltCoeff(SampRate);
			end
		end
		if ch_HighPassFilt 
			if isempty(HighA) | isempty(HighB) | isempty(HighFreqsMat)
				[ch_LowPassFilt,NoUse,NoUse,ch_HighPassFilt,HighB,HighA,NoUse,HighFreqsMat] = ...
				GetHighLowFiltCoeff(SampRate);
			end
		end
		if ch_LowPassFilt | ch_HighPassFilt
			[ChanToFilt]=ChooseChanInd(NChan,ch_ChooseChannelInd,1,['Do you want to choose'],['special channels to filter ?'],0);
			NChanToFilt=length(ChanToFilt);
		else
			return;
		end
	end

	%==================================================================
	if FileStatus==1
		NoPointsInCell=NPoints.*NTrials.*NChan;
		NoPointsInAllCells=sum(NoPointsInCell);
		NoBytesInFile=NoPointsInAllCells*2+LHeader;
	end	
	
	
	if ch_HighPassFilt & ~ch_LowPassFilt
		SaveFilePath=[FilePath '.h.' num2str(HighFreqsMat(2))];
	elseif ~ch_HighPassFilt & ch_LowPassFilt
		SaveFilePath=[FilePath '.l.' num2str(LowFreqsMat(2))];
	elseif ch_HighPassFilt & ch_LowPassFilt
		SaveFilePath=[FilePath '.h.' num2str(HighFreqsMat(2)) '.l.' num2str(LowFreqsMat(2))];
	end
	if FileStatus==1
		FilterFid=fopen(SaveFilePath,'w','b');
		fwrite(FilterFid,Header,'int8');
		clear Header
	end
	TotalTrialIndex=0;
	NFrontPoints=[]; NBackPoints=[]; NZPoints=[]; 
	if FileIndex==1
		[PlotStatus]=IfEmptyInputBo('Do you want to plot the un- and filtered data ?',[],PlotStatus,0,1);
		if PlotStatus
			if FileStatus==1 
				[PlotChan]=IfEmptyInputValInt('Please choose a special channel to plot:',[],[],1,1,NChan,1);
				[MinPlotPoint]=IfEmptyInputValInt('Please insert the minimal value to be plotted:',[],[],1,1,NPoints-1,1);
				if NPoints-MinPlotPoint>1000; DefPlotPoints=1000; else; DefPlotPoints=NPoints; end
				[MaxPlotPoint]=IfEmptyInputValInt('Please insert the maximal value to be plotted:',[],[],MinPlotPoint+DefPlotPoints-1,MinPlotPoint+1,NPoints,1);
			elseif TestStatus
				if NChan>1
					[PlotChan]=IfEmptyInputValInt('Please choose a special channel to plot:',[],[],1,1,NChan,1);
				else
					PlotChan=1;
				end
				MinPlotPoint=1;
				MaxPlotPoint=NPoints;
			end
		end
	end
	for CellIndex=1:NCells
	
		[TimeToFilt,NoTimeInd,TimeStatusVec,NInt,IntStartVec,IntEndVec,ch_ChooseTimeInd]=ChooseTimeInd(TimeStatusVec,ch_ChooseTimeInd,NPoints(CellIndex));

		if NInt>1; ch_CorrIntAmp=1; end
		for TrialIndex=1:NTrials(CellIndex);
			if FileStatus==1
				OrigDataMat=fread(OriginalFid,[NChan,NPoints(CellIndex)],'int16');
			end
		
			FiltDataMat=OrigDataMat;
		
			%====================Gleiche Einzelintervalle an=====================================
			if NInt>1
				%Falls der erste Wert den Status 1 hat ist der erste Punkt das Ende des ersten Intervalles
				if IntStartVec(1)==1; IntStartVec=IntStartVec(2:NInt); end
				if IntEndVec(NInt)==NPoints(CellIndex); IntEndVec=IntEndVec(1:NInt-1); end
				IntStartVec=IntStartVec
				IntEndVec=IntEndVec
				FiltAddSubMat=FiltDataMat(:,IntStartVec)-FiltDataMat(:,IntEndVec)
				for IntIndex=1:NInt-1
					for TimeIndex=1:IntEndVec(IntIndex)
						FiltDataMat(:,TimeIndex)=FiltDataMat(:,TimeIndex)+FiltAddSubMat(:,IntIndex);
					end
				end
			end
		
			%=====================================================================================
	
					
			if (PlotStatus & FileStatus==1) | (PlotStatus & TestStatus)
				if TrialIndex==1
					hfig=figure('Units','normal','NumberTitle','off','Position',[0,.5,.5,.5],'Name','Original / Filter (--w):');			
				end		
				plot(OrigDataMat(PlotChan,MinPlotPoint:MaxPlotPoint)); 
				hold on
			end
			if NTrials(CellIndex)>1
				TotalTrialIndex=TotalTrialIndex+1;
				fprintf(1,'Filter trial %g of %g total trials !\n',TotalTrialIndex,NTrialsInFile);
			end
			%=============================================
			[TmpFiltDataMat,NFrontPoints,NBackPoints,ch_AddSubZeros,NZPoints] = AddSubZeros(1,FiltDataMat(ChanToFilt,TimeToFilt),NFrontPoints,NBackPoints,ch_AddSubZeros,NZPoints);				
			% ==========  Filter the data  ===================
			if ch_HighPassFilt
				for i=1:NChanToFilt
					TmpFiltDataMat(i,:)=(filtfilt(HighB,HighA,TmpFiltDataMat(i,:)'))';
				end
			end
			if ch_LowPassFilt
				for i=1:NChanToFilt
					TmpFiltDataMat(i,:)=(filtfilt(LowB,LowA,TmpFiltDataMat(i,:)'))';
				end
			end
			% =================================================
			[FiltDataMat(ChanToFilt,TimeToFilt),NFrontPoints,NBackPoints,ch_AddSubZeros,NZPoints] = AddSubZeros(-1,TmpFiltDataMat,NFrontPoints,NBackPoints,ch_AddSubZeros,NZPoints);
			%=============================================
		
			%====================Gleiche Einzelintervalle an=====================================
			if NInt>1
				for IntIndex=1:NInt-1
					for TimeIndex=1:IntEndVec(IntIndex)
						FiltDataMat(:,TimeIndex)=FiltDataMat(:,TimeIndex)-FiltAddSubMat(:,IntIndex);
					end
				end
			end
			%=====================================================================================
		
			if (PlotStatus & FileStatus==1) | (PlotStatus & TestStatus)
				plot(FiltDataMat(PlotChan,MinPlotPoint:MaxPlotPoint),'--w'); 
				hold off
				title(['Channel:   ' int2str(PlotChan)])
				xlabel('Time'); ylabel('Amplitude')
				axis([MinPlotPoint MaxPlotPoint .9.*min(FiltDataMat(PlotChan,MinPlotPoint:MaxPlotPoint)) 1.1.*max(FiltDataMat(PlotChan,MinPlotPoint:MaxPlotPoint))])
			end
			if FileStatus==1
				fwrite(FilterFid,FiltDataMat,'int16');
			end
		end
	end

	if FileStatus==1; fclose(OriginalFid); end

	if FileStatus==2
		SaveAvgFile(SaveFilePath,FiltDataMat, NTrialAvgVec, StdChanTimeMat,SampRate);
		if PlotStatus
			hfig=figure('Units','normal','NumberTitle','off','Position',[0,.5,.5,.5],'Name','Original / Filter (--w):');			
			for ChanInd=1:NChan
				plot(OrigDataMat(ChanInd,:)); 
				hold on
				plot(FiltDataMat(ChanInd,:),'--w'); 
				hold off
				title(['Channel:   ' int2str(ChanInd)])
				xlabel('Time'); ylabel('Amplitude')
				%axis([1 NPoints -10 10])
				pause
				clf;
			end
		end
	elseif ~TestStatus
		fclose(FilterFid);
	end
end
return;
		
