%==================================================================
%
%	FilterStopNetFile.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 [FilterFileMat]=FilterStopNetFile(NetFileMat,StopN,StopFrequ,ch_ChooseChannelInd,ChanStatusVec,PlotStatus, ...
                                              ch_AddSubZeros,NZPoints,TRStatus,MaxMemory,NFrontPoints,NBackPoints);
%
%	

%==================================================================
if nargin<12; NBackPoints=[]; end
if nargin<11; NFrontPoints=[]; end
if nargin<10; MaxMemory=60; end		%Memory RAM in MB
if nargin<9; TRStatus=[]; end
if nargin<8; NZPoints=[]; end
if nargin<7; ch_AddSubZeros=[]; end
if nargin<6; PlotStatus=[]; end
if nargin<5; ChanStatusVec=[]; end
if nargin<4; ch_ChooseChannelInd=[]; end
if nargin<3; StopFrequ=[49 51]; end
if nargin<2; StopN=6; end
if nargin<1; NetFileMat=[]; end
 
fprintf('\n');clc;
[TRStatus]=IfEmptyInputSpecVal(TRStatus,[-1 1],1,...
'Please insert the file status:\n',...
'RAW                   [ 1]',...
'TAW                   [-1]');

if TRStatus==1; Mask='*.RAW'; Ext='.RAW'; elseif TRStatus==-1; Mask='*.TAW'; Ext='.TAW'; end
[NNetFiles,NetFileMat]=ReadFileNames(NetFileMat,Mask,'Choose NetStation file:');
if NNetFiles==0; return; end

if isempty(ChanStatusVec) | ch_ChooseChannelInd
	[ChanIndVec,NoChanIndVec,ChanStatusVec]=ChooseChanInd(NChan,[],1,['Do you want to choose'],['special channels to filter ?'],0);
end

[FilterFileMat]=FileNameExt2Mat(NetFileMat,'fs',3);
[FilterFileMat]=FileNameExt2Mat(FilterFileMat,Ext);

%========================================================================================================

for FileIndex=1:NNetFiles
	
	[NetFile,NetPath,NetFilePath]=GetFileNameOfMat(NetFileMat,FileIndex);
	
	[FilterFile,FilterPath,FilterFilePath]=GetFileNameOfMat(FilterFileMat,FileIndex);
	fprintf('Save filtered data to file:\n\n')
	disp(FilterFilePath)
	pause(.5);

	%=========NetStation Header===============

	[FilePath,NetHeader,LNetHeader,Version, ...
	Year, Month, Day, Hour, Minute, Sec, MSec, ...
	SampRate, NChan, Gain, Bits, Range, NPoints, NEvents, EventCodes] = ...
	ReadRAWHeader(NetFilePath);
	NTotChan=NChan+NEvents;
    
    Nyquist=SampRate./2;
    [BStop,AStop]=butter(StopN,StopFrequ./Nyquist,'stop');

	%==================================================================
	NetFid=fopen(NetFilePath,'r','b');
	fseek(NetFid,LNetHeader,-1);
	FilterFid=fopen(FilterFilePath,'w','b');
	fwrite(FilterFid,NetHeader,'int8');

	clear NetHeader

	[PlotStatus]=IfEmptyInputBo('Do you want to plot the un- and filtered data ?',[],PlotStatus,0,1);
	if PlotStatus
		[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);
		hfig1=figure('Units','normal','NumberTitle','off','Position',[0,.5,.5,.5],'Name','Original / Filter (--w):');					
	end
	if TRStatus==1 %RAW
		MaxNPointToRead=floor(MaxMemory.*1000000./(8.*NTotChan));
		if MaxNPointToRead>=NPoints; 
			ReadAllStatus=1; 
			DataMat=fread(NetFid,[NTotChan,NPoints],'int16');	
		else; 
			ReadAllStatus=0; 
		end
	end
	for ChanIndex=1:NTotChan
		if TRStatus==1 %RAW
			if ReadAllStatus
				DataVec=(DataMat(ChanIndex,:))';
			else
				fseek(NetFid,LNetHeader+(ChanIndex-1).*2,-1);
				DataVec=zeros(NPoints,1);
				for TimeInd=1:NPoints-1
					DataVec(TimeInd)=fread(NetFid,1,'int16');
					fseek(NetFid,(NTotChan-1).*2,0);
				end
			end
		elseif TRStatus==-1  %TAW
			DataVec=fread(NetFid,NPoints,'int16');	
		end
		if 	ChanIndex<=NChan
			if ChanStatusVec(ChanIndex)
				[DataVec,NFrontPoints,NBackPoints,ch_AddSubZeros,NZPoints] = AddSubZeros(1,DataVec',NFrontPoints,NBackPoints,ch_AddSubZeros,NZPoints);
				if ch_AddSubZeros; DataVec=DataVec'; end
				if PlotStatus
					figure(hfig1)
					plot(DataVec(MinPlotPoint:MaxPlotPoint)); 
					hold on
				else
					fprintf(1,'Filter channel %g of %g total channels !\n',ChanIndex,NChan);		
				end
				% ==========  Filter the data  ===================
				DataVec=(filtfilt(BStop,AStop,DataVec));
				% =================================================
				if PlotStatus
					plot(DataVec(MinPlotPoint:MaxPlotPoint),'--w'); 
					hold off
					title(['Channel:   ' int2str(ChanIndex)])
					xlabel('Time'); ylabel('Amplitude')
				end
				[DataVec,NFrontPoints,NBackPoints,ch_AddSubZeros,NZPoints] = AddSubZeros(-1,DataVec',NFrontPoints,NBackPoints,ch_AddSubZeros,NZPoints);	
				if ch_AddSubZeros; DataVec=DataVec'; end
			end
		end
		fwrite(FilterFid,DataVec,'int16');
	end
	fclose(NetFid);
	fclose(FilterFid);
end
return;
		
