% newtimef() - Return estimates and plots of mean event-related (log) spectral
%           perturbation (ERSP) and inter-trial coherence (ITC) events
%           across event-related trials (epochs) of a single input time series.
%
%         * Also can compute and statistically compare transforms for two time
%           series. Use to compare ERSP and ITC means in two conditions.
%
%         * Uses either fixed-window, zero-padded FFTs (fastest), wavelet
%           0-padded DFTs. FFT uses Hanning tapers; wavelets use similar Morlet
%           tapers.
%
%         * For the wavelet and FFT methods, output frequency spacing
%           is the lowest frequency ('srate'/'winsize') divided by 'padratio'.
%           NaN input values (such as returned by eventlock()) are ignored.
%
%         * If 'alpha' is given (see below), permutation statistics are computed
%           (from a distribution of 'naccu' surrogate data trials) and
%           non-significant features of the output plots are zeroed out
%           (and plotted in green).
%
%         * Given a 'topovec' topo vector and 'elocs' electrode location file,
%           the figure also shows a topoplot() view of the specified scalp map.
%
%         * Note: Left-click on subplots to view and zoom in separate windows.
%
% Usage with single dataset:
%        >> [ersp,itc,powbase,times,freqs,erspboot,itcboot] = ...
%                  newtimef(data, frames, tlimits, srate, cycles,...
%                       'key1',value1, 'key2',value2, ... );
%
% Example to compare two condition (channel 1 EEG versus ALLEEG(2)):
%        >> [ersp,itc,powbase,times,freqs,erspboot,itcboot] = ...
%                  newtimef({EEG.data(1,:,:) ALLEEG(2).data(1,:,:)}, ,,,
%                       EEG.pnts, [EEG.xmin EEG.xmax]*1000, EEG.srate, cycles);
% NOTE:
%        >> timef details  % presents more detailed argument information
%                          % Note: timef() also computes multitaper transforms
%
% Required inputs:    Value                                 {default}
%       data        = Single-channel data vector (1,frames*ntrials), else 
%                     2-D array (frames,trials) or 3-D array (1,frames,trials).
%                     To compare two conditions (data1 and data2), in place of 
%                     a single data matrix enter a cell array {data1 data2}
%       frames      = Frames per trial. Ignored if data is 2-D or 3-D.  {750}
%       tlimits     = [mintime maxtime] (ms). Epoch time limits default is 
%                     [-1000 2000]. Note that this is the time limit of the 
%                     data itself NOT A SUB-WINDOW TO EXTRACT FROM EXISTING
%                     EPOCHS as this is the case in pop_newtimef.
%       srate       = data sampling rate (Hz)  {default: from icadefs.m or 250}
%       '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.
%
%    Optional inter-trial coherence (ITC) Type:
%       'itctype'   = ['coher'|'phasecoher'|'phasecoher2'] Compute either linear
%                     coherence ('coher') or phase coherence ('phasecoher'),
%                     first called 'phase-locking factor' {default: 'phasecoher'}
%
%    Optional detrending:
%       'detrend'   = ['on'|'off'], Linearly detrend each data epoch   {'off'}
%       'rmerp'     = ['on'|'off'], Remove epoch mean from data epochs {'off'}
%
%    Optional FFT/DFT parameters:
%       'winsize'   = If cycles==0: data subwindow length (fastest, 2^n<frames);
%                     If cycles >0: *longest* window length to use. This
%                      determines the lowest output frequency. Note that this
%                     parameter is overwritten if the minimum frequency requires
%                     a longer time window {default: ~frames/8}
%       'timesout'  = Number of output times (int<frames-winframes). Enter a
%                     negative value [-S] to subsample original times by S.
%                     Enter an array to obtain spectral decomposition at
%                     specific times (Note: algorithm finds the closest time
%                     point in data; this could give a slightly unevenly spaced
%                     time array                                    {default: 200}
%       'padratio'  = FFT-length/winframes (2^k)                    {default: 2}
%                     Multiplies the number of output frequencies by dividing
%                     their spacing (standard FFT padding). When cycles~=0,
%                     frequency spacing is divided by padratio.
%       'maxfreq'   = Maximum frequency (Hz) to plot (& to output, if cycles>0)
%                     If cycles==0, all FFT frequencies are output. {default: 50}
%                     DEPRECATED, use 'freqs' instead,and never both.
%       'freqs'     = [min max] frequency limits. {default [minfreq 50],
%                     minfreq being determined by the number of data points,
%                     cycles and sampling frequency.
%       '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: 'linear'}
%                     Note that for obtaining 'log' spaced freqs using FFT,
%                     closest correspondant frequencies in the 'linear' space
%                     are returned.
%       'lowmem'    = ['on'|'off'] compute frequency by frequency to save
%                     memory. {default: 'off'}
%       'verbose'   = ['on'|'off'] print text {'on'}
%       'subitc'    = ['on'|'off'] subtract stimulus locked Inter-Trial Coherence
%                     (ITC) from x and y. This computes an 'intrinsic' coherence
%                     of x and y not arising directly from common phase locking 
%                     to experimental events. See notes.    {default: 'off'}
%       'wletmethod' = ['dftfilt'|'dftfilt2'|'dftfilt3'] Wavelet type to use.
%                     'dftfilt2' -> Morlet-variant wavelets, or Hanning DFT.
%                     'dftfilt3' -> Morlet wavelets.
%                     See the timefreq function for more detials {default: 'dftfilt3'}
%       'cycleinc'    ['linear'|'log'] increase mode if [min max] cycles is
%                     provided in 'cycle' parameter. Applies only to 
%                     'wletmethod','dftfilt'  {default: 'linear'}
%       
%   Optional baseline parameters:
%       'baseline'  = Spectral baseline end-time (in ms). NaN imply that no
%                     baseline is used. A range [min max] may also be entered
%                     You may also enter one row per region for baseline
%                     e.g. [0 100; 300 400] considers the window 0 to 100 ms and
%                     300 to 400 ms {default: 0 -> all negative time values}. 
%                     This parameter is valid to define all baseline types below.
%                     Once more "NaN" prevents baseline subtraction.
%       'powbase'   = Baseline spectrum to log-subtract {default|NaN -> from data}
%       'commonbase' = ['on'|'off'] use common baseline when comparing two 
%                     conditions (default is 'on').
%       'basenorm'  = ['on'|'off'] 'on' normalize baseline in the power spectral
%                     average instead of 'off' dividing by the average power across 
%                     trials at each frequency (gain model). Default: 'off'.
%       'trialbase' = ['on'|'off'] perform baseline (normalization or division 
%                     above in single trial instead of the trial average. Default
%                     if 'off'. 
%
%    Optional time warping parameter: 
%       'timewarp'  = [eventms matrix] Time-warp amplitude and phase time-
%                     courses(following time/freq transform but before 
%                     smoothing across trials). 'eventms' is a matrix 
%                     of size (all_trials,epoch_events) whose columns
%                     specify the epoch latencies (in ms) at which the
%                     same series of successive events occur in each 
%                     trial. If two data conditions, eventms should be 
%                     [eventms1;eventms2] = all trials stacked vertically.
%      'timewarpms' = [warpms] optional vector of event latencies (in ms) 
%                     to which the series of events should be warped.
%                     (Note: Epoch start and end should not be declared
%                     as eventms or warpms}. If 'warpms' is absent or [], 
%                     the median of each 'eventms' column will be used;
%                     If two datasets, the grand medians of the two.
%     'timewarpidx' = [plotidx] is an vector of indices telling which of 
%                     the time-warped 'eventms' columns (above) to plot with 
%                     vertical lines. If undefined, all columns are plotted. 
%                     Overwrites the 'vert' argument (below) if any.
%
%    Optional permutation parameters:
%       'alpha'     = If non-0, compute two-tailed permutation significance 
%                      probability level. Show non-signif. output values 
%                      as green.                                     {0}
%       'mcorrect'  = ['none'|'fdr'] correction for multiple comparison
%                     'fdr' uses false detection rate (see function fdr()).
%                     Default is 'none'. Not available for condition 
%                     comparisons.
%       'pcontour'  = ['on'|'off'] draw contour around significant regions
%                     instead of masking them. Default is 'off'. Not 
%                     available for condition comparisons.
%       'naccu'     = Number of permutation replications to accumulate {200}
%       'baseboot'  = permutation baseline subtract (1 -> use 'baseline';
%                                                  0 -> use whole trial
%                                          [min max] -> use time range) {1}
%                     You may also enter one row per region for baseline
%                     e.g. [0 100; 300 400] considers the window 0 to 100 ms 
%                     and 300 to 400 ms.
%       'boottype'  = ['shuffle'|'rand'|'randall'] 'shuffle' -> shuffle times 
%                     and trials; 'rand' -> invert polarity of spectral data 
%                     (for ERSP) or randomize phase (for ITC); 'randall' -> 
%                     compute significances by accumulating random-polarity 
%                     inversions for each time/frequency point (slow!). Note
%                     that in previous revision of this function, we use to 
%                     call this method bootstrap while it is really permutation
%                     per say {'shuffle'}
%       'condboot'  = ['abs'|'angle'|'complex'] for comparing 2 conditions,
%                     either subtract ITC absolute vales ('abs'), angles
%                     ('angles') or complex values ('complex').   {'abs'}
%       'pboot'     = permutation power limits (e.g., from newtimef()) {from data}
%       'rboot'     = permutation ITC limits (e.g., from newtimef()). 
%                     Note: Both pboot and rboot must be provided to avoid 
%                     recomputing the surrogate data!           {from data}
%
%    Optional Scalp Map:
%       'topovec'   = Scalp topography (map) to plot              {none}
%       'elocs'     = Electrode location file for scalp map       {none}
%                     Value should be a string array containing the path
%                     and name of the file.  For file format, see
%                         >> topoplot example
%       'chaninfo'    Passed to topoplot, if called.
%                     [struct] optional structure containing fields 
%                     'nosedir', 'plotrad', and/or 'chantype'. See these 
%                     field definitions above, below.
%                     {default: nosedir +X, plotrad 0.5, all channels}
%
%     Optional Plotting Parameters:
%       'scale'     = ['log'|'abs'] visualize power in log scale (dB) or absolute
%                     scale. Default is 'log'.
%       'plottype'  = ['image'|'curve'] plot time frequency images or
%                     curves (one curve per frequency). {default: 'image'}
%       'plotmean'  = ['on'|'off'] For 'curve' plots only. Average all
%                     frequencies given as input. {default: 'on'}
%       'highlightmode'  = ['background'|'bottom'] For 'curve' plots only,
%                     display significant time regions either in the plot background
%                     or underneath the curve.
%       'plotersp'  = ['on'|'off'] Plot power spectral perturbations    {'on'}
%       'plotitc'   = ['on'|'off'] Plot inter trial coherence           {'on'}
%       'plotphasesign' = ['on'|'off'] Plot phase sign in the inter trial coherence {'on'}
%       'plotphaseonly' = ['on'|'off'] Plot ITC phase instead of ITC amplitude {'off'}
%       'erspmax'   = [real] set the ERSP max. For the color scale (min= -max) {auto}
%       'itcmax'    = [real] set the ITC image maximum for the color scale {auto}
%       'hzdir'     = ['up' or 'normal'|'down' or 'reverse'] Direction of
%                     the frequency axes {as in icadefs.m, or 'up'}
%       'ydir'      = ['up' or 'normal'|'down' or 'reverse'] Direction of
%                     the ERP axis plotted below the ITC {as in icadefs.m, or 'up'}
%       'erplim'    = [min max] ERP limits for ITC (below ITC image)       {auto}
%       'itcavglim' = [min max] average ITC limits for all freq. (left of ITC) {auto}
%       'speclim'   = [min max] average spectrum limits (left of ERSP image)   {auto}
%       'erspmarglim' = [min max] average marginal ERSP limits (below ERSP image) {auto}
%       'title'     = Optional figure or (brief) title {none}. For multiple conditions
%                     this must contain a cell array of 2 or 3 title strings.
%       'marktimes' = Non-0 times to mark with a dotted vertical line (ms)     {none}
%       'linewidth' = Line width for 'marktimes' traces (thick=2, thin=1)      {2}
%       'axesfont'  = Axes text font size                                      {10}
%       'titlefont' = Title text font size                                     {8}
%       'vert'      = [times_vector] -> plot vertical dashed lines at specified times
%                     in ms. {default: none}
%       'newfig'    = ['on'|'off'] Create new figure for difference plots {'on'}
%       'outputformat' = ['old'|'plot'] for compatibility with script that used the 
%                        old output format, set to 'old' (mbase in absolute amplitude (not
%                        dB) and real itc instead of complex itc). 'plot' returns
%                        the plotted result {default: 'plot'}
% Outputs:
%            ersp   = (nfreqs,timesout) matrix of log spectral diffs from baseline
%                     (in dB log scale or absolute scale). Use the 'plot' output format
%                     above to output the ERSP as shown on the plot.
%            itc    = (nfreqs,timesout) matrix of complex inter-trial coherencies.
%                     itc is complex -- ITC magnitude is abs(itc); ITC phase in radians
%                     is angle(itc), or in deg phase(itc)*180/pi.
%          powbase  = baseline power spectrum. Note that even, when selecting the 
%                     the 'trialbase' option, the average power spectrum is
%                     returned (not trial based). To obtain the baseline of
%                     each trial, recompute it manually using the tfdata
%                     output described below.
%            times  = vector of output times (spectral time window centers) (in ms).
%            freqs  = vector of frequency bin centers (in Hz).
%         erspboot  = (nfreqs,2) matrix of [lower upper] ERSP significance.
%          itcboot  = (nfreqs) matrix of [upper] abs(itc) threshold.
%           tfdata  = optional (nfreqs,timesout,trials) time/frequency decomposition 
%                      of the single data trials. Values are complex.
%
% Plot description:
%   Assuming both 'plotersp' and 'plotitc' options are 'on' (= default). 
%   The upper panel presents the data ERSP (Event-Related Spectral Perturbation) 
%   in dB, with mean baseline spectral activity (in dB) subtracted. Use 
%   "'baseline', NaN" to prevent timef() from removing the baseline. 
%   The lower panel presents the data ITC (Inter-Trial Coherence). 
%   Click on any plot axes to pop up a new window (using 'axcopy()')
%   -- Upper left marginal panel presents the mean spectrum during the baseline 
%      period (blue), and when significance is set, the significance threshold 
%      at each frequency (dotted green-black trace).
%   -- The marginal panel under the ERSP image shows the maximum (green) and 
%      minimum (blue) ERSP values relative to baseline power at each frequency.
%   -- The lower left marginal panel shows mean ITC across the imaged time range 
%      (blue), and when significance is set, the significance threshold (dotted 
%      green-black).  
%   -- The marginal panel under the ITC image shows the ERP (which is produced by 
%      ITC across the data spectral pass band).
%
% Authors: Arnaud Delorme, Sigurd Enghoff, Jean Hausser, && Scott Makeig
%          CNL / Salk Institute 1998- | SCCN/INC, UCSD 2002-
%
% See also: timefreq(), condstat(), newcrossf(), tftopo()

%    Deprecated Multitaper Parameters: [not included here]
%       'mtaper'    = If [N W], performs multitaper decomposition.
%                      (N is the time resolution and W the frequency resolution;
%                      maximum taper number is 2NW-1). Overwrites 'winsize' and 'padratio'.
%                     If [N W K], forces the use of K Slepian tapers (if possible).
%                      Phase is calculated using standard methods.
%                      The use of mutitaper with wavelets (cycles>0) is not
%                      recommended (as multiwavelets are not implemented).
%                      Uses Matlab functions DPSS, PMTM.      {no multitaper}

%    Deprecated time warp keywords (working?)
%      'timewarpfr' = {{[events], [warpfr], [plotidx]}} Time warp amplitude and phase
%                     time-courses (after time/freq transform but before smoothingtimefreqfunc
%                     across trials). 'events' is a matrix whose columns specify the
%                     epoch frames [1 ... end] at which a series of successive events
%                     occur in each trial. 'warpfr' is an optional vector of event
%                     frames to which the series of events should be time locked.
%                     (Note: Epoch start and end should not be declared as events or
%                     warpfr}. If 'warpfr' is absent or [], the median of each 'events'
%                     column will be used. [plotidx] is an optional vector of indices
%                     telling which of the warpfr to plot with vertical lines. If
%                     undefined, all marks are plotted. Overwrites 'vert' argument,
%                     if any. [Note: In future releases, 'timewarpfr' will be deprecated
%                     in favor of 'timewarp' using latencies in ms instead of frames].

%    Deprecated original time warp keywords (working?)
%       'timeStretchMarks' = [(marks,trials) matrix] Each trial data will be
%                     linearly warped (after time/freq. transform) so that the
%                     event marks are time locked to the reference frames
%                     (see timeStretchRefs). Marks must be specified in frames
%       'timeStretchRefs' = [1 x marks] Common reference frames to all trials.
%                     If empty or undefined, median latency for each mark will be used.boottype
%       'timeStretchPlot' = [vector] Indicates the indices of the reference frames
%                     (in StretchRefs) should be overplotted on the ERSP and ITC.
%
%
% Copyright (C) 1998- Arnaud Delorme, Scott Makeig
% first built as timef.m at CNL / Salk Institute 8/1/98-8/28/01 by
% Sigurd Enghoff and Scott Makeig, edited by Arnaud Delorme
% SCCN/INC/UCSD/ reprogrammed as newtimef -Arnaud Delorme 2002-
% SCCN/INC/UCSD/ added time warping capabilities -Jean Hausser 2005
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 10-19-98 avoided division by zero (using MIN_ABS) -sm
% 10-19-98 improved usage message and commandline info printing -sm
% 10-19-98 made valid [] values for tvec and g.elocs -sm
% 04-01-99 added missing freq in freqs and plots, fixed log scaling bug -se && -tpj
% 06-29-99 fixed frequency indexing for constant-Q -se
% 08-24-99 reworked to handle NaN input values -sm
% 12-07-99 adjusted ERPtimes to plot ERP under ITC -sm
% 12-22-99 debugged ERPtimes, added BASE_BOOT -sm
% 01-10-00 debugged BASE_BOOT=0 -sm
% 02-28-00 added NOTE on formula derivation below -sm
% 03-16-00 added axcopy() feature -sm && tpj
% 04-16-00 added multiple marktimes loop -sm
% 04-20-00 fixed ITC cbar limits when spcified in input -sm
% 07-29-00 changed frequencies displayed msg -sm
% 10-12-00 fixed bug in freqs when cycles>0 -sm
% 02-07-01 fixed inconsistency in BASE_BOOT use -sm
% 08-28-01 matlab 'key' value arguments -ad
% 08-28-01 multitaper decomposition -ad
% 01-25-02 reformated help && license -ad
% 03-08-02 debug && compare to old timef function -ad
% 03-16-02 timeout automatically adjusted if too high -ad
% 04-02-02 added 'coher' option -ad

function [P,R,mbase,timesout,freqs,Pboot,Rboot,alltfX,PA] = newtimef( data, frames, tlimits, Fs, varwin, varargin);

% Note: Above, PA is output of 'phsamp','on'

% For future 'timewarp' keyword help: 'timewarp' 3rd element {colors} contains a
%               list of Matlab linestyles to use for vertical lines marking the occurence
%               of the time warped events. If '', no line will be drawn for this event
%               column. If fewer colors than event columns, cycles through the given color
%               labels.  Note: Not compatible with 'vert' (below).

%varwin,winsize,g.timesout,g.padratio,g.maxfreq,g.topovec,g.elocs,g.alpha,g.marktimes,g.powbase,g.pboot,g.rboot)

% ITC:   Normally, R = |Sum(Pxy)| / (Sum(|Pxx|)*Sum(|Pyy|)) is coherence.
%        But here, we consider    Phase(Pyy) = 0 and |Pyy| = 1 -> Pxy = Pxx
%        Giving, R = |Sum(Pxx)|/Sum(|Pxx|), the inter-trial coherence (ITC)
%        Also called 'phase-locking factor' by Tallon-Baudry et al. (1996)

if nargin < 1
    help newtimef;
    return;
end;

% Read system (or directory) constants and preferences:
% ------------------------------------------------------
icadefs % read local EEGLAB constants: HZDIR, YDIR, DEFAULT_SRATE, DEFAULT_TIMLIM

if ~exist('HZDIR'), HZDIR = 'up'; end; % ascending freqs
if ~exist('YDIR'), YDIR = 'up'; end;   % positive up

if YDIR == 1, YDIR = 'up'; end;        % convert from [-1|1] as set in icadefs.m  
if YDIR == -1, YDIR = 'down'; end;     % and read by other plotting functions

if ~exist('DEFAULT_SRATE'), DEFAULT_SRATE = 250; end;            % 250 Hz
if ~exist('DEFAULT_TIMLIM'), DEFAULT_TIMLIM = [-1000 2000]; end; % [-1 2] s epochs

% Constants set here:
% ------------------
ERSP_CAXIS_LIMIT = 0;           % 0 -> use data limits; else positive value
% giving symmetric +/- caxis limits.
ITC_CAXIS_LIMIT  = 0;           % 0 -> use data limits; else positive value
% giving symmetric +/- caxis limits.
MIN_ABS          = 1e-8;        % avoid division by ~zero

% Command line argument defaults:
% ------------------------------
DEFAULT_NWIN	= 200;		% Number of windows = horizontal resolution
DEFAULT_VARWIN	= 0;		% Fixed window length or fixed number of cycles.
% =0: fix window length to that determined by nwin
% >0: set window length equal to varwin cycles
%     Bounded above by winsize, which determines
%     the min. freq. to be computed.

DEFAULT_OVERSMP	= 2;		% Number of times to oversample frequencies
DEFAULT_MAXFREQ = 50;		% Maximum frequency to display (Hz)
DEFAULT_TITLE	= '';		% Figure title (no default)
DEFAULT_ELOC    = 'chan.locs';	% Channel location file
DEFAULT_ALPHA   = NaN;		% Percentile of bins to keep
DEFAULT_MARKTIME= NaN;

% Font sizes:
AXES_FONT       = 10;           % axes text FontSize
TITLE_FONT      =  8;

if (nargin < 2)
    frames = floor((DEFAULT_TIMLIN(2)-DEFAULT_TIMLIM(1))/DEFAULT_SRATE);
elseif (~isnumeric(frames) | length(frames)~=1 | frames~=round(frames))
    error('Value of frames must be an integer.');
elseif (frames <= 0)
    error('Value of frames must be positive.');
end;

DEFAULT_WINSIZE = max(pow2(nextpow2(frames)-3),4);
DEFAULT_PAD = max(pow2(nextpow2(DEFAULT_WINSIZE)),4);

if (nargin < 1)
    help newtimef
    return
end

if isstr(data) && strcmp(data,'details')
    more on
    help timefdetails
    more off
    return
end
if ~iscell(data)
    data = reshape_data(data, frames);
    trials = size(data,ndims(data));
else
    if ndims(data) == 3 && size(data,1) == 1
        error('Cannot process multiple channel component in compare mode');
    end;
    [data{1}, frames] = reshape_data(data{1}, frames);
    [data{2}, frames] = reshape_data(data{2}, frames);
    trials = size(data{1},2);
end;

if (nargin < 3)
    tlimits = DEFAULT_TIMLIM;
elseif (~isnumeric(tlimits) | sum(size(tlimits))~=3)
    error('Value of tlimits must be a vector containing two numbers.');
elseif (tlimits(1) >= tlimits(2))
    error('tlimits interval must be ascending.');
end

if (nargin < 4)
    Fs = DEFAULT_SRATE;
elseif (~isnumeric(Fs) | length(Fs)~=1)
    error('Value of srate must be a number.');
elseif (Fs <= 0)
    error('Value of srate must be positive.');
end

if (nargin < 5)
    varwin = DEFAULT_VARWIN;
elseif ~isnumeric(varwin) && strcmpi(varwin, 'cycles')
    varwin = varargin{1};
    varargin(1) = [];
elseif (varwin < 0)
    error('Value of cycles must be zero or positive.');
end

% build a structure for keyword arguments
% --------------------------------------
if ~isempty(varargin)
    [tmp indices] = unique(varargin(1:2:end));
    varargin = varargin(sort(union(indices*2-1, indices*2))); % these 2 lines remove duplicate arguments
    try, g = struct(varargin{:});
    catch, error('Argument error in the {''param'', value} sequence'); end;
end
%}
[ g timefreqopts ] = finputcheck(varargin, ...
    {'boottype'      'string'    {'shuffle','rand','randall'}    'shuffle'; ...
    'condboot'      'string'    {'abs','angle','complex'}       'abs'; ...
    'title'         { 'string' 'cell' }   { [] [] }         DEFAULT_TITLE; ...
    'title2'        'string'    []          DEFAULT_TITLE; ...
    'winsize'       'integer'      [0 Inf]  DEFAULT_WINSIZE; ...
    'pad'           'real'      []          DEFAULT_PAD; ...
    'timesout'      'integer'   []          DEFAULT_NWIN; ...
    'padratio'      'integer'   [0 Inf]     DEFAULT_OVERSMP; ...
    'topovec'       'real'      []          []; ...
    'elocs'         {'string','struct'} []  DEFAULT_ELOC; ...
    'alpha'         'real'      [0 0.5]     DEFAULT_ALPHA; ...
    'marktimes'     'real'      []          DEFAULT_MARKTIME; ...
    'powbase'       'real'      []          NaN; ...
    'pboot'         'real'      []          NaN; ...
    'rboot'         'real'      []          NaN; ...
    'plotersp'      'string'    {'on','off'} 'on'; ...
    'plotamp'       'string'    {'on','off'} 'on'; ...
    'plotitc'       'string'    {'on','off'} 'on'; ...
    'detrend'       'string'    {'on','off'} 'off'; ...
    'rmerp'         'string'    {'on','off'} 'off'; ...
    'basenorm'      'string'    {'on','off'} 'off'; ...
    'commonbase'    'string'    {'on','off'} 'on'; ...
    'baseline'      'real'      []           0; ...
    'baseboot'      'real'      []           1; ...
    'linewidth'     'integer'   [1 2]        2; ...
    'naccu'         'integer'   [1 Inf]      200; ...
    'mtaper'        'real'      []           []; ...
    'maxfreq'       'real'      [0 Inf]      DEFAULT_MAXFREQ; ...
    'freqs'         'real'      [0 Inf]      [0 DEFAULT_MAXFREQ]; ...
    'cycles'        'integer'   []           []; ...
    'nfreqs'        'integer'   []           []; ...
    'freqscale'     'string'    []           'linear'; ...
    'vert'          'real'      []           [];  ...
    'newfig'        'string'    {'on','off'} 'on'; ...
    'type'          'string'    {'coher','phasecoher','phasecoher2'}  'phasecoher'; ...
    'itctype'       'string'    {'coher','phasecoher','phasecoher2'}  'phasecoher'; ...
    'phsamp'        'string'    {'on','off'} 'off'; ...  % phsamp not completed - Toby 9.28.2006
    'plotphaseonly' 'string'    {'on','off'} 'off'; ...
    'plotphasesign' 'string'    {'on','off'} 'on'; ...
    'plotphase'     'string'    {'on','off'} 'on'; ... % same as above for backward compatibility
    'pcontour'      'string'    {'on','off'} 'off'; ... 
    'outputformat'  'string'    {'old','new' 'plot' } 'plot'; ...
    'itcmax'        'real'      []           []; ...
    'erspmax'       'real'      []           []; ...
    'lowmem'        'string'    {'on','off'} 'off'; ...
    'verbose'       'string'    {'on','off'} 'on'; ...
    'plottype'      'string'    {'image','curve'}   'image'; ...
    'mcorrect'      'string'    {'fdr','none'}      'none'; ...
    'plotmean'      'string'    {'on','off'} 'on'; ...
    'highlightmode' 'string'    {'background','bottom'}     'background'; ...
    'chaninfo'      'struct'    []           struct([]); ...
    'erspmarglim'   'real'      []           []; ...
    'itcavglim'     'real'      []           []; ...
    'erplim'        'real'      []           []; ...
    'speclim'       'real'      []           []; ...
    'ntimesout'     'real'      []           []; ...
    'scale'         'string'    { 'log' 'abs'} 'log'; ...
    'timewarp'      'real'      []           []; ...
    'precomputed'   'struct'    []           struct([]); ...
    'timewarpms'    'real'      []           []; ...
    'timewarpfr'    'real'      []           []; ...
    'timewarpidx'   'real'      []           []; ...
    'timewarpidx'   'real'      []           []; ...
    'timeStretchMarks'  'real'  []           []; ...
    'timeStretchRefs'   'real'  []           []; ...
    'timeStretchPlot'   'real'  []           []; ...
    'trialbase'     'string'    {'on','off'} 'off'; 
    'hzdir'         'string'    {'up','down','normal','reverse'}   HZDIR; ...
    'ydir'          'string'    {'up','down','normal','reverse'}   YDIR; ...
    'cycleinc'      'string'   {'linear' 'log'}        'linear'
    }, 'newtimef', 'ignore');
if isstr(g), error(g); end;
if strcmpi(g.plotamp, 'off'), g.plotersp = 'off'; end;    
if strcmpi(g.basenorm, 'on'), g.scale = 'abs'; end;
if ~strcmpi(g.itctype , 'phasecoher'), g.type = g.itctype; end;

g.tlimits = tlimits;
g.frames  = frames;
g.srate   = Fs;
if isempty(g.cycles)
    g.cycles  = varwin;
end;
g.AXES_FONT        = AXES_FONT;      % axes text FontSize
g.TITLE_FONT       = TITLE_FONT;
g.ERSP_CAXIS_LIMIT = ERSP_CAXIS_LIMIT;
g.ITC_CAXIS_LIMIT  = ITC_CAXIS_LIMIT;
if ~strcmpi(g.plotphase, 'on'), g.plotphasesign = g.plotphase; end;

% unpack 'timewarp' (and undocumented 'timewarpfr') arguments
%------------------------------------------------------------
if isfield(g,'timewarpfr')
    if iscell(g.timewarpfr) && length(g.timewarpfr) > 3
        error('undocumented ''timewarpfr'' cell array may have at most 3 elements');
    end
end

if isempty(g.nfreqs)
    verboseprintf(g.verbose, 'Warning: ''nfreqs'' input overwrite ''padratio''\n');
end;
if strcmpi(g.basenorm, 'on')
    verboseprintf(g.verbose, 'Baseline normalization is on (results will be shown as z-scores)\n');
end;

if isfield(g,'timewarp') && ~isempty(g.timewarp)
    if ndims(data) == 3
        error('Cannot perform time warping on 3-D data input');
    end;
    if ~isempty(g.timewarp) % convert timewarp ms to timewarpfr frames -sm
        fprintf('\n')
        if iscell(g.timewarp)
           error('timewarp argument must be a (total_trials,epoch_events) matrix');
        end
        evntms = g.timewarp;
        warpfr = round((evntms - g.tlimits(1))/1000*g.srate)+1;
        g.timewarpfr{1} = warpfr';

        if isfield(g,'timewarpms')
           refms = g.timewarpms;
           reffr = round((refms - g.tlimits(1))/1000*g.srate)+1;
           g.timewarpfr{2} = reffr';
        end
        if isfield(g,'timewarpidx')
           g.timewarpfr{3} = g.timewarpidx;
        end
    end

    % convert again to timeStretch parameters
    % ---------------------------------------
    if ~isempty(g.timewarpfr)
        g.timeStretchMarks = g.timewarpfr{1};
        if length(g.timewarpfr) > 1
            g.timeStretchRefs = g.timewarpfr{2};
        end

        if length(g.timewarpfr) > 2
          if isempty(g.timewarpfr{3})
            stretchevents = size(g.timeStretchMarks,1);
            g.timeStretchPlot = [1:stretchevents]; % default to plotting all lines
          else
            g.timeStretchPlot = g.timewarpfr{3};
          end
        end

        if max(max(g.timeStretchMarks)) > frames-2 | min(min(g.timeStretchMarks)) < 3
            error('Time warping events must be inside the epochs.');
        end
        if ~isempty(g.timeStretchRefs)
            if max(g.timeStretchRefs) > frames-2 | min(g.timeStretchRefs) < 3
                error('Time warping reference latencies must be within the epochs.');
            end
        end
    end
end

% test argument consistency
% --------------------------
if g.tlimits(2)-g.tlimits(1) < 30
    verboseprintf(g.verbose, 'newtimef(): WARNING: Specified time range is very small (< 30 ms)???\n');
    verboseprintf(g.verbose, '                     Epoch time limits should be in msec, not seconds!\n');
end

if (g.winsize > g.frames)
    error('Value of winsize must be smaller than epoch frames.');
end

if length(g.timesout) == 1 && g.timesout > 0
    if g.timesout > g.frames-g.winsize
        g.timesout = g.frames-g.winsize;
        disp(['Value of timesout must be <= frames-winsize, timeout adjusted to ' int2str(g.timesout) ]);
    end
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

if (g.maxfreq > Fs/2)
    verboseprintf(g.verbose, ['Warning: value of maxfreq reduced to Nyquist rate' ...
        ' (%3.2f)\n\n'], Fs/2);
    g.maxfreq = Fs/2;
end
if g.maxfreq ~= DEFAULT_MAXFREQ, g.freqs(2) = g.maxfreq; end;

if isempty(g.topovec)
    g.topovec = [];
    if isempty(g.elocs)
        error('Channel location file must be specified.');
    end;
end

if (round(g.naccu*g.alpha) < 2)
    verboseprintf(g.verbose, 'Value of alpha is outside its normal range [%g,0.5]\n',2/g.naccu);
    g.naccu = round(2/g.alpha);
    verboseprintf(g.verbose, '  Increasing the number of iterations to %d\n',g.naccu);
end

if ~isnan(g.alpha)
    if length(g.baseboot) == 2
        verboseprintf(g.verbose, 'Permutation analysis will use data from %3.2g to %3.2g ms.\n', ...
            g.baseboot(1),  g.baseboot(2))
    elseif g.baseboot > 0
        verboseprintf(g.verbose, 'Permutation analysis will use data in (pre-0) baseline subwindows only.\n')
    else
        verboseprintf(g.verbose, 'Permutation analysis will use data in all subwindows.\n')
    end
end

if ~isempty(g.timeStretchMarks) % timeStretch code by Jean Hauser
    if isempty(g.timeStretchRefs)
        verboseprintf(g.verbose, ['Using median event latencies as reference event times for time warping.\n']);
        g.timeStretchRefs = median(g.timeStretchMarks,2); 
                                          % Note: Uses (grand) median latencies for two conditions
    else
        verboseprintf(g.verbose, ['Using supplied latencies as reference event times for time warping.\n']);
    end
    if isempty(g.timeStretchPlot)
        verboseprintf(g.verbose, 'Will not overplot the reference event times on the ERSP.\n');
    elseif length(g.timeStretchPlot) > 0
        g.vert = ((g.timeStretchRefs(g.timeStretchPlot)-1) ...
            /g.srate+g.tlimits(1)/1000)*1000;
        fprintf('Plotting timewarp markers at ')
           for li = 1:length(g.vert), fprintf('%d ',g.vert(li)); end
        fprintf(' ms.\n')
    end
end 

if min(g.vert) < g.tlimits(1) | max(g.vert) > g.tlimits(2)
    error('vertical line (''vert'') latency outside of epoch boundaries');
end

if strcmp(g.hzdir,'up')| strcmp(g.hzdir,'normal')
    g.hzdir = 'normal'; % convert to Matlab graphics constants
elseif strcmp(g.hzdir,'down') | strcmp(g.hzdir,'reverse')| g.hzdir==-1
    g.hzdir = 'reverse';
else
    error('unknown ''hzdir'' argument'); 
end

if strcmp(g.ydir,'up')| strcmp(g.ydir,'normal')
    g.ydir = 'normal'; % convert to Matlab graphics constants
elseif strcmp(g.ydir,'down') | strcmp(g.ydir,'reverse')
    g.ydir = 'reverse';
else
    error('unknown ''ydir'' argument'); 
end

% -----------------
% ERSP scaling unit
% -----------------
if strcmpi(g.scale, 'log')
    if strcmpi(g.basenorm, 'on')
        g.unitpower = '10*log(std.)'; % impossible
    elseif isnan(g.baseline)
        g.unitpower = '10*log10(\muV^{2}/Hz)';
    else
        g.unitpower = 'dB';
    end;
else
    if strcmpi(g.basenorm, 'on')
        g.unitpower = 'std.';
    elseif isnan(g.baseline)
        g.unitpower = '\muV^{2}/Hz';
    else
        g.unitpower = '% of baseline';
    end;
end;

% Multitaper - used in timef
% --------------------------
if ~isempty(g.mtaper) % multitaper, inspired from a Bijan Pesaran matlab function
    if length(g.mtaper) < 3
        %error('mtaper arguement must be [N W] or [N W K]');

        if g.mtaper(1) * g.mtaper(2) < 1
            error('mtaper 2 first arguments'' product must be larger than 1');
        end;
        if length(g.mtaper) == 2
            g.mtaper(3) = floor( 2*g.mtaper(2)*g.mtaper(1) - 1);
        end
        if length(g.mtaper) == 3
            if g.mtaper(3) > 2 * g.mtaper(1) * g.mtaper(2) -1
                error('mtaper number too high (maximum (2*N*W-1))');
            end;
        end
        disp(['Using ' num2str(g.mtaper(3)) ' tapers.']);
        NW = g.mtaper(1)*g.mtaper(2);   % product NW
        N  = g.mtaper(1)*g.srate;
        [e,v] = dpss(N, NW, 'calc');
        e=e(:,1:g.mtaper(3));
        g.alltapers = e;
    else
        g.alltapers = g.mtaper;
        disp('mtaper argument not [N W] or [N W K]; considering raw taper matrix');
    end;
    g.winsize = size(g.alltapers, 1);
    g.pad = max(pow2(nextpow2(g.winsize)),256); % pad*nextpow
    nfk = floor([0 g.maxfreq]./g.srate.*g.pad);
    g.padratio = 2*nfk(2)/g.winsize;

    %compute number of frequencies
    %nf = max(256, g.pad*2^nextpow2(g.winsize+1));
    %nfk = floor([0 g.maxfreq]./g.srate.*nf);

    %freqs = linspace( 0, g.maxfreq, diff(nfk)); % this also works in the case of a FFT
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute frequency by frequency if low memory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmpi(g.lowmem, 'on') && length(data) ~= g.frames && isempty(g.nfreqs) && ~iscell(data)

    % compute for first 2 trials to get freqsout
    XX = reshape(data, 1, frames, prod(size(data))/g.frames);
    [P,R,mbase,timesout,freqsout] = newtimef(XX(1,:,1), frames, tlimits, Fs, varwin, 'plotitc', 'off', 'plotamp', 'off',varargin{:}, 'lowmem', 'off');

    % scan all frequencies
    for index = 1:length(freqsout)
        if nargout < 8
            [P(index,:),R(index,:),mbase(index),timesout,tmpfreqs(index),Pboot(index,:),Rboot(index,:)] = ...
                newtimef(data, frames, tlimits, Fs, varwin, ...
                          'freqs', [freqsout(index) freqsout(index)], 'nfreqs', 1, ...
                             'plotamp', 'off', 'plotphasesign', 'off',varargin{:}, ...
                                  'lowmem', 'off', 'timesout', timesout);
        else
            [P(index,:),R(index,:),mbase(index),timesout,tmpfreqs(index),Pboot(index,:),Rboot(index,:), ...
                alltfX(index,:,:)] = ...
                newtimef(data, frames, tlimits, Fs, varwin, ...
                            'freqs', [freqsout(index) freqsout(index)], 'nfreqs', 1, ...
                                  'plotamp', 'off', 'plotphasesign', 'off',varargin{:}, ...
                                          'lowmem', 'off', 'timesout', timesout);
        end;
    end;

    % compute trial-average ERP 
    % -------------------------
    ERP = mean(data,2);

    % plot results 
    %-------------
    plottimef(P, R, Pboot, Rboot, ERP, freqsout, timesout, mbase, [], [], g);

    return; % finished
end;


%%%%%%%%%%%%%%%%%%%%%%%
% compare 2 conditions 
%%%%%%%%%%%%%%%%%%%%%%%
if iscell(data)
    Pboot = [];
    Rboot = [];
    if ~strcmpi(g.mcorrect, 'none')
        error('Correction for multiple comparison not implemented for comparing conditions');
    end;
    
    vararginori = varargin;
    if length(data) ~= 2
        error('newtimef: to compare two conditions, data must be a length-2 cell array');
    end;

    % deal with titles
    % ----------------
    for index = 1:2:length(vararginori)
        if index<=length(vararginori) % needed if elements are deleted

          %  if      strcmp(vararginori{index}, 'title') | ... % Added by Jean Hauser
          %          strcmp(vararginori{index}, 'title2') | ...
                 if strcmp(vararginori{index}, 'timeStretchMarks') | ...
                    strcmp(vararginori{index}, 'timeStretchRefs') | ...
                    strcmp(vararginori{index}, 'timeStretchPlots')
                vararginori(index:index+1) = [];
            end;
        end;
    end;
    if iscell(g.title) && length(g.title) >= 2 % Changed that part because providing titles
                                              % as cells caused the function to crash (why?) 
                                              % at line 704 (g.tlimits = tlimits) -Jean
        if length(g.title) == 2,
             g.title{3} = [ g.title{1} ' - '  g.title{2} ];
        end;
    else
        disp('Warning: title must be a cell array');
        g.title = { 'Condition 1' 'Condition 2' 'Condition 1 minus Condition 2' };
    end;

    verboseprintf(g.verbose, '\nRunning newtimef() on Condition 1 **********************\n\n');

    verboseprintf(g.verbose, 'Note: If an out-of-memory error occurs, try reducing the\n');
    verboseprintf(g.verbose, '      the number of time points or number of frequencies\n');
    verboseprintf(g.verbose, '(''coher'' options take 3 times the memory of other options)\n\n');

    cond_1_epochs = size(data{1},2);

    if ~isempty(g.timeStretchMarks)
        [P1,R1,mbase1,timesout,freqs,Pboot1,Rboot1,alltfX1] = ...
        newtimef( data{1}, frames, tlimits, Fs, varwin, 'plotitc', 'off', ...
          'plotersp', 'off', vararginori{:}, 'lowmem', 'off', ...
            'timeStretchMarks', g.timeStretchMarks(:,1:cond_1_epochs), ... 
              'timeStretchRefs', g.timeStretchRefs);
    else
        [P1,R1,mbase1,timesout,freqs,Pboot1,Rboot1,alltfX1] = ...
        newtimef( data{1}, frames, tlimits, Fs, varwin, 'plotitc', 'off', ...
          'plotersp', 'off', vararginori{:}, 'lowmem', 'off');
    end

    verboseprintf(g.verbose,'\nRunning newtimef() on Condition 2 **********************\n\n');

    [P2,R2,mbase2,timesout,freqs,Pboot2,Rboot2,alltfX2] = ...
        newtimef( data{2}, frames, tlimits, Fs, varwin, 'plotitc', 'off', ...
          'plotersp', 'off', vararginori{:}, 'lowmem', 'off', ...
            'timeStretchMarks', g.timeStretchMarks(:,cond_1_epochs+1:end), ... 
              'timeStretchRefs', g.timeStretchRefs);

    verboseprintf(g.verbose,'\nComputing difference **********************\n\n');

    % recompute power baselines 
    % -------------------------
    if ~isnan( g.baseline(1) ) && ~isnan( mbase1(1) ) && isnan(g.powbase(1)) && strcmpi(g.commonbase, 'on')
        disp('Recomputing baseline power: using the grand mean of both conditions ...');
        mbase = (mbase1 + mbase2)/2;
        P1 = P1 + repmat(mbase1(1:size(P1,1))',[1 size(P1,2)]);
        P2 = P2 + repmat(mbase2(1:size(P1,1))',[1 size(P1,2)]);
        P1 = P1 - repmat(mbase (1:size(P1,1))',[1 size(P1,2)]);
        P2 = P2 - repmat(mbase (1:size(P1,1))',[1 size(P1,2)]);
        if ~isnan(g.alpha)
            Pboot1 = Pboot1 + repmat(mbase1(1:size(Pboot1,1))',[1 size(Pboot1,2) size(Pboot1,3)]);
            Pboot2 = Pboot2 + repmat(mbase2(1:size(Pboot1,1))',[1 size(Pboot1,2) size(Pboot1,3)]);
            Pboot1 = Pboot1 - repmat(mbase (1:size(Pboot1,1))',[1 size(Pboot1,2) size(Pboot1,3)]);
            Pboot2 = Pboot2 - repmat(mbase (1:size(Pboot1,1))',[1 size(Pboot1,2) size(Pboot1,3)]);
        end;
        verboseprintf(g.verbose, '\nSubtracting the common power baseline ...\n');
        meanmbase = mbase;
        mbase = { mbase mbase };
    elseif strcmpi(g.commonbase, 'on')
        mbase = { NaN NaN };
    else 
        meanmbase = (mbase1 + mbase2)/2;
        mbase = { mbase1 mbase2 };
    end;

    % plotting
    % --------
    if strcmpi(g.plotersp, 'on') | strcmpi(g.plotitc, 'on')
        g.titleall = g.title;
        if strcmpi(g.newfig, 'on'), figure; end; % declare a new figure

        % using same color scale
        % ----------------------
        if ~isfield(g, 'erspmax')
            g.erspmax = max( max(max(abs(Pboot1))), max(max(abs(Pboot2))) );
        end;
        if ~isfield(g, 'itcmax')
            g.itcmax  = max( max(max(abs(Rboot1))), max(max(abs(Rboot2))) );
        end;

        subplot(1,3,1); % plot Condition 1
        g.title = g.titleall{1};
        g = plottimef(P1, R1, Pboot1, Rboot1, mean(data{1},2), freqs, timesout, mbase{1}, [], [], g);
        g.itcavglim = [];

        subplot(1,3,2); % plot Condition 2
        g.title = g.titleall{2};
        plottimef(P2, R2, Pboot2, Rboot2, mean(data{2},2), freqs, timesout, mbase{2}, [], [], g);

        subplot(1,3,3); % plot Condition 1 - Condition 2
        g.title =  g.titleall{3};
    end;

    if isnan(g.alpha)
        switch(g.condboot)
            case 'abs',  Rdiff = abs(R1)-abs(R2);
            case 'angle',  Rdiff = angle(R1)-angle(R2);
            case 'complex',  Rdiff = R1-R2;
        end;
        if strcmpi(g.plotersp, 'on') | strcmpi(g.plotitc, 'on')
            plottimef(P1-P2, Rdiff, [], [], mean(data{1},2)-mean(data{2},2), freqs, timesout, meanmbase, [], [], g);
        end;
    else
        % preprocess data and run compstat() function
        % -------------------------------------------
        alltfX1power = alltfX1.*conj(alltfX1);
        alltfX2power = alltfX2.*conj(alltfX2);

        if ~isnan(mbase{1}(1))
            mbase1 = 10.^(mbase{1}(1:size(alltfX1,1))'/20);
            mbase2 = 10.^(mbase{2}(1:size(alltfX1,1))'/20);
            alltfX1 = alltfX1./repmat(mbase1/2,[1 size(alltfX1,2) size(alltfX1,3)]);
            alltfX2 = alltfX2./repmat(mbase2/2,[1 size(alltfX2,2) size(alltfX2,3)]);
            alltfX1power = alltfX1power./repmat(mbase1,[1 size(alltfX1power,2) size(alltfX1power,3)]);
            alltfX2power = alltfX2power./repmat(mbase2,[1 size(alltfX2power,2) size(alltfX2power,3)]);
        end;
        
        %formula = {'log10(mean(arg1,3))'};              % toby 10.02.2006
        %formula = {'log10(mean(arg1(:,:,data),3))'};

        formula = {'log10(mean(arg1(:,:,X),3))'};
        switch g.type
            case 'coher', % take the square of alltfx and alltfy first to speed up
                formula = { formula{1} ['sum(arg2(:,:,data),3)./sqrt(sum(arg1(:,:,data),3)*length(data) )'] };
                if strcmpi(g.lowmem, 'on')
                    for ind = 1:2:size(alltfX1power,1)
                        if ind == size(alltfX1,1), indarr = ind; else indarr = [ind:ind+1]; end;
                        [resdifftmp resimagestmp res1tmp res2tmp] = ...
                            condstat(formula, g.naccu, g.alpha, {'both' 'upper'}, { '' g.condboot}, ...
                            { alltfX1power(indarr,:,:) alltfX2power(indarr,:,:) }, {alltfX1(indarr,:,:) alltfX2(indarr,:,:)});
                        resdiff{1}(indarr,:)     = resdifftmp{1};   resdiff{2}(indarr,:)     = resdifftmp{2};
                        resimages{1}(indarr,:,:) = resimagestmp{1}; resimages{2}(indarr,:,:) = resimagestmp{2};
                        res1{1}(indarr,:)        = res1tmp{1};      res1{2}(indarr,:)        = res1tmp{2};
                        res2{1}(indarr,:)        = res2tmp{1};      res2{2}(indarr,:)        = res2tmp{2};
                    end;
                else
                    alltfXpower = { alltfX1power alltfX2power };
                    alltfX      = { alltfX1 alltfX2 };
                    alltfXabs   = { alltfX1abs alltfX2abs };
                    [resdiff resimages res1 res2] = condstat(formula, g.naccu, g.alpha, {'both' 'upper'}, { '' g.condboot}, alltfXpower, alltfX, alltfXabs);
                end;
            case 'phasecoher2', % normalize first to speed up

                %formula = { formula{1} ['sum(arg2(:,:,data),3)./sum(arg3(:,:,data),3)'] }; 
                % toby 10/3/2006

                formula = { formula{1} ['sum(arg2(:,:,X),3)./sum(arg3(:,:,X),3)'] };
                alltfX1abs = sqrt(alltfX1power); % these 2 lines can be suppressed
                alltfX2abs = sqrt(alltfX2power); % by inserting sqrt(arg1(:,:,data)) instead of arg3(:,:,data))
                if strcmpi(g.lowmem, 'on')
                    for ind = 1:2:size(alltfX1abs,1)
                        if ind == size(alltfX1,1), indarr = ind; else indarr = [ind:ind+1]; end;
                        [resdifftmp resimagestmp res1tmp res2tmp] = ...
                            condstat(formula, g.naccu, g.alpha, {'both' 'upper'}, { '' g.condboot}, ...
                            { alltfX1power(indarr,:,:) alltfX2power(indarr,:,:) }, {alltfX1(indarr,:,:) ...
                            alltfX2(indarr,:,:)}, { alltfX1abs(indarr,:,:) alltfX2abs(indarr,:,:) });
                        resdiff{1}(indarr,:)     = resdifftmp{1};   resdiff{2}(indarr,:)     = resdifftmp{2};
                        resimages{1}(indarr,:,:) = resimagestmp{1}; resimages{2}(indarr,:,:) = resimagestmp{2};
                        res1{1}(indarr,:)        = res1tmp{1};      res1{2}(indarr,:)        = res1tmp{2};
                        res2{1}(indarr,:)        = res2tmp{1};      res2{2}(indarr,:)        = res2tmp{2};
                    end;
                else
                    alltfXpower = { alltfX1power alltfX2power };
                    alltfX      = { alltfX1 alltfX2 };
                    alltfXabs   = { alltfX1abs alltfX2abs };
                    [resdiff resimages res1 res2] = condstat(formula, g.naccu, g.alpha, {'both' 'upper'}, { '' g.condboot}, alltfXpower, alltfX, alltfXabs);
                end;
            case 'phasecoher',

                %formula = { formula{1} ['mean(arg2,3)'] };              % toby 10.02.2006
                %formula = { formula{1} ['mean(arg2(:,:,data),3)'] };

                formula = { formula{1} ['mean(arg2(:,:,X),3)'] };
                if strcmpi(g.lowmem, 'on')
                    for ind = 1:2:size(alltfX1,1)
                        if ind == size(alltfX1,1), indarr = ind; else indarr = [ind:ind+1]; end;
                        alltfX1norm = alltfX1(indarr,:,:)./sqrt(alltfX1(indarr,:,:).*conj(alltfX1(indarr,:,:)));
                        alltfX2norm = alltfX2(indarr,:,:)./sqrt(alltfX2(indarr,:,:).*conj(alltfX2(indarr,:,:)));
                        alltfXpower = { alltfX1power(indarr,:,:) alltfX2power(indarr,:,:) };
                        alltfXnorm  = { alltfX1norm alltfX2norm };
                        [resdifftmp resimagestmp res1tmp res2tmp] = ...
                            condstat(formula, g.naccu, g.alpha, {'both' 'both'}, { '' g.condboot}, ...
                            alltfXpower, alltfXnorm);
                        resdiff{1}(indarr,:)     = resdifftmp{1};   resdiff{2}(indarr,:)     = resdifftmp{2};
                        resimages{1}(indarr,:,:) = resimagestmp{1}; resimages{2}(indarr,:,:) = resimagestmp{2};
                        res1{1}(indarr,:)        = res1tmp{1};      res1{2}(indarr,:)        = res1tmp{2};
                        res2{1}(indarr,:)        = res2tmp{1};      res2{2}(indarr,:)        = res2tmp{2};
                    end;
                else
                    alltfX1norm = alltfX1./sqrt(alltfX1.*conj(alltfX1));
                    alltfX2norm = alltfX2./sqrt(alltfX2.*conj(alltfX2)); % maybe have to suppress preprocessing -> lot of memory
                    alltfXpower = { alltfX1power alltfX2power };
                    alltfXnorm  = { alltfX1norm alltfX2norm };
                    [resdiff resimages res1 res2] = condstat(formula, g.naccu, g.alpha, {'both' 'both'}, { '' g.condboot}, ...
                        alltfXpower, alltfXnorm);
		end;
        end;

        % same as below: plottimef(P1-P2, R2-R1, 10*resimages{1}, resimages{2}, mean(data{1},2)-mean(data{2},2), freqs, times, mbase, g);
        if strcmpi(g.plotersp, 'on') | strcmpi(g.plotitc, 'on')
            g.erspmax = []; % auto scale
            g.itcmax  = []; % auto scale
            plottimef(10*resdiff{1}, resdiff{2}, 10*resimages{1}, resimages{2}, ...
                mean(data{1},2)-mean(data{2},2), freqs, timesout, meanmbase, [], [], g);
        end;
        R1 = res1{2};
        R2 = res2{2};
        Rdiff = resdiff{2};
        Pboot = { Pboot1 Pboot2 10*resimages{1} };
        Rboot = { Rboot1 Rboot2 resimages{2} };
    end;
    P = { P1 P2 P1-P2 };
    R = { R1 R2 Rdiff };

    if nargout >= 8, alltfX = { alltfX1 alltfX2 }; end;

    return; % ********************************** END FOR MULTIPLE CONDITIONS
end;

%%%%%%%%%%%%%%%%%%%%%%
% display text to user (computation perfomed only for display)
%%%%%%%%%%%%%%%%%%%%%%
verboseprintf(g.verbose, 'Computing Event-Related Spectral Perturbation (ERSP) and\n');
switch g.type
    case 'phasecoher',  verboseprintf(g.verbose, '  Inter-Trial Phase Coherence (ITC) images based on %d trials\n',trials);
    case 'phasecoher2', verboseprintf(g.verbose, '  Inter-Trial Phase Coherence 2 (ITC) images based on %d trials\n',trials);
    case 'coher',       verboseprintf(g.verbose, '  Linear Inter-Trial Coherence (ITC) images based on %d trials\n',trials);
end;
verboseprintf(g.verbose, '  of %d frames sampled at %g Hz.\n',g.frames,g.srate);
verboseprintf(g.verbose, 'Each trial contains samples from %1.0f ms before to\n',g.tlimits(1));
verboseprintf(g.verbose, '  %1.0f ms after the timelocking event.\n',g.tlimits(2));
if ~isnan(g.alpha)
    verboseprintf(g.verbose, 'Only significant values (permutation statistics p<%g) will be colored;\n',g.alpha)
    verboseprintf(g.verbose, '  non-significant values will be plotted in green\n');
end
verboseprintf(g.verbose,'  Image frequency direction: %s\n',g.hzdir);

if isempty(g.precomputed)
    % -----------------------------------------
    % detrend over epochs (trials) if requested
    % -----------------------------------------
    if strcmpi(g.rmerp, 'on')
        if ndims(data) == 2
             data = data - mean(data,2)*ones(1, length(data(:))/g.frames);
        else data = data - repmat(mean(data,3), [1 1 trials]);
        end;
    end;

    % ----------------------------------------------------
    % compute time frequency decompositions, power and ITC
    % ----------------------------------------------------
    if length(g.timesout) > 1,   tmioutopt = { 'timesout' , g.timesout };
    elseif ~isempty(g.ntimesout) tmioutopt = { 'ntimesout', g.ntimesout };
    else                         tmioutopt = { 'ntimesout', g.timesout };
    end;

    [alltfX freqs timesout R] = timefreq(data, g.srate, tmioutopt{:}, ...
        'winsize', g.winsize, 'tlimits', g.tlimits, 'detrend', g.detrend, ...
        'itctype', g.type, 'wavelet', g.cycles, 'verbose', g.verbose, ...
        'padratio', g.padratio, 'freqs', g.freqs, 'freqscale', g.freqscale, ...
        'nfreqs', g.nfreqs, 'timestretch', {g.timeStretchMarks', g.timeStretchRefs}, timefreqopts{:});
else
    alltfX   = g.precomputed.tfdata;
    timesout = g.precomputed.times;
    freqs    = g.precomputed.freqs;
    if strcmpi(g.precomputed.recompute, 'ersp')
        R = [];
    else
        switch g.itctype
            case 'coher',       R = alltfX ./ repmat(sqrt(sum(alltfX .* conj(alltfX),3) * size(alltfX,3)), [1 1 size(alltfX,3)]);
            case 'phasecoher2', R = alltfX ./ repmat(sum(sqrt(alltfX .* conj(alltfX)),3), [1 1 size(alltfX,3)]);
            case 'phasecoher',  R = alltfX ./ sqrt(alltfX .* conj(alltfX));
        end;
        P = []; mbase = []; return;
    end;
end;

if g.cycles(1) == 0
    alltfX = 2*0.375*alltfX/g.winsize; % TF and MC (12/11/2006): normalization, divide by g.winsize
    P  = alltfX.*conj(alltfX); % power    
    % TF and MC (12/14/2006): multiply by 2 account for negative frequencies,
    % and ounteract the reduction by a factor 0.375 that occurs as a result of 
    % cosine (Hann) tapering. Refer to Bug 446
else 
    P  = alltfX.*conj(alltfX); % power for wavelets
end;

% ---------------
% baseline length
% ---------------
if size(g.baseline,2) == 2
    baseln = [];
    for index = 1:size(g.baseline,1)
        tmptime   = find(timesout >= g.baseline(index,1) & timesout <= g.baseline(index,2));
        baseln = union(baseln, tmptime);
    end;
    if length(baseln)==0
        error('No point found in baseline');
    end
else
    if ~isempty(find(timesout < g.baseline))
         baseln = find(timesout < g.baseline); % subtract means of pre-0 (centered) windows
    else baseln = 1:length(timesout); % use all times as baseline
    end
end;

if ~isnan(g.alpha) && length(baseln)==0
    verboseprintf(g.verbose, 'timef(): no window centers in baseline (times<%g) - shorten (max) window length.\n', g.baseline)
    return
end

% -----------------------------------------
% remove baseline on a trial by trial basis
% -----------------------------------------
if ndims(P) == 4
    Pori = mean(P, 4);
    if strcmpi(g.trialbase, 'on') && isnan( g.powbase(1) )
        mbase = mean(P(:,:,baseln,:),3);
        if strcmpi(g.basenorm, 'on')
             mstd = std(P(:,:,baseln,:),[],3);
             P = bsxfun(@rdivide, bsxfun(@minus, P, mbase), mstd);
        else P = bsxfun(@rdivide, P, mbase);
        end;
    end;
else
    Pori = mean(P, 3);
    if strcmpi(g.trialbase, 'on') && isnan( g.powbase(1) )
        mbase = mean(P(:,baseln,:),2);
        if strcmpi(g.basenorm, 'on')
            mstd = std(P(:,baseln,:),[],2);
            P = (P-repmat(mbase,[1 size(P,2) 1]))./repmat(mstd,[1 size(P,2) 1]); % convert to log then back to normal
        else
            P = P./repmat(mbase,[1 size(P,2) 1]); 
            %P = 10 .^ (log10(P) - repmat(log10(mbase),[1 size(P,2) 1])); % same as above
        end;
    end;
end;
if ~isempty(g.precomputed)
    return; % return single trial power
end;

% ----------------
% phase amp option
% ----------------
if strcmpi(g.phsamp, 'on')
    disp( 'phsamp option is deprecated');
    %  switch g.phsamp
    %  case 'on'
    %PA = zeros(size(P,1),size(P,1),g.timesout); % NB: (freqs,freqs,times)
    % $$$ end                                             %       phs   amp
    %PA (freq x freq x time)
    %PA(:,:,j) = PA(:,:,j)  + (tmpX ./ abs(tmpX)) * ((P(:,j)))';
    % x-product: unit phase column
    % times amplitude row

    %tmpcx(1,:,:) = cumulX; % allow ./ below
    %for jj=1:g.timesout
    %    PA(:,:,jj) = PA(:,:,jj) ./ repmat(P(:,jj)', [size(P,1) 1]);
    %end
end

% ---------
% bootstrap
% --------- % this ensures that if bootstrap limits provided that no
% 'alpha' won't prevent application of the provided limits
if ~isnan(g.alpha) | ~isempty(find(~isnan(g.pboot))) | ~isempty(find(~isnan(g.rboot)))% if bootstrap analysis requested . . .
    
    % ERSP bootstrap
    % --------------
    if ~isempty(find(~isnan(g.pboot))) % if ERSP bootstrap limits provided already
        Pboot = g.pboot(:);
    else
        if size(g.baseboot,2) == 1
            if g.baseboot == 0, baselntmp = [];
            elseif ~isnan(g.baseline(1))
                baselntmp = baseln;
            else baselntmp = find(timesout <= 0); % if it is empty use whole epoch
            end;
        else
            baselntmp = [];
            for index = 1:size(g.baseboot,1)
                tmptime   = find(timesout >= g.baseboot(index,1) && timesout <= g.baseboot(index,2));
                if isempty(tmptime),
                    fprintf('Warning: empty baseline interval [%3.2f %3.2f]\n', g.baseboot(index,1), g.baseboot(index,2));
                end;
                baselntmp = union(baselntmp, tmptime);
            end;
        end;
        if prod(size(g.baseboot)) > 2
            fprintf('Permutation statistics will use data in multiple selected windows.\n');
        elseif size(g.baseboot,2) == 2
            fprintf('Permutation statistics will use data in range %3.2g-%3.2g ms.\n', g.baseboot(1),  g.baseboot(2));
        elseif g.baseboot
            fprintf('   %d permutation statistics windows in baseline (times<%g).\n', length(baselntmp), g.baseboot)
        end;
        
        % power significance
        % ------------------
        formula = 'mean(arg1,3);';
        [ Pboot Pboottrialstmp Pboottrials] = bootstat(P, formula, 'boottype', 'shuffle', ...
            'label', 'ERSP', 'bootside', 'both', 'naccu', g.naccu, ...
            'basevect', baselntmp, 'alpha', g.alpha, 'dimaccu', 2 );
        clear Pboottrialstmp;
        
        if size(Pboot,2) == 1, Pboot = Pboot'; end;
    end;
    
    % ITC bootstrap
    % -------------
    if ~isempty(find(~isnan(g.rboot))) % if itc bootstrap provided
        Rboot = g.rboot;
    else
        if ~isempty(find(~isnan(g.pboot))) % if ERSP limits were provided (but ITC not)
            if size(g.baseboot,2) == 1
                if g.baseboot == 0, baselntmp = [];
                elseif ~isnan(g.baseline(1))
                    baselntmp = baseln;
                else baselntmp = find(timesout <= 0); % if it is empty use whole epoch
                end;
            else
                baselntmp = [];
                for index = 1:size(g.baseboot,1)
                    tmptime   = find(timesout >= g.baseboot(index,1) && timesout <= g.baseboot(index,2));
                    if isempty(tmptime),
                        fprintf('Warning: empty baseline interval [%3.2f %3.2f]\n', g.baseboot(index,1), g.baseboot(index,2));
                    end;
                    baselntmp = union(baselntmp, tmptime);
                end;
            end;
            if prod(size(g.baseboot)) > 2
                fprintf('Permutation statistics will use data in multiple selected windows.\n');
            elseif size(g.baseboot,2) == 2
                fprintf('Permutation statistics will use data in range %3.2g-%3.2g ms.\n', g.baseboot(1),  g.baseboot(2));
            elseif g.baseboot
                fprintf('   %d permutation statistics windows in baseline (times<%g).\n', length(baselntmp), g.baseboot)
            end;
        end;        
        % ITC significance
        % ----------------
        inputdata = alltfX;
        switch g.type
            case 'coher',       formula = [ 'sum(arg1,3)./sqrt(sum(arg1.*conj(arg1),3))/ sqrt(' int2str(size(alltfX,3)) ');' ];
            case 'phasecoher',  formula = [ 'mean(arg1,3);' ]; inputdata = alltfX./sqrt(alltfX.*conj(alltfX));
            case 'phasecoher2', formula = [ 'sum(arg1,3)./sum(sqrt(arg1.*conj(arg1)),3);' ];
        end;
        if strcmpi(g.boottype, 'randall'), dimaccu = []; g.boottype = 'rand';
        else										 dimaccu = 2;
        end;
        [Rboot Rboottmp Rboottrials] = bootstat(inputdata, formula, 'boottype', g.boottype, ...
            'label', 'ITC', 'bootside', 'upper', 'naccu', g.naccu, ...
            'basevect', baselntmp, 'alpha', g.alpha, 'dimaccu', 2 );
        fprintf('\n');
        clear Rboottmp;        
    end;
else
    Pboot = []; Rboot = [];
end

% average the power
% -----------------
if ndims(P) == 4, P = mean(P, 4);
else              P = mean(P, 3);
end;

% correction for multiple comparisons
% -----------------------------------
maskersp = [];
maskitc  = []; 
if ~isnan(g.alpha)
    if isempty(find(~isnan(g.pboot))) % if ERSP lims not provided
        exactp_ersp = compute_pvals(P, Pboottrials');
        if strcmpi(g.mcorrect, 'fdr')
            alphafdr = fdr(exactp_ersp, g.alpha);
            if alphafdr ~= 0
                fprintf('ERSP correction for multiple comparisons using FDR, alpha_fdr = %3.6f\n', alphafdr);
            else fprintf('ERSP correction for multiple comparisons using FDR, nothing significant\n', alphafdr);
            end;
            maskersp = exactp_ersp <= alphafdr;
        else
            maskersp = exactp_ersp <= g.alpha;
        end;
    end;    
    if isempty(find(~isnan(g.rboot))) % if ITC lims not provided
        exactp_itc  = compute_pvals(abs(R), abs(Rboottrials'));        
        if strcmpi(g.mcorrect, 'fdr')
            alphafdr = fdr(exactp_itc, g.alpha);
            if alphafdr ~= 0
                fprintf('ITC  correction for multiple comparisons using FDR, alpha_fdr = %3.6f\n', alphafdr);
            else fprintf('ITC  correction for multiple comparisons using FDR, nothing significant\n', alphafdr);
            end;
            maskitc = exactp_itc <= alphafdr;
        else
            maskitc = exactp_itc  <= g.alpha;
        end
    end;
end;

% computer baseline values
% ------------------------
if isnan(g.powbase(1))
    verboseprintf(g.verbose, 'Computing the mean baseline spectrum\n');
    if ndims(P) == 3
         mbase    = mean(P(:,:,baseln),3);
         mbaseori = mean(Pori(:,:,baseln),3);
    else mbase    = mean(P(:,baseln),2);
         mbaseori = mean(Pori(:,baseln),2);
    end;
else
    verboseprintf(g.verbose, 'Using the input baseline spectrum\n');
    mbase    = g.powbase; 
    if strcmpi(g.scale, 'log'), mbase = 10.^(mbase/10); end; 
    if size(mbase,1) == 1 % if input was a row vector, flip to be a column
        mbase = mbase';
    end;
    mbaseori = mbase;
end
baselength = length(baseln);

% ---------------
% remove baseline
% ---------------
% original ERSP baseline removal
if ~isnan( g.baseline(1) ) && any(~isnan( mbase(1) )) && strcmpi(g.trialbase, 'off') && strcmpi(g.basenorm, 'off')
    
    P = bsxfun(@rdivide, P, mbase);
    if ~isempty(Pboot)
         Pboot = bsxfun(@rdivide, Pboot, mbase);
    end;   
% ERSP baseline normalized
elseif ~isnan( g.baseline(1) ) && ~isnan( mbase(1) ) && strcmpi(g.basenorm, 'on') && strcmpi(g.trialbase, 'off')

    if ndims(P) == 3, 
         mstd = std(P(:,:,baseln),[],3);
    else mstd = std(P(:,baseln),[],2);
    end;
    P = bsxfun(@rdivide, bsxfun(@minus, P, mbase), mstd);
    if ~isempty(Pboot) && isnan(g.pboot)
        Pboot = bsxfun(@rdivide, bsxfun(@minus, Pboot, mbase), mstd);
    end;
end;
    
% convert to log if necessary
% ---------------------------
if strcmpi(g.scale, 'log')
    if ~isnan( g.baseline(1) ) && ~isnan( mbase(1) ) && strcmpi(g.trialbase, 'off'), mbase = log10(mbase)*10; end;
    P = 10 * log10(P);
    if ~isempty(Pboot)
        Pboot = 10 * log10(Pboot);
    end;
end;

% auto scalling
% -------------
if isempty(g.erspmax)
    g.erspmax = [max(max(abs(P)))/2];
    %g.erspmax = [-g.erspmax g.erspmax]+1;
end;

% --------
% plotting
% --------
if strcmpi(g.plotersp, 'on') || strcmpi(g.plotitc, 'on')
    if ndims(P) == 3
        P = squeeze(P(2,:,:,:));
        R = squeeze(R(2,:,:,:));
        mbase = squeeze(mbase(2,:));
        ERP = mean(squeeze(data(1,:,:)),2);
    else      
        ERP = mean(data,2);
    end;
    if strcmpi(g.plottype, 'image')
        plottimef(P, R, Pboot, Rboot, ERP, freqs, timesout, mbase, maskersp, maskitc, g);
    else
        plotallcurves(P, R, Pboot, Rboot, ERP, freqs, timesout, mbase, g);
    end;
end;

% --------------
% format outputs
% --------------
if strcmpi(g.outputformat, 'old')
    R = abs(R); % convert coherence vector to magnitude
    if strcmpi(g.scale, 'log'), mbase = 10^(mbase/10); end;
end;
if strcmpi(g.verbose, 'on')
    disp('Note: Add output variables to command line call in history to');
    disp('      retrieve results and use the tftopo function to replot them');
end;
mbase = mbase';
return;

% -----------------
% plotting function
% -----------------
function g = plottimef(P, R, Pboot, Rboot, ERP, freqs, times, mbase, maskersp, maskitc, g);
%
% compute ERP
%
ERPtimes = [g.tlimits(1):(g.tlimits(2)-g.tlimits(1))/(g.frames-1):g.tlimits(2)+0.000001];
ERPindices = zeros(1, length(times));
for ti=1:length(times)
    [tmp ERPindices(ti)] = min(abs(ERPtimes-times(ti)));
end
ERPtimes = ERPtimes(ERPindices); % subset of ERP frames on t/f window centers
ERP = ERP(ERPindices);

if ~isreal(R)
    Rangle = angle(R);
    Rsign = sign(imag(R));
    R = abs(R); % convert coherence vector to magnitude
    setylim = 1;
else
    Rsign = ones(size(R));
    setylim = 0;
end;
switch lower(g.plotitc)
    case 'on',
        switch lower(g.plotersp),
            case 'on', ordinate1 = 0.67; ordinate2 = 0.1; height = 0.33; g.plot = 1;
            case 'off', ordinate2 = 0.1; height = 0.9; g.plot = 1;
        end;
    case 'off', ordinate1 = 0.1; height = 0.9;
        switch lower(g.plotersp),
            case 'on', ordinate1 = 0.1; height = 0.9;  g.plot = 1;
            case 'off', g.plot = 0;
        end;
end;

if g.plot
    % verboseprintf(g.verbose, '\nNow plotting...\n');
    set(gcf,'DefaultAxesFontSize',g.AXES_FONT)
    colormap(jet(256));
    pos = get(gca,'position');
    q = [pos(1) pos(2) 0 0];
    s = [pos(3) pos(4) pos(3) pos(4)];
    axis off;
end;

switch lower(g.plotersp)
    case 'on'
        %
        %%%%%%% image the ERSP %%%%%%%%%%%%%%%%%%%%%%%%%%
        %

        h(1) = axes('Position',[.1 ordinate1 .9 height].*s+q);
        set(h(1), 'tag', 'ersp');

        PP = P;
        if strcmpi(g.scale, 'abs') && strcmpi(g.basenorm, 'off')
             baseval = 1;
        else baseval = 0;
        end;
        if ~isnan(g.alpha)
            if strcmpi(g.pcontour, 'off') && ~isempty(maskersp) % zero out nonsignif. power differences
                PP(~maskersp) = baseval;
                %PP = PP .* maskersp;
            elseif isempty(maskersp)
                if size(PP,1) == size(Pboot,1) && size(PP,2) == size(Pboot,2)
                    PP(find(PP > Pboot(:,:,1) & (PP < Pboot(:,:,2)))) = baseval;
                    Pboot = squeeze(mean(Pboot,2));
                    if size(Pboot,2) == 1, Pboot = Pboot'; end;
                else
                    PP(find((PP > repmat(Pboot(:,1),[1 length(times)])) ...
                        & (PP < repmat(Pboot(:,2),[1 length(times)])))) = baseval;
                end
            end;
        end;
 
        % find color limits
        % -----------------
        if isempty(g.erspmax)
            if g.ERSP_CAXIS_LIMIT == 0
                g.erspmax = [-1 1]*1.1*max(max(abs(P(:,:))));
            else
                g.erspmax = g.ERSP_CAXIS_LIMIT*[-1 1];
            end
        elseif length(g.erspmax) == 1
            g.erspmax = [ -g.erspmax g.erspmax];
        end
        if isnan( g.baseline(1) ) && g.erspmax(1) < 0
            g.erspmax = [ min(min(P(:,:))) max(max(P(:,:)))];
        end;

        % plot image
        % ----------
        if ~strcmpi(g.freqscale, 'log')
            imagesc(times,freqs,PP(:,:), g.erspmax);
        else
            imagesclogy(times,freqs,PP(:,:),g.erspmax);
        end;
        set(gca,'ydir',g.hzdir);  % make frequency ascend or descend

        % put contour for multiple comparison masking
        if ~isempty(maskersp) && strcmpi(g.pcontour, 'on')
            hold on; [tmpc tmph] = contour(times, freqs, maskersp);
            set(tmph, 'linecolor', 'k', 'linewidth', 0.25)
        end;
        
        hold on
        plot([0 0],[0 freqs(end)],'--m','LineWidth',g.linewidth); % plot time 0
        if ~isnan(g.marktimes) % plot marked time
            for mt = g.marktimes(:)'
                plot([mt mt],[0 freqs(end)],'--k','LineWidth',g.linewidth);
            end
        end
        hold off
        set(h(1),'YTickLabel',[],'YTick',[])
        set(h(1),'XTickLabel',[],'XTick',[])
        if ~isempty(g.vert)
            for index = 1:length(g.vert)
                line([g.vert(index), g.vert(index)], [min(freqs) max(freqs)], 'linewidth', 1, 'color', 'm');
            end;
        end;

        h(2) = gca;
        h(3) = cbar('vert'); % ERSP colorbar axes
        set(h(2),'Position',[.1 ordinate1 .8 height].*s+q)
        set(h(3),'Position',[.95 ordinate1 .05 height].*s+q)
        title([ 'ERSP(' g.unitpower ')' ])

        %
        %%%%% plot marginal ERSP mean below ERSP image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %

        h(4) = axes('Position',[.1 ordinate1-0.1 .8 .1].*s+q);

        E = [min(P(:,:),[],1);max(P(:,:),[],1)];

        % plotting limits
        if isempty(g.erspmarglim)
            g.erspmarglim = [min(E(1,:))-max(max(abs(E)))/3 max(E(2,:))+max(max(abs(E)))/3];
        end;

        plot(times,E,[0 0],g.erspmarglim, '--m','LineWidth',g.linewidth)
        xlim([min(times) max(times)])
        ylim(g.erspmarglim)

        tick = get(h(4),'YTick');
        set(h(4),'YTick',[tick(1) ; tick(end)])
        set(h(4),'YAxisLocation','right')
        set(h(4),'TickLength',[0.020 0.025]);
        xlabel('Time (ms)')
        ylabel(g.unitpower)

        %
        %%%%% plot mean spectrum to left of ERSP image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %

        h(5) = axes('Position',[0 ordinate1 .1 height].*s+q);

        E = mbase;

        if ~isnan(E(1))

            % plotting limits
            if isempty(g.speclim)
                g.speclim = [min(E)-max(abs(E))/3 max(E)+max(abs(E))/3];
            end;

            % plot curves
            if ~strcmpi(g.freqscale, 'log')
                plot(freqs,E,'LineWidth',g.linewidth); hold on;
                if ~isnan(g.alpha)
                    try
                        plot(freqs,Pboot(:,:)'+[E;E], 'g', 'LineWidth',g.linewidth)
                        plot(freqs,Pboot(:,:)'+[E;E], 'k:','LineWidth',g.linewidth)
                    catch
                        plot(freqs,Pboot(:,:)+[E E], 'g', 'LineWidth',g.linewidth)
                        plot(freqs,Pboot(:,:)+[E E], 'k:','LineWidth',g.linewidth)
                    end;
                end
                if freqs(1) ~= freqs(end), xlim([freqs(1) freqs(end)]); end;
                ylim(g.speclim)

            else % 'log'
                semilogx(freqs,E,'LineWidth',g.linewidth); hold on;
                if ~isnan(g.alpha)
                    try
                        semilogx(freqs,Pboot(:,:)'+[E;E],'g', 'LineWidth',g.linewidth)
                        semilogx(freqs,Pboot(:,:)'+[E;E],'k:','LineWidth',g.linewidth)
                    catch
                        semilogx(freqs,Pboot(:,:)+[E E],'g', 'LineWidth',g.linewidth)
                        semilogx(freqs,Pboot(:,:)+[E E],'k:','LineWidth',g.linewidth)
                    end;
                end
                if freqs(1) ~= freqs(end), xlim([freqs(1) freqs(end)]); end;
                ylim(g.speclim)
                set(h(5),'View',[90 90])
                divs = linspace(log(freqs(1)), log(freqs(end)), 10);
                set(gca, 'xtickmode', 'manual');
                divs = ceil(exp(divs)); divs = unique(divs); % ceil is critical here, round might misalign
                set(gca, 'xtick', divs);
            end;
            set(h(5),'TickLength',[0.020 0.025]);
            set(h(5),'View',[90 90])
            xlabel('Frequency (Hz)')
            if strcmp(g.hzdir,'normal')
                set(gca,'xdir','reverse');
            else
                set(gca,'xdir','normal');
            end
            ylabel(g.unitpower)
            tick = get(h(5),'YTick');
            if (length(tick)>2)
                set(h(5),'YTick',[tick(1) ; tick(end-1)])
            end
        end;
end;

switch lower(g.plotitc)
    case 'on'
        %
        %%%%%%%%%%%% Image the ITC %%%%%%%%%%%%%%%%%%
        %
        h(6) = axes('Position',[.1 ordinate2 .9 height].*s+q); % ITC image
        set(h(1), 'tag', 'itc');

        if abs(R(1,1)-1) < 0.0001, g.plotphaseonly = 'on'; end;
        if strcmpi(g.plotphaseonly, 'on')
            RR = Rangle/pi*180;
        else
            RR = R;
        end;
        if ~isnan(g.alpha)
            if ~isempty(maskitc) && strcmpi(g.pcontour, 'off')
                RR = RR .* maskitc;
            elseif isempty(maskitc)
                if size(RR,1) == size(Rboot,1) && size(RR,2) == size(Rboot,2)
                    tmp = gcf;
                    if size(Rboot,3) == 2	 RR(find(RR > Rboot(:,:,1) & RR < Rboot(:,:,2))) = 0;
                    else                   RR(find(RR < Rboot)) = 0;
                    end;
                    Rboot = mean(Rboot(:,:,end),2);
                else
                    RR(find(RR < repmat(Rboot(:),[1 length(times)]))) = 0;
                end;
            end;
        end

        if g.ITC_CAXIS_LIMIT == 0
            coh_caxis = min(max(max(R(:,:))),1)*[-1 1]; % 1 WAS 0.4 !
        else
            coh_caxis = g.ITC_CAXIS_LIMIT*[-1 1];
        end

        if strcmpi(g.plotphaseonly, 'on')
            if ~strcmpi(g.freqscale, 'log')
                imagesc(times,freqs,RR(:,:)); % <---
            else
                imagesclogy(times,freqs,RR(:,:)); % <---
            end;
            g.itcmax = [-180 180];
            setylim = 0;
        else
            if max(coh_caxis) == 0,              % toby 10.02.2006
                coh_caxis = [-1 1];
            end
            if ~strcmpi(g.freqscale, 'log')
                if exist('Rsign') && strcmp(g.plotphasesign, 'on')
                    imagesc(times,freqs,Rsign(:,:).*RR(:,:),coh_caxis); % <---
                else
                    imagesc(times,freqs,RR(:,:),coh_caxis); % <---
                end
            else
                if exist('Rsign') && strcmp(g.plotphasesign, 'on')
                    imagesclogy(times,freqs,Rsign(:,:).*RR(:,:),coh_caxis); % <---
                else
                    imagesclogy(times,freqs,RR(:,:),coh_caxis); % <---
                end
            end;
        end;
        set(gca,'ydir',g.hzdir);  % make frequency ascend or descend

        % plot contour if necessary
        if ~isempty(maskitc) && strcmpi(g.pcontour, 'on')
            hold on; [tmpc tmph] = contour(times, freqs, maskitc);
            set(tmph, 'linecolor', 'k', 'linewidth', 0.25)
        end;

        if isempty(g.itcmax)
            g.itcmax = caxis;
        elseif length(g.itcmax) == 1
            g.itcmax = [ -g.itcmax g.itcmax ];
        end;
        caxis(g.itcmax);

        hold on
        plot([0 0],[0 freqs(end)],'--m','LineWidth',g.linewidth);
        if ~isnan(g.marktimes)
            for mt = g.marktimes(:)'
                plot([mt mt],[0 freqs(end)],'--k','LineWidth',g.linewidth);
            end
        end
        hold off
        set(h(6),'YTickLabel',[],'YTick',[])
        set(h(6),'XTickLabel',[],'XTick',[])
        if ~isempty(g.vert)
            for index = 1:length(g.vert)
                line([g.vert(index), g.vert(index)], [min(freqs) max(freqs)], 'linewidth', 1, 'color', 'm');
            end;
        end;

        h(7) = gca;
        h(8) = cbar('vert');
        %h(9) = get(h(8),'Children'); % make the function crash
        set(h(7),'Position',[.1 ordinate2 .8 height].*s+q)
        set(h(8),'Position',[.95 ordinate2 .05 height].*s+q)
        if setylim
            set(h(8),'YLim',[0 g.itcmax(2)]);
        end;
        if strcmpi(g.plotphaseonly, 'on')
            title('ITC phase')
        else
            title('ITC')
        end;

        %
        %%%%% plot the ERP below the ITC image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %

        h(10) = axes('Position',[.1 ordinate2-0.1 .8 .1].*s+q); % ERP

        if isempty(g.erplim)
            ERPmax = max(ERP);
            ERPmin = min(ERP);
            g.erplim = [ ERPmin - 0.1*(ERPmax-ERPmin) ERPmax + 0.1*(ERPmax-ERPmin) ];
        end;

        plot(ERPtimes,ERP, [0 0],g.erplim,'--m','LineWidth',g.linewidth);
        hold on;
        plot([times(1) times(length(times))],[0 0], 'k');
        xlim([min(ERPtimes) max(ERPtimes)]);
        ylim(g.erplim)
        set(gca,'ydir',g.ydir);

        tick = get(h(10),'YTick');
        set(h(10),'YTick',[tick(1) ; tick(end)])
        set(h(10),'TickLength',[0.02 0.025]);
        set(h(10),'YAxisLocation','right')
        xlabel('Time (ms)')
        ylabel('\muV')
        if (~isempty(g.topovec))
            if length(g.topovec) ~= 1, ylabel(''); end; % ICA component
        end;
        E = nan_mean(R(:,:)'); % don't let a few NaN's crash this

        %
        %%%%% plot the marginal mean left of the ITC image %%%%%%%%%%%%%%%%%%%%%
        %

        h(11) = axes('Position',[0 ordinate2 .1 height].*s+q); % plot the marginal mean
        % ITC left of the ITC image
        % set plotting limits
        if isempty(g.itcavglim)
            if ~isnan(g.alpha)
                g.itcavglim = [ min(E)-max(E)/3 max(Rboot)+max(Rboot)/3];
            else
                g.itcavglim = [ min(E)-max(E)/3 max(E)+max(E)/3];
            end;
        end;
        if max(g.itcavglim) == 0        % toby 10.02.2006
            g.itcavglim = [-1 1];
        end
        
        % plot marginal ITC
        if ~strcmpi(g.freqscale, 'log')
            plot(freqs,E,'LineWidth',g.linewidth); hold on;
            if ~isnan(g.alpha)
                plot(freqs,Rboot,'g', 'LineWidth',g.linewidth)
                plot(freqs,Rboot,'k:','LineWidth',g.linewidth)
            end
            if freqs(1) ~= freqs(end), xlim([freqs(1) freqs(end)]); end
            ylim(g.itcavglim)
        else
            semilogx(freqs,E,'LineWidth',g.linewidth); hold on;
            if ~isnan(g.alpha)
                semilogx(freqs,Rboot(:),'g', 'LineWidth',g.linewidth)
                semilogx(freqs,Rboot(:),'k:','LineWidth',g.linewidth)
            end
            if freqs(1) ~= freqs(end), xlim([freqs(1) freqs(end)]); end;
            ylim(g.itcavglim)
            divs = linspace(log(freqs(1)), log(freqs(end)), 10);
            set(gca, 'xtickmode', 'manual');
            divs = ceil(exp(divs)); divs = unique(divs); % ceil is critical here, round might misalign
            set(gca, 'xtick', divs);
         end;

        % ITC plot details
        tick = get(h(11),'YTick');
        if length(tick) > 1
            set(h(11),'YTick',[tick(1) ; tick(length(tick))])
        end;
        set(h(11),'View',[90 90])
        %set(h(11),'TickLength',[0.020 0.025]);
        xlabel('Frequency (Hz)')
        if strcmp(g.hzdir,'normal')
            set(gca,'xdir','reverse');
        else
            set(gca,'xdir','normal');
        end
        ylabel('ERP')

end; %switch

%
%%%%%%%%%%%%%%% plot a topoplot() %%%%%%%%%%%%%%%%%%%%%%%
%
if (~isempty(g.topovec)) && strcmpi(g.plotitc, 'on') && strcmpi(g.plotersp, 'on')
    
    if strcmp(g.plotersp,'off')
        h(12) = axes('Position',[-.207 .95 .2 .14].*s+q); % place the scalp map at top-left
    else
        h(12) = axes('Position',[-.1 .43 .2 .14].*s+q);   % place the scalp map at middle-left
    end;
    if length(g.topovec) == 1
        topoplot(g.topovec,g.elocs,'electrodes','off', ...
                 'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', g.chaninfo);
    else
        topoplot(g.topovec,g.elocs,'electrodes','off', 'chaninfo', g.chaninfo);
    end;
    axis('square')
end

if g.plot
    try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
    if (length(g.title) > 0) && ~iscell(g.title)
        axes('Position',pos,'Visible','Off');
        h(13) = text(-.05,1.01,g.title);
        set(h(13),'VerticalAlignment','bottom')
        set(h(13),'HorizontalAlignment','left')
        set(h(13),'FontSize',g.TITLE_FONT);
    end

    try, axcopy(gcf); catch, end;
end;

% ---------------
% Plotting curves
% ---------------
function plotallcurves(P, R, Pboot, Rboot, ERP, freqs, times, mbase, g);

if ~isreal(R)
    Rangle = angle(R);
    R = abs(R); % convert coherence vector to magnitude
    setylim = 1;
else
    Rsign = ones(size(R));
    setylim = 0;
end;

if strcmpi(g.plotitc, 'on') | strcmpi(g.plotersp, 'on')
    verboseprintf(g.verbose, '\nNow plotting...\n');
    pos = get(gca,'position');
    q = [pos(1) pos(2) 0 0];
    s = [pos(3) pos(4) pos(3) pos(4)];
end;

% time unit
% ---------
if times(end) > 10000
    times = times/1000;
    timeunit = 's';
else
    timeunit = 'ms';
end;

if strcmpi(g.plotersp, 'on')
    %
    %%%%%%% image the ERSP %%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    if strcmpi(g.plotitc, 'on'), subplot(2,1,1); end;
    set(gca, 'tag', 'ersp');
    alllegend = {};

    for index = 1:length(freqs)
        alllegend{index} = [ num2str(freqs(index)) 'Hz baseline ' num2str(mbase(index)) ' dB' ];
    end;
    if strcmpi(g.plotmean, 'on') && freqs(1) ~= freqs(end)
        alllegend = { alllegend{:} [ num2str(freqs(1)) '-' num2str(freqs(end)) ...
            'Hz mean baseline ' num2str(mean(mbase)) ' dB' ] };
    end;
    plotcurve(times, P, 'maskarray', Pboot, 'title', 'ERSP', ...
        'xlabel', [ 'Time (' timeunit ')' ], 'ylabel', 'dB', 'ylim', [-g.erspmax g.erspmax], ...
        'vert', g.vert, 'marktimes', g.marktimes, 'legend', alllegend, ...
        'linewidth', g.linewidth, 'highlightmode', g.highlightmode, 'plotmean', g.plotmean);
end;

if strcmpi(g.plotitc, 'on')
    %
    %%%%%%%%%%%% Image the ITC %%%%%%%%%%%%%%%%%%
    %
    if strcmpi(g.plotersp, 'on'), subplot(2,1,2); end;
    set(gca, 'tag', 'itc');
    if abs(R(1,1)-1) < 0.0001, g.plotphaseonly = 'on'; end;
    if strcmpi(g.plotphaseonly, 'on') % plot ITC phase instead of amplitude (e.g. for continuous data)
        RR = Rangle/pi*180;
    else RR = R;
    end;

    % find regions of significance
    % ----------------------------
    alllegend = {};
    for index = 1:length(freqs)
        alllegend{index} = [ num2str(freqs(index)) 'Hz baseline ' num2str(mbase(index)) ' dB' ];
    end;
    if strcmpi(g.plotmean, 'on') && freqs(1) ~= freqs(end)
        alllegend = { alllegend{:} [ num2str(freqs(1)) '-' num2str(freqs(end)) ...
            'Hz mean baseline ' num2str(mean(mbase)) ' dB' ] };
    end;
    plotcurve(times, RR, 'maskarray', Rboot, 'val2mask', R, 'title', 'ITC', ...
        'xlabel', [ 'Time (' timeunit ')' ], 'ylabel', 'dB', 'ylim', g.itcmax, ...
        'vert', g.vert, 'marktimes', g.marktimes, 'legend', alllegend, ...
        'linewidth', g.linewidth, 'highlightmode', g.highlightmode, 'plotmean', g.plotmean);
end;

if strcmpi(g.plotitc, 'on') | strcmpi(g.plotersp, 'on')
    %
    %%%%%%%%%%%%%%% plot a topoplot() %%%%%%%%%%%%%%%%%%%%%%%
    %
    if (~isempty(g.topovec))
        h(12) = axes('Position',[-.1 .43 .2 .14].*s+q);
        if length(g.topovec) == 1
            topoplot(g.topovec,g.elocs,'electrodes','off', ...
                'style', 'blank', 'emarkersize1chan', 10);
        else
            topoplot(g.topovec,g.elocs,'electrodes','off');
        end;
        axis('square')
    end

    try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end;
    if (length(g.title) > 0) && ~iscell(g.title)
        axes('Position',pos,'Visible','Off');
        h(13) = text(-.05,1.01,g.title);
        set(h(13),'VerticalAlignment','bottom')
        set(h(13),'HorizontalAlignment','left')
        set(h(13),'FontSize',g.TITLE_FONT);
    end

    try, axcopy(gcf); catch, end;
end;

%
%%%%%%%%%%%%%%%%%%%%%%% Highlight regions %%%%%%%%%%%%%%%%%%%%%%%%%%
%
function highlight(ax, times, regions, highlightmode);
color1 = [0.75 0.75 0.75];
color2 = [0 0 0];
yl  = ylim;

if ~strcmpi(highlightmode, 'background')
    yl2 = [ yl(1)-(yl(2)-yl(1))*0.15   yl(1)-(yl(2)-yl(1))*0.1 ];
    tmph = patch([times(1) times(end) times(end) times(1)], ...
        [yl2(1) yl2(1) yl2(2) yl2(2)], [1 1 1]); hold on;
    ylim([ yl2(1) yl(2)]);
    set(tmph, 'edgecolor', [1 1 1]);
end;

if ~isempty(regions)
    axes(ax);
    in_a_region = 0;
    for index=1:length(regions)
        if regions(index) && ~in_a_region
            tmpreg(1) = times(index);
            in_a_region = 1;
        end;
        if ~regions(index) && in_a_region
            tmpreg(2) = times(index);
            in_a_region = 0;
            if strcmpi(highlightmode, 'background')
                tmph = patch([tmpreg(1) tmpreg(2) tmpreg(2) tmpreg(1)], ...
                    [yl(1) yl(1) yl(2) yl(2)], color1); hold on;
                set(tmph, 'edgecolor', color1);
            else
                tmph = patch([tmpreg(1) tmpreg(2) tmpreg(2) tmpreg(1)], ...
                    [yl2(1) yl2(1) yl2(2) yl2(2)], color2); hold on;
                set(tmph, 'edgecolor', color2);
            end;
        end;
    end;
end;

% reshaping data
% -----------
function [data, frames] = reshape_data(data, frames)
data = squeeze(data);
if min(size(data)) == 1
    if (rem(length(data),frames) ~= 0)
        error('Length of data vector must be divisible by frames.');
    end
    data = reshape(data, frames, length(data)/frames);
else
    frames = size(data,1);
end

function verboseprintf(verbose, varargin)
if strcmpi(verbose, 'on')
    fprintf(varargin{:});
end;

% reshaping data
% -----------
function pvals = compute_pvals(oridat, surrog, tail)
    
    if nargin < 3
        tail = 'both';
    end;
    
    if myndims(oridat) > 1        
        if size(oridat,2) ~= size(surrog, 2) | myndims(surrog) == 2
            if size(oridat,1) == size(surrog, 1)
                surrog = repmat( reshape(surrog, [size(surrog,1) 1 size(surrog,2)]), [1 size(oridat,2) 1]);
            elseif size(oridat,2) == size(surrog, 1)
                surrog = repmat( reshape(surrog, [1 size(surrog,1) size(surrog,2)]), [size(oridat,1) 1 1]);
            else
                error('Permutation statistics array size error');
            end;
        end;
    end;

    surrog = sort(surrog, myndims(surrog)); % sort last dimension
    
    if myndims(surrog) == 1    
        surrog(end+1) = oridat;        
    elseif myndims(surrog) == 2
        surrog(:,end+1) = oridat;        
    elseif myndims(surrog) == 3
        surrog(:,:,end+1) = oridat;
    else
        surrog(:,:,:,end+1) = oridat;
    end;

    [tmp idx] = sort( surrog, myndims(surrog) );
    [tmp mx]  = max( idx,[], myndims(surrog));        
                
    len = size(surrog,  myndims(surrog) );
    pvals = 1-(mx-0.5)/len;
    if strcmpi(tail, 'both')
        pvals = min(pvals, 1-pvals);
        pvals = 2*pvals;
    end;    
    
function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end;
    end; 


