/
ft_preproc_highpassfilter.m
311 lines (277 loc) · 11 KB
/
ft_preproc_highpassfilter.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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
function [filt, B, A] = ft_preproc_highpassfilter(dat, Fs, Fhp, N, type, dir, instabilityfix, df, wintype, dev, plotfiltresp, usefftfilt)
% FT_PREPROC_HIGHPASSFILTER applies a high-pass filter to the data and thereby removes
% the low frequency components in the data
%
% Use as
% [filt] = ft_preproc_highpassfilter(dat, Fsample, Fhp, N, type, dir, instabilityfix)
% where
% dat data matrix (Nchans X Ntime)
% Fs sampling frequency in Hz
% Fhp filter frequency in Hz
% order optional filter order, default is 6 (but) or dependent on frequency band and data length (fir/firls)
% type optional filter type, can be
% 'but' Butterworth IIR filter (default)
% 'firws' FIR filter with windowed sinc
% 'fir' FIR filter using MATLAB fir1 function
% 'firls' FIR filter using MATLAB firls function (requires MATLAB Signal Processing Toolbox)
% 'brickwall' frequency-domain filter using forward and inverse FFT
% dir optional filter direction, can be
% 'onepass' forward filter only
% 'onepass-reverse' reverse filter only, i.e. backward in time
% 'onepass-zerophase' zero-phase forward filter with delay compensation (default for firws, linear-phase symmetric FIR only)
% 'onepass-reverse-zerophase' zero-phase reverse filter with delay compensation
% 'onepass-minphase' minimum-phase converted forward filter (non-linear, only for firws)
% 'twopass' zero-phase forward and reverse filter (default, except for firws)
% 'twopass-reverse' zero-phase reverse and forward filter
% 'twopass-average' average of the twopass and the twopass-reverse
% instabilityfix optional method to deal with filter instabilities
% 'no' only detect and give error (default)
% 'reduce' reduce the filter order
% 'split' split the filter in two lower-order filters, apply sequentially
% df optional transition width (firws)
% wintype optional window type (firws), can be
% 'hamming' (default) maximum passband deviation 0.0022 [0.22%], stopband attenuation -53dB
% 'hann' maximum passband deviation 0.0063 [0.63%], stopband attenuation -44dB
% 'blackman' maximum passband deviation 0.0002 [0.02%], stopband attenuation -74dB
% 'kaiser'
% dev optional max passband deviation/stopband attenuation (only for firws with kaiser window, default = 0.001 [0.1%, -60 dB])
% plotfiltresp optional, 'yes' or 'no', plot filter responses (only for firws, default = 'no')
% usefftfilt optional, 'yes' or 'no', use fftfilt instead of filter (only for firws, default = 'no')
%
% Note that a one- or two-pass filter has consequences for the strength of the filter,
% i.e. a two-pass filter with the same filter order will attenuate the signal twice as
% strong.
%
% Further note that the filter type 'brickwall' filters in the frequency domain,
% but may have severe issues. For instance, it has the implication that the time
% domain signal is periodic. Another issue pertains to that frequencies are
% not well defined over short time intervals; particularly for low frequencies.
%
% See also PREPROC
% Copyright (c) 2003-2022, Robert Oostenveld, Arjen Stolk, Andreas Widmann,
% Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
% FieldTrip is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% FieldTrip is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$
% determine the size of the data
[nchans, nsamples] = size(dat);
% set the default filter order later
if nargin<4 || isempty(N)
N = [];
end
% set the default filter type
if nargin<5 || isempty(type)
type = 'but';
end
% set the default filter direction
if nargin<6 || isempty(dir)
if strcmp(type, 'firws')
dir = 'onepass-zerophase';
else
dir = 'twopass';
end
end
% set the default method to deal with filter instabilities
if nargin<7|| isempty(instabilityfix)
instabilityfix = 'no';
end
% Set default transition width later
if nargin < 8 || isempty(df)
df = [];
end
% Set default window type
if nargin < 9 || isempty(wintype)
wintype = 'hamming';
end
% Set default passband deviation/stopband attenuation for Kaiser window
if nargin < 10 || isempty(dev)
if strcmp(wintype, 'kaiser')
dev = 0.001;
else
dev = [];
end
end
% Set default passband deviation/stopband attenuation for Kaiser window
if nargin < 11 || isempty(plotfiltresp)
plotfiltresp = 'no';
end
% Set default filter function
if nargin < 12 || isempty(usefftfilt)
usefftfilt = false;
else
% convert to boolean value
usefftfilt = istrue(usefftfilt);
end
% Filtering does not work on integer data
if ~isa(dat, 'double') && ~isa(dat, 'single')
dat = cast(dat, 'double');
end
% preprocessing fails on channels that contain NaN
if any(isnan(dat(:)))
ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
end
% Nyquist frequency
Fn = Fs/2;
% demean the data before filtering
meandat = nanmean(dat,2);
dat = bsxfun(@minus, dat, meandat);
% compute filter coefficients
switch type
case 'but'
if isempty(N)
N = 6;
end
[B, A] = butter(N, max(Fhp)/Fn, 'high');
case 'firws'
% Input arguments
if length(Fhp) ~= 1
ft_error('One cutoff frequency required.')
end
% Filter order AND transition width set?
if ~isempty(N) && ~isempty(df)
ft_warning('firws:dfOverridesN', 'Filter order AND transition width set - transition width setting will override filter order.')
elseif isempty(N) && isempty(df) % Default transition width heuristic
df = fir_df(Fhp, Fs);
end
% Compute filter order from transition width
[foo, maxDf] = fir_df(Fhp, Fs); %#ok<ASGLU>
isOrderLow = false;
if ~isempty(df)
if df > maxDf
ft_error('Transition band too wide. Maximum transition width is %.2f Hz.', maxDf)
end
[N, dev] = firwsord(wintype, Fs, df, dev);
else % Check filter order otherwise
[df, dev] = invfirwsord(wintype, Fs, N, dev);
if df > maxDf
nOpt = firwsord(wintype, Fs, maxDf, dev);
ft_warning('firws:filterOrderLow', 'Filter order too low. For better results a minimum filter order of %d is recommended. Effective cutoff frequency might deviate from requested cutoff frequency.', nOpt)
isOrderLow = true;
end
end
% Window
if strcmp(wintype, 'kaiser')
beta = kaiserbeta(dev);
win = windows('kaiser', N + 1, beta);
else
win = windows(wintype, N + 1);
end
% Impulse response
B = firws(N, Fhp / Fn, 'high', win);
A = 1;
% Convert to minimum phase
if strcmp(dir, 'onepass-minphase')
B = minphaserceps(B);
end
% Twopass filtering
if strncmp(dir, 'twopass', 7)
pbDev = (dev + 1)^2 - 1;
sbAtt = 20 * log10(dev^2);
order = 2 * (length(B) - 1);
isTwopass = true;
else
pbDev = dev;
sbAtt = 20 * log10(dev);
order = length(B) - 1;
isTwopass = false;
end
% Reporting
ft_info once
ft_info('Highpass filtering data: %s, order %d, %s-windowed sinc FIR\n', dir, order, wintype);
if ~isTwopass && ~isOrderLow % Do not report shifted cutoffs
ft_info(' cutoff (-6 dB) %g Hz\n', Fhp);
tb = [max([Fhp - df / 2 0]), min([Fhp + df / 2 Fn])]; % Transition band edges
ft_info(' transition width %.1f Hz, stopband 0-%.1f Hz, passband %.1f-%.0f Hz\n', df, tb, Fn);
end
if ~isOrderLow
ft_info(' maximum passband deviation %.4f (%.2f%%), stopband attenuation %.0f dB\n', pbDev, pbDev * 100, sbAtt);
end
case 'fir'
if isempty(N)
N = 3*fix(Fs / Fhp);
if rem(N,2)==1, N=N+1; end
end
if N > floor( (size(dat,2) - 1) / 3)
N=floor(size(dat,2)/3) - 2;
if rem(N,2)==1, N=N+1; end
end
B = fir1(N, max(Fhp)/Fn, 'high');
A = 1;
case 'firls' % from NUTMEG's implementation
% Deprecated: see bug 2453
ft_warning('The filter type you requested is not recommended for neural signals, only proceed if you know what you are doing.')
if isempty(N)
N = 3*fix(Fs / Fhp);
if rem(N,2)==1, N=N+1; end
end
if N > floor( (size(dat,2) - 1) / 3)
N=floor(size(dat,2)/3) - 2;
if rem(N,2)==1, N=N+1; end
end
f = 0:0.001:1;
if rem(length(f),2)~=0
f(end)=[];
end
z = zeros(1,length(f));
[val,pos1] = min(abs(Fs*f/2 - Fhp));
pos2 = length(f);
z(pos1:pos2) = 1;
A = 1;
B = firls(N,f,z); % requires MATLAB signal processing toolbox
case 'brickwall'
ax = (0:(size(dat,2)-1))./(Fs/size(dat,2)); % frequency axis
a = ones(1, size(dat,2));
fbin1 = nearest(ax, Fhp);
fbin2 = nearest(ax, Fs-Fhp);
a(1:(fbin1-1)) = 0;
a((fbin2+1):end) = 0;
f = fft(dat,[],2); % FFT
f = f.*a(ones(size(dat,1),1),:); % brickwall
filt = real(ifft(f,[],2)); % iFFT
otherwise
ft_error('unsupported filter type "%s"', type);
end
% Plot filter responses
if strcmp(plotfiltresp, 'yes')
plotfresp(B, A, [], Fs, dir)
end
if ~isequal(type, 'brickwall')
try
filt = filter_with_correction(B,A,dat,dir,usefftfilt);
catch
switch instabilityfix
case 'no'
rethrow(lasterror);
case 'reduce'
ft_warning('off','backtrace');
ft_warning('filter instability detected - reducing the %dth order filter to an %dth order filter', N, N-1);
ft_warning('on','backtrace');
filt = ft_preproc_highpassfilter(dat,Fs,Fhp,N-1,type,dir,instabilityfix);
case 'split'
N1 = ceil(N/2);
N2 = floor(N/2);
ft_warning('off','backtrace');
ft_warning('filter instability detected - splitting the %dth order filter in a sequential %dth and a %dth order filter', N, N1, N2);
ft_warning('on','backtrace');
filt = ft_preproc_highpassfilter(dat ,Fs,Fhp,N1,type,dir,instabilityfix);
filt = ft_preproc_highpassfilter(filt,Fs,Fhp,N2,type,dir,instabilityfix);
otherwise
ft_error('incorrect specification of instabilityfix');
end % switch
end
end