function ActualBandFrequency=nic_BandPassFilter(Path,Dimension,ScanNum,TR,FreqBand,Filefilter,PrefixOUT,Retrend)
% Bandpass filter on the post-processed images obtained from the directory
% (Path). The directory may contain a large number of subdirectory.
% However, the subdirectory and images cann't be put in the main directory
% (Path) together.
%
% FORMAT [ActualHighpass ActualLowpass]=nic_BandPassFilter(Path,Dimension,ScanNum,TR,FreqBand,Filefilter,PrefixOUT,Retrend)
%
% Input:
% 1.Path - the main file directory, also allow to contain subdirectory
% 2.Dimension - the image Dimensionension e.g., [61 73 61]
% 3.ScanNum - the number of scan (time points)
% 4.TR - acquisition time. e.g., [2] pointing to 2s
% 5.FreqBand - filtered frequency band. e.g.,[0.01 0.08]
% 6.Filefilter - file filter set. e.g., ['w'] or ['sw'] meaning that only
%               process the image involving the 'w' or 'sw' as the prefix
% 7.PrefixOUT - arranged as the prefix of the filtered images
% 8.Retrend - whether to add linear trend to filtered time series. 
%            'YES': added; 'NO': not added
%
% Output:
% 1.ActualHighpass - the actually adopted high pass frequency due to
%                    computered error
% 2.ActualLowpass - the actually adopted low pass frequency due to
%                    computered error
%_______________________________________________________________________
% Copyright (C) 2007 Neuroimage Computing Group, State Key Laboratory of
% Cognitive Neuroscience and Learning
%
% Liang Wang, wanglbit@gmail.com
% @(#)nic_BandPassFilter.m  ver 2.0, 07/11/29
CurrentPath=pwd;
cd(Path);
spmPath=nic_getSPMPath('nifti');
addpath(spmPath);
Datafolder=dir;
Datafolder=Datafolder(3:end);
% Determine whether the directory is the deepest catalog
if isdir(Datafolder(1).name)
    nsubjects=size(Datafolder,1);
    Seconddir=1;
else
    Seconddir=0;
    nsubjects=1;
end

%Arrange the border of filtered frequency band
LowPass=FreqBand(2);
HighPass=FreqBand(1);
sampleFreq 	 = 1/TR;
sampleLength = ScanNum;
paddedLength = nic_rest_nextpow2_one35(sampleLength); 
freqPrecision= sampleFreq/paddedLength;

%-------------------------------------------------------------
%| Generate low- and high- pass filter mask (refer to
%| rest_bandpass.m written by Song Xiaowei)
%-------------------------------------------------------------
% Generate low-pass filter mask
MaskData=ones(Dimension);
MaskData=logical(MaskData);
maskLowPass =repmat(MaskData, [1, 1, 1, paddedLength]);
maskHighPass=maskLowPass;
if (LowPass>=(sampleFreq/2))||(LowPass==0)
    maskLowPass(:,:,:,:)=1;	%All pass
    ActualLowpass = sampleFreq/2;
elseif (LowPass>0)&&(LowPass< freqPrecision)
    maskLowPass(:,:,:,:)=0;	% All stop
    ActualLowpass = 0;
else
    % Low pass, e.g., freq < 0.08 Hz
    idxCutoff	=round(LowPass *paddedLength *TR);
    ActualLowpass = idxCutoff./(TR*paddedLength);
    idxCutoff2	=paddedLength+2 -idxCutoff;	
    maskLowPass(:,:,:,idxCutoff+1:idxCutoff2-1)=0; 
end
   
% Generate high-pass filter mask
if (HighPass < freqPrecision)
    maskHighPass(:,:,:,:)=1;	%All pass
    ActualHighpass=0;
elseif (HighPass >= (sampleFreq/2))
    maskHighPass(:,:,:,:)=0;	% All stop
    ActualHighpass = sampleFreq/2;
else
    % high pass, e.g., freq > 0.01 Hz
    idxCutoff	=round(HighPass *paddedLength *TR);
    ActualHighpass = idxCutoff./(TR*paddedLength);
    idxCutoff2	=paddedLength+2 -idxCutoff;	
    maskHighPass(:,:,:,1:idxCutoff-1)=0;
    maskHighPass(:,:,:,idxCutoff2+1:paddedLength)=0;
end
% filtering the images each subject
for iFolder=1:nsubjects
    if Seconddir ==1
    workpath=[Path '\' Datafolder(iFolder).name];
        cd(workpath);
    end
    % get file identifiers
    FileString=nic_spm_select('List',pwd, ['^' Filefilter '.*\.img$']);
    FileFid=nic_spm_vol(FileString);
    nimgo = size(FileString,1);
    nslices= FileFid(1).dim(3);
    nslices_t= FileFid(1).dim(3);
    if ( nslices_t ~= nslices )
        fprintf('Number of slices differ! %d %\n', nimg);
    else
        % create new header files for filter
        Vout 	= FileFid;
        for k=1:nimgo,
            [pth,nm,xt,vr] = fileparts(deblank(FileFid(k).fname));
            Vout(k).fname  = fullfile(workpath,[PrefixOUT '_' nm xt vr]);
            if isfield(Vout(k),'descrip'),
                desc = [Vout(k).descrip ' '];
            else,
                desc = '';
            end;
        end;        
        Vout = nic_spm_create_vol(Vout);
        slices = zeros([Vout(1).dim(1:2) nimgo]);
        for k = 1:nslices
            B  = nic_spm_matrix([0 0 k]);
            for m=1:nimgo
                slices(:,:,m) = nic_spm_slice_vol(FileFid(m),B,FileFid(1).dim(1:2),1);
            end
            %Save the linear trend
            theTrend_Intercept=slices(:,:,1);
            theTrend_Slope= (slices(:,:,end) -theTrend_Intercept)/(sampleLength-1);
            for y=1:sampleLength
                %remove the linear trend first
                slices(:,:,y)=slices(:,:, y) -(theTrend_Intercept + y*theTrend_Slope);
            end
            %FFT
            slicesfreq =fft(slices, paddedLength, 3);
            %Mask redundant frequency components
            FilterMask= squeeze(maskLowPass(:,:,k,:));
            slicesfreq(~FilterMask)=0;
            FilterMask=squeeze(maskHighPass(:,:,k,:));
            slicesfreq(~FilterMask)=0;
            %inverse FFT
            FilteredSlices =ifft(slicesfreq, paddedLength, 3);
            FilteredSlices =FilteredSlices(:,:,1:sampleLength);%remove the padded parts
            %retrend the time course
            if strcmpi(Retrend, 'yes')
                for y=1:nimgo
                    %add the linear trend after filter
                    FilteredSlices(:,:, y)=FilteredSlices(:,:, y) + (theTrend_Intercept+y*theTrend_Slope);
                end
            end
            % write out the slice for all volumes
            for p = 1:nimgo,
                Vout(p) = nic_spm_write_plane(Vout(p),FilteredSlices(:,:,p),k);
            end;
        end
        Vout = nic_spm_close_vol(Vout);
    end
end
cd(CurrentPath);
ActualBandFrequency=[ActualHighpass ActualLowpass];
return



