source: misclib/fmri_filter.m @ 94

Last change on this file since 94 was 80, checked in by tjniskan, 10 years ago
  • Changed the historical "an2_" prefix to "aedes_" in all files. NOTE:

Any script or function relying to Aedes functions will be broken
because of this. Just do a search/replace from "an2_" to "aedes_" in
your files and all should be well...

  • Changed the name of an2_readtab.m to a more informative

aedes_readphasetable.m

File size: 2.5 KB
Line 
1function filtered_data_out = fmri_filter(data,TR,varargin)
2
3% Defaults
4hipass = [];
5lowpass = [];
6detrending = false;
7smoothpr = false;
8reserve_mean = true;
9out_file = '';
10
11% Parse varargin
12for 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
37end
38
39if isstruct(data)
40  data = data.FTDATA;
41elseif ischar(data)
42  data = aedes_data_read(data);
43  data = data.FTDATA;
44end
45filtered_data = double(data);
46
47if isempty(hipass) && isempty(lowpass) && ~detrend
48  return
49end
50
51% Construct hi-pass filter
52if ~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);
57end
58
59% Construct low-pass filter
60if ~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);
65end
66
67
68% Filter image data -----------------------------
69if ~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)
117end
118
119if nargout>0
120  filtered_data_out = filtered_data;
121end
122
123% Write output to file
124if ~isempty(out_file)
125  aedes_write_nifti(filtered_data,out_file);
126end
127
128
129
130
131
Note: See TracBrowser for help on using the repository browser.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008