-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathf_alff.m
executable file
·150 lines (126 loc) · 6.5 KB
/
f_alff.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
function [fALFFBrain, Header] = f_alff(AllVolume,ASamplePeriod, HighCutoff, LowCutoff, AMaskFilename,AResultFilename, TemporalMask, ScrubbingMethod, Header, CUTNUMBER)
% Use fALFF method to compute the brain and return a fALFF brain map.
% 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.
% FORMAT y_f_alff(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
% 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:
% fALFFBrain - The fALFF results
% Header - The NIfTI Header
% AResultFilename the filename of ALFF result
%-----------------------------------------------------------
% Algorithm (ALFF) originally Written by Xiao-Wei Song ([email protected]).
% Revised to calculate (fALFF) by CHENG Wen-Lian.
% Algorithm Re-Written by YAN Chao-Gan ([email protected]) on 120328.
% Note: the fALFF generated by the new version is slightly different from
% the orignial version. When calculate the sum of full band amplitude,
% the new version treat the component at Nyquist Frequency as the same as
% other frequency component (2*abs(FFT_AtNyquistFrequency)/N), while the
% old version calculate as sqrt(abs(FFT_AtNyquistFrequency)^2/N).
% http://restfmri.net
if ~exist('CUTNUMBER','var')
CUTNUMBER = 10;
end
theElapsedTime =cputime;
fprintf('\nComputing fALFF...');
if ~isnumeric(AllVolume)
[AllVolume,VoxelSize,theImgFileList, Header] =rest_to4d(AllVolume);
end
[nDim1 nDim2 nDim3 nDimTimePoints]=size(AllVolume);
BrainSize = [nDim1 nDim2 nDim3];
VoxelSize = sqrt(sum(Header.mat(1:3,1:3).^2));
fprintf('\n\t Load mask "%s".', AMaskFilename);
MaskData = rest_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
% Get the frequency index
sampleFreq = 1/ASamplePeriod;
sampleLength = nDimTimePoints;
paddedLength = rest_nextpow2_one35(sampleLength); %2^nextpow2(sampleLength);
if (LowCutoff >= sampleFreq/2) % All high included
idx_LowCutoff = paddedLength/2 + 1;
else % high cut off, such as freq > 0.01 Hz
idx_LowCutoff = ceil(LowCutoff * paddedLength * ASamplePeriod + 1);
% Change from round to ceil: idx_LowCutoff = round(LowCutoff *paddedLength *ASamplePeriod + 1);
end
if (HighCutoff>=sampleFreq/2)||(HighCutoff==0) % All low pass
idx_HighCutoff = paddedLength/2 + 1;
else % Low pass, such as freq < 0.08 Hz
idx_HighCutoff = fix(HighCutoff *paddedLength *ASamplePeriod + 1);
% Change from round to fix: idx_HighCutoff =round(HighCutoff *paddedLength *ASamplePeriod + 1);
end
% Detrend before fft as did in the previous f_alff.m
%AllVolume=detrend(AllVolume);
% 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));
end
% Zero Padding
AllVolume = [AllVolume;zeros(paddedLength -sampleLength,size(AllVolume,2))]; %padded with zero
fprintf('\n\t Performing FFT ...');
%AllVolume = 2*abs(fft(AllVolume))/sampleLength;
% 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) = 2*abs(fft(AllVolume(:,Segment)))/sampleLength;
fprintf('.');
end
% Generate fALFF
fALFF_2D = sum(AllVolume(idx_LowCutoff:idx_HighCutoff,:)) ./ sum(AllVolume(2:(paddedLength/2 + 1),:));
fALFF_2D(~isfinite(fALFF_2D))=0;
% Get the 3D brain back
fALFFBrain = zeros(size(MaskDataOneDim));
fALFFBrain(1,find(MaskDataOneDim)) = fALFF_2D;
fALFFBrain = reshape(fALFFBrain,nDim1, nDim2, nDim3);
Header.pinfo = [1;0;0];
Header.dt =[16,0];
%Save fALFF image to disk
fprintf('\n\t Saving fALFF map.\tWait...');
rest_writefile(single(fALFFBrain), ...
AResultFilename, ...
BrainSize,VoxelSize,Header, 'single');
theElapsedTime = cputime - theElapsedTime;
fprintf('\n\t fALFF compution over, elapsed time: %g seconds.\n', theElapsedTime);