Skip to content

Commit

Permalink
ENH - added documentation for Hilbert method, cleaned up cfg for supe…
Browse files Browse the repository at this point in the history 
…rlets
  • Loading branch information
@schoffelen
schoffelen committed Sep 16, 2021
1 parent 68c25fe commit 75a170d
Show file tree
Hide file tree
Showing 3 changed files with 117 additions and 63 deletions.
91 changes: 58 additions & 33 deletions ft_freqanalysis.m
View file
Expand Up @@ -127,36 +127,54 @@
% 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
% The standard deviation in the temporal domain (st) at frequency f0 is
% 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
% defined as: sf = f0/width
% 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.
Expand Down Expand Up @@ -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'
Expand All @@ -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'))));
Expand Down Expand Up @@ -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;
Expand All @@ -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

Expand Down Expand Up @@ -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));
Expand Down
81 changes: 51 additions & 30 deletions specest/ft_specest_hilbert.m
View file
Expand Up @@ -52,33 +52,70 @@
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)
fbopt.i = 1;
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')
Expand All @@ -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');
Expand All @@ -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);
Expand All @@ -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;
Expand Down Expand Up @@ -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))));
Expand Down
8 changes: 8 additions & 0 deletions utilities/ft_checkconfig.m
View file
Expand Up @@ -311,6 +311,14 @@
'normalizeparam'
'weight'
};

case {'superlet'}
fieldname = {
'basewidth'
'gwidth'
'combine'
'order'
};

otherwise
ft_error('unexpected name of the subfunction');
Expand Down

0 comments on commit 75a170d

Please sign in to comment.