1 | function filtered_data_out = fmri_filter(data,TR,varargin) |
---|
2 | |
---|
3 | % Defaults |
---|
4 | hipass = []; |
---|
5 | lowpass = []; |
---|
6 | detrending = false; |
---|
7 | smoothpr = false; |
---|
8 | reserve_mean = true; |
---|
9 | out_file = ''; |
---|
10 | |
---|
11 | % Parse varargin |
---|
12 | for ii=1:2:length(varargin) |
---|
13 | param = varargin{ii}; |
---|
14 | value = varargin{ii+1}; |
---|
15 | switch lower(param) |
---|
16 | case 'hipass' |
---|
17 | hipass = value; |
---|
18 | case 'lowpass' |
---|
19 | lowpass = value; |
---|
20 | case 'detrending' |
---|
21 | if strcmpi(value,'on') |
---|
22 | detrending = true; |
---|
23 | end |
---|
24 | case 'smoothpr' |
---|
25 | if strcmpi(value,'on') |
---|
26 | smoothpr = true; |
---|
27 | end |
---|
28 | case 'reservemean' |
---|
29 | if strcmpi(value,'off') |
---|
30 | reserve_mean = false; |
---|
31 | end |
---|
32 | case 'filename' |
---|
33 | out_file = value; |
---|
34 | otherwise |
---|
35 | error('Unknown parameter "%s"',param) |
---|
36 | end |
---|
37 | end |
---|
38 | |
---|
39 | if isstruct(data) |
---|
40 | data = data.FTDATA; |
---|
41 | elseif ischar(data) |
---|
42 | data = aedes_data_read(data); |
---|
43 | data = data.FTDATA; |
---|
44 | end |
---|
45 | filtered_data = double(data); |
---|
46 | |
---|
47 | if isempty(hipass) && isempty(lowpass) && ~detrending |
---|
48 | return |
---|
49 | end |
---|
50 | |
---|
51 | % Construct hi-pass filter |
---|
52 | if ~isempty(hipass) |
---|
53 | [B_hipass,A_hipass]=butter(5,hipass/((1/TR)/2),'high'); |
---|
54 | %[bz_h,bp_h,bk_h]=butter(15,hipass/((1/TR)/2),'high'); |
---|
55 | %[sos_var_high,g_high] = zp2sos(bz_h, bp_h, bk_h); |
---|
56 | %Hd_high = dfilt.df2sos(sos_var_high, g_high); |
---|
57 | end |
---|
58 | |
---|
59 | % Construct low-pass filter |
---|
60 | if ~isempty(lowpass) |
---|
61 | [B_lowpass,A_lowpass]=butter(10,lowpass/((1/TR)/2),'low'); |
---|
62 | %[bz_l,bp_l,bk_l]=butter(15,hipass/((1/TR)/2),'low'); |
---|
63 | %[sos_var_low,g_low] = zp2sos(bz_l, bp_l, bk_l); |
---|
64 | %Hd_low = dfilt.df2sos(sos_var_low, g_low); |
---|
65 | end |
---|
66 | |
---|
67 | |
---|
68 | % Filter image data ----------------------------- |
---|
69 | if ~isempty(lowpass) || ~isempty(hipass) || detrending |
---|
70 | |
---|
71 | do_lp = true; |
---|
72 | do_hp = true; |
---|
73 | if isempty(lowpass) |
---|
74 | do_lp = false; |
---|
75 | end |
---|
76 | |
---|
77 | if isempty(hipass) |
---|
78 | do_hp = false; |
---|
79 | end |
---|
80 | |
---|
81 | % Waitbar |
---|
82 | wbh=aedes_wbar(0,'Filtering image data...'); |
---|
83 | |
---|
84 | counter=1; |
---|
85 | sz = ones(1,4); |
---|
86 | sz(1:length(size(data)))=size(data); |
---|
87 | nVox = prod(sz(1:3)); |
---|
88 | for ii=1:size(data,1) |
---|
89 | for kk=1:size(data,2) |
---|
90 | for zz=1:size(data,3) |
---|
91 | tmp = squeeze(double(data(ii,kk,zz,:))); |
---|
92 | mean_tmp = mean(tmp); |
---|
93 | if do_hp |
---|
94 | tmp = filtfilt(B_hipass,A_hipass,tmp); |
---|
95 | %tmp = filter(Hd_high,tmp); |
---|
96 | if reserve_mean |
---|
97 | tmp = tmp+mean_tmp; |
---|
98 | end |
---|
99 | end |
---|
100 | if do_lp |
---|
101 | tmp = filtfilt(B_lowpass,A_lowpass,tmp); |
---|
102 | %tmp = filter(Hd_low,tmp); |
---|
103 | end |
---|
104 | if detrending |
---|
105 | tmp = detrend(tmp,'linear'); |
---|
106 | if reserve_mean |
---|
107 | tmp = tmp+mean_tmp; |
---|
108 | end |
---|
109 | end |
---|
110 | filtered_data(ii,kk,zz,:) = tmp; |
---|
111 | aedes_wbar(counter/nVox,wbh) |
---|
112 | counter=counter+1; |
---|
113 | end |
---|
114 | end |
---|
115 | end |
---|
116 | close(wbh) |
---|
117 | end |
---|
118 | |
---|
119 | if nargout>0 |
---|
120 | filtered_data_out = filtered_data; |
---|
121 | end |
---|
122 | |
---|
123 | % Write output to file |
---|
124 | if ~isempty(out_file) |
---|
125 | aedes_write_nifti(filtered_data,out_file); |
---|
126 | end |
---|
127 | |
---|