f_alff

 Use fALFF method to compute the brain and return a fALFF brain map.
 FORMAT    function [] = f_alff(ADataDir,ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, AMaskFilename,AResultFilename)
 Ref: Zou QH, Zhu CZ, Yang Y, Zuo XN, Long XY, Cao QJ, Wang YF, Zang YF (2008) An improved approach to detection of amplitude of low-frequency fluctuation (ALFF) for resting-state fMRI: fractional ALFF. Journal of neuroscience methods 172:137-141.
 Input:
     ADataDir            where the 3d+time dataset stay, and there should be 3d EPI functional image files. It must not contain / or \ at the end.
     ASamplePeriod        TR, or like the variable name
     AHighPass_LowCutoff            the low edge of the pass band
     ALowPass_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
 Output:
    AResultFilename    the filename of fALFF result
-----------------------------------------------------------
    Copyright(c) 2007~2010
    State Key Laboratory of Cognitive Neuroscience and Learning in Beijing Normal University
    Written by Xiao-Wei Song 
    http://resting-fmri.sourceforge.net
-----------------------------------------------------------
     Mail to Authors:  <a href="Dawnwei.Song@gmail.com">Xiaowei Song</a>; <a href="ycg.yan@gmail.com">Chaogan Yan</a> 
    Version=1.3;
    Release=20090321;
   Revised by YAN Chao-Gan, 080610. NIFTI compatible
   Last Revised by YAN Chao-Gan, 090321. (1) Added the fALFF module, thank Dr. CHENG Wen-Lian for the helpful work. (2) Result data will be saved in the format 'single'.

0001 function [] = f_alff(ADataDir,ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, AMaskFilename,AResultFilename)
0002 % Use fALFF method to compute the brain and return a fALFF brain map.
0003 % FORMAT    function [] = f_alff(ADataDir,ASamplePeriod, ALowPass_HighCutoff, AHighPass_LowCutoff, AMaskFilename,AResultFilename)
0004 % Ref: Zou QH, Zhu CZ, Yang Y, Zuo XN, Long XY, Cao QJ, Wang YF, Zang YF (2008) An improved approach to detection of amplitude of low-frequency fluctuation (ALFF) for resting-state fMRI: fractional ALFF. Journal of neuroscience methods 172:137-141.
0005 % Input:
0006 %     ADataDir            where the 3d+time dataset stay, and there should be 3d EPI functional image files. It must not contain / or \ at the end.
0007 %     ASamplePeriod        TR, or like the variable name
0008 %     AHighPass_LowCutoff            the low edge of the pass band
0009 %     ALowPass_HighCutoff            the High edge of the pass band
0010 %     AMaskFilename        the mask file name, I only compute the point within the mask
0011 %    AResultFilename        the output filename
0012 % Output:
0013 %    AResultFilename    the filename of fALFF result
0014 %-----------------------------------------------------------
0015 %    Copyright(c) 2007~2010
0016 %    State Key Laboratory of Cognitive Neuroscience and Learning in Beijing Normal University
0017 %    Written by Xiao-Wei Song
0018 %    http://resting-fmri.sourceforge.net
0019 %-----------------------------------------------------------
0020 %     Mail to Authors:  <a href="Dawnwei.Song@gmail.com">Xiaowei Song</a>; <a href="ycg.yan@gmail.com">Chaogan Yan</a>
0021 %    Version=1.3;
0022 %    Release=20090321;
0023 %   Revised by YAN Chao-Gan, 080610. NIFTI compatible
0024 %   Last Revised by YAN Chao-Gan, 090321. (1) Added the fALFF module, thank Dr. CHENG Wen-Lian for the helpful work. (2) Result data will be saved in the format 'single'.
0025 
0026 
0027 
0028     if nargin~=6
0029         error(' Error using ==> f_alff. 6 arguments wanted.'); 
0030     end
0031         
0032     theElapsedTime =cputime;
0033     fprintf('\nComputing fALFF with:\t"%s"', ADataDir);
0034     [AllVolume,vsize,theImgFileList, Header] =rest_to4d(ADataDir);
0035     
0036     % examin the dimensions of the functional images and set mask
0037     [nDim1 nDim2 nDim3 nDim4]=size(AllVolume);
0038     %nDim1 = size(AllVolume,1); nDim2 = size(AllVolume,2); nDim3 = size(AllVolume,3);
0039     isize = [nDim1 nDim2 nDim3]; 
0040     
0041     %20070512    Saving a big 3D+time Dataset to small pieces by its first dimension to make this process run at least
0042     % put pieces of 4D dataset to the temp dir determined by the current time
0043     theTempDatasetDirName =sprintf('fALFF_%d', fix((1e4) *rem(now, 1) ));    
0044     theTempDatasetDir =[tempdir theTempDatasetDirName] ;
0045     ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0046     mkdir(tempdir, theTempDatasetDirName);    %Matlab 6.5 compatible
0047     
0048     Save1stDimPieces(theTempDatasetDir, AllVolume, 'dim1_');
0049     clear AllVolume;%Free large memory
0050 %mask selection, added by Xiaowei Song, 20070421
0051     fprintf('\n\t Load mask "%s".', AMaskFilename);    
0052     mask=rest_loadmask(nDim1, nDim2, nDim3, AMaskFilename);
0053     
0054     fprintf('\n\t Build fALFF filtered mask.\tWait...');
0055     sampleFreq      = 1/ASamplePeriod; 
0056     sampleLength = size(theImgFileList,1);
0057     paddedLength = rest_nextpow2_one35(sampleLength); %2^nextpow2(sampleLength);
0058     freqPrecision= sampleFreq/paddedLength;
0059         
0060     mask =logical(mask);%Revise the mask to ensure that it contain only 0 and 1
0061     maskLowPass =    repmat(mask, [1, 1, 1, paddedLength]);
0062     maskHighPass=    maskLowPass;
0063     clear mask;
0064     
0065     %20071226 fALFF parameters
0066     idxLowPass_HighCutoff    =round(ALowPass_HighCutoff *paddedLength *ASamplePeriod);
0067     %% GENERATE LOW PASS WINDOW    20070514, reference: fourior_filter.c in AFNI
0068     if (ALowPass_HighCutoff>=sampleFreq/2)||(ALowPass_HighCutoff==0)        
0069         maskLowPass(:,:,:,:)=1;    %All pass
0070     elseif (ALowPass_HighCutoff>0)&&(ALowPass_HighCutoff< freqPrecision)        
0071         maskLowPass(:,:,:,:)=0;    % All stop
0072     else
0073         % Low pass, such as freq < 0.08 Hz
0074         idxLowPass_HighCutoff    =round(ALowPass_HighCutoff *paddedLength *ASamplePeriod);
0075         idxLowPass_HighCutoff2    =paddedLength+2 -idxLowPass_HighCutoff;                %Center Index =(paddedLength/2 +1)
0076         %maskLowPass(:,:,:,1:idxCutoff)=1;            %Low pass, contain DC
0077         maskLowPass(:,:,:,idxLowPass_HighCutoff+1:idxLowPass_HighCutoff2-1)=0; %High eliminate
0078         %maskLowPass(:,:,:,idxCutoff2:paddedLength)=1;    %Low pass
0079     end
0080     
0081     %fALFF parameters 20071226
0082     idxHighPass_LowCutoff    =round(AHighPass_LowCutoff *paddedLength *ASamplePeriod);
0083     %%GENERATE HIGH PASS WINDOW
0084     if (AHighPass_LowCutoff < freqPrecision)    
0085         maskHighPass(:,:,:,:)=1;    %All pass
0086     elseif (AHighPass_LowCutoff >= sampleFreq/2)
0087         maskHighPass(:,:,:,:)=0;    % All stop
0088     else
0089         % high pass, such as freq > 0.01 Hz
0090         idxHighPass_LowCutoff    =round(AHighPass_LowCutoff *paddedLength *ASamplePeriod);
0091         idxHighPass_LowCutoff2    =paddedLength+2 -idxHighPass_LowCutoff;                %Center Index =(paddedLength/2 +1)
0092         maskHighPass(:,:,:,1:idxHighPass_LowCutoff-1)=0;    %Low eliminate
0093         %maskHighPass(:,:,:,idxCutoff:idxCutoff2)=1;    %High Pass
0094         maskHighPass(:,:,:,idxHighPass_LowCutoff2+1:paddedLength)=0;    %Low eliminate
0095     end    
0096     %Combine the low pass mask and the high pass mask
0097     %(~maskHighPass)=0;    %Don't combine because filter will not work when I only want low-pass or high-pass after combination, 20070517
0098     %Save mask pieces to disk to make this program at least run
0099     Save1stDimPieces(theTempDatasetDir, maskLowPass, 'fmLow_');    
0100     Save1stDimPieces(theTempDatasetDir, maskHighPass, 'fmHigh_');    
0101     clear maskLowPass maskHighPass; %Free large memory
0102 
0103     %20070513    remove trend --> FFT --> filter --> inverse FFT --> retrend
0104     if rest_misc( 'GetMatlabVersion')>=7.3
0105         fftw('dwisdom');
0106     end    
0107     fprintf('\n\t fALFF computing.\tWait...');        
0108     NumPieces_Dim1 =4;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0109     NumComputingCount =floor(nDim1/NumPieces_Dim1);
0110     if NumComputingCount< (nDim1/NumPieces_Dim1),
0111         NumComputingCount =NumComputingCount +1;
0112     else
0113     end
0114     for x=1:(NumComputingCount),
0115 %    for x=1:(floor(nDim1/NumPieces_Dim1) +1)
0116         rest_waitbar(x/(floor(nDim1/NumPieces_Dim1) +1), ...
0117                     'Computing fALFF. Please wait...', ...
0118                     'REST working','Child','NeedCancelBtn');
0119                     
0120         %%Remove the linear trend first, ref: fourier_filter.c in AFNI, 20070509
0121         %Get every slope and intercept within the mask
0122         theFilename =fullfile(theTempDatasetDir, sprintf('dim1_%.8d', x));
0123         theDim1Volume4D =Load1stDimVolume(theFilename);
0124         theDim1Volume4D =double(theDim1Volume4D);
0125                 
0126         %Save the linear trend
0127         % theTrend_Intercept=theDim1Volume4D(:,:,:, 1);
0128         % theTrend_Slope= (theDim1Volume4D(:,:,:, end) -theTrend_Intercept) /double(sampleLength-1);
0129         % for y=1:sampleLength % Does fALFF need this routine?
0130             % remove the linear trend first
0131             % theDim1Volume4D(:,:,:, y)=theDim1Volume4D(:,:,:, y) -(theTrend_Intercept + y*theTrend_Slope);
0132         % end
0133         
0134         %I must Detrend it first, 20070703
0135         % for xx=1:size(theDim1Volume4D,1), for yy=1:size(theDim1Volume4D,2),
0136             % dim3PlusTimeCourse =squeeze( theDim1Volume4D(xx, yy, :, :) );
0137             % detrend only support 2-D operations at most, so I have to do like this
0138             % dim3PlusTimeCourse =detrend(dim3PlusTimeCourse');
0139             % I didn't add the mean back this time, 20070703
0140             % theDim1Volume4D(xx, yy, :, :) =dim3PlusTimeCourse';
0141         % end;end;
0142         %20071110
0143         for xx=1:size(theDim1Volume4D,1),
0144             oneAxialSlice =double(theDim1Volume4D(xx, :, :, :));
0145             oneAxialSlice =reshape(oneAxialSlice, 1*nDim2*nDim3, nDim4)';
0146             oneAxialSlice =detrend(oneAxialSlice);% +repmat(mean(oneAxialSlice), [size(oneAxialSlice,1), 1]);
0147             oneAxialSlice =reshape(oneAxialSlice', 1,nDim2,nDim3, nDim4);
0148             theDim1Volume4D(xx, :, :, :) =(oneAxialSlice);
0149         end;
0150                 
0151         theDim1Volume4D =cat(4,theDim1Volume4D,zeros(size(theDim1Volume4D,1),nDim2,nDim3,paddedLength -sampleLength));    %padded with zero
0152         
0153         %FFT
0154         theDim1Volume4D =fft(theDim1Volume4D, [], 4);
0155         theDim1Volume4D_without_filter=theDim1Volume4D;% Save the FFT (without filter)  Revised by YAN Chao-Gan, 20090321. Thank Dr. CHENG Wen-Lian for the helpful work.
0156         %Low-pass Filter mask
0157         theFilename =fullfile(theTempDatasetDir, sprintf('fmLow_%.8d', x));
0158         theDim1FilterMask4D =Load1stDimVolume(theFilename);            
0159         %Apply the filter Low Pass
0160         theDim1Volume4D(~theDim1FilterMask4D)=0;
0161         
0162         %High-pass Filter mask
0163         theFilename =fullfile(theTempDatasetDir, sprintf('fmHigh_%.8d', x));
0164         theDim1FilterMask4D =Load1stDimVolume(theFilename);            
0165         %Apply the filter High Pass
0166         theDim1Volume4D(~theDim1FilterMask4D)=0;
0167         
0168         %Get the amplitude only in one of the symmetric sides after FFT
0169         theDim1Volume4D =theDim1Volume4D(:, :, :, 2:(paddedLength/2+1));
0170         theDim1Volume4D =abs(theDim1Volume4D);
0171         %Get the Power Spectrum, double it because I only want one side of both sides
0172         theDim1Volume4D =2*(theDim1Volume4D .* theDim1Volume4D) /sampleLength;
0173         theDim1Volume4D(:,:,:, end) =theDim1Volume4D(:,:,:, end) /2;
0174         %The DC component didn't double because it didn't have its symetric side
0175         %theDim1Volume4D(:,:,:,1) =theDim1Volume4D(:,:,:,1)/2;
0176         %Get the Square root of the power spectrum between 0.01 and 0.08, i.e., fALFF
0177         theDim1Volume4D =sqrt(theDim1Volume4D);
0178         theDim1Volume4D =sum(theDim1Volume4D,4);
0179         
0180         %Compute the fALFF. Revised by YAN Chao-Gan, 20090321. Thank Dr. CHENG Wen-Lian for the helpful work.
0181         %The square root was calculated at each frequency of the power spectrum. The sum of amplitude across 0.01¨C0.08 Hz was divided by that across the entire frequency range, i.e., 0¨C0.25 Hz.
0182         %Get the amplitude only in one of the symmetric sides after FFT
0183         theDim1Volume4D_without_filter =theDim1Volume4D_without_filter(:, :, :, 2:(paddedLength/2+1));
0184         theDim1Volume4D_without_filter =abs(theDim1Volume4D_without_filter);
0185         %Get the Power Spectrum, double it because I only want one side of both sides
0186         theDim1Volume4D_without_filter =2*(theDim1Volume4D_without_filter .* theDim1Volume4D_without_filter) /sampleLength;
0187         theDim1Volume4D_without_filter(:,:,:, end) =theDim1Volume4D_without_filter(:,:,:, end) /2;
0188         %The DC component didn't double because it didn't have its symetric side
0189         %theDim1Volume4D_without_filter(:,:,:,1) =theDim1Volume4D_without_filter(:,:,:,1)/2;
0190         %Get the Square root of the power spectrum between 0.01 and 0.08, i.e., fALFF
0191         theDim1Volume4D_without_filter =sqrt(theDim1Volume4D_without_filter);
0192         theDim1Volume4D_without_filter =sum(theDim1Volume4D_without_filter,4);
0193         
0194         temp_pos=find(theDim1Volume4D_without_filter);
0195         theDim1Volume4D_fALFF=zeros(size(theDim1Volume4D_without_filter));
0196         theDim1Volume4D_fALFF(temp_pos) =theDim1Volume4D(temp_pos) ./ theDim1Volume4D_without_filter(temp_pos);
0197         theDim1Volume4D=theDim1Volume4D_fALFF;                
0198         
0199         %Save to file
0200         theFilename =fullfile(theTempDatasetDir, sprintf('result_%.8d', x));        
0201         save(theFilename, 'theDim1Volume4D');         
0202         fprintf('.');
0203     end
0204     clear theDim1Volume4D theTrend_Intercept theTrend_Slope theDim1FilterMask4D oneAxialSlice theDim1Volume4D_fALFF theDim1Volume4D_without_filter;
0205     
0206     %Construct the 3D+time Dataset from files again
0207     fprintf('\n\t ReConstructing 3D Dataset fALFF.\tWait...');
0208     theDataset3D=zeros(nDim1, nDim2, nDim3);
0209     for x=1:(NumComputingCount)
0210         rest_waitbar(x/(floor(nDim1/NumPieces_Dim1)+1), ...
0211                     'fALFF 3D Brain reconstructing. Please wait...', ...
0212                     'REST working','Child','NeedCancelBtn');
0213         
0214         theFilename =fullfile(theTempDatasetDir, sprintf('result_%.8d', x));
0215         %fprintf('\t%d',x);% Just for debugging
0216         if x~=(floor(nDim1/NumPieces_Dim1)+1)
0217             theDataset3D(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1),:,:)=Load1stDimVolume(theFilename);
0218         else
0219             theDataset3D(((x-1)*NumPieces_Dim1+1):end,:,:)=Load1stDimVolume(theFilename);
0220         end        
0221         fprintf('.');
0222     end
0223     
0224     %Save fALFF image to disk
0225     fprintf('\n\t Saving fALFF map.\tWait...');    
0226     rest_writefile(single(theDataset3D), ...
0227         AResultFilename, ...
0228         isize,vsize,Header, 'single'); %Revised by YAN Chao-Gan, 090321. Result data will be stored in 'single' format. %'double');
0229 
0230     theElapsedTime =cputime - theElapsedTime;
0231     fprintf('\n\t fALFF compution over, elapsed time: %g seconds.\n', theElapsedTime);
0232 
0233     %After Band pass filter, remove the temporary files
0234     ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
0235 %end
0236 
0237 %Save the 1st dimension of the 4D dataset to files
0238 function Save1stDimPieces(ATempDir, A4DVolume, AFilenamePrefix)
0239     NumPieces_Dim1 =4;    %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
0240     NumComputingCount =floor(size(A4DVolume,1)/NumPieces_Dim1);
0241     if NumComputingCount< (size(A4DVolume,1)/NumPieces_Dim1),
0242         NumComputingCount =NumComputingCount +1;
0243     else
0244     end
0245     for x = 1:(NumComputingCount),
0246     %for x = 1:(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
0247         rest_waitbar(x/(floor(size(A4DVolume,1)/NumPieces_Dim1)+1), ...
0248                     'Cut one Big 3D+time Dataset into pieces of 3D+time Dataset Before fALFF. Please wait...', ...
0249                     'REST working','Child','NeedCancelBtn');
0250                     
0251         theFilename =fullfile(ATempDir, sprintf('%s%.8d',AFilenamePrefix, x));
0252         if x~=(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
0253             the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1), :,:,:);
0254         else
0255             the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):end, :,:,:);
0256         end
0257         save(theFilename, 'the1stDim');         
0258     end    
0259 
0260 %Load the 1st dimension of the 4D dataset from files, return a Matrix not a struct
0261 function Result=Load1stDimVolume(AFilename)    
0262     Result =load(AFilename);
0263     theFieldnames=fieldnames(Result);    
0264     % Result =eval(sprintf('Result.%s',the1stField));%remove the struct variable to any named variable with a matrix
0265     Result = Result.(theFieldnames{1});