123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253 |
- function [Dtf, Dtph] = spm_eeg_tf(S)
- % Compute instantaneous power and phase in peri-stimulus time and frequency
- % FORMAT [Dtf, Dtph] = spm_eeg_tf(S)
- %
- % S - input structure
- %
- % fields of S:
- % S.D - MEEG object or filename of M/EEG mat-file with
- %
- % S.channels - cell array of channel names. Can include generic
- % wildcards: 'All', 'EEG', 'MEG' etc.
- %
- % S.frequencies - vector of frequencies of interest
- %
- % S.timewin - time window of interest in PST in ms.
- %
- % S.method - name for the spectral estimation to use. This
- % corresponds to the name of a plug-in function that comes
- % after 'spm_eeg_specest_' prefix.
- % S.settings - plug-in specific settings
- %
- % S.phase - also save phase dataset (1) or not (0)
- % phase dataset cannot be computed for some
- % spectral estimation methods
- % S.prefix - prefix added before the standard prefix (tf_ or tph_)
- %
- % Output:
- % Dtf - M/EEG object with power (also written on disk)
- % Dtph - M/EEG object with phase (also written on disk)
- %__________________________________________________________________________
- % This is a modular function for which plugins can be developed implementing
- % specific spectral estimation methods. There are 3 basic plugins presently
- % implemented and they can be used as templates for new plugins.
- % The name of a plugin function should start with 'spm_eeg_specest_'
- %
- % morlet (spm_eeg_specest_morlet) - Morlet wavelet transform
- %
- % hilbert (spm_eeg_specest_hilbert) - filtering + Hilbert transform
- %
- % ft_mtmconvol (spm_eeg_specest_mtmconvol) - Fieldtrip implementation
- % of multi-taper spectral
- % analysis
- %__________________________________________________________________________
- % Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging
- % Vladimir Litvak
- % $Id: spm_eeg_tf.m 7449 2018-10-16 13:52:04Z vladimir $
- SVNrev = '$Rev: 7449 $';
- %-Startup
- %--------------------------------------------------------------------------
- spm('FnBanner', mfilename, SVNrev);
- spm('FigName','M/EEG Time-Frequency'); spm('Pointer','Watch');
- %-Configure the analysis
- %--------------------------------------------------------------------------
- if ~isfield(S, 'channels'), S.channels = 'all'; end
- if ~isfield(S, 'timewin'), S.timewin = [-Inf Inf]; end
- if ~isfield(S, 'phase'), S.phase = 0; end
- if ~isfield(S, 'prefix'), S.prefix = ''; end
- if ~isfield(S, 'frequencies') || isempty(S.frequencies)
- S.frequencies = 1:48;
- end
- if ~isfield(S, 'method')
- S.method = 'morlet';
- S.settings = [];
- end
- D = spm_eeg_load(S.D);
- chanind = D.selectchannels(S.channels);
- if isempty(chanind)
- error('No channels selected.');
- end
- if isfield(S, 'settings')
- S1 = S.settings;
- else
- S1 = [];
- end
- timeind = D.indsample(1e-3*min(S.timewin)):D.indsample(1e-3*max(S.timewin));
- % remove uninetended non-uniformities in the frequency axis
- frequencies = S.frequencies;
- if length(frequencies)>1
- df = unique(diff(frequencies));
- if length(df) == 1 || (max(diff(df))/mean(df))<0.1
- df = mean(diff(frequencies));
- frequencies = (0:(length(frequencies)-1))*df + frequencies(1);
- end
- end
- S1.frequencies = frequencies;
- if ~isequal(D.type, 'continuous')
-
- %-Run the analysis on all trials
- %--------------------------------------------------------------------------
- spm_progress_bar('Init', D.ntrials, 'trials done');
- if D.ntrials > 100, Ibar = floor(linspace(1, D.ntrials, 100));
- else Ibar = 1:D.ntrials; end
-
- for k = 1:D.ntrials
- spm_progress_bar('Set','ylabel','reading and computing...');
- trial = feval(['spm_eeg_specest_' S.method], S1, D(chanind, timeind, k), D.time(timeind));
-
- if k == 1
- spm_progress_bar('Set','ylabel','initialising output...');
- if isfield(trial, 'fourier')
- outdata = trial.fourier;
- else
- outdata = trial.pow;
- end
-
- Nchannels = size(outdata, 1);
- Nfrequencies = size(outdata, 2);
- Nsamples = size(outdata, 3);
-
- %-Generate output datasets
- %--------------------------------------------------------------------------
- Dtf = clone(D, [S.prefix 'tf_' D.fname], [Nchannels Nfrequencies Nsamples D.ntrials], 0, 1);
- Dtf = Dtf.frequencies(':', trial.freq);
- Dtf = timeonset(Dtf, trial.time(1));
- Dtf = fsample(Dtf, 1/diff(trial.time(1:2)));
- Dtf = transformtype(Dtf, 'TF');
-
- Dtf = chanlabels(Dtf, 1:Nchannels, D.chanlabels(chanind));
- Dtf = badchannels(Dtf, 1:Nchannels, D.badchannels(chanind));
- Dtf = chantype(Dtf, 1:Nchannels, D.chantype(chanind));
- Dtf = coor2D(Dtf, 1:Nchannels, coor2D(D,chanind));
-
- if S.phase && isfield(trial, 'fourier')
- Dtph = clone(Dtf, [S.prefix 'tph_' D.fname]);
- Dtph = transformtype(Dtph, 'TFphase');
- else
- if ~isfield(trial, 'fourier')
- warning('Phase cannot be estimated with the requested method. Estimating power only.');
- end
-
- Dtph = [];
- end
-
- end
-
- spm_progress_bar('Set','ylabel','writing...');
- if isfield(trial, 'fourier')
- Dtf(:, :, :, k) = trial.fourier.*conj(trial.fourier);
-
- if S.phase
- Dtph(:, :, :, k) = angle(trial.fourier);
- end
- elseif isfield(trial, 'pow')
- Dtf(:, :, :, k) = trial.pow;
- else
- error('The plug-in returned unexpected output');
- end
-
- if ismember(k, Ibar), spm_progress_bar('Set', k); end
- end
-
- spm_progress_bar('Clear');
-
- else % by channel for continuous data
- Nchannels = length(chanind);
-
- spm_progress_bar('Init', Nchannels , 'channels done');
- if Nchannels > 100, Ibar = floor(linspace(1, Nchannels, 100));
- else Ibar = 1:Nchannels; end
-
- for k = 1:Nchannels
- spm_progress_bar('Set','ylabel','reading and computing...');
- trial = feval(['spm_eeg_specest_' S.method], S1, D(chanind(k), timeind), D.time(timeind));
-
- if k == 1
- spm_progress_bar('Set','ylabel','initialising output...');
- if isfield(trial, 'fourier')
- outdata = trial.fourier;
- else
- outdata = trial.pow;
- end
-
- Nfrequencies = size(outdata, 2);
- Nsamples = size(outdata, 3);
-
- %-Generate output datasets
- %--------------------------------------------------------------------------
- Dtf = clone(D, [S.prefix 'tf_' D.fname], [Nchannels Nfrequencies Nsamples 1]);
- Dtf = Dtf.frequencies(':', trial.freq);
- Dtf = timeonset(Dtf, trial.time(1));
- Dtf = fsample(Dtf, 1/diff(trial.time(1:2)));
- Dtf = transformtype(Dtf, 'TF');
-
- Dtf = chanlabels(Dtf, 1:Nchannels, D.chanlabels(chanind));
- Dtf = badchannels(Dtf, 1:Nchannels, D.badchannels(chanind));
- Dtf = chantype(Dtf, 1:Nchannels, D.chantype(chanind));
- Dtf = coor2D(Dtf, 1:Nchannels, coor2D(D,chanind));
-
- ev = Dtf.events;
-
- if ~isempty(ev)
- ev = ev([ev.time]>=Dtf.time(1) & [ev.time]<=Dtf.time(end));
- Dtf = events(Dtf, 1, ev);
- end
-
- if S.phase && isfield(trial, 'fourier')
- Dtph = clone(Dtf, [S.prefix 'tph_' D.fname]);
- Dtph = transformtype(Dtph, 'TFphase');
- else
- if ~isfield(trial, 'fourier')
- warning('Phase cannot be estimated with the requested method. Estimating power only.');
- end
-
- Dtph = [];
- end
- end
-
- spm_progress_bar('Set','ylabel','writing...');
- if isfield(trial, 'fourier')
- Dtf(k, :, :, 1) = trial.fourier.*conj(trial.fourier);
-
- if S.phase
- Dtph(k, :, :, 1) = angle(trial.fourier);
- end
- elseif isfield(trial, 'pow')
- Dtf(k, :, :, 1) = trial.pow;
- else
- error('The plug-in returned unexpected output');
- end
-
- if ismember(k, Ibar), spm_progress_bar('Set', k); end
- end
-
- spm_progress_bar('Clear');
-
- end
- %-Save new M/EEG dataset(s)
- %--------------------------------------------------------------------------
- Dtf = Dtf.history('spm_eeg_tf', S);
- save(Dtf);
- if ~isempty(Dtph)
- Dtph = Dtph.history('spm_eeg_tf', S);
- save(Dtph);
- end
- %-Cleanup
- %--------------------------------------------------------------------------
- fprintf('%-40s: %30s\n','Completed',spm('time')); %-#
- spm('FigName','M/EEG Time Frequency: done'); spm('Pointer','Arrow');
|