function infodr_WaveletALFF(AllVolume,outfile_WaveletALFF,infile_Mask,Parameter)
%-----------------------------------------------------------
%   Copyright(c) 2015
%	Center for Cognition and Brain Disorders, Hangzhou Normal University, Hangzhou 310015, China
%	Written by JIA Xi-Ze 201412
%	http://www.restfmri.net/
% 	Mail to Authors: jxz.rest@gmail.com, jiaxize@foxmail.com
%   Wavelet-ALFF was added by SUN Jia-Wei   20_04_20

ASamplePeriod=Parameter.SamplePeriod;
ALowPass_HighCutoff=Parameter.LowPass_HighCutoff;
AHighPass_LowCutoff=Parameter.HighPass_LowCutoff;
CUTNUMBER=Parameter.CutNumber;

rp_WaveletALFF(AllVolume,...
    ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, ...
    infile_Mask,outfile_WaveletALFF,...
    '', '', '', CUTNUMBER);

end

function [WaveletALFFBrain,Header] = rp_WaveletALFF(AllVolume,ASamplePeriod, HighCutoff, LowCutoff, AMaskFilename, AResultFilename, TemporalMask, ScrubbingMethod, Header, CUTNUMBER)
% Use continuous wavelet method to compute and return a WaveletALFF brain map which reflects the "energy" of the voxels' BOLD signal
% Ref: Luo, F.F., Wang, J.B., Yuan, L.X., Zhou, Z.W., Xu, H., Ma, S.H., Zang, Y.F., Zhang, M.,2020. Higher sensitivity and reproducibility of wavelet-based amplitude of resting-state fMRI. Front. Neurosci., doi: 10.3389/fnins.2020.00224 .
% FORMAT    [WaveletALFFBrain,Header] = WaveletALFF(AllVolume,ASamplePeriod, HighCutoff, LowCutoff, AMaskFilename, AResultFilename, TemporalMask, ScrubbingMethod, Header, CUTNUMBER)
% Input:
% 	AllVolume		-	4D data matrix (DimX*DimY*DimZ*DimTimePoints) or the directory of 3D image data file or the filename of one 4D data file
% 	ASamplePeriod		TR, or like the variable name
% 	LowCutoff			the low edge of the pass band
% 	HighCutoff			the High edge of the pass band
% 	AMaskFilename		the mask file name, I only compute the point within the mask
%	AResultFilename		the output filename for WaveletALFF
%                       or   string: name for WaveletALFF.
%   TemporalMask    -   Temporal mask for scrubbing (DimTimePoints*1)
%                   -   Empty (blank: '' or []) means do not need scrube. Then ScrubbingMethod can leave blank
%   ScrubbingMethod -   The methods for scrubbing.
%                       -1. 'cut': discarding the timepoints with TemporalMask == 0
%                       -2. 'nearest': interpolate the timepoints with TemporalMask == 0 by Nearest neighbor interpolation 
%                       -3. 'linear': interpolate the timepoints with TemporalMask == 0 by Linear interpolation
%                       -4. 'spline': interpolate the timepoints with TemporalMask == 0 by Cubic spline interpolation
%                       -5. 'pchip': interpolate the timepoints with TemporalMask == 0 by Piecewise cubic Hermite interpolation
%   Header          -   If AllVolume is given as a 4D Brain matrix, then Header should be designated.
%   CUTNUMBER           Cut the data into pieces if small RAM memory e.g. 4GB is available on PC. It can be set to 1 on server with big memory (e.g., 50GB).
%                       default: 10
% Output:
%	WaveletALFFBrain       -   The WaveletALFF results
%   Header          -   The NIfTI Header
%	AResultFilename	- the filename of WaveletALFF results.
%-----------------------------------------------------------
%   Algorithm originally Written by LUO Fei-Fei (luofeifei@stu.xjtu.edu.cn) on 20200312 based on the FFT-ALFF re-written by YAN Chao-Gan (ycg.yan@gmail.com) on 120328 .
%	http://restfmri.net 

if ~exist('CUTNUMBER','var')
    CUTNUMBER = 10;
end

theElapsedTime =cputime;

fprintf('\nComputing WaveletALFF...');

if ~isnumeric(AllVolume)
    [AllVolume,VoxelSize,theImgFileList, Header] =rp_to4d(AllVolume);
end

[nDim1 nDim2 nDim3 nDimTimePoints]=size(AllVolume);
BrainSize = [nDim1 nDim2 nDim3];
VoxelSize = sqrt(sum(Header.mat(1:3,1:3).^2));

fprintf('\nLoad mask "%s".\n', AMaskFilename);
MaskData = rp_loadmask(nDim1, nDim2, nDim3, AMaskFilename);
MaskData = logical(MaskData);

% Convert into 2D
AllVolume=reshape(AllVolume,[],nDimTimePoints)';

MaskDataOneDim=reshape(MaskData,1,[]);
AllVolume=AllVolume(:,find(MaskDataOneDim));

% Scrubbing
if exist('TemporalMask','var') && ~isempty(TemporalMask)
    if ~all(TemporalMask)
        fprintf('\n\t Scrubbing...');
        AllVolume = AllVolume(find(TemporalMask),:); %'cut'
        if ~strcmpi(ScrubbingMethod,'cut')
            xi=1:length(TemporalMask);
            x=xi(find(TemporalMask));
            AllVolume = interp1(x,AllVolume,xi,ScrubbingMethod);
        end
        nDimTimePoints = size(AllVolume,1);
    end
end


%Wavelet parameter setting
wavename='db2';   %the optimal mother wavelet proved in this paper
totalscale=64;    % the commonly used parameter 
Fc=centfrq(wavename);  %wavelet center frequency
c=2*Fc*totalscale;
scales=c./(1:totalscale);  %scale sequence
% f=scal2frq(scales,wavename,1/fs); % corresponding frequency sequence % fs-sampling frequency
%64 frequency bins in the current study, between 0 and fs/2 Hz at an interval of (fs/2)/64 Hz

% Get the frequency index
sampleFreq = 1/ASamplePeriod;
if LowCutoff<HighCutoff && HighCutoff<=sampleFreq/2 && LowCutoff>=0
    if (LowCutoff == 0)
        idx_LowCutoff = 1;
    else % high cut off, such as freq > 0.01 Hz
        idx_LowCutoff = ceil(LowCutoff * totalscale * ASamplePeriod *2);
        % Change from round to ceil: idx_LowCutoff = round(LowCutoff *paddedLength *ASamplePeriod + 1);
    end
    if (HighCutoff==sampleFreq/2)
        idx_HighCutoff = totalscale;
    else % Low pass, such as freq < 0.08 Hz
        idx_HighCutoff = fix(HighCutoff * totalscale * ASamplePeriod * 2);
        % Change from round to fix: idx_HighCutoff	=round(HighCutoff *paddedLength *ASamplePeriod + 1);
    end
else
    error('Frequency band range input error.');
end

% Detrend before wavelet transform as did in the previous alff.m for fft
% Cut to be friendly with the RAM Memory
SegmentLength = ceil(size(AllVolume,2) / CUTNUMBER);
for iCut=1:CUTNUMBER
    if iCut~=CUTNUMBER
        Segment = (iCut-1)*SegmentLength+1 : iCut*SegmentLength;
    else
        Segment = (iCut-1)*SegmentLength+1 : size(AllVolume,2);
    end
    AllVolume(:,Segment) = detrend(AllVolume(:,Segment));  % Detrend before wavelet transform 

end


% Perform wavelet transform and calculate WaveletALFF
fprintf('\n\t Performing WaveletALFF ...');
for SNumber=1:size(AllVolume,2)
    AllVolume(1:totalscale,SNumber) = sum(abs(cwt(AllVolume(:,SNumber),scales,wavename)),2);
%     fprintf('.'); % cut off by SUN Jia-Wei 200506
end
WaveletALFF_2D = mean(AllVolume(idx_LowCutoff:idx_HighCutoff,:));


% Get the 3D brain back
WaveletALFFBrain = zeros(size(MaskDataOneDim));
WaveletALFFBrain(1,find(MaskDataOneDim)) = WaveletALFF_2D;
WaveletALFFBrain = reshape(WaveletALFFBrain,nDim1, nDim2, nDim3);


Header.pinfo = [1;0;0];
Header.dt    =[16,0];

%Save WaveletALFF image to disk
fprintf('\n\t Saving WaveletALFF map.\tWait...');
rp_writefile(single(WaveletALFFBrain), ...
    AResultFilename, ...
    BrainSize,VoxelSize,Header, 'single');


theElapsedTime = cputime - theElapsedTime;
fprintf('\n\t WaveletALFF compution over, elapsed time: %g seconds.\n', theElapsedTime);

end





