timefreq

 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'  = [0] Use FFTs and Hanning window tapering.  {default: 0}
                   or [real positive scalar] Number of cycles in each Morlet
                   wavelet, constant across frequencies.
                   or [cycles(1) cycles(2)] wavelet cycles increase with frequency
                   starting at cycles(1) and, 
                   if cycles(2) > 1, increasing to cycles(2) at
                   the upper frequency,
                   if cycles(2) < 1, increasing by a factor of cycles(2),
                   or if cycles(2) = 1, not increasing (same as giving
                   only one value for 'cycles').
       '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 correspondant 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 betweeen 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 innacuracy, bug 874 by WuQiang

 See also: timef(), newtimef(), crossf(), newcrossf()

0001 % timefreq() - compute time/frequency decomposition of data trials. This
0002 %              function is a compute-only function called by
0003 %              the more complete time/frequency functions newtimef()
0004 %              and newcrossf() which also plot timefreq() results.
0005 %
0006 % Usage:
0007 %     >> [tf, freqs, times]          = timefreq(data, srate);
0008 %     >> [tf, freqs, times, itcvals] = timefreq(data, srate, ...
0009 %                                        'key1', 'val1', 'key2', 'val2' ...)
0010 % Inputs:
0011 %         data    = [float array] 2-D data array of size (times,trials)
0012 %         srate   = sampling rate
0013 %
0014 % Optional inputs:
0015 %       'cycles'  = [0] Use FFTs and Hanning window tapering.  {default: 0}
0016 %                   or [real positive scalar] Number of cycles in each Morlet
0017 %                   wavelet, constant across frequencies.
0018 %                   or [cycles(1) cycles(2)] wavelet cycles increase with frequency
0019 %                   starting at cycles(1) and,
0020 %                   if cycles(2) > 1, increasing to cycles(2) at
0021 %                   the upper frequency,
0022 %                   if cycles(2) < 1, increasing by a factor of cycles(2),
0023 %                   or if cycles(2) = 1, not increasing (same as giving
0024 %                   only one value for 'cycles').
0025 %       'wavelet' = DEPRECATED, please use 'cycles'. This function does not
0026 %                   support multitaper. For multitaper, use timef().
0027 %       'wletmethod' = ['dftfilt2'|'dftfilt3'] Wavelet method/program to use.
0028 %                   {default: 'dftfilt3'}
0029 %                   'dftfilt'  DEPRECATED. Method used in regular timef()
0030 %                              program. Not available any more.
0031 %                   'dftfilt2' Morlet-variant or Hanning DFT (calls dftfilt2()
0032 %                              to generate wavelets).
0033 %                   'dftfilt3' Morlet wavelet or Hanning DFT (exact Tallon
0034 %                              Baudry). Calls dftfilt3().
0035 %       'ffttaper' = ['none'|'hanning'|'hamming'|'blackmanharris'] FFT tapering
0036 %                   function. Default is 'hanning'. Note that 'hamming' and
0037 %                   'blackmanharris' require the signal processing toolbox.
0038 % Optional ITC type:
0039 %        'type'   = ['coher'|'phasecoher'] Compute either linear coherence
0040 %                   ('coher') or phase coherence ('phasecoher') also known
0041 %                   as phase coupling factor' {default: 'phasecoher'}.
0042 %        'subitc' = ['on'|'off'] subtract stimulus locked Inter-Trial Coherence
0043 %                   (ITC) from x and y. This computes the  'intrinsic' coherence
0044 %                   x and y not arising from common synchronization to
0045 %                   experimental events. See notes. {default: 'off'}
0046 %
0047 % Optional detrending:
0048 %       'detrend' = ['on'|'off'], Linearly detrend each data epoch  {'off'}
0049 %
0050 % Optional FFT/DFT parameters:
0051 %      'tlimits'  = [min max] time limits in ms.
0052 %      'winsize'  = If cycles==0 (FFT, see 'wavelet' input): data subwindow
0053 %                   length (fastest, 2^n<frames);
0054 %                   if cycles >0: *longest* window length to use. This
0055 %                   determines the lowest output frequency  {~frames/8}
0056 %      'ntimesout' = Number of output times (int<frames-winsize). Enter a
0057 %                   negative value [-S] to subsample original time by S.
0058 %      'timesout' = Enter an array to obtain spectral decomposition at
0059 %                   specific time values (note: algorithm find closest time
0060 %                   point in data and this might result in an unevenly spaced
0061 %                   time array). Overwrite 'ntimesout'. {def: automatic}
0062 %      'freqs'    = [min max] frequency limits. Default [minfreq srate/2],
0063 %                   minfreq being determined by the number of data points,
0064 %                   cycles and sampling frequency. Enter a single value
0065 %                   to compute spectral decompisition at a single frequency
0066 %                   (note: for FFT the closest frequency will be estimated).
0067 %                   For wavelet, reducing the max frequency reduce
0068 %                   the computation load.
0069 %      'padratio' = FFTlength/winsize (2^k)                     {def: 2}
0070 %                   Multiplies the number of output frequencies by
0071 %                   dividing their spacing. When cycles==0, frequency
0072 %                   spacing is (low_frequency/padratio).
0073 %      'nfreqs'   = number of output frequencies. For FFT, closest computed
0074 %                   frequency will be returned. Overwrite 'padratio' effects
0075 %                   for wavelets. Default: use 'padratio'.
0076 %     'freqscale' = ['log'|'linear'] frequency scale. Default is 'linear'.
0077 %                   Note that for obtaining 'log' spaced freqs using FFT,
0078 %                   closest correspondant frequencies in the 'linear' space
0079 %                   are returned.
0080 %     'wletmethod'= ['dftfilt2'|'dftfilt3'] Wavelet method/program to use.
0081 %                   Default is 'dftfilt3'
0082 %                   'dftfilt3' Morlet wavelet or Hanning DFT
0083 %                   'dftfilt2' Morlet-variant or Hanning DFT.
0084 %                   Note that there are differences betweeen the Hanning
0085 %                   DFTs in the two programs.
0086 %       'causal'  = ['on'|'off'] apply FFT or time-frequency in a causal
0087 %                   way where only data before any given latency can
0088 %                   influence the spectral decomposition. (default: 'off'}
0089 %
0090 % Optional time warping:
0091 %   'timestretch' = {[Refmarks], [Refframes]}
0092 %                   Stretch amplitude and phase time-course before smoothing.
0093 %                   Refmarks is a (trials,eventframes) matrix in which rows are
0094 %                   event marks to time-lock to for a given trial. Each trial
0095 %                   will be stretched so that its marked events occur at frames
0096 %                   Refframes. If Refframes is [], the median frame, across trials,
0097 %                   of the Refmarks latencies will be used. Both Refmarks and
0098 %                   Refframes are given in frames in this version - will be
0099 %                   changed to ms in future.
0100 % Outputs:
0101 %         tf      = complex time frequency array for all trials (freqs,
0102 %                   times, trials)
0103 %         freqs   = vector of computed frequencies (Hz)
0104 %         times   = vector of computed time points (ms)
0105 %         itcvals = time frequency "average" for all trials (freqs, times).
0106 %                   In the coherence case, it is the real mean of the time
0107 %                   frequency decomposition, but in the phase coherence case
0108 %                   (see 'type' input'), this is the mean of the normalized
0109 %                   spectral estimate.
0110 %
0111 % Authors: Arnaud Delorme, Jean Hausser & Scott Makeig
0112 %          CNL/Salk Institute 1998-2001; SCCN/INC/UCSD, La Jolla, 2002-
0113 %          Fix FFT frequency innacuracy, bug 874 by WuQiang
0114 %
0115 % See also: timef(), newtimef(), crossf(), newcrossf()
0116 
0117 % Note: it is not advised to use a FFT decomposition in a log scale. Output
0118 %       value are accurate but plotting might not be because of the non-uniform
0119 %       frequency output values in log-space. If you have to do it, use a
0120 %       padratio as large as possible, or interpolate time-freq image at
0121 %       exact log scale values before plotting.
0122 
0123 % Copyright (C) 8/1/98  Arnaud Delorme, Sigurd Enghoff & Scott Makeig, SCCN/INC/UCSD
0124 %
0125 % This program is free software; you can redistribute it and/or modify
0126 % it under the terms of the GNU General Public License as published by
0127 % the Free Software Foundation; either version 2 of the License, or
0128 % (at your option) any later version.
0129 %
0130 % This program is distributed in the hope that it will be useful,
0131 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0132 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0133 % GNU General Public License for more details.
0134 %
0135 % You should have received a copy of the GNU General Public License
0136 % along with this program; if not, write to the Free Software
0137 % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
0138 
0139 function [tmpall, freqs, timesout, itcvals] = timefreq(data, srate, varargin)
0140 
0141 if nargin < 2
0142     help timefreq;
0143     return;
0144 end;
0145 
0146 [chan frame trials]= size(data);
0147 if trials == 1 && chan ~= 1
0148     trials = frame;
0149     frame  = chan;
0150     chan   = 1;
0151 end;
0152 g = finputcheck(varargin, ...
0153     { 'ntimesout'     'integer'  []                     []; ...
0154     'timesout'      'real'     []                       []; ...
0155     'winsize'       'integer'  [0 Inf]                  []; ...
0156     'tlimits'       'real'     []                       []; ...
0157     'detrend'       'string'   {'on' 'off'}             'off'; ...
0158     'causal'        'string'   {'on' 'off'}             'off'; ...
0159     'verbose'       'string'   {'on' 'off'}             'on'; ...
0160     'freqs'         'real'     [0 Inf]                  []; ...
0161     'nfreqs'        'integer'  [0 Inf]                  []; ...
0162     'freqscale'     'string'   { 'linear' 'log' '' }    'linear'; ...
0163     'ffttaper'      'string'   { 'hanning' 'hamming' 'blackmanharris' 'none' }  'hanning';
0164     'wavelet'       'real'     [0 Inf]                  0; ...
0165     'cycles'        {'real','integer'}    [0 Inf]       0; ...
0166     'padratio'      'integer'  [1 Inf]                  2; ...
0167     'itctype'       'string'   {'phasecoher' 'phasecoher2' 'coher'}  'phasecoher'; ...
0168     'subitc'        'string'   {'on' 'off'}             'off'; ...
0169     'timestretch'   'cell'     []                       {}; ...
0170     'wletmethod'    'string'   {'dftfilt2' 'dftfilt3'}    'dftfilt3'; ...
0171     });
0172 if isstr(g), error(g); end;
0173 if isempty(g.freqscale), g.freqscale = 'linear'; end;
0174 if isempty(g.winsize),   g.winsize   = max(pow2(nextpow2(frame)-3),4); end;
0175 if isempty(g.ntimesout), g.ntimesout = 200; end;
0176 if isempty(g.freqs),     g.freqs     = [0 srate/2]; end;
0177 if isempty(g.tlimits),   g.tlimits   = [0 frame/srate*1000]; end;
0178 
0179 % checkin parameters
0180 % ------------------
0181 
0182 % Use 'wavelet' if 'cycles' undefined for backwards compatibility
0183 if g.cycles == 0
0184     g.cycles = g.wavelet;
0185 end
0186 
0187 if (g.winsize > frame)
0188     error('Value of winsize must be less than frame length.');
0189 end
0190 if (pow2(nextpow2(g.padratio)) ~= g.padratio)
0191     error('Value of padratio must be an integer power of two [1,2,4,8,16,...]');
0192 end
0193 
0194 % finding frequency limits
0195 % ------------------------
0196 if g.cycles(1) ~= 0 & g.freqs(1) == 0, g.freqs(1) = srate*g.cycles(1)/g.winsize; end;
0197 
0198 % finding frequencies
0199 % -------------------
0200 if length(g.freqs) == 2
0201 
0202     % min and max
0203     % -----------
0204     if g.freqs(1) == 0 & g.cycles(1) ~= 0
0205         g.freqs(1) = srate*g.cycles(1)/g.winsize;
0206     end;
0207 
0208     % default number of freqs using padratio
0209     % --------------------------------------
0210     if isempty(g.nfreqs)
0211         g.nfreqs = g.winsize/2*g.padratio+1;
0212         % adjust nfreqs depending on frequency range
0213         tmpfreqs = linspace(0, srate/2, g.nfreqs);
0214         tmpfreqs = tmpfreqs(2:end);  % remove DC (match the output of PSD)
0215 
0216         % adjust limits for FFT (only linear scale)
0217         if g.cycles(1) == 0 & ~strcmpi(g.freqscale, 'log')
0218             if ~any(tmpfreqs == g.freqs(1))
0219                 [tmp minind] = min(abs(tmpfreqs-g.freqs(1)));
0220                 g.freqs(1)   = tmpfreqs(minind);
0221                 verboseprintf(g.verbose, 'Adjust min freq. to %3.2f Hz to match FFT output frequencies\n', g.freqs(1));
0222             end;
0223             if ~any(tmpfreqs == g.freqs(2))
0224                 [tmp minind] = min(abs(tmpfreqs-g.freqs(2)));
0225                 g.freqs(2)   = tmpfreqs(minind);
0226                 verboseprintf(g.verbose, 'Adjust max freq. to %3.2f Hz to match FFT output frequencies\n', g.freqs(2));
0227             end;
0228         end;
0229 
0230         % find number of frequencies
0231         % --------------------------
0232         g.nfreqs = length(tmpfreqs( intersect( find(tmpfreqs >= g.freqs(1)), find(tmpfreqs <= g.freqs(2)))));
0233         if g.freqs(1)==g.freqs(2), g.nfreqs = 1; end;
0234     end;
0235 
0236     % find closest freqs for FFT
0237     % --------------------------
0238     if strcmpi(g.freqscale, 'log')
0239         g.freqs = linspace(log(g.freqs(1)), log(g.freqs(end)), g.nfreqs);
0240         g.freqs = exp(g.freqs);
0241     else
0242         g.freqs = linspace(g.freqs(1), g.freqs(2), g.nfreqs); % this should be OK for FFT
0243         % because of the limit adjustment
0244     end;
0245 end;
0246 g.nfreqs = length(g.freqs);
0247 
0248 % function for time freq initialisation
0249 % -------------------------------------
0250 if (g.cycles(1) == 0) %%%%%%%%%%%%%% constant window-length FFTs %%%%%%%%%%%%%%%%
0251     freqs = linspace(0, srate/2, g.winsize*g.padratio/2+1);
0252     freqs = freqs(2:end); % remove DC (match the output of PSD)
0253     %srate/g.winsize*[1:2/g.padratio:g.winsize]/2
0254     verboseprintf(g.verbose, 'Using %s FFT tapering\n', g.ffttaper);
0255     switch g.ffttaper
0256         case 'hanning',        g.win   = hanning(g.winsize);
0257         case 'hamming',        g.win   = hamming(g.winsize);
0258         case 'blackmanharris', g.win   = blackmanharris(g.winsize);
0259         case 'none',           g.win   = ones(g.winsize,1);
0260     end;
0261 else % %%%%%%%%%%%%%%%%%% Constant-Q (wavelet) DFTs %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0262     %freqs = srate*g.cycles/g.winsize*[2:2/g.padratio:g.winsize]/2;
0263     %g.win = dftfilt(g.winsize,g.freqs(2)/srate,g.cycles,g.padratio,g.cyclesfact);
0264 
0265     freqs = g.freqs;
0266     if length(g.cycles) == 2
0267         if g.cycles(2) < 1
0268             g.cycles = [ g.cycles(1) g.cycles(1)*g.freqs(end)/g.freqs(1)*(1-g.cycles(2))];
0269         end
0270         verboseprintf(g.verbose, 'Using %g cycles at lowest frequency to %g at highest.\n', g.cycles(1), g.cycles(2));
0271     elseif length(g.cycles) == 1
0272         verboseprintf(g.verbose, 'Using %d cycles at all frequencies.\n',g.cycles);
0273     else
0274         verboseprintf(g.verbose, 'Using user-defined cycle for each frequency\n');
0275     end
0276     if strcmp(g.wletmethod, 'dftfilt2')
0277         g.win    = dftfilt2(g.freqs,g.cycles,srate, g.freqscale); % uses Morlet taper by default
0278     elseif strcmp(g.wletmethod, 'dftfilt3')     % Default
0279         g.win    = dftfilt3(g.freqs,g.cycles,srate, 'cycleinc', g.freqscale); % uses Morlet taper by default
0280     else return
0281     end
0282     g.winsize = 0;
0283     for index = 1:length(g.win)
0284         g.winsize = max(g.winsize,length(g.win{index}));
0285     end;
0286 end;
0287 
0288 % compute time vector
0289 % -------------------
0290 [ g.timesout g.indexout ] = gettimes(frame, g.tlimits, g.timesout, g.winsize, g.ntimesout, g.causal, g.verbose);
0291 
0292 % -------------------------------
0293 % compute time freq decomposition
0294 % -------------------------------
0295 verboseprintf(g.verbose, 'The window size used is %d samples (%g ms) wide.\n',g.winsize, 1000/srate*g.winsize);
0296 if strcmpi(g.freqscale, 'log') % fastif was having strange "function not available" messages
0297     scaletoprint = 'log';
0298 else scaletoprint = 'linear';
0299 end
0300 verboseprintf(g.verbose, 'Estimating %d %s-spaced frequencies from %2.1f Hz to %3.1f Hz.\n', length(g.freqs), ...
0301     scaletoprint, g.freqs(1), g.freqs(end));
0302 %verboseprintf(g.verbose, 'Estimating %d %s-spaced frequencies from %2.1f Hz to %3.1f Hz.\n', length(g.freqs), ...
0303 %    fastif(strcmpi(g.freqscale, 'log'), 'log', 'linear'), g.freqs(1), g.freqs(end));
0304 
0305 if g.cycles(1) == 0
0306     if 1
0307         % build large matrix to compute FFT
0308         % ---------------------------------
0309         indices = repmat([-g.winsize/2+1:g.winsize/2]', [1 length(g.indexout) trials]);
0310         indices = indices + repmat(g.indexout, [size(indices,1) 1 trials]);
0311         indices = indices + repmat(reshape(([1:trials]-1)*frame,1,1,trials), [size(indices,1) length(g.indexout) 1]);
0312         if chan > 1
0313             tmpall = repmat(nan,[chan length(freqs) length(g.timesout) trials]);
0314             tmpX = reshape(data(:,indices), [ size(data,1) size(indices)]);
0315             tmpX = bsxfun(@minus, tmpX, mean( tmpX, 2)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
0316             tmpX = bsxfun(@times, tmpX, g.win');
0317             tmpX = fft(tmpX,g.padratio*g.winsize,2);
0318             tmpall = squeeze(tmpX(:,2:g.padratio*g.winsize/2+1,:,:));
0319         else
0320             tmpall = repmat(nan,[length(freqs) length(g.timesout) trials]);
0321             tmpX    = data(indices);
0322             tmpX = bsxfun(@minus, tmpX, mean( tmpX, 1)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
0323             tmpX = bsxfun(@times, tmpX, g.win);
0324             %tmpX = fft(tmpX,2^ceil(log2(g.padratio*g.winsize)));
0325             %tmpall = tmpX(2:g.padratio*g.winsize/2+1,:,:);
0326             tmpX = fft(tmpX,g.padratio*g.winsize);
0327             tmpall = tmpX(2:g.padratio*g.winsize/2+1,:,:);
0328         end;
0329     else % old iterative computation
0330         tmpall = repmat(nan,[length(freqs) length(g.timesout) trials]);
0331         verboseprintf(g.verbose, 'Processing trial (of %d):',trials);
0332         for trial = 1:trials
0333             if rem(trial,10) == 0,  verboseprintf(g.verbose, ' %d',trial); end
0334             if rem(trial,120) == 0, verboseprintf(g.verbose, '\n'); end
0335             for index = 1:length(g.indexout)
0336                 if strcmpi(g.causal, 'off')
0337                     tmpX = data([-g.winsize/2+1:g.winsize/2]+g.indexout(index)+(trial-1)*frame); % 1 point imprecision
0338                 else
0339                     tmpX = data([-g.winsize+1:0]+g.indexout(index)+(trial-1)*frame); % 1 point imprecision
0340                 end;
0341 
0342                 tmpX = tmpX - mean(tmpX);
0343                 if strcmpi(g.detrend, 'on'),
0344                     tmpX = detrend(tmpX);
0345                 end;
0346 
0347                 tmpX = g.win .* tmpX(:);
0348                 tmpX = fft(tmpX,g.padratio*g.winsize);
0349                 tmpX = tmpX(2:g.padratio*g.winsize/2+1);
0350                 tmpall(:,index, trial) = tmpX(:);
0351             end;
0352         end;
0353     end;
0354 else % wavelet
0355     if chan > 1
0356         % wavelets are processed in groups of the same size
0357         % to speed up computation. Wavelet of groups of different size
0358         % can be processed together but at a cost of a lot of RAM and
0359         % a lot of extra computation -> not efficient
0360         tmpall = repmat(nan,[chan length(freqs) length(g.timesout) trials]);
0361         wt = [ 1 find(diff(cellfun(@length,g.win)))+1 length(g.win)+1];
0362         verboseprintf(g.verbose, 'Computing of %d:', length(wt));
0363         for ind = 1:length(wt)-1
0364             verboseprintf(g.verbose, '.');
0365             wavarray = reshape([ g.win{wt(ind):wt(ind+1)-1} ], [ length(g.win{wt(ind)}) wt(ind+1)-wt(ind) ]);
0366             sizewav = size(wavarray,1)-1;
0367             indices = repmat([-sizewav/2:sizewav/2]', [1 size(wavarray,2) length(g.indexout) trials]);
0368             indices = indices + repmat(reshape(g.indexout, 1,1,length(g.indexout)), [size(indices,1) size(indices,2) 1 trials]);
0369             indices = indices + repmat(reshape(([1:trials]-1)*frame,1,1,1,trials),  [size(indices,1) size(indices,2) size(indices,3) 1]);
0370             szfreqdata = [ size(data,1) size(indices) ];
0371             tmpX    = reshape(data(:,indices), szfreqdata);
0372             tmpX    = bsxfun(@minus, tmpX, mean( tmpX, 2)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
0373             wavarray = reshape(wavarray, [1 size(wavarray,1) size(wavarray,2)]);
0374             tmpall(:,wt(ind):wt(ind+1)-1,:,:,:) = reshape(sum(bsxfun(@times, tmpX, wavarray),2), [szfreqdata(1) szfreqdata(3:end)]);
0375         end;
0376         verboseprintf(g.verbose, '\n');
0377         %tmpall = squeeze(tmpall(1,:,:,:));
0378     elseif 0
0379         tmpall = repmat(nan,[length(freqs) length(g.timesout) trials]);
0380         % wavelets are processed in groups of the same size
0381         % to speed up computation. Wavelet of groups of different size
0382         % can be processed together but at a cost of a lot of RAM and
0383         % a lot of extra computation -> not faster than the regular
0384         % iterative method
0385         wt = [ 1 find(diff(cellfun(@length,g.win)))+1 length(g.win)+1];
0386         for ind = 1:length(wt)-1
0387             wavarray = reshape([ g.win{wt(ind):wt(ind+1)-1} ], [ length(g.win{wt(ind)}) wt(ind+1)-wt(ind) ]);
0388             sizewav = size(wavarray,1)-1;
0389             indices = repmat([-sizewav/2:sizewav/2]', [1 size(wavarray,2) length(g.indexout) trials]);
0390             indices = indices + repmat(reshape(g.indexout, 1,1,length(g.indexout)), [size(indices,1) size(indices,2) 1 trials]);
0391             indices = indices + repmat(reshape(([1:trials]-1)*frame,1,1,1,trials), [size(indices,1) size(indices,2) size(indices,3) 1]);
0392             tmpX    = data(indices);
0393             tmpX    = bsxfun(@minus, tmpX, mean( tmpX, 1)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
0394             tmpall(wt(ind):wt(ind+1)-1,:,:)  = squeeze(sum(bsxfun(@times, tmpX, wavarray),1));
0395         end;
0396     elseif 0
0397         % wavelets are processed one by one but all windows simultaneously
0398         % -> not faster than the regular iterative method
0399         tmpall  = repmat(nan,[length(freqs) length(g.timesout) trials]);
0400         sizewav = length(g.win{1})-1; % max window size
0401         mainc   = sizewav/2;
0402         indices = repmat([-sizewav/2:sizewav/2]', [1 length(g.indexout) trials]);
0403         indices = indices + repmat(g.indexout, [size(indices,1) 1 trials]);
0404         indices = indices + repmat(reshape(([1:trials]-1)*frame,1,1,trials), [size(indices,1) length(g.indexout) 1]);
0405         
0406         for freqind = 1:length(g.win)
0407             winc = (length(g.win{freqind})-1)/2;
0408             wins = length(g.win{freqind})-1;
0409             wini = [-wins/2:wins/2]+winc+mainc-winc+1;
0410             tmpX = data(indices(wini,:,:));
0411             tmpX = bsxfun(@minus, tmpX, mean( tmpX, 1)); % avoids repmat - faster than tmpX = tmpX - repmat(mean(tmpX), [size(tmpX,1) 1 1]);
0412             tmpX = sum(bsxfun(@times, tmpX, g.win{freqind}'),1);
0413             tmpall(freqind,:,:) = tmpX;
0414         end;
0415     else
0416         % prepare wavelet filters
0417         % -----------------------
0418         for index = 1:length(g.win)
0419             g.win{index} = transpose(repmat(g.win{index}, [trials 1]));
0420         end;
0421 
0422         % apply filters
0423         % -------------
0424         verboseprintf(g.verbose, 'Processing time point (of %d):',length(g.timesout));
0425         for index = 1:length(g.indexout)
0426             if rem(index,10) == 0,  verboseprintf(g.verbose, ' %d',index); end
0427             if rem(index,120) == 0, verboseprintf(g.verbose, '\n'); end
0428             for freqind = 1:length(g.win)
0429                 wav = g.win{freqind}; 
0430                 sizewav = size(wav,1)-1;
0431                 %g.indexout(index), size(wav,1), g.freqs(freqind)
0432                 if strcmpi(g.causal, 'off')
0433                     tmpX = data([-sizewav/2:sizewav/2]+g.indexout(index),:);
0434                 else
0435                     tmpX = data([-sizewav:0]+g.indexout(index),:);
0436                 end;
0437 
0438                 tmpX = tmpX - ones(size(tmpX,1),1)*mean(tmpX);
0439                 if strcmpi(g.detrend, 'on'),
0440                     for trial = 1:trials
0441                         tmpX(:,trial) = detrend(tmpX(:,trial));
0442                     end;
0443                 end;
0444 
0445                 tmpX = sum(wav .* tmpX);
0446                 tmpall( freqind, index, :) = tmpX;
0447             end;
0448         end;
0449     end;
0450 end;
0451 verboseprintf(g.verbose, '\n');
0452 
0453 % time-warp code begins -Jean
0454 % ---------------------------
0455 if ~isempty(g.timestretch) && length(g.timestretch{1}) > 0
0456 
0457     timemarks = g.timestretch{1}';
0458     if isempty(g.timestretch{2}) | length(g.timestretch{2}) == 0
0459         timerefs = median(g.timestretch{1}',2);
0460     else
0461         timerefs = g.timestretch{2};
0462     end
0463     trials = size(tmpall,3);
0464 
0465     % convert timerefs to subsampled ERSP space
0466     % -----------------------------------------
0467 
0468     [dummy refsPos] = min(transpose(abs( ...
0469         repmat(timerefs, [1 length(g.indexout)]) - repmat(g.indexout, [length(timerefs) 1]))));
0470     refsPos(end+1) = 1;
0471     refsPos(end+1) = length(g.indexout);
0472     refsPos = sort(refsPos);
0473 
0474     for t=1:trials
0475 
0476         % convert timemarks to subsampled ERSP space
0477         % ------------------------------------------
0478 
0479         %[dummy pos]=min(abs(repmat(timemarks(2:7,1), [1 length(g.indexout)])-repmat(g.indexout,[6 1])));
0480 
0481         outOfTimeRangeTimeWarpMarkers = find(timemarks(:,t) < min(g.indexout) | timemarks(:,t) > max(g.indexout));
0482         
0483 %         if ~isempty(outOfTimeRangeTimeWarpMarkers)
0484 %             verboseprintf(g.verbose, 'Timefreq warning: time-warp latencies in epoch %d are out of time range defined for calculation of ERSP.\n', t);
0485 %         end;
0486         
0487         [dummy marksPos] = min(transpose( ...
0488             abs( ...
0489             repmat(timemarks(:,t), [1 length(g.indexout)]) ...
0490             - repmat(g.indexout, [size(timemarks,1) 1]) ...
0491             ) ...
0492             ));
0493          
0494 
0495         marksPos(end+1) = 1;
0496         marksPos(end+1) = length(g.indexout);
0497         marksPos = sort(marksPos);
0498 
0499         %now warp tmpall
0500         mytmpall = tmpall(:,:,t);
0501         r = sqrt(mytmpall.*conj(mytmpall));
0502         theta = angle(mytmpall);
0503 
0504         % So mytmpall is almost equal to r.*exp(i*theta)
0505         % whos marksPos refsPos
0506 
0507         M = timeWarp(marksPos, refsPos);
0508         
0509         TSr = transpose(M*r');
0510         TStheta = zeros(size(theta,1), size(theta,2));
0511 
0512         for freqInd=1:size(TStheta,1)
0513             TStheta(freqInd, :) = angTimeWarp(marksPos, refsPos, theta(freqInd, :));            
0514         end
0515         TStmpall = TSr.*exp(i*TStheta);
0516 
0517         % $$$     keyboard;
0518 
0519         tmpall(:,:,t) =  TStmpall;
0520     end
0521 end
0522 %time-warp ends
0523 
0524 % compute and subtract ITC
0525 % ------------------------
0526 if nargout > 3 || strcmpi(g.subitc, 'on')
0527     itcvals = tfitc(tmpall, g.itctype);
0528 end;
0529 if strcmpi(g.subitc, 'on')
0530     %a = gcf; figure; imagesc(abs(itcvals)); cbar; figure(a);
0531     if ndims(tmpall) <= 3
0532          tmpall = (tmpall - abs(tmpall) .* repmat(itcvals, [1 1 trials])) ./ abs(tmpall);
0533     else tmpall = (tmpall - abs(tmpall) .* repmat(itcvals, [1 1 1 trials])) ./ abs(tmpall);
0534     end;
0535 end;
0536 
0537 % find closest output frequencies
0538 % -------------------------------
0539 if length(g.freqs) ~= length(freqs) || any(g.freqs ~= freqs)
0540     allindices = zeros(1,length(g.freqs));
0541     for index = 1:length(g.freqs)
0542         [dum ind] = min(abs(freqs-g.freqs(index)));
0543         allindices(index) = ind;
0544     end;
0545     verboseprintf(g.verbose, 'finding closest frequencies: %d freqs removed\n', length(freqs)-length(allindices));
0546     freqs = freqs(allindices);
0547     if ndims(tmpall) <= 3
0548          tmpall = tmpall(allindices,:,:);
0549     else tmpall = tmpall(:,allindices,:,:);
0550     end;
0551     if nargout > 3 | strcmpi(g.subitc, 'on')
0552         if ndims(tmpall) <= 3
0553              itcvals = itcvals(allindices,:,:);
0554         else itcvals = itcvals(:,allindices,:,:);
0555         end;
0556     end;
0557 end;
0558 
0559 timesout = g.timesout;
0560 
0561 %figure; imagesc(abs(sum(itcvals,3))); cbar;
0562 return;
0563 
0564 % function for itc
0565 % ----------------
0566 function [itcvals] = tfitc(tfdecomp, itctype);
0567 % first dimension are trials
0568 nd = max(3,ndims(tfdecomp));
0569 switch itctype
0570     case 'coher',
0571         try,
0572             itcvals = sum(tfdecomp,nd) ./ sqrt(sum(tfdecomp .* conj(tfdecomp),nd) * size(tfdecomp,nd));
0573         catch, % scan rows if out of memory
0574             for index =1:size(tfdecomp,1)
0575                 itcvals(index,:,:) = sum(tfdecomp(index,:,:,:),nd) ./ sqrt(sum(tfdecomp(index,:,:,:) .* conj(tfdecomp(index,:,:,:)),nd) * size(tfdecomp,nd));
0576             end;
0577         end;
0578     case 'phasecoher2',
0579         try,
0580             itcvals = sum(tfdecomp,nd) ./ sum(sqrt(tfdecomp .* conj(tfdecomp)),nd);
0581         catch, % scan rows if out of memory
0582             for index =1:size(tfdecomp,1)
0583                 itcvals(index,:,:) = sum(tfdecomp(index,:,:,:),nd) ./ sum(sqrt(tfdecomp(index,:,:,:) .* conj(tfdecomp(index,:,:,:))),nd);
0584             end;
0585         end;
0586     case 'phasecoher',
0587         try,
0588             itcvals = sum(tfdecomp ./ sqrt(tfdecomp .* conj(tfdecomp)) ,nd) / size(tfdecomp,nd);
0589         catch, % scan rows if out of memory
0590             for index =1:size(tfdecomp,1)
0591                 itcvals(index,:,:) = sum(tfdecomp(index,:,:,:) ./ sqrt(tfdecomp(index,:,:,:) .* conj(tfdecomp(index,:,:,:))) ,nd) / size(tfdecomp,nd);
0592             end;
0593         end;
0594 end % ~any(isnan())
0595 return;
0596 
0597 function w = hanning(n)
0598 if ~rem(n,2)
0599     w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
0600     w = [w; w(end:-1:1)];
0601 else
0602     w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
0603     w = [w; w(end-1:-1:1)];
0604 end
0605 
0606 % get time points
0607 % ---------------
0608 function [ timevals, timeindices ] = gettimes(frames, tlimits, timevar, winsize, ntimevar, causal, verbose);
0609 timevect = linspace(tlimits(1), tlimits(2), frames);
0610 srate    = 1000*(frames-1)/(tlimits(2)-tlimits(1));
0611 
0612 if isempty(timevar) % no pre-defined time points
0613     if ntimevar(1) > 0
0614         % generate linearly space vector
0615         % ------------------------------
0616         if (ntimevar > frames-winsize)
0617             ntimevar = frames-winsize;
0618             if ntimevar < 0
0619                 error('Not enough data points, reduce the window size or lowest frequency');
0620             end;
0621             verboseprintf(verbose, ['Value of ''timesout'' must be <= frame-winsize, ''timesout'' adjusted to ' int2str(ntimevar) '\n']);
0622         end
0623         npoints = ntimevar(1);
0624         wintime = 500*winsize/srate;
0625         if strcmpi(causal, 'on')
0626              timevals = linspace(tlimits(1)+2*wintime, tlimits(2), npoints);
0627         else timevals = linspace(tlimits(1)+wintime, tlimits(2)-wintime, npoints);
0628         end;
0629         verboseprintf(verbose, 'Generating %d time points (%1.1f to %1.1f ms)\n', npoints, min(timevals), max(timevals));
0630     else
0631         % subsample data
0632         % --------------
0633         nsub     = -ntimevar(1);
0634         if strcmpi(causal, 'on')
0635              timeindices = [ceil(winsize+nsub):nsub:length(timevect)];
0636         else timeindices = [ceil(winsize/2+nsub/2):nsub:length(timevect)-ceil(winsize/2)-1];
0637         end;
0638         timevals    = timevect( timeindices ); % the conversion at line 741 leaves timeindices unchanged
0639         verboseprintf(verbose, 'Subsampling by %d (%1.1f to %1.1f ms)\n', nsub, min(timevals), max(timevals));
0640     end;
0641 else
0642     timevals = timevar;
0643     % check boundaries
0644     % ----------------
0645     wintime = 500*winsize/srate;
0646     if strcmpi(causal, 'on')
0647          tmpind  = find( (timevals >= tlimits(1)+2*wintime-0.0001) & (timevals <= tlimits(2)) ); 
0648     else tmpind  = find( (timevals >= tlimits(1)+wintime-0.0001) & (timevals <= tlimits(2)-wintime+0.0001) ); 
0649     end;
0650     % 0.0001 account for numerical innacuracies on opteron computers
0651     if isempty(tmpind)
0652         error('No time points. Reduce time window or minimum frequency.');
0653     end;
0654     if  length(timevals) ~= length(tmpind)
0655         verboseprintf(verbose, 'Warning: %d out of %d time values were removed (now %3.2f to %3.2f ms) so the lowest\n', ...
0656             length(timevals)-length(tmpind), length(timevals), timevals(tmpind(1)), timevals(tmpind(end)));
0657         verboseprintf(verbose, '         frequency could be computed with the requested accuracy\n');
0658     end;
0659     timevals = timevals(tmpind);
0660 end;
0661 
0662 % find closet points in data
0663 % --------------------------
0664 timeindices = round(eeg_lat2point(timevals, 1, srate, tlimits, 1E-3));
0665 if length(timeindices) < length(unique(timeindices))
0666     timeindices = unique(timeindices)
0667     verboseprintf(verbose, 'Warning: duplicate times, reduce the number of output times\n');
0668 end;
0669 if length(unique(timeindices(2:end)-timeindices(1:end-1))) > 1
0670     verboseprintf(verbose, 'Finding closest points for time variable');
0671     verboseprintf(verbose, 'Time values for time/freq decomposition is not perfectly uniformly distributed\n');
0672 else
0673     verboseprintf(verbose, 'Distribution of data point for time/freq decomposition is perfectly uniform\n');
0674 end;
0675 timevals    = timevect(timeindices);
0676 
0677 % DEPRECATED, FOR C INTERFACE
0678 function nofunction()
0679 % C PART %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0680 filename = [ 'tmpcrossf' num2str(round(rand(1)*1000)) ];
0681 f = fopen([ filename '.in'], 'w');
0682 fwrite(f, tmpsaveall, 'int32');
0683 fwrite(f, g.detret, 'int32');
0684 fwrite(f, g.srate, 'int32');
0685 fwrite(f, g.maxfreq, 'int32');
0686 fwrite(f, g.padratio, 'int32');
0687 fwrite(f, g.cycles, 'int32');
0688 fwrite(f, g.winsize, 'int32');
0689 fwrite(f, g.timesout, 'int32');
0690 fwrite(f, g.subitc, 'int32');
0691 fwrite(f, g.type, 'int32');
0692 fwrite(f, trials, 'int32');
0693 fwrite(f, g.naccu, 'int32');
0694 fwrite(f, length(X), 'int32');
0695 fwrite(f, X, 'double');
0696 fwrite(f, Y, 'double');
0697 fclose(f);
0698 
0699 command = [ '!cppcrosff ' filename '.in ' filename '.out' ];
0700 eval(command);
0701 
0702 f = fopen([ filename '.out'], 'r');
0703 size1 = fread(f, 'int32', 1);
0704 size2 = fread(f, 'int32', 1);
0705 Rreal = fread(f, 'double', [size1 size2]);
0706 Rimg  = fread(f, 'double', [size1 size2]);
0707 Coher.R = Rreal + j*Rimg;
0708 Boot.Coherboot.R = [];
0709 Boot.Rsignif = [];
0710 
0711 function verboseprintf(verbose, varargin)
0712 if strcmpi(verbose, 'on')
0713     fprintf(varargin{:});
0714 end;
0715