% TIMEFREQ - compute time/frequency decomposition of data trials. This
% function is a compute-only function called by
% the more complete time/frequency functions NEWTIMEF
% and NEWCROSSF which also plot TIMEFREQ results.
%
% Usage:
% >> [tf, freqs, times] = timefreq(data, srate);
% >> [tf, freqs, times, itcvals] = timefreq(data, srate, ...
% 'key1', 'val1', 'key2', 'val2' ...)
% Inputs:
% data = [float array] 2-D data array of size (times,trials)
% srate = sampling rate
%
% Optional inputs:
% 'cycles' = [real] indicates the number of cycles for the
% time-frequency decomposition {default: 0}
% if 0, use FFTs and Hanning window tapering.
% or [real positive scalar] Number of cycles in each Morlet
% wavelet, constant across frequencies.
% or [cycles cycles(2)] wavelet cycles increase with
% frequency starting at cycles(1) and,
% if cycles(2) > 1, increasing to cycles(2) at
% the upper frequency,
% or if cycles(2) = 0, same window size at all
% frequencies (similar to FFT if cycles(1) = 1)
% or if cycles(2) = 1, not increasing (same as giving
% only one value for 'cycles'). This corresponds to pure
% wavelet with the same number of cycles at each frequencies
% if 0 < cycles(2) < 1, linear variation in between pure
% wavelets (1) and FFT (0). The exact number of cycles
% at the highest frequency is indicated on the command line.
% 'wavelet' = DEPRECATED, please use 'cycles'. This function does not
% support multitaper. For multitaper, use TIMEF.
% 'wletmethod' = ['dftfilt2'|'dftfilt3'] Wavelet method/program to use.
% {default: 'dftfilt3'}
% 'dftfilt' DEPRECATED. Method used in regular TIMEF
% program. Not available any more.
% 'dftfilt2' Morlet-variant or Hanning DFT (calls DFTFILT2
% to generate wavelets).
% 'dftfilt3' Morlet wavelet or Hanning DFT (exact Tallon
% Baudry). Calls DFTFILT3.
% 'ffttaper' = ['none'|'hanning'|'hamming'|'blackmanharris'] FFT tapering
% function. Default is 'hanning'. Note that 'hamming' and
% 'blackmanharris' require the signal processing toolbox.
% Optional ITC type:
% 'type' = ['coher'|'phasecoher'] Compute either linear coherence
% ('coher') or phase coherence ('phasecoher') also known
% as phase coupling factor' {default: 'phasecoher'}.
% 'subitc' = ['on'|'off'] subtract stimulus locked Inter-Trial Coherence
% (ITC) from x and y. This computes the 'intrinsic' coherence
% x and y not arising from common synchronization to
% experimental events. See notes. {default: 'off'}
%
% Optional detrending:
% 'detrend' = ['on'|'off'], Linearly detrend each data epoch {'off'}
%
% Optional FFT/DFT parameters:
% 'tlimits' = [min max] time limits in ms.
% 'winsize' = If cycles==0 (FFT, see 'wavelet' input): data subwindow
% length (fastest, 2^n<frames);
% if cycles >0: *longest* window length to use. This
% determines the lowest output frequency {~frames/8}
% 'ntimesout' = Number of output times (int<frames-winsize). Enter a
% negative value [-S] to subsample original time by S.
% 'timesout' = Enter an array to obtain spectral decomposition at
% specific time values (note: algorithm find closest time
% point in data and this might result in an unevenly spaced
% time array). Overwrite 'ntimesout'. {def: automatic}
% 'freqs' = [min max] frequency limits. Default [minfreq srate/2],
% minfreq being determined by the number of data points,
% cycles and sampling frequency. Enter a single value
% to compute spectral decompisition at a single frequency
% (note: for FFT the closest frequency will be estimated).
% For wavelet, reducing the max frequency reduce
% the computation load.
% 'padratio' = FFTlength/winsize (2^k) {def: 2}
% Multiplies the number of output frequencies by
% dividing their spacing. When cycles==0, frequency
% spacing is (low_frequency/padratio).
% 'nfreqs' = number of output frequencies. For FFT, closest computed
% frequency will be returned. Overwrite 'padratio' effects
% for wavelets. Default: use 'padratio'.
% 'freqscale' = ['log'|'linear'] frequency scale. Default is 'linear'.
% Note that for obtaining 'log' spaced freqs using FFT,
% closest correspondent frequencies in the 'linear' space
% are returned.
% 'wletmethod'= ['dftfilt2'|'dftfilt3'] Wavelet method/program to use.
% Default is 'dftfilt3'
% 'dftfilt3' Morlet wavelet or Hanning DFT
% 'dftfilt2' Morlet-variant or Hanning DFT.
% Note that there are differences between the Hanning
% DFTs in the two programs.
% 'causal' = ['on'|'off'] apply FFT or time-frequency in a causal
% way where only data before any given latency can
% influence the spectral decomposition. (default: 'off'}
%
% Optional time warping:
% 'timestretch' = {[Refmarks], [Refframes]}
% Stretch amplitude and phase time-course before smoothing.
% Refmarks is a (trials,eventframes) matrix in which rows are
% event marks to time-lock to for a given trial. Each trial
% will be stretched so that its marked events occur at frames
% Refframes. If Refframes is [], the median frame, across trials,
% of the Refmarks latencies will be used. Both Refmarks and
% Refframes are given in frames in this version - will be
% changed to ms in future.
% Outputs:
% tf = complex time frequency array for all trials (freqs,
% times, trials)
% freqs = vector of computed frequencies (Hz)
% times = vector of computed time points (ms)
% itcvals = time frequency "average" for all trials (freqs, times).
% In the coherence case, it is the real mean of the time
% frequency decomposition, but in the phase coherence case
% (see 'type' input'), this is the mean of the normalized
% spectral estimate.
%
% Authors: Arnaud Delorme, Jean Hausser & Scott Makeig
% CNL/Salk Institute 1998-2001; SCCN/INC/UCSD, La Jolla, 2002-
% Fix FFT frequency inaccuracy, bug 874 by WuQiang
%
% See also: TIMEF, NEWTIMEF, CROSSF, NEWCROSSF
% Note: it is not advised to use a FFT decomposition in a log scale. Output
% value are accurate but plotting might not be because of the non-uniform
% frequency output values in log-space. If you have to do it, use a
% padratio as large as possible, or interpolate time-freq image at
% exact log scale values before plotting.
% Copyright (C) 8/1/98 Arnaud Delorme, Sigurd Enghoff & Scott Makeig, SCCN/INC/UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
function [tmpall, freqs, timesout, itcvals] = timefreq(data, srate, varargin)
if nargin < 2
help timefreq;
return;
end
[chan frame trials]= size(data);
if trials == 1 && chan ~= 1
trials = frame;
frame = chan;
chan = 1;
end
g = finputcheck(varargin, ...
{ 'ntimesout' 'integer' [] []; ...
'timesout' 'real' [] []; ...
'winsize' 'integer' [0 Inf] []; ...
'tlimits' 'real' [] []; ...
'detrend' 'string' {'on','off'} 'off'; ...
'causal' 'string' {'on','off'} 'off'; ...
'verbose' 'string' {'on','off'} 'on'; ...
'freqs' 'real' [0 Inf] []; ...
'nfreqs' 'integer' [0 Inf] []; ...
'freqscale' 'string' { 'linear','log','' } 'linear'; ...
'ffttaper' 'string' { 'hanning','hamming','blackmanharris','none' } 'hanning';
'wavelet' 'real' [0 Inf] 0; ...
'cycles' {'real','integer'} [0 Inf] 0; ...
'padratio' 'integer' [1 Inf] 2; ...
'itctype' 'string' {'phasecoher','phasecoher2','coher'} 'phasecoher'; ...
'subitc' 'string' {'on','off'} 'off'; ...
'timestretch' 'cell' [] {}; ...
'wletmethod' 'string' {'dftfilt2','dftfilt3'} 'dftfilt3'; ...
});
if ischar(g), error(g); end
if isempty(g.freqscale), g.freqscale = 'linear'; end
if isempty(g.winsize), g.winsize = max(pow2(nextpow2(frame)-3),4); end
if isempty(g.ntimesout), g.ntimesout = 200; end
if isempty(g.freqs), g.freqs = [0 srate/2]; end
if isempty(g.tlimits), g.tlimits = [0 frame/srate*1000]; end
% checkin parameters
% ------------------
% Use 'wavelet' if 'cycles' undefined for backwards compatibility
if g.cycles == 0
g.cycles = g.wavelet;
end
if (g.winsize > frame)
error('Value of winsize must be less than frame length.');
end
if (pow2(nextpow2(g.padratio)) ~= g.padratio)
error('Value of padratio must be an integer power of two [1,2,4,8,16,...]');
end
% finding frequency limits
% ------------------------
if g.cycles(1) ~= 0 && g.freqs(1) == 0, g.freqs(1) = srate*g.cycles(1)/g.winsize; end
% finding frequencies
% -------------------
if length(g.freqs) == 2
% min and max
% -----------
if g.freqs(1) == 0 && g.cycles(1) ~= 0
g.freqs(1) = srate*g.cycles(1)/g.winsize;
end
% default number of freqs using padratio
% --------------------------------------
if isempty(g.nfreqs)
g.nfreqs = g.winsize/2*g.padratio+1;
% adjust nfreqs depending on frequency range
tmpfreqs = linspace(0, srate/2, g.nfreqs);
tmpfreqs = tmpfreqs(2:end); % remove DC (match the output of PSD)
% adjust limits for FFT (only linear scale)
if g.cycles(1) == 0 && ~strcmpi(g.freqscale, 'log')
if ~any(tmpfreqs == g.freqs(1))
[tmp minind] = min(abs(tmpfreqs-g.freqs(1)));
g.freqs(1) = tmpfreqs(minind);
verboseprintf(g.verbose, 'Adjust min freq. to %3.2f Hz to match FFT output frequencies\n', g.freqs(1));
end
if ~any(tmpfreqs == g.freqs(2))
[tmp minind] = min(abs(tmpfreqs-g.freqs(2)));
g.freqs(2) = tmpfreqs(minind);
verboseprintf(g.verbose, 'Adjust max freq. to %3.2f Hz to match FFT output frequencies\n', g.freqs(2));
end
end
% find number of frequencies
% --------------------------
g.nfreqs = length(tmpfreqs( intersect( find(tmpfreqs >= g.freqs(1)), find(tmpfreqs <= g.freqs(2)))));
if g.freqs(1)==g.freqs(2), g.nfreqs = 1; end
end
% find closest freqs for FFT
% --------------------------
if strcmpi(g.freqscale, 'log')
g.freqs = linspace(log(g.freqs(1)), log(g.freqs(end)), g.nfreqs);
g.freqs = exp(g.freqs);
else
g.freqs = linspace(g.freqs(1), g.freqs(2), g.nfreqs); % this should be OK for FFT
% because of the limit adjustment
end
end
g.nfreqs = length(g.freqs);
% function for time freq initialisation
% -------------------------------------
if (g.cycles(1) == 0) %%%%%%%%%%%%%% constant window-length FFTs %%%%%%%%%%%%%%%%
freqs = linspace(0, srate/2, g.winsize*g.padratio/2+1);
freqs = freqs(2:end); % remove DC (match the output of PSD)
%srate/g.winsize*[1:2/g.padratio:g.winsize]/2
verboseprintf(g.verbose, 'Using %s FFT tapering\n', g.ffttaper);
switch g.ffttaper
case 'hanning', g.win = hanning(g.winsize);
case 'hamming', g.win = hamming(g.winsize);
case 'blackmanharris', g.win = blackmanharris(g.winsize);
case 'none', g.win = ones(g.winsize,1);
end
else % %%%%%%%%%%%%%%%%%% Constant-Q (wavelet) DFTs %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%freqs = srate*g.cycles/g.winsize*[2:2/g.padratio:g.winsize]/2;
%g.win = dftfilt(g.winsize,g.freqs(2)/srate,g.cycles,g.padratio,g.cyclesfact);
freqs = g.freqs;
if length(g.cycles) == 2
if g.cycles(2) < 1
g.cycles = [ g.cycles(1) g.cycles(1)*g.freqs(end)/g.freqs(1)*(1-g.cycles(2))];
end
verboseprintf(g.verbose, 'Using %g cycles at lowest frequency to %g at highest.\n', g.cycles(1), g.cycles(2));
elseif length(g.cycles) == 1
verboseprintf(g.verbose, 'Using %d cycles at all frequencies.\n',g.cycles);
else
verboseprintf(g.verbose, 'Using user-defined cycle for each frequency\n');
end
if strcmp(g.wletmethod, 'dftfilt2')
g.win = dftfilt2(g.freqs,g.cycles,srate, g.freqscale); % uses Morlet taper by default
elseif strcmp(g.wletmethod, 'dftfilt3') % Default
g.win = dftfilt3(g.freqs,g.cycles,srate, 'cycleinc', g.freqscale); % uses Morlet taper by default
else return
end
g.winsize = 0;
for index = 1:length(g.win)
g.winsize = max(g.winsize,length(g.win{index}));
end
end
% compute time vector
% -------------------
[ g.timesout g.indexout ] = gettimes(frame, g.tlimits, g.timesout, g.winsize, g.ntimesout, g.causal, g.verbose);
% -------------------------------
% compute time freq decomposition
% -------------------------------
verboseprintf(g.verbose, 'The window size used is %d samples (%g ms) wide.\n',g.winsize, 1000/srate*g.winsize);
if strcmpi(g.freqscale, 'log') % fastif was having strange "function not available" messages
scaletoprint = 'log';
else scaletoprint = 'linear';
end
verboseprintf(g.verbose, 'Estimating %d %s-spaced frequencies from %2.1f Hz to %3.1f Hz.\n', length(g.freqs), ...
scaletoprint, g.freqs(1), g.freqs(end));
%verboseprintf(g.verbose, 'Estimating %d %s-spaced frequencies from %2.1f Hz to %3.1f Hz.\n', length(g.freqs), ...
% fastif(strcmpi(g.freqscale, 'log'), 'log', 'linear'), g.freqs(1), g.freqs(end));
if g.cycles(1) == 0
if 1
% build large matrix to compute FFT
% ---------------------------------
indices = repmat([-g.winsize/2+1:g.winsize/2]', [1 length(g.indexout) trials]);
indices = indices + repmat(g.indexout, [size(indices,1) 1 trials]);
indices = indices + repmat(reshape(([1:trials]-1)*frame,1,1,trials), [size(indices,1) length(g.indexout) 1]);
if chan > 1
tmpall = repmat(nan,[chan length(freqs) length(g.timesout) trials]);
tmpX = reshape(data(:,indices), [ size(data,1) size(indices)]);
tmpX = bsxfun(@minus, tmpX, mean( tmpX, 2)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
tmpX = bsxfun(@times, tmpX, g.win');
tmpX = fft(tmpX,g.padratio*g.winsize,2);
tmpall = squeeze(tmpX(:,2:g.padratio*g.winsize/2+1,:,:));
else
tmpall = repmat(nan,[length(freqs) length(g.timesout) trials]);
tmpX = data(indices);
tmpX = bsxfun(@minus, tmpX, mean( tmpX, 1)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
tmpX = bsxfun(@times, tmpX, g.win);
%tmpX = fft(tmpX,2^ceil(log2(g.padratio*g.winsize)));
%tmpall = tmpX(2:g.padratio*g.winsize/2+1,:,:);
tmpX = fft(tmpX,g.padratio*g.winsize);
tmpall = tmpX(2:g.padratio*g.winsize/2+1,:,:);
end
else % old iterative computation
tmpall = repmat(nan,[length(freqs) length(g.timesout) trials]);
verboseprintf(g.verbose, 'Processing trial (of %d):',trials);
for trial = 1:trials
if rem(trial,10) == 0, verboseprintf(g.verbose, ' %d',trial); end
if rem(trial,120) == 0, verboseprintf(g.verbose, '\n'); end
for index = 1:length(g.indexout)
if strcmpi(g.causal, 'off')
tmpX = data([-g.winsize/2+1:g.winsize/2]+g.indexout(index)+(trial-1)*frame); % 1 point imprecision
else
tmpX = data([-g.winsize+1:0]+g.indexout(index)+(trial-1)*frame); % 1 point imprecision
end
tmpX = tmpX - mean(tmpX);
if strcmpi(g.detrend, 'on'),
tmpX = detrend(tmpX);
end
tmpX = g.win .* tmpX(:);
tmpX = fft(tmpX,g.padratio*g.winsize);
tmpX = tmpX(2:g.padratio*g.winsize/2+1);
tmpall(:,index, trial) = tmpX(:);
end
end
end
else % wavelet
if chan > 1
% wavelets are processed in groups of the same size
% to speed up computation. Wavelet of groups of different size
% can be processed together but at a cost of a lot of RAM and
% a lot of extra computation -> not efficient
tmpall = repmat(nan,[chan length(freqs) length(g.timesout) trials]);
wt = [ 1 find(diff(cellfun(@length,g.win)))+1 length(g.win)+1];
verboseprintf(g.verbose, 'Computing of %d:', length(wt));
for ind = 1:length(wt)-1
verboseprintf(g.verbose, '.');
wavarray = reshape([ g.win{wt(ind):wt(ind+1)-1} ], [ length(g.win{wt(ind)}) wt(ind+1)-wt(ind) ]);
sizewav = size(wavarray,1)-1;
indices = repmat([-sizewav/2:sizewav/2]', [1 size(wavarray,2) length(g.indexout) trials]);
indices = indices + repmat(reshape(g.indexout, 1,1,length(g.indexout)), [size(indices,1) size(indices,2) 1 trials]);
indices = indices + repmat(reshape(([1:trials]-1)*frame,1,1,1,trials), [size(indices,1) size(indices,2) size(indices,3) 1]);
szfreqdata = [ size(data,1) size(indices) ];
tmpX = reshape(data(:,indices), szfreqdata);
tmpX = bsxfun(@minus, tmpX, mean( tmpX, 2)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
wavarray = reshape(wavarray, [1 size(wavarray,1) size(wavarray,2)]);
tmpall(:,wt(ind):wt(ind+1)-1,:,:,:) = reshape(sum(bsxfun(@times, tmpX, wavarray),2), [szfreqdata(1) szfreqdata(3:end)]);
end
verboseprintf(g.verbose, '\n');
%tmpall = squeeze(tmpall(1,:,:,:));
elseif 0
tmpall = repmat(nan,[length(freqs) length(g.timesout) trials]);
% wavelets are processed in groups of the same size
% to speed up computation. Wavelet of groups of different size
% can be processed together but at a cost of a lot of RAM and
% a lot of extra computation -> not faster than the regular
% iterative method
wt = [ 1 find(diff(cellfun(@length,g.win)))+1 length(g.win)+1];
for ind = 1:length(wt)-1
wavarray = reshape([ g.win{wt(ind):wt(ind+1)-1} ], [ length(g.win{wt(ind)}) wt(ind+1)-wt(ind) ]);
sizewav = size(wavarray,1)-1;
indices = repmat([-sizewav/2:sizewav/2]', [1 size(wavarray,2) length(g.indexout) trials]);
indices = indices + repmat(reshape(g.indexout, 1,1,length(g.indexout)), [size(indices,1) size(indices,2) 1 trials]);
indices = indices + repmat(reshape(([1:trials]-1)*frame,1,1,1,trials), [size(indices,1) size(indices,2) size(indices,3) 1]);
tmpX = data(indices);
tmpX = bsxfun(@minus, tmpX, mean( tmpX, 1)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
tmpall(wt(ind):wt(ind+1)-1,:,:) = squeeze(sum(bsxfun(@times, tmpX, wavarray),1));
end
elseif 0
% wavelets are processed one by one but all windows simultaneously
% -> not faster than the regular iterative method
tmpall = repmat(nan,[length(freqs) length(g.timesout) trials]);
sizewav = length(g.win{1})-1; % max window size
mainc = sizewav/2;
indices = repmat([-sizewav/2:sizewav/2]', [1 length(g.indexout) trials]);
indices = indices + repmat(g.indexout, [size(indices,1) 1 trials]);
indices = indices + repmat(reshape(([1:trials]-1)*frame,1,1,trials), [size(indices,1) length(g.indexout) 1]);
for freqind = 1:length(g.win)
winc = (length(g.win{freqind})-1)/2;
wins = length(g.win{freqind})-1;
wini = [-wins/2:wins/2]+winc+mainc-winc+1;
tmpX = data(indices(wini,:,:));
tmpX = bsxfun(@minus, tmpX, mean( tmpX, 1)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
tmpX = sum(bsxfun(@times, tmpX, g.win{freqind}'),1);
tmpall(freqind,:,:) = tmpX;
end
else
% prepare wavelet filters
% -----------------------
for index = 1:length(g.win)
g.win{index} = transpose(repmat(g.win{index}, [trials 1]));
end
% apply filters
% -------------
verboseprintf(g.verbose, 'Processing time point (of %d):',length(g.timesout));
tmpall = zeros(length(g.win), length(g.indexout), size(data,2));
for index = 1:length(g.indexout)
if rem(index,10) == 0, verboseprintf(g.verbose, ' %d',index); end
if rem(index,120) == 0, verboseprintf(g.verbose, '\n'); end
for freqind = 1:length(g.win)
wav = g.win{freqind};
sizewav = size(wav,1)-1;
%g.indexout(index), size(wav,1), g.freqs(freqind)
if strcmpi(g.causal, 'off')
tmpX = data([-sizewav/2:sizewav/2]+g.indexout(index),:);
else
tmpX = data([-sizewav:0]+g.indexout(index),:);
end
tmpX = tmpX - ones(size(tmpX,1),1)*mean(tmpX);
if strcmpi(g.detrend, 'on'),
for trial = 1:trials
tmpX(:,trial) = detrend(tmpX(:,trial));
end
end
tmpX = sum(wav .* tmpX);
tmpall( freqind, index, :) = tmpX;
end
end
end
end
verboseprintf(g.verbose, '\n');
% time-warp code begins -Jean
% ---------------------------
if ~isempty(g.timestretch) && length(g.timestretch{1}) > 0
timemarks = g.timestretch{1}';
if isempty(g.timestretch{2}) || length(g.timestretch{2}) == 0
timerefs = median(g.timestretch{1}',2);
else
timerefs = g.timestretch{2};
end
trials = size(tmpall,3);
% convert timerefs to subsampled ERSP space
% -----------------------------------------
[dummy refsPos] = min(transpose(abs( ...
repmat(timerefs, [1 length(g.indexout)]) - repmat(g.indexout, [length(timerefs) 1]))));
refsPos(end+1) = 1;
refsPos(end+1) = length(g.indexout);
refsPos = sort(refsPos);
for t=1:trials
% convert timemarks to subsampled ERSP space
% ------------------------------------------
%[dummy pos]=min(abs(repmat(timemarks(2:7,1), [1 length(g.indexout)])-repmat(g.indexout,[6 1])));
outOfTimeRangeTimeWarpMarkers = find(timemarks(:,t) < min(g.indexout) | timemarks(:,t) > max(g.indexout));
% if ~isempty(outOfTimeRangeTimeWarpMarkers)
% verboseprintf(g.verbose, 'Timefreq warning: time-warp latencies in epoch %d are out of time range defined for calculation of ERSP.\n', t);
% end
[dummy marksPos] = min(transpose( ...
abs( ...
repmat(timemarks(:,t), [1 length(g.indexout)]) ...
- repmat(g.indexout, [size(timemarks,1) 1]) ...
) ...
));
marksPos(end+1) = 1;
marksPos(end+1) = length(g.indexout);
marksPos = sort(marksPos);
%now warp tmpall
mytmpall = tmpall(:,:,t);
r = sqrt(mytmpall.*conj(mytmpall));
theta = angle(mytmpall);
% So mytmpall is almost equal to r.*exp(i*theta)
% whos marksPos refsPos
M = timewarp(marksPos, refsPos);
TSr = transpose(M*r');
TStheta = zeros(size(theta,1), size(theta,2));
for freqInd=1:size(TStheta,1)
TStheta(freqInd, :) = angtimewarp(marksPos, refsPos, theta(freqInd, :));
end
TStmpall = TSr.*exp(i*TStheta);
% $$$ keyboard;
tmpall(:,:,t) = TStmpall;
end
end
%time-warp ends
zerovals = tmpall == 0;
if any(reshape(zerovals, 1, prod(size(zerovals))))
tmpall(zerovals) = Inf;
minval = min(tmpall(:)); % remove bug
tmpall(zerovals) = minval;
end
% compute and subtract ITC
% ------------------------
if nargout > 3 || strcmpi(g.subitc, 'on')
itcvals = newtimefitc(tmpall, g.itctype);
end
if strcmpi(g.subitc, 'on')
%a = gcf; figure; imagesc(abs(itcvals)); cbar; figure(a);
if ndims(tmpall) <= 3
tmpall = (tmpall - abs(tmpall) .* repmat(itcvals, [1 1 trials])) ./ abs(tmpall);
else tmpall = (tmpall - abs(tmpall) .* repmat(itcvals, [1 1 1 trials])) ./ abs(tmpall);
end
end
% find closest output frequencies
% -------------------------------
if length(g.freqs) ~= length(freqs) || any(g.freqs ~= freqs)
allindices = zeros(1,length(g.freqs));
for index = 1:length(g.freqs)
[dum ind] = min(abs(freqs-g.freqs(index)));
allindices(index) = ind;
end
verboseprintf(g.verbose, 'finding closest frequencies: %d freqs removed\n', length(freqs)-length(allindices));
freqs = freqs(allindices);
if ndims(tmpall) <= 3
tmpall = tmpall(allindices,:,:);
else tmpall = tmpall(:,allindices,:,:);
end
if nargout > 3 || strcmpi(g.subitc, 'on')
if ndims(tmpall) <= 3
itcvals = itcvals(allindices,:,:);
else itcvals = itcvals(:,allindices,:,:);
end
end
end
timesout = g.timesout;
%figure; imagesc(abs(sum(itcvals,3))); cbar;
return;
function w = hanning(n)
if ~rem(n,2)
w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
w = [w; w(end:-1:1)];
else
w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
w = [w; w(end-1:-1:1)];
end
% get time points
% ---------------
function [ timevals, timeindices ] = gettimes(frames, tlimits, timevar, winsize, ntimevar, causal, verbose);
timevect = linspace(tlimits(1), tlimits(2), frames);
srate = 1000*(frames-1)/(tlimits(2)-tlimits(1));
if isempty(timevar) % no pre-defined time points
if ntimevar(1) > 0
% generate linearly space vector
% ------------------------------
if (ntimevar > frames-winsize)
ntimevar = frames-winsize;
if ntimevar < 0
error('Not enough data points, reduce the window size or lowest frequency');
end
verboseprintf(verbose, ['Value of ''timesout'' must be <= frame-winsize, ''timesout'' adjusted to ' int2str(ntimevar) '\n']);
end
npoints = ntimevar(1);
wintime = 500*winsize/srate;
if strcmpi(causal, 'on')
timevals = linspace(tlimits(1)+2*wintime, tlimits(2), npoints);
else timevals = linspace(tlimits(1)+wintime, tlimits(2)-wintime, npoints);
end
verboseprintf(verbose, 'Generating %d time points (%1.1f to %1.1f ms)\n', npoints, min(timevals), max(timevals));
else
% subsample data
% --------------
nsub = -ntimevar(1);
if strcmpi(causal, 'on')
timeindices = [ceil(winsize+nsub):nsub:length(timevect)];
else timeindices = [ceil(winsize/2+nsub/2):nsub:length(timevect)-ceil(winsize/2)-1];
end
timevals = timevect( timeindices ); % the conversion at line 741 leaves timeindices unchanged
verboseprintf(verbose, 'Subsampling by %d (%1.1f to %1.1f ms)\n', nsub, min(timevals), max(timevals));
end
else
timevals = timevar;
% check boundaries
% ----------------
wintime = 500*winsize/srate;
if strcmpi(causal, 'on')
tmpind = find( (timevals >= tlimits(1)+2*wintime-0.0001) & (timevals <= tlimits(2)) );
else tmpind = find( (timevals >= tlimits(1)+wintime-0.0001) & (timevals <= tlimits(2)-wintime+0.0001) );
end
% 0.0001 account for numerical inaccuracies on opteron computers
if isempty(tmpind)
error('No time points. Reduce time window or minimum frequency.');
end
if length(timevals) ~= length(tmpind)
verboseprintf(verbose, 'Warning: %d out of %d time values were removed (now %3.2f to %3.2f ms) so the lowest\n', ...
length(timevals)-length(tmpind), length(timevals), timevals(tmpind(1)), timevals(tmpind(end)));
verboseprintf(verbose, ' frequency could be computed with the requested accuracy\n');
end
timevals = timevals(tmpind);
end
% find closet points in data
% --------------------------
timeindices = round(eeg_lat2point(timevals, 1, srate, tlimits, 1E-3));
if length(timeindices) < length(unique(timeindices))
timeindices = unique_bc(timeindices)
verboseprintf(verbose, 'Warning: duplicate times, reduce the number of output times\n');
end
if length(unique(timeindices(2:end)-timeindices(1:end-1))) > 1
verboseprintf(verbose, 'Finding closest points for time variable\n');
verboseprintf(verbose, 'Time values for time/freq decomposition is not perfectly uniformly distributed\n');
else
verboseprintf(verbose, 'Distribution of data point for time/freq decomposition is perfectly uniform\n');
end
timevals = timevect(timeindices);
% DEPRECATED, FOR C INTERFACE
function nofunction()
% C PART %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename = [ 'tmpcrossf' num2str(round(rand(1)*1000)) ];
f = fopen([ filename '.in'], 'w');
fwrite(f, tmpsaveall, 'int32');
fwrite(f, g.detret, 'int32');
fwrite(f, g.srate, 'int32');
fwrite(f, g.maxfreq, 'int32');
fwrite(f, g.padratio, 'int32');
fwrite(f, g.cycles, 'int32');
fwrite(f, g.winsize, 'int32');
fwrite(f, g.timesout, 'int32');
fwrite(f, g.subitc, 'int32');
fwrite(f, g.type, 'int32');
fwrite(f, trials, 'int32');
fwrite(f, g.naccu, 'int32');
fwrite(f, length(X), 'int32');
fwrite(f, X, 'double');
fwrite(f, Y, 'double');
fclose(f);
command = [ '!cppcrosff ' filename '.in ' filename '.out' ];
eval(command);
f = fopen([ filename '.out'], 'r');
size1 = fread(f, 'int32', 1);
size2 = fread(f, 'int32', 1);
Rreal = fread(f, 'double', [size1 size2]);
Rimg = fread(f, 'double', [size1 size2]);
Coher.R = Rreal + j*Rimg;
Boot.Coherboot.R = [];
Boot.Rsignif = [];
function verboseprintf(verbose, varargin)
if strcmpi(verbose, 'on')
fprintf(varargin{:});
end
Datasets
Models
