From 75a170d1bad96ca9004ea8a657cac70d70c5577d Mon Sep 17 00:00:00 2001 From: Jan Mathijs Schoffelen Date: Thu, 16 Sep 2021 12:23:28 +0200 Subject: [PATCH] ENH - added documentation for Hilbert method, cleaned up cfg for superlets --- ft_freqanalysis.m | 91 +++++++++++++++++++++++------------- specest/ft_specest_hilbert.m | 81 ++++++++++++++++++++------------ utilities/ft_checkconfig.m | 8 ++++ 3 files changed, 117 insertions(+), 63 deletions(-) diff --git a/ft_freqanalysis.m b/ft_freqanalysis.m index 4f5558c82b..586e160105 100644 --- a/ft_freqanalysis.m +++ b/ft_freqanalysis.m @@ -127,7 +127,7 @@ % cfg.width = 'width', or number of cycles, of the wavelet (default = 7) % cfg.gwidth = determines the length of the used wavelets in standard % deviations of the implicit Gaussian kernel and should -% be choosen >= 3; (default = 3) +% be chosen >= 3; (default = 3) % % The standard deviation in the frequency domain (sf) at frequency f0 is % defined as: sf = f0/width @@ -135,21 +135,21 @@ % defined as: st = 1/(2*pi*sf) % % SUPERLET performs time-frequency analysis on any time series trial data using the -% 'wavelet method' based on a frequency-wise combination of Morlet wavelets of varying cycle +% 'superlet method' based on a frequency-wise combination of Morlet wavelets of varying cycle % widths (see Moca et al. 2019, https://doi.org/10.1101/583732). % cfg.foi = vector 1 x numfoi, frequencies of interest % OR % cfg.foilim = [begin end], frequency band of interest % cfg.toi = vector 1 x numtoi, the times on which the analysis % windows should be centered (in seconds) -% cfg.superlet.basewidth = 'width', or number of cycles, of the base wavelet (default = 3) -% cfg.superlet.gwidth = determines the length of the used wavelets in standard +% cfg.width = 'width', or number of cycles, of the base wavelet (default = 3) +% cfg.gwidth = determines the length of the used wavelets in standard % deviations of the implicit Gaussian kernel and should -% be choosen >= 3; (default = 3) -% cfg.superlet.combine = 'additive', 'multiplicative' (default = 'additive') +% be chosen >= 3; (default = 3) +% cfg.combine = 'additive', 'multiplicative' (default = 'additive') % determines if cycle numbers of wavelets comprising a superlet % are chosen additively or multiplicatively -% cfg.superlet.order = vector 1 x numfoi, superlet order, i.e. number of combined +% cfg.order = vector 1 x numfoi, superlet order, i.e. number of combined % wavelets, for individual frequencies of interest. % % The standard deviation in the frequency domain (sf) at frequency f0 is @@ -157,6 +157,24 @@ % The standard deviation in the temporal domain (st) at frequency f0 is % defined as: st = 1/(2*pi*sf) % +% HILBERT performs time-frequency analysis on any time series data using a frequency specific +% bandpass filter, followed by the Hilbert transform. +% cfg.foi = vector 1 x numfoi, frequencies of interest +% cfg.toi = vector 1 x numtoi, the time points for which the estimates will be returned (in seconds) +% cfg.width = scalar, or vector (default: 1), specifying the half bandwidth of the filter; +% cfg.edgartnan = 'no' (default) or 'yes', replace filter edges with nans, works only for finite impulse response (FIR) filters, and +% requires a user specification of the filter order +% +% For the bandpass filtering the following options can be specified, the default values are as in FT_PREPROC_BANDPASSFILTER, for more +% information see the help of FT_PREPROCESSING +% cfg.bpfilttype +% cfg.bpfiltord = (optional) scalar, or vector 1 x numfoi; +% cfg.bpfiltdir +% cfg.bpinstabilityfix +% cfg.bpfiltdf +% cfg.bpfiltwintype +% cfg.bpfiltdev +% % TFR performs time-frequency analysis on any time series trial data using the % 'wavelet method' based on Morlet wavelets. Using convolution in the time domain % instead of multiplication in the frequency domain. @@ -349,13 +367,19 @@ cfg.gwidth = ft_getopt(cfg, 'gwidth', 3); case 'superlet' - cfg.superlet = ft_getopt(cfg, 'superlet'); - cfg.superlet.basewidth = ft_getopt(cfg.superlet, 'basewidth', 3); - cfg.superlet.gwidth = ft_getopt(cfg.superlet, 'gwidth', 3); - cfg.superlet.combine = ft_getopt(cfg.superlet, 'combine', 'additive'); - cfg.superlet.order = ft_getopt(cfg.superlet, 'order', ones(1, numel(cfg.foi))); - if size(cfg.superlet.order) ~= size(cfg.foi) - ft_error('cfg.foi and cfg.superlet.order must be the same size'); + % reorganize the cfg, a nested cfg is not consistent with the othe methods + cfg = ft_checkconfig(cfg, 'createtopcfg', 'superlet'); + cfg = removefields(cfg, 'superlet'); + cfg = ft_checkconfig(cfg, 'renamed', {'basewidth', 'width'}); + + cfg.width = ft_getopt(cfg, 'width', 3); + cfg.gwidth = ft_getopt(cfg, 'gwidth', 3); + cfg.combine = ft_getopt(cfg, 'combine', 'additive'); + cfg.order = ft_getopt(cfg, 'order', ones(1, numel(cfg.foi))); + if numel(cfg.order) == 1 + cfg.order = cfg.order.*length(cfg.foi); + elseif numel(cfg.order)~= numel(cfg.foi) + ft_error('cfg.foi must have the same number of elements as cfg.foi, or must be a scalar'); end case 'tfr' @@ -365,18 +389,19 @@ cfg.gwidth = ft_getopt(cfg, 'gwidth', 3); case 'hilbert' - ft_warning('method = hilbert requires user action to deal with filtering-artifacts') + ft_warning('method = hilbert may require user action to deal with filtering-artifacts') cfg = ft_checkconfig(cfg, 'renamed', {'filttype', 'bpfilttype'}); cfg = ft_checkconfig(cfg, 'renamed', {'filtorder', 'bpfiltord'}); cfg = ft_checkconfig(cfg, 'renamed', {'filtdir', 'bpfiltdir'}); - cfg.bpfilttype = ft_getopt(cfg, 'bpfilttype', 'but'); - cfg.bpfiltord = ft_getopt(cfg, 'bpfiltord', 4); - cfg.bpfiltdir = ft_getopt(cfg, 'bpfiltdir', 'twopass'); - cfg.bpinstabilityfix = ft_getopt(cfg, 'bpinstabilityfix', 'no'); + cfg.bpfilttype = ft_getopt(cfg, 'bpfilttype'); + cfg.bpfiltord = ft_getopt(cfg, 'bpfiltord'); + cfg.bpfiltdir = ft_getopt(cfg, 'bpfiltdir'); + cfg.bpinstabilityfix = ft_getopt(cfg, 'bpinstabilityfix'); cfg.bpfiltdf = ft_getopt(cfg, 'bpfiltdf'); cfg.bpfiltwintype = ft_getopt(cfg, 'bpfiltwintype'); cfg.bpfiltdev = ft_getopt(cfg, 'bpfiltdev'); cfg.width = ft_getopt(cfg, 'width', 1); + cfg.edgartnan = istrue(ft_getopt(cfg, 'edgeartnan', 'no')); fn = fieldnames(cfg); bpfiltoptions = ft_cfg2keyval(keepfields(cfg, fn(startsWith(fn, 'bp')))); @@ -623,18 +648,18 @@ case 'superlet' % calculate number of wavelets and respective cycle width dependent on superlet order % equivalent one-liners: - % multiplicative: cycles = arrayfun(@(order) arrayfun(@(wl_num) cfg.superlet.basewidth*wl_num, 1:order), cfg.superlet.order,'uni',0) - % additive: cycles = arrayfun(@(order) arrayfun(@(wl_num) cfg.superlet.basewidth+wl_num-1, 1:order), cfg.superlet.order,'uni',0) + % multiplicative: cycles = arrayfun(@(order) arrayfun(@(wl_num) cfg.width*wl_num, 1:order), cfg.order,'uni',0) + % additive: cycles = arrayfun(@(order) arrayfun(@(wl_num) cfg.width+wl_num-1, 1:order), cfg.order,'uni',0) cycles = cell(length(cfg.foi),1); for i_f = 1:length(cfg.foi) - frq_cyc = NaN(1,cfg.superlet.order(i_f)); - if strcmp(cfg.superlet.combine, 'multiplicative') - for i_wl = 1:cfg.superlet.order(i_f) - frq_cyc(i_wl) = cfg.superlet.basewidth*i_wl; + frq_cyc = NaN(1,cfg.order(i_f)); + if strcmp(cfg.combine, 'multiplicative') + for i_wl = 1:cfg.order(i_f) + frq_cyc(i_wl) = cfg.width*i_wl; end - elseif strcmp(cfg.superlet.combine, 'additive') - for i_wl = 1:cfg.superlet.order(i_f) - frq_cyc(i_wl) = cfg.superlet.basewidth+i_wl-1; + elseif strcmp(cfg.combine, 'additive') + for i_wl = 1:cfg.order(i_f) + frq_cyc(i_wl) = cfg.width+i_wl-1; end end cycles{i_f} = frq_cyc; @@ -647,15 +672,15 @@ foi = options{idx_freqoi}; for i_f = 1:length(cfg.foi) % collext individual wavelets' responses per frequency - spec_f = NaN(cfg.superlet.order(i_f), nchan, length(cfg.toi)); + spec_f = NaN(cfg.order(i_f), nchan, length(cfg.toi)); opt = options; opt{idx_freqoi} = cfg.foi(i_f); % compute responses for individual wavelets - for i_wl = 1:cfg.superlet.order(i_f) - [spec_f(i_wl,:,:), dum, toi] = ft_specest_wavelet(dat, time, 'timeoi', cfg.toi, 'width', cycles{i_f}(i_wl), 'gwidth', cfg.superlet.gwidth, opt{:}, 'feedback', fbopt); + for i_wl = 1:cfg.order(i_f) + [spec_f(i_wl,:,:), dum, toi] = ft_specest_wavelet(dat, time, 'timeoi', cfg.toi, 'width', cycles{i_f}(i_wl), 'gwidth', cfg.gwidth, opt{:}, 'feedback', fbopt); end % geometric mean across individual wavelets - spectrum(:,i_f,:) = prod(spec_f, 1).^(1/cfg.superlet.order(i_f)); + spectrum(:,i_f,:) = prod(spec_f, 1).^(1/cfg.order(i_f)); end clear spec_f @@ -685,7 +710,7 @@ spectrum = reshape(spectrum,[1 nchan numel(foi) numel(toi)]); case 'hilbert' - [spectrum,foi,toi] = ft_specest_hilbert(dat, time, 'timeoi', cfg.toi, 'width', cfg.width, bpfiltoptions{:}, options{:}, 'feedback', fbopt); + [spectrum,foi,toi] = ft_specest_hilbert(dat, time, 'timeoi', cfg.toi, 'width', cfg.width, bpfiltoptions{:}, options{:}, 'feedback', fbopt, 'edgeartnan', cfg.edgeartnan); hastime = true; % create FAKE ntaper (this requires very minimal code change below for compatibility with the other specest functions) ntaper = ones(1,numel(foi)); diff --git a/specest/ft_specest_hilbert.m b/specest/ft_specest_hilbert.m index 2f3ac184d0..ad5522e44b 100644 --- a/specest/ft_specest_hilbert.m +++ b/specest/ft_specest_hilbert.m @@ -52,19 +52,19 @@ freqoi = ft_getopt(varargin, 'freqoi'); timeoi = ft_getopt(varargin, 'timeoi', 'all'); width = ft_getopt(varargin, 'width', 1); -bpfilttype = ft_getopt(varargin, 'bpfilttype'); if isempty(bpfilttype), ft_error('you need to specify filter type'), end -bpfiltord = ft_getopt(varargin, 'bpfiltord'); if isempty(bpfiltord), ft_error('you need to specify filter order'), end -bpfiltdir = ft_getopt(varargin, 'bpfiltdir'); if isempty(bpfiltdir), ft_error('you need to specify filter direction'), end +bpfilttype = ft_getopt(varargin, 'bpfilttype'); +bpfiltord = ft_getopt(varargin, 'bpfiltord'); +bpfiltdir = ft_getopt(varargin, 'bpfiltdir'); bpinstabilityfix = ft_getopt(varargin, 'bpinstabilityfix', 'no'); bpfiltdf = ft_getopt(varargin, 'bpfiltdf'); bpfiltwintype = ft_getopt(varargin, 'bpfiltwintype'); bpfiltdev = ft_getopt(varargin, 'bpfiltdev'); pad = ft_getopt(varargin, 'pad'); -padtype = ft_getopt(varargin, 'padtype', 'zero'); -polyorder = ft_getopt(varargin, 'polyorder', 0); +padtype = ft_getopt(varargin, 'padtype', 'zero'); +polyorder = ft_getopt(varargin, 'polyorder', 0); fbopt = ft_getopt(varargin, 'feedback'); -verbose = ft_getopt(varargin, 'verbose', true); +verbose = ft_getopt(varargin, 'verbose', true); edgeartnan = ft_getopt(varargin, 'edgeartnan', false); if isempty(fbopt) @@ -72,13 +72,50 @@ fbopt.n = 1; end +% set the default filter type +if isempty(bpfilttype) + bpfilttype = 'but'; +end + +% set the default filter direction +if isempty(bpfiltdir) + if strcmp(bpfilttype, 'firws') + bpfiltdir = 'onepass-zerophase'; + else + bpfiltdir = 'twopass'; + end +end + % edge artifacts from filters are only exactly defined for FIR filters if edgeartnan && ~any(strcmp(bpfilttype, {'firws','fir','firls'})) - ft_error('edge artifacts are only exactly defined, and can only be replaced by NaNs, for FIR filters') + ft_error('edge artifacts are only exactly defined, and can only be replaced by NaNs, for FIR filters'); +elseif edgeartnan + if isempty(bpfiltord) + ft_error('the filter order for the FIR filters needs to be specified for edge artifact replacement'); + end end % Set n's -[nchan,ndatsample] = size(dat); +[nchan, ndatsample] = size(dat); + +% Determine fsample and set total time-length of data +fsample = 1./mean(diff(time)); +dattime = ndatsample./fsample; % total time in seconds of input data + +% Set a default sampling for the frequencies-of-interest if not defined +if isempty(freqoi) + freqoi = linspace(2*width, (fsample/3), 50); +end +% check for freqoi = 0 and remove it +if any(freqoi==0) + freqoi(freqoi==0) = []; +end +nfreqoi = length(freqoi); + +% expand width to array if constant width +if numel(width) == 1 + width = ones(1,nfreqoi) * width; +end % This does not work on integer data if ~isa(dat, 'double') && ~isa(dat, 'single') @@ -90,11 +127,6 @@ dat = ft_preproc_polyremoval(dat, polyorder, 1, ndatsample); end -% Determine fsample and set total time-length of data -fsample = 1./mean(diff(time)); -dattime = ndatsample / fsample; % total time in seconds of input data - - % Zero padding if round(pad * fsample) < ndatsample ft_error('the padding that you specified is shorter than the data'); @@ -105,17 +137,6 @@ prepad = floor(((pad - dattime) * fsample) ./ 2); postpad = ceil(((pad - dattime) * fsample) ./ 2); - -% set a default sampling for the frequencies-of-interest -if isempty(freqoi) - freqoi = linspace(2*width, (fsample/3), 50); -end -% check for freqoi = 0 and remove it -if any(freqoi==0) - freqoi(freqoi==0) = []; -end -nfreqoi = length(freqoi); - % Set timeboi and timeoi timeoiinput = timeoi; offset = round(time(1)*fsample); @@ -140,11 +161,6 @@ end end -% expand width to array if constant width -if numel(width) == 1 - width = ones(1,nfreqoi) * width; -end - % expand filter order to array if constant filterorder if numel(bpfiltord) == 1 bpfiltord = ones(1,nfreqoi) * bpfiltord; @@ -179,7 +195,12 @@ end % filter - flt = ft_preproc_bandpassfilter(dat, fsample, filtfreq(ifreqoi,:), bpfiltord(ifreqoi), bpfilttype, bpfiltdir, bpinstabilityfix, bpfiltdf, bpfiltwintype, bpfiltdev); + if isempty(bpfiltord) + ibpfiltord = []; + else + ibpfiltord = bpfiltord(ifreqoi); + end + flt = ft_preproc_bandpassfilter(dat, fsample, filtfreq(ifreqoi,:), ibpfiltord, bpfilttype, bpfiltdir, bpinstabilityfix, bpfiltdf, bpfiltwintype, bpfiltdev); % transform dum = transpose(hilbert(transpose(ft_preproc_padding(flt, padtype, prepad, postpad)))); diff --git a/utilities/ft_checkconfig.m b/utilities/ft_checkconfig.m index 72e01a96d4..7a6777f7fe 100644 --- a/utilities/ft_checkconfig.m +++ b/utilities/ft_checkconfig.m @@ -311,6 +311,14 @@ 'normalizeparam' 'weight' }; + + case {'superlet'} + fieldname = { + 'basewidth' + 'gwidth' + 'combine' + 'order' + }; otherwise ft_error('unexpected name of the subfunction');