text:
function [data] = ft_freqsimulation(cfg)

% FT_FREQSIMULATION simulates channel-level time-series data . The data is built up
% from different frequencies and can contain a signal in which the different
% frequencies interact (i.e. cross-frequency coherent). Different methods are
% possible to make data with specific properties.
%
% Use as
%   [data] = ft_freqsimulation(cfg)
% which will return a raw data structure that resembles the output of
% FT_PREPROCESSING.
%
% The configuration options can include
%   cfg.method     = The methods are explained in more detail below, but they can be
%                     'superimposed'    simply add the contribution of the different frequencies
%                     'broadband'       create a single broadband signal component
%                     'phalow_amphigh'  phase of low freq correlated with amplitude of high freq
%                     'amplow_amphigh'  amplitude of low freq correlated with amplithude of high freq
%                     'phalow_freqhigh' phase of low freq correlated with frequency of high signal
%                     'asymmetric'      single signal component with asymmetric positive/negative deflections
%   cfg.output     = which channels should be in the output data, can be 'mixed' or 'all' (default = 'all')
%   cfg.randomseed = 'yes' or a number or vector with the seed value (default = 'yes')
%
% The number of trials and the time axes of the trials can be specified by
%   cfg.fsample    = simulated sample frequency (default = 1200)
%   cfg.trllen     = length of simulated trials in seconds (default = 1)
%   cfg.numtrl     = number of simulated trials (default = 1)
%   cfg.baseline   = number (default = 0)
% or by
%   cfg.time       = cell-array with one time axis per trial, which are for example obtained from an existing dataset
%
% For each of the methods default parameters are configured to generate
% example data, including noise. To get full control over the generated
% data you should explicitely set all parameters involved in the method
% of your choise. The interpretation of the following signal components
% depends on the specified method:
%
% cfg.s1.freq     = frequency of signal 1
% cfg.s1.phase    = phase (in rad) relative to cosine of signal 1  (default depends on method)
%                 = number or 'random'
% cfg.s1.ampl     = amplitude of signal 1
% cfg.s2.freq     = frequency of signal 2
% cfg.s2.phase    = phase (in rad) relative to cosine of signal 1  (default depends on method)
%                 = number or 'random'
% cfg.s2.ampl     = amplitude of signal 2
% cfg.s3.freq     = frequency of signal 3
% cfg.s3.phase    = phase (in rad) relative to cosine of signal 1  (default depends on method)
%                 = number or 'random'
% cfg.s3.ampl     = amplitude of signal 3
% cfg.s4.freq     = frequency of signal 4
% cfg.s4.phase    = phase (in rad) relative to cosine of signal 1  (default depends on method)
%                 = number or 'random'
% cfg.s4.ampl     = amplitude of signal 4
%
% cfg.n1.ampl     = root-mean-square amplitude of wide-band signal prior to filtering
% cfg.n1.bpfreq   = [Flow Fhigh]
% cfg.n2.ampl     = root-mean-square amplitude of wide-band signal prior to filtering
% cfg.n2.bpfreq   = [Flow Fhigh]
%
% cfg.asymmetry   = amount of asymmetry (default = 0, which is none)
% cfg.noise.ampl  = amplitude of noise
%
%
% In the method 'superimposed' the signal contains just the sum of the different frequency contributions:
%     s1: first frequency
%     s2: second frequency
%     s3: third frequency
% and the output consists of the following channels:
%     1st channel: mixed signal = s1 + s2 + s3 + noise
%     2nd channel: s1
%     3rd channel: s2
%     4th channel: s3
%     5th channel: noise
%
% In the method 'broadband' the signal contains a the superposition of two
% broadband signal components, which are created by bandpass filtering a
% Gaussian noise signal:
%     n1: first broadband signal
%     n2: second broadband signal
% and the output consists of the following channels:
%     1st channel: mixed signal = n1 + n2 + noise
%     2nd channel: n1
%     3rd channel: n2
%     4th channel: noise
%
% In the method 'phalow_amphigh' the signal is build up of 4 components; s1, s2, s3 and noise:
%     s1: amplitude modulation (AM), frequency of this signal should be lower than s2
%     s2: second frequency, frequncy that becomes amplitude modulated
%     s3: DC shift of s1, should have frequency of 0
% and the output consists of the following channels:
%     1st channel: mixed signal = (s1 + s3)*s2 + noise,
%     2nd channel: s1
%     3rd channel: s2
%     4th channel: s3
%     5th channel: noise
%
% In the method 'amplow_amphigh' the signal is build up of 5 components; s1, s2, s3, s4 and noise.
%     s1: first frequency
%     s2: second frequency
%     s3: DC shift of s1 and s2, should have frequency of 0
%     s4: amplitude modulation (AM), frequency of this signal should be lower than s1 and s2
% and the output consists of the following channels:
%     1st channel: mixed signal = (s4 + s3)*s1 + (s4 + s3)*s2 + noise,
%     2nd channel: s1
%     3rd channel: s2
%     4th channel: s3
%     5th channel: noise
%     6th channel: s4
%     7th channel: mixed part 1: (s4 + s3)*s1
%     8th channel: mixed part 2: (s4 + s3)*s2
%
% In the method 'phalow_freqhigh' a frequency modulated signal is created.
%   signal is build up of 3 components; s1, s2 and noise.
%     s1: represents the base signal that will be modulated
%     s2: signal that will be used for the frequency modulation
% and the output consists of the following channels:
%     1st channel: mixed signal = s1.ampl * cos(ins_pha) + noise
%     2nd channel: s1
%     3rd channel: s2
%     4th channel: noise
%     5th channel: inst_pha_base   instantaneous phase of the high (=base) frequency signal s1
%     6th channel: inst_pha_mod    low frequency phase modulation, this is equal to s2
%     7th channel: inst_pha        instantaneous phase, i.e. inst_pha_base + inst_pha_mod
%
% In the method 'asymmetric' there is only one periodic signal, but that
% signal is more peaked for the positive than for the negative deflections.
% The average of the signal over time is zero.
%     s1: represents the frequency of the base signal
% and the output consists of the following channels:
%     1st channel: mixed signal = asymmetric signal + noise
%     2nd channel: sine wave with base frequency and phase, i.e. s1
%     3rd channel: asymmetric signal
%     4th channel: noise
%
% See also FT_FREQANALYSIS, FT_TIMELOCKSIMULATION, FT_DIPOLESIMULATION,
% FT_CONNECTIVITYSIMULATION

% Copyright (C) 2007-2008, Ingrid Nieuwenhuis & Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble provenance
ft_preamble randomseed

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% return immediately after distributed execution
if ~isempty(ft_getopt(cfg, 'distribute'))
  return
end

% set defaults
if ~isfield(cfg, 'method'),        cfg.method = 'phalow_amphigh';         end
if ~isfield(cfg, 'output'),        cfg.output = 'all';                    end
if ~isfield(cfg, 'time'),          cfg.time = [];                         end

if isempty(cfg.time)
  cfg.fsample   = ft_getopt(cfg, 'fsample', 1200);
  cfg.trllen    = ft_getopt(cfg, 'trllen', 1);
  cfg.numtrl    = ft_getopt(cfg, 'numtrl', 1);
  cfg.baseline  = ft_getopt(cfg, 'baseline', 0);
else
  cfg.trllen  = [];                         % can be variable
  cfg.fsample = 1/mean(diff(cfg.time{1}));  % determine from time-axis
  cfg.numtrl  = length(cfg.time);
end

if strcmp(cfg.method, 'superimposed')
  if ~isfield(cfg, 's1'),           cfg.s1 = [];                          end
  if ~isfield(cfg.s1, 'freq'),      cfg.s1.freq  = 10;                    end
  if ~isfield(cfg.s1, 'phase'),     cfg.s1.phase = 0;                     end
  if ~isfield(cfg.s1, 'ampl'),      cfg.s1.ampl  = 1;                     end
  if ~isfield(cfg, 's2'),           cfg.s2 = [];                          end
  if ~isfield(cfg.s2, 'freq'),      cfg.s2.freq  = 20;                    end
  if ~isfield(cfg.s2, 'phase'),     cfg.s2.phase = 0;                     end
  if ~isfield(cfg.s2, 'ampl'),      cfg.s2.ampl  = 0;                     end
  if ~isfield(cfg, 's3'),           cfg.s3 = [];                          end
  if ~isfield(cfg.s3, 'freq'),      cfg.s3.freq  = 30;                    end
  if ~isfield(cfg.s3, 'phase'),     cfg.s3.phase = 0;                     end
  if ~isfield(cfg.s3, 'ampl'),      cfg.s3.ampl  = 0;                     end
end

if strcmp(cfg.method, 'broadband')
  if ~isfield(cfg, 'n1'),           cfg.n1 = [];                          end
  if ~isfield(cfg.n1, 'ampl'),      cfg.n1.ampl  = 1;                     end
  if ~isfield(cfg.n1, 'bpfreq'),    cfg.n1.bpfreq  = [30 50];             end
  if ~isfield(cfg, 'n2'),           cfg.n2 = [];                          end
  if ~isfield(cfg.n2, 'ampl'),      cfg.n2.ampl  = 1;                     end
  if ~isfield(cfg.n2, 'bpfreq'),    cfg.n2.bpfreq  = [80 120];            end
end

if strcmp(cfg.method, 'phalow_amphigh')
  if ~isfield(cfg, 's1'),           cfg.s1 = [];                          end
  if ~isfield(cfg.s1, 'freq'),      cfg.s1.freq = 3;                      end
  if ~isfield(cfg.s1, 'phase'),     cfg.s1.phase = -1*pi;                 end
  if ~isfield(cfg.s1, 'ampl'),      cfg.s1.ampl = 1;                      end
  if ~isfield(cfg, 's2'),           cfg.s2 = [];                          end
  if ~isfield(cfg.s2, 'freq'),      cfg.s2.freq = 20;                     end
  if ~isfield(cfg.s2, 'phase'),     cfg.s2.phase = 0;                     end
  if ~isfield(cfg.s2, 'ampl'),      cfg.s2.ampl = 1;                      end
  if ~isfield(cfg, 's3'),           cfg.s3 = [];                          end
  if ~isfield(cfg.s3, 'freq'),      cfg.s3.freq = 0;                      end
  if ~isfield(cfg.s3, 'phase'),     cfg.s3.phase = 0;                     end
  if ~isfield(cfg.s3, 'ampl'),      cfg.s3.ampl = cfg.s1.ampl;            end
end

if strcmp(cfg.method, 'amplow_amphigh')
  if ~isfield(cfg, 's1'),           cfg.s1 = [];                          end
  if ~isfield(cfg.s1, 'freq'),      cfg.s1.freq = 6;                      end
  if ~isfield(cfg.s1, 'phase'),     cfg.s1.phase = 0;                     end
  if ~isfield(cfg.s1, 'ampl'),      cfg.s1.ampl = 1;                      end
  if ~isfield(cfg, 's2'),           cfg.s2 = [];                          end
  if ~isfield(cfg.s2, 'freq'),      cfg.s2.freq = 20;                     end
  if ~isfield(cfg.s2, 'phase'),     cfg.s2.phase = 0;                     end
  if ~isfield(cfg.s2, 'ampl'),      cfg.s2.ampl = 1;                      end
  if ~isfield(cfg, 's4'),           cfg.s4 = [];                          end
  if ~isfield(cfg.s4, 'freq'),      cfg.s4.freq = 1;                      end
  if ~isfield(cfg.s4, 'phase'),     cfg.s4.phase = -1*pi;                 end
  if ~isfield(cfg.s4, 'ampl'),      cfg.s4.ampl = 1;                      end
  if ~isfield(cfg, 's3'),           cfg.s3 = [];                          end
  if ~isfield(cfg.s3, 'freq'),      cfg.s3.freq = 0;                      end
  if ~isfield(cfg.s3, 'phase'),     cfg.s3.phase = 0;                     end
  if ~isfield(cfg.s3, 'ampl'),      cfg.s3.ampl = cfg.s4.ampl;            end
end

if strcmp(cfg.method, 'phalow_freqhigh')
  if ~isfield(cfg, 's1'),           cfg.s1 = [];                          end
  if ~isfield(cfg.s1, 'freq'),      cfg.s1.freq = 20;                     end
  if ~isfield(cfg.s1, 'phase'),     cfg.s1.phase = 0;                     end
  if ~isfield(cfg.s1, 'ampl'),      cfg.s1.ampl = 1;                      end
  if ~isfield(cfg, 's2'),           cfg.s2 = [];                          end
  if ~isfield(cfg.s2, 'freq'),      cfg.s2.freq = 2;                      end
  if ~isfield(cfg.s2, 'phase'),     cfg.s2.phase = -0.5 * pi;             end %then base freq at t=0
  if ~isfield(cfg.s2, 'ampl'),      cfg.s2.ampl = pi;                     end
end

if strcmp(cfg.method, 'asymmetric')
  if ~isfield(cfg, 's1'),           cfg.s1 = [];                          end
  if ~isfield(cfg.s1, 'freq'),      cfg.s1.freq = 6;                      end
  if ~isfield(cfg.s1, 'phase'),     cfg.s1.phase = 0;                     end
  if ~isfield(cfg.s1, 'ampl'),      cfg.s1.ampl = 1;                      end
  if ~isfield(cfg, 'noise'),        cfg.noise = [];                       end
  if ~isfield(cfg.noise, 'ampl'),   cfg.noise.ampl = 0.1;                 end % default should not be too high
end

if ~isfield(cfg, 'noise'),         cfg.noise = [];                        end
if ~isfield(cfg.noise, 'ampl'),    cfg.noise.ampl = 1;                    end

if ~isempty(cfg.time)
  % use the user-supplied time vectors
  timevec = cfg.time;
else
  % give the user some feedback
  ft_debug('using %f as samping frequency', cfg.fsample);
  ft_debug('using %d trials of %f seconds long', cfg.numtrl, cfg.trllen);
  nsample = round(cfg.trllen*cfg.fsample);
  timevec = cell(1, cfg.numtrl);
  for iTr = 1:cfg.numtrl
    timevec{iTr} = (((1:nsample)-1)/cfg.fsample) - cfg.baseline;
  end
end

% give the user some feedback
ft_info('simulating data using %s method', cfg.method);

%%%%%%% SUPERIMPOSED, SIMPLY ADD THE SIGNALS %%%%%%%%%
if strcmp(cfg.method, 'superimposed')

  % make data
  for iTr = 1:length(timevec)
    if ischar(cfg.s1.phase); phase_s1 = rand * 2 * pi; else phase_s1 = cfg.s1.phase; end
    if ischar(cfg.s2.phase); phase_s2 = rand * 2 * pi; else phase_s2 = cfg.s2.phase; end
    if ischar(cfg.s3.phase); phase_s3 = rand * 2 * pi; else phase_s3 = cfg.s3.phase; end

    s1    = cfg.s1.ampl*cos(2*pi*cfg.s1.freq*timevec{iTr} + phase_s1);
    s2    = cfg.s2.ampl*cos(2*pi*cfg.s2.freq*timevec{iTr} + phase_s2);
    s3    = cfg.s3.ampl*cos(2*pi*cfg.s3.freq*timevec{iTr} + phase_s3);
    noise = cfg.noise.ampl*randn(size(timevec{iTr}));
    mix   = s1 + s2 + s3 + noise;

    data.trial{iTr}(1,:) = mix;
    if strcmp(cfg.output, 'all')
      data.trial{iTr}(2,:) = s1;
      data.trial{iTr}(3,:) = s2;
      data.trial{iTr}(4,:) = s3;
      data.trial{iTr}(5,:) = noise;
    end
    data.time{iTr} = timevec{iTr};
  end % for iTr

  data.label{1} = 'mix';
  if strcmp(cfg.output, 'all')
    data.label{2} = 's1';
    data.label{3} = 's2';
    data.label{4} = 's3';
    data.label{5} = 'noise';
  end
  data.fsample = cfg.fsample;

  %%%%%%% SUPERIMPOSED BROADBAND SIGNAL %%%%%%%%%
elseif strcmp(cfg.method, 'broadband')

  % make data
  for iTr = 1:length(timevec)
    n1    = ft_preproc_bandpassfilter(cfg.n1.ampl*randn(size(timevec{iTr})), cfg.fsample, cfg.n1.bpfreq);
    n2    = ft_preproc_bandpassfilter(cfg.n2.ampl*randn(size(timevec{iTr})), cfg.fsample, cfg.n2.bpfreq);
    noise = cfg.noise.ampl*randn(size(timevec{iTr}));
    mix   = n1 + n2 + noise;

    data.trial{iTr}(1,:) = mix;
    if strcmp(cfg.output, 'all')
      data.trial{iTr}(2,:) = n1;
      data.trial{iTr}(3,:) = n2;
      data.trial{iTr}(4,:) = noise;
    end
    data.time{iTr} = timevec{iTr};
  end % for iTr

  data.label{1} = 'mix';
  if strcmp(cfg.output, 'all')
    data.label{2} = 'n1';
    data.label{3} = 'n2';
    data.label{4} = 'noise';
  end
  data.fsample = cfg.fsample;

  %%%%%%% PHASE TO AMPLITUDE CORRELATION %%%%%%%%%
elseif strcmp(cfg.method, 'phalow_amphigh')

  % sanity checks
  if cfg.s2.freq < cfg.s1.freq
    ft_error('with method is phalow_amphigh freq s2 should be higher than freq s1')
  end
  if cfg.s2.freq > cfg.fsample/2
    ft_error('you cannot have a frequency higher than the sample frequency/2')
  end
  if cfg.s3.freq ~= 0 || cfg.s3.phase ~= 0
    ft_warning('for method phalow_amphigh s3 is DC and therefore expect freq and phase to be zero but they are not')
  end
  if cfg.s3.ampl < cfg.s1.ampl
    ft_warning('expect amplitude s3 (=DC) not to be smaller than amplitude s1 (=low frequency)')
  end

  % make data
  for iTr = 1:length(timevec)

    if ischar(cfg.s1.phase); phase_AM   = rand * 2 * pi; else phase_AM   = cfg.s1.phase;  end
    if ischar(cfg.s2.phase); phase_high = rand * 2 * pi; else phase_high = cfg.s2.phase; end
    if ischar(cfg.s3.phase); phase_DC   = rand * 2 * pi; else phase_DC   = cfg.s3.phase;   end
    high  = cfg.s2.ampl*cos(2*pi*cfg.s2.freq*timevec{iTr} + phase_high);
    AM    = cfg.s1.ampl*cos(2*pi*cfg.s1.freq*timevec{iTr} + phase_AM);
    DC    = cfg.s3.ampl*cos(2*pi*0*timevec{iTr} + phase_DC);
    noise = cfg.noise.ampl*randn(size(timevec{iTr}));
    mix   = ((AM + DC) .* high) + noise;

    data.trial{iTr}(1,:) = mix;
    if strcmp(cfg.output, 'all')
      data.trial{iTr}(2,:) = AM;
      data.trial{iTr}(3,:) = high;
      data.trial{iTr}(4,:) = DC;
      data.trial{iTr}(5,:) = noise;
    end
    data.time{iTr} = timevec{iTr};
  end % for iTr

  data.label{1} = 'mix';
  if strcmp(cfg.output, 'all')
    data.label{2} = 's1 (AM)';
    data.label{3} = 's2 (high)';
    data.label{4} = 's3 (DC)';
    data.label{5} = 'noise';
  end
  data.fsample = cfg.fsample;

  %%%%%%% POWER TO POWER CORRELATION %%%%%%%%%
elseif strcmp(cfg.method, 'amplow_amphigh')

  % sanity checks
  if cfg.s2.freq < cfg.s1.freq || cfg.s1.freq < cfg.s4.freq
    ft_error('with method is powlow_powhigh freq s4 < s1 < s2')
  end
  if cfg.s2.freq > cfg.fsample/2
    ft_error('you cannot have a frequency higher than the sample frequency/2')
  end
  if cfg.s3.freq ~= 0 || cfg.s3.phase ~= 0
    ft_warning('for method powlow_powhigh s3 is DC and therefore expect freq and phase to be zero but they are not')
  end
  if cfg.s3.ampl < cfg.s4.ampl
    ft_warning('expect amplitude s3 (=DC) not to be smaller than amplitude s4 (= AM frequency)')
  end

  % make data
  for iTr = 1:length(timevec)

    if ischar(cfg.s1.phase); phase_low  = rand * 2 * pi; else phase_low = cfg.s1.phase;    end
    if ischar(cfg.s2.phase); phase_high = rand * 2 * pi; else phase_high = cfg.s2.phase;   end
    if ischar(cfg.s3.phase); phase_DC   = rand * 2 * pi; else phase_DC = cfg.s3.phase;     end
    if ischar(cfg.s4.phase); phase_AM   = rand * 2 * pi; else phase_AM = cfg.s4.phase; end
    high     = cfg.s2.ampl*cos(2*pi*cfg.s2.freq*timevec{iTr} + phase_high);
    low      = cfg.s1.ampl*cos(2*pi*cfg.s1.freq*timevec{iTr} + phase_low);
    AM       = cfg.s4.ampl*cos(2*pi*cfg.s4.freq*timevec{iTr} + phase_AM);
    DC       = cfg.s3.ampl*cos(2*pi*0*timevec{iTr} + phase_DC);
    noise    = cfg.noise.ampl*randn(size(timevec{iTr}));
    lowmix  = ((AM + DC) .* low);
    highmix = ((AM + DC) .* high);
    mix     = lowmix + highmix + noise;

    data.trial{iTr}(1,:) = mix;
    if strcmp(cfg.output, 'all')
      data.trial{iTr}(2,:) = low;
      data.trial{iTr}(3,:) = high;
      data.trial{iTr}(4,:) = DC;
      data.trial{iTr}(5,:) = noise;
      data.trial{iTr}(6,:) = AM;
      data.trial{iTr}(7,:) = lowmix;
      data.trial{iTr}(8,:) = highmix;
    end
    data.time{iTr} = timevec{iTr};
  end % for iTr

  data.label{1} = 'mix';
  if strcmp(cfg.output, 'all')
    data.label{2} = 's1 (low)';
    data.label{3} = 's2 (high)';
    data.label{4} = 's3 (DC)';
    data.label{5} = 'noise';
    data.label{6} = 's4 (AM)';
    data.label{7} = 'mixlow';
    data.label{8} = 'mixhigh';
  end
  data.fsample = cfg.fsample;

  %%%%%%% PHASE TO FREQUENCY CORRELATION %%%%%%%%%
elseif strcmp(cfg.method, 'phalow_freqhigh')

  % sanity checks
  if cfg.s1.freq > cfg.fsample/2 || cfg.s2.freq > cfg.fsample/2
    ft_error('you cannot have a frequency higher than the sample frequency/2')
  end

  % make data
  for iTr = 1:length(timevec)

    if ischar(cfg.s1.phase); phase_s1 = rand * 2 * pi; else phase_s1 = cfg.s1.phase;    end
    if ischar(cfg.s2.phase); phase_s2 = rand * 2 * pi; else phase_s2= cfg.s2.phase;    end
    s1            = cfg.s1.ampl .* cos(2*pi*cfg.s1.freq * timevec{iTr} + phase_s1); % to be modulated signal
    s2            = cfg.s2.ampl .* cos(2*pi*cfg.s2.freq * timevec{iTr} + phase_s2); % modulation of instantaneous phase
    inst_pha_base = 2*pi*cfg.s1.freq * timevec{iTr} + phase_s1; % unmodulated instantaneous phase s1 (linear)
    inst_pha_mod  = s2;                                    % modulation of instantaneous phase
    inst_pha      = inst_pha_base + inst_pha_mod;
    noise         = cfg.noise.ampl*randn(size(timevec{iTr}));
    mix           = cfg.s1.ampl .* cos(inst_pha) + noise;

    data.trial{iTr}(1,:) = mix;
    if strcmp(cfg.output, 'all')
      data.trial{iTr}(2,:) = s1;
      data.trial{iTr}(3,:) = s2;
      data.trial{iTr}(4,:) = noise;
      data.trial{iTr}(5,:) = inst_pha_base;
      data.trial{iTr}(6,:) = inst_pha_mod;
      data.trial{iTr}(7,:) = inst_pha;
    end
    data.time{iTr} = timevec{iTr};
  end % for iTr

  data.label{1} = 'mix';
  if strcmp(cfg.output, 'all')
    data.label{2} = 's1';
    data.label{3} = 's2';
    data.label{4} = 'noise';
    data.label{5} = 'inst phase base';
    data.label{6} = 'inst phase modulation (=s2)';
    data.label{7} = 'inst phase';
  end
  data.fsample = cfg.fsample;

  %%%%%%% ASYMETRIC POSITIVE AND NEGATIVE PEAKS %%%%%%%%%
elseif strcmp(cfg.method, 'asymmetric')

  % make data
  for iTr = 1:length(timevec)
    if ischar(cfg.s1.phase); phase_s1 = rand * 2 *pi; else phase_s1 = cfg.s1.phase; end

    s1    = cfg.s1.ampl*cos(2*pi*cfg.s1.freq*timevec{iTr} + phase_s1);
    tmp   = cos(2*pi*cfg.s1.freq*timevec{iTr} + phase_s1);  % same signal but with unit amplitude
    tmp   = (tmp+1)/2;                                 % scaled and shifted between 0 and 1
    tmp   = tmp.^(cfg.asymmetry+1);                    % made asymmetric
    tmp   = (tmp - mean(tmp))*2*cfg.s1.ampl;           % rescale
    s2    = tmp;
    noise = cfg.noise.ampl*randn(size(timevec{iTr}));
    mix   = s2 + noise;

    data.trial{iTr}(1,:) = mix;
    if strcmp(cfg.output, 'all')
      data.trial{iTr}(2,:) = s1;
      data.trial{iTr}(3,:) = s2;
      data.trial{iTr}(4,:) = noise;
    end
    data.time{iTr} = timevec{iTr};
  end % for iTr

  data.label{1} = 'mix';
  if strcmp(cfg.output, 'all')
    data.label{2} = 's1';
    data.label{3} = 's2';
    data.label{4} = 'noise';
  end
  data.fsample = cfg.fsample;

else
  ft_error('unknown method specified')
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble randomseed
ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function r8vec_print ( n, a, title )

%*****************************************************************************80
%
%% R8VEC_PRINT prints an R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, real A(N), the vector to be printed.
%
%    Input, string TITLE, a title to be printed first.
%    TITLE may be blank.
%
  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '%6d  %12f\n', i, a(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% IVP with respect to Alfa
function dy = JacobianAlfa (t, y)
    global SOL;
    global ETA_INF;
    global BETA0;
    global BETA;
    global EX;
    method = 'spline';
    % y(1) = df\dalfa
    % y(2) = du\dalfa
    % y(3) = dv\dalfa
    f = interp1(EX, SOL(:,1), t, method); % Interpolate the data set for f
    u = interp1(EX, SOL(:,2), t, method); % Interpolate the data set for u
    v = interp1(EX, SOL(:,3), t, method); % Interpolate the data set for v
    dy = [ETA_INF*y(2) ETA_INF*y(3), -ETA_INF*(BETA0*(y(1)*v+f*y(3))-2*BETA*u*y(2))]';

author: Sable

repo: mcbench-benchmarks


text:
function lvPerfTable = getLVPerfTable(lvDef, liftoffUt)
    lvPerfTable = [];
    lvDefCumData = getStageData(lvDef);
    startUt = liftoffUt;
    for(i=1:size(lvDefCumData,1)) %#ok<*NO4LP>
        row = lvDefCumData(i,:);
        coastDur = row(8);
        coastStart = startUt + row(1);
        endUt = startUt + row(1) + coastDur;
        thrust = row(2);
        mdot = row(3);
        massStart = row(4);
        massEnd = row(5);
        dryMassCum = row(6);
        fuMassCum = row(7);
        cDa = row(9);

        lvPerfTable(i,:) = [startUt, endUt, thrust, mdot, massStart, massEnd, dryMassCum, fuMassCum, coastStart, cDa]; %#ok<AGROW>
        startUt = endUt;
    end

    lvPerfTable(end+1,:) = lvPerfTable(end,:);
    lvPerfTable(end, 1) = lvPerfTable(end-1, 2);
    lvPerfTable(end, 2) = Inf;
    lvPerfTable(end, 3) = 0.0;
    lvPerfTable(end, 4) = 0.0;
    lvPerfTable(end, 5) = lvPerfTable(end-1, 6);
    lvPerfTable(end, 6) = lvPerfTable(end, 5);
    lvPerfTable(end, 7) = lvPerfTable(end-1, 7);
    lvPerfTable(end, 8) = 0.0;
    lvPerfTable(end, 9) = Inf;
    lvPerfTable(end, 10) = 0.0;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% op_alignAverages.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [out,fs,phs]=op_alignAverages(in,tmax,med,ref);
%
% DESCRIPTION:
% Perform spectral registration in the time domain to correct frequency and
% phase drifts.  As described in Near et al.  Frequency and phase drift
% correction of magnetic resonance spectroscopy data by spectral
% registration in the time domain. Magn Reson Med 2015; 73(1):44-50.
%
% June 15th 2017:  Made the tmax and med arguments optional.  If tmax is
% not specified, the value is determined automatically by finding the time
% at which the SNR of the FID drops permanently below 5.  This idea
% was suggested by Mark Mikkelsen.  Thanks Mark!!
%
% INPUTS:
% in        = Input data structure.
% tmax      = Maximum time (s) in time domain to use for alignment.
%             (Optional.  Default is the time at which SNR drops below 5)
% med       = Align averages to the median of the averages? ('y','n', 'a', or
%             'r').  (Optional.  Default = 'n').  If you select 'n', all
%             averages will be aligned to a single average.  The average
%             chosen as the reference average will be the one with the
%             lowest 'unlikeness' metric (see 'op_rmbadaverages.m').  If
%             select 'y', all averages will be alinged to the median of
%             the averages.  If you select 'a', all averages will be
%             aligned to the average of the averages.  If you select 'r',
%             all averages will be aligned to an externally provided
%             reference spectrum.
% ref       = An externally provided reference spectrum that you would like
%             to align everything to (Required only if med = 'r').
%
% OUTPUTS:
% out       = Output following alignment of averages.
% fs        = Vector of frequency shifts (in Hz) used for alignment.
% phs       = Vector of phase shifts (in degrees) used for alignment.

function [out,fs,phs]=op_alignAverages(in,tmax,med,initPars)

if ~in.flags.addedrcvrs
    error('ERROR:  I think it only makes sense to do this after you have combined the channels using op_addrcvrs.  ABORTING!!');
end

if in.dims.averages==0
    %DO NOTHING
    disp('WARNING:  No averages found.  Returning input without modification!');
    out=in;
    fs=0;
    phs=0;

else

    parsFit=[0,0];

    if nargin<3
        med='n'
        if nargin<2
            %if tmax is not specified, find the time at which the SNR
            %drops below 5
            disp('tmax not supplied.  Calculating tmax....');
            sig=abs(in.fids);
            noise=std(real(in.fids(ceil(0.75*end):end,:,:)),[]);
            noise=mean(mean(mean(noise,2),3),4);
            snr=sig/noise;

            for n=1:(numel(snr)/size(snr,1))
                N=find(snr(:,n)>5);
                tmax_est(n)=in.t(N(end));
            end
            tmax=median(tmax_est);
            disp(['tmax = ' num2str(tmax*1000) 'ms.']);
        end
    end

    if (strcmp(med,'r') || strcmp(med,'R'))
        if nargin<4
            error('ERROR:  If using the ''r'' option for input variable ''med'', then a 4th input argument must be provided');
        end
    else
        if nargin<4
            ref=struct();
        end
    end

    if in.dims.subSpecs==0
        B=1;
    else
        B=in.sz(in.dims.subSpecs);
    end

    fs=zeros(in.sz(in.dims.averages),B);
    phs=zeros(in.sz(in.dims.averages),B);
    fids=zeros(in.sz(in.dims.t),1,B);
    for m=1:B
        if med=='y' || med=='Y'
            disp('Aligning all averages to the median of the averages.');
            base=op_median(in);
            base=[real(base.fids( in.t>=0 & in.t<tmax ,m));imag(base.fids( in.t>=0 & in.t<tmax ,m))];
            ind_min=0;
        elseif med=='a' || med=='a'
            disp('Aligning all averages to the average of the averages.');
            base=op_averaging(in);
            base=[real(base.fids( in.t>=0 & in.t<tmax ,m));imag(base.fids( in.t>=0 & in.t<tmax ,m))];
            ind_min=0;
        elseif med=='n' || med=='N'
            %First find the average that is most similar to the total average:
            inavg=op_median(in);
            for k=1:in.sz(in.dims.averages)
                for l=1:B
                    metric(k,l)=sum((real(in.fids(in.t>=0 & in.t<=tmax,k,l))-(real(inavg.fids(inavg.t>=0 & inavg.t<=tmax,l)))).^2);
                end
            end
            [temp,ind_min]=min(metric(:,m));

            %Now set the base function using the index of the most similar
            %average:
            disp(['Aligning all averages to average number ' num2str(ind_min) '.']);
            base=[real(in.fids(in.t>=0 & in.t<tmax,ind_min,m));imag(in.fids(in.t>=0 & in.t<tmax,ind_min,m))];
            fids(:,ind_min,m)=in.fids(:,ind_min,m);
        elseif med=='r' || med=='R'
            disp('Aligning all averages to an externally provided reference spectrum.');
            base=ref;
            base=[real(base.fids( in.t>=0 & in.t<tmax ,m));imag(base.fids( in.t>=0 & in.t<tmax ,m))];
            ind_min=0;
        end
        for n=1:in.sz(in.dims.averages)
            if n~=ind_min
                parsGuess=parsFit;
                %disp(['fitting subspec number ' num2str(m) ' and average number ' num2str(n)]);
                parsFit=nlinfit(in.fids(in.t>=0 & in.t<tmax,n,m),base,@op_freqPhaseShiftComplexNest,parsGuess);
                fids(:,n,m)=op_freqPhaseShiftNest(parsFit,in.fids(:,n,m));
                fs(n,m)=parsFit(1);
                phs(n,m)=parsFit(2);
                %plot(in.ppm,fftshift(ifft(fids(:,1,m))),in.ppm,fftshift(ifft(fids(:,n,m))));
            end
        end
    end

    %re-calculate Specs using fft
    specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

    %FILLING IN DATA STRUCTURE
    out=in;
    out.fids=fids;
    out.specs=specs;

    %FILLING IN THE FLAGS
    out.flags=in.flags;
    out.flags.writtentostruct=1;
    out.flags.freqcorrected=1;

end

    function y=op_freqPhaseShiftComplexNest(pars,input)
        f=pars(1);     %Frequency Shift [Hz]
        p=pars(2);     %Phase Shift [deg]

        dwelltime=in.dwelltime;
        t=0:dwelltime:(length(input)-1)*dwelltime;
        fid=input(:);

        shifted=addphase(fid.*exp(1i*t'*f*2*pi),p);

        y=[real(shifted);imag(shifted)];
        %y=real(fid.*exp(-1i*t'*f*2*pi));

    end

    function y=op_freqPhaseShiftNest(pars,input)
        f=pars(1);     %Frequency Shift [Hz]
        p=pars(2);     %Phase Shift [deg]

        dwelltime=in.dwelltime;
        t=0:dwelltime:(length(input)-1)*dwelltime;
        fid=input(:);

        y=addphase(fid.*exp(1i*t'*f*2*pi),p);
        %y=real(fid.*exp(-1i*t'*f*2*pi));

    end

end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = sigmoid(z) .* (1 - sigmoid(z));

% =============================================================

end

repo: machine-learning


text:
function test06 ( sample_num )

%*****************************************************************************80
%
%% TEST06 times R4_UNI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06\n' );
  fprintf ( 1, '  Measure the time it takes R4_UNI to generate\n' );
  fprintf ( 1, '  %d uniform deviates.\n', sample_num );

  seed = uint32 ( 123456789 );

  time1 = cputime;

  for sample = 1 : sample_num
    [ value, seed ] = r4_uni ( seed );
  end

  time2 = cputime;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  %f seconds.\n', time2 - time1 );

  return
end


text:
%%% Extraction of the final intrinsic and extrinsic paramaters:

check_active_images;

if ~exist('solution_error')
   solution_error = zeros(6*n_ima + 15,1);
end;

fc = solution(1:2);%***
cc = solution(3:4);%***
alpha_c = solution(5);%***
kc = solution(6:10);%***

fc_error = solution_error(1:2);
cc_error = solution_error(3:4);
alpha_c_error = solution_error(5);
kc_error = solution_error(6:10);

% Calibration matrix:

KK = [fc(1) fc(1)*alpha_c cc(1);0 fc(2) cc(2); 0 0 1];
inv_KK = inv(KK);

% Extract the extrinsic paramters, and recomputer the collineations

for kk = 1:n_ima,

   if active_images(kk),

      omckk = solution(15+6*(kk-1) + 1:15+6*(kk-1) + 3);%***
      Tckk = solution(15+6*(kk-1) + 4:15+6*(kk-1) + 6);%***

      omckk_error = solution_error(15+6*(kk-1) + 1:15+6*(kk-1) + 3);
      Tckk_error = solution_error(15+6*(kk-1) + 4:15+6*(kk-1) + 6);

   	Rckk = rodrigues(omckk);

   	Hkk = KK * [Rckk(:,1) Rckk(:,2) Tckk];

   	Hkk = Hkk / Hkk(3,3);

   else

      omckk = NaN*ones(3,1);
      Tckk = NaN*ones(3,1);
      Rckk = NaN*ones(3,3);
      Hkk = NaN*ones(3,3);
      omckk_error = NaN*ones(3,1);
      Tckk_error = NaN*ones(3,1);

   end;

   eval(['omc_' num2str(kk) ' = omckk;']);
   eval(['Rc_' num2str(kk) ' = Rckk;']);
   eval(['Tc_' num2str(kk) ' = Tckk;']);
   eval(['H_' num2str(kk) '= Hkk;']);
   eval(['omc_error_' num2str(kk) ' = omckk_error;']);
   eval(['Tc_error_' num2str(kk) ' = Tckk_error;']);

end;

author: strawlab

repo: MultiCamSelfCal


text:
function anim = generateSFromLSBasis( anim )
% Generate 3D data from the basis of an Animation
%
% USAGE
%  anim = anim.generateSFromLSBasis()
%
% INPUTS
%  anim     - Animation object (help Animation for details)
%
% OUTPUTS
%  anim     - modified Animation anime
%
% EXAMPLE
%
% See also GENERATETOYANIMATION
%
% Vincent's Structure From Motion Toolbox      Version 3.0
% Copyright (C) 2009 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

if size(anim.l,1)>=(size(anim.SBasis,3)-1) && ~isempty(anim.SBasis)
  anim.nPoint=size(anim.SBasis,2);
  switch size(anim.l,1)
    case anim.nBasis,
      anim.S=reshape(sum(bsxfun(@times,reshape(anim.l,1,1,anim.nBasis,...
        anim.nFrame),anim.SBasis),3),3,anim.nPoint,anim.nFrame);
    case anim.nBasis-1,
      % add the first basis with a coeff of 1
      anim.S=bsxfun(@plus,anim.SBasis(:,:,1), reshape(sum(bsxfun(@times,...
        reshape(anim.l,1,1,anim.nBasis-1,anim.nFrame),...
        anim.SBasis(:,:,2:end)),3),3,anim.nPoint,anim.nFrame));
    otherwise,
      error('l and nBasis do not have matching dimensions');
  end
end

author: vrabaud


text:
function b = PVsubsasgn_2idx(a,L,RHS)
% function B =  PVsubsasgn_2idx(A,L,RHS)
%
% DESCRIPTION
%   (INTERNAL FUNCTION)
%   (:)-subsassign for polynomial objects.
%
% INPUTS
%   A: polynomial
%   L: a 1x1 structure with the fields:
%    type -- string containing '()'
%    subs -- 1x2 cell array containing indices
%   RHS: Value to be assigned.
%
% OUTPUTS
%   B: object after subsassignment
%
% SYNTAX
%   B =  subsasgn(A,L,RHS)

% 6/9/2002: PJS  Initial Coding
% 11/11/2008: PJS Mods to reduce computation
% 1/13/2011: PJS  Fixed bug in sub2ind conversion

% Check Indices
ridx = L(1).subs{1};
lr = length(ridx);
cidx = L(1).subs{2};
lc = length(cidx);
if min(ridx)<1 || min(cidx)<1
    error('Index into matrix is negative or zero.');
end

% Get info about polynomials
b = a;
[nrb,ncb]=size(b);
ntb = size(b.degmat,1);

temp = RHS;
[nrt,nct]=size(temp);
ntt = size(temp.degmat,1);

% Fan out scalar RHS to correct dimension
if all([nrt nct]==[1 1])
    temp.coefficient = repmat(temp.coefficient,1,lr*lc);
    nrt = lr;
    nct = lc;
    temp.matdim = [nrt nct];
end

% Check dimensions
if lr*lc~=nrt*nct
    error(['In an assignment  A(I) = B, the number of' ...
        ' elements in B and I must be the same.']);
end

% Determine proper dimensions for poly after the assignment
nr = max([nrb; ridx(:)]);
nc = max([ncb; cidx(:)]);

% Convert from row/col indices to single index
if isempty(ridx) || isempty(cidx)
    % No substitution
    b = a;
    return;
end
% XXX PJS If ridx and cidx are N-by-1 then the commented code produces
% subsidx of size N-by-1. However, it should be N^2-by-1 (all combos)
% elseif length(cidx)==1
%     cidx = repmat(cidx,length(ridx),1);
% elseif length(ridx)==1
%     ridx = repmat(ridx,length(cidx),1);
% end
% subsidx = sub2ind([nr nc],ridx(:),cidx(:));
tmp = reshape(1:nr*nc,[nr,nc]);
subsidx = tmp(ridx,cidx);
subsidx = subsidx(:);

% Perform the subsassignment
if isempty(b)
    tempcoef = sparse(ntt,nr*nc);
    tempcoef(:,subsidx) = temp.coefficient;
    temp.coefficient = tempcoef;
    temp.matdim = [nr nc];

    b = temp;
else
    % Place coefs of b in proper location, zero out terms of b to be
    % assigned, stack in new terms, and then combine
    coef = sparse(ntb,nr*nc);
    [ii,jj]=ind2sub([nrb ncb],1:(nrb*ncb));
    bidx = sub2ind([nr nc],ii,jj);
    coef(:,bidx) = b.coefficient;

    coef(:,subsidx) = 0;
    coef(end+1:end+ntt,subsidx) = temp.coefficient;

    degmat = blkdiag(b.degmat,temp.degmat);
    varname = [b.varname; temp.varname];
    matdim = [nr nc];
    chkval = 0; % skip validity check
    b = polynomial(coef,degmat,varname,matdim,chkval);
    b = combine(b);
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% openbdf() - Opens an BDF File (European Data Format for Biosignals) in MATLAB (R)
%
% Usage:
%   >> EDF=openedf(FILENAME)
%
% Note: About EDF -> www.biosemi.com/faq/file_format.htm
%
% Author: Alois Schloegl, 5.Nov.1998
%
% See also: readedf()

%   Copyright (C) 1997-1998 by Alois Schloegl
%   a.schloegl@ieee.org
%   Ver 2.20    18.Aug.1998
%   Ver 2.21    10.Oct.1998
%   Ver 2.30     5.Nov.1998
%
%   For use under Octave define the following function
% function s=upper(s); s=toupper(s); end

% V2.12    Warning for missing Header information
% V2.20    EDF.AS.* changed
% V2.30    EDF.T0 made Y2K compatible until Year 2090

% 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.

% Name changed for bdf files Sept 6,2002  T.S. Lorig
% Header updated for EEGLAB format (update web link too) - Arnaud Delorme 14 Oct 2002

function [DAT,H1]=openbdf(FILENAME)

SLASH='/';      % defines Seperator for Subdirectories
BSLASH=char(92);

cname=computer;
if cname(1:2)=='PC' SLASH=BSLASH; end

fid=fopen_or_error(FILENAME,'r','ieee-le');

EDF.FILE.FID=fid;
EDF.FILE.OPEN = 1;
EDF.FileName = FILENAME;

PPos=min([max(find(FILENAME=='.')) length(FILENAME)+1]);
SPos=max([0 find((FILENAME=='/') | (FILENAME==BSLASH))]);
EDF.FILE.Ext = FILENAME(PPos+1:length(FILENAME));
EDF.FILE.Name = FILENAME(SPos+1:PPos-1);
if SPos==0
    EDF.FILE.Path = pwd;
else
    EDF.FILE.Path = FILENAME(1:SPos-1);
end
EDF.FileName = [EDF.FILE.Path SLASH EDF.FILE.Name '.' EDF.FILE.Ext];

H1=char(fread(EDF.FILE.FID,256,'char')');     %
EDF.VERSION=H1(1:8);                     % 8 Byte  Versionsnummer
%if 0 fprintf(2,'LOADEDF: WARNING  Version EDF Format %i',ver); end
EDF.PID = deblank(H1(9:88));                  % 80 Byte local patient identification
EDF.RID = deblank(H1(89:168));                % 80 Byte local recording identification
%EDF.H.StartDate = H1(169:176);         % 8 Byte
%EDF.H.StartTime = H1(177:184);         % 8 Byte
EDF.T0=[str2num(H1(168+[7 8])) str2num(H1(168+[4 5])) str2num(H1(168+[1 2])) str2num(H1(168+[9 10])) str2num(H1(168+[12 13])) str2num(H1(168+[15 16])) ];

% Y2K compatibility until year 2090
if EDF.VERSION(1)=='0'
        if EDF.T0(1) < 91
                EDF.T0(1)=2000+EDF.T0(1);
        else
                EDF.T0(1)=1900+EDF.T0(1);
        end
else
        % in a future version, this is hopefully not needed
end

EDF.HeadLen = str2num(H1(185:192));  % 8 Byte  Length of Header
% reserved = H1(193:236);            % 44 Byte
EDF.NRec = str2num(H1(237:244));     % 8 Byte  # of data records
EDF.Dur = str2num(H1(245:252));      % 8 Byte  # duration of data record in sec
EDF.NS = str2num(H1(253:256));       % 8 Byte  # of signals

EDF.Label = char(fread(EDF.FILE.FID,[16,EDF.NS],'char')');
EDF.Transducer = char(fread(EDF.FILE.FID,[80,EDF.NS],'char')');
EDF.PhysDim = char(fread(EDF.FILE.FID,[8,EDF.NS],'char')');

EDF.PhysMin= str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));
EDF.PhysMax= str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')'));
EDF.DigMin = str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')')); %
EDF.DigMax = str2num(char(fread(EDF.FILE.FID,[8,EDF.NS],'char')')); %

% check validity of DigMin and DigMax
if (length(EDF.DigMin) ~= EDF.NS)
        fprintf(2,'Warning OPENEDF: Failing Digital Minimum\n');
        EDF.DigMin = -(2^15)*ones(EDF.NS,1);
end
if (length(EDF.DigMax) ~= EDF.NS)
        fprintf(2,'Warning OPENEDF: Failing Digital Maximum\n');
        EDF.DigMax = (2^15-1)*ones(EDF.NS,1);
end
if (any(EDF.DigMin >= EDF.DigMax))
        fprintf(2,'Warning OPENEDF: Digital Minimum larger than Maximum\n');
end
% check validity of PhysMin and PhysMax
if (length(EDF.PhysMin) ~= EDF.NS)
        fprintf(2,'Warning OPENEDF: Failing Physical Minimum\n');
        EDF.PhysMin = EDF.DigMin;
end
if (length(EDF.PhysMax) ~= EDF.NS)
        fprintf(2,'Warning OPENEDF: Failing Physical Maximum\n');
        EDF.PhysMax = EDF.DigMax;
end
if (any(EDF.PhysMin >= EDF.PhysMax))
        fprintf(2,'Warning OPENEDF: Physical Minimum larger than Maximum\n');
        EDF.PhysMin = EDF.DigMin;
        EDF.PhysMax = EDF.DigMax;
end
EDF.PreFilt= char(fread(EDF.FILE.FID,[80,EDF.NS],'char')'); %
tmp = fread(EDF.FILE.FID,[8,EDF.NS],'char')';   %   samples per data record
EDF.SPR = str2num(char(tmp));   %   samples per data record

fseek(EDF.FILE.FID,32*EDF.NS,0);

EDF.Cal = (EDF.PhysMax-EDF.PhysMin)./ ...
    (EDF.DigMax-EDF.DigMin);
EDF.Off = EDF.PhysMin - EDF.Cal .* EDF.DigMin;
tmp = find(EDF.Cal < 0);
EDF.Cal(tmp) = ones(size(tmp));
EDF.Off(tmp) = zeros(size(tmp));

EDF.Calib=[EDF.Off';(diag(EDF.Cal))];
%EDF.Calib=sparse(diag([1; EDF.Cal]));
%EDF.Calib(1,2:EDF.NS+1)=EDF.Off';

EDF.SampleRate = EDF.SPR / EDF.Dur;

EDF.FILE.POS = ftell(EDF.FILE.FID);
if EDF.NRec == -1   % unknown record size, determine correct NRec
  fseek(EDF.FILE.FID, 0, 'eof');
  endpos = ftell(EDF.FILE.FID);
  EDF.NRec = floor((endpos - EDF.FILE.POS) / (sum(EDF.SPR) * 2));
  fseek(EDF.FILE.FID, EDF.FILE.POS, 'bof');
  H1(237:244)=sprintf('%-8i',EDF.NRec);      % write number of records
end

EDF.Chan_Select=(EDF.SPR==max(EDF.SPR));
for k=1:EDF.NS
    if EDF.Chan_Select(k)
        EDF.ChanTyp(k)='N';
    else
        EDF.ChanTyp(k)=' ';
    end
    if contains(upper(EDF.Label(k,:)),'ECG')
        EDF.ChanTyp(k)='C';
    elseif contains(upper(EDF.Label(k,:)),'EKG')
        EDF.ChanTyp(k)='C';
    elseif contains(upper(EDF.Label(k,:)),'EEG')
        EDF.ChanTyp(k)='E';
    elseif contains(upper(EDF.Label(k,:)),'EOG')
        EDF.ChanTyp(k)='O';
    elseif contains(upper(EDF.Label(k,:)),'EMG')
        EDF.ChanTyp(k)='M';
    end
end

EDF.AS.spb = sum(EDF.SPR);  % Samples per Block
bi=[0;cumsum(EDF.SPR)];

idx=[];idx2=[];
for k=1:EDF.NS,
    idx2=[idx2, (k-1)*max(EDF.SPR)+(1:EDF.SPR(k))];
end
maxspr=max(EDF.SPR);
idx3=zeros(EDF.NS*maxspr,1);
for k=1:EDF.NS, idx3(maxspr*(k-1)+(1:maxspr))=bi(k)+ceil((1:maxspr)'/maxspr*EDF.SPR(k));end;

%EDF.AS.bi=bi;
EDF.AS.IDX2=idx2;
%EDF.AS.IDX3=idx3;

DAT.Head=EDF;
DAT.MX.ReRef=1;

%DAT.MX=feval('loadxcm',EDF);

return


text:
function [C,theta] = x2majoraxis(A,B)
%X2MAJORAXIS Aligns coordinate x with the major axis of a region.
%  [C,THETA] = X2MAJORAXIS(A,B) aligns the x-coordinate axis with the
%  major axis of a region or boundary. The y-axis is perpendicular to
%  the x-axis.  The rows of 2-by-2 matrix A are the coordinates of the
%  two end points of the major axis, in the form A = [x1 y1; x2 y2].
%  Input B is either a an image of class logical containing a single
%  region, or it is an np-by-2 set of points representing a (connected)
%  boundary. In the latter case, the first column of B must represent
%  x-coordinates and the second column must represent the corresponding
%  y-coordinates. Output C contains the same data as the input, but
%  aligned with the major axis. If the input is an image, so is the
%  output; similarly the output is a sequence of coordinates if the
%  input is such a sequence. Parameter THETA is the initial angle
%  between the major axis and the x-axis. The origin of the xy-axis
%  system is at the bottom left; the x-axis is the horizontal axis and
%  the y-axis is the vertical.
%
%  Keep in mind that rotations can introduce round-off errors when the
%  data are converted to integer (pixel) coordinates, which typically is
%  a requirement.  Thus, postprocessing (e.g., with bwmorph) of the
%  output may be required to reconnect a boundary.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Preliminaries.
if islogical(B)
  type = 'region';
elseif size(B,2) == 2
  type = 'boundary';
  [M,N] = size(B);
  if M < N
    error('B is boundary. It must be of size np-by-2; np > 2.')
  end
  % Compute centroid for later use. c is a 1-by-2 vector.
  % Its 1st component is the mean of the boundary in the x-direction.
  % The second is the mean in the y-direction.
  c(1) = round((min(B(:,1)) + max(B(:,1)))/2);
  c(2) = round((min(B(:,2)) + max(B(:,2)))/2);

  % It is possible for a connected boundary to develop small breaks
  % after rotation. To prevent this, the input boundary is filled,
  % processed as a region, and then the boundary is re-extracted.
  % This guarantees that the output will be a connected boundary.
  m = max(max(B));
  % The following image is of size m-by-m to make sure that there
  % there will be no size truncation after rotation.
  B = bound2im(B,m,m);
  B = imfill(B,'holes');
else
  error('Input must be a boundary or a binary image.')
end

% Major axis in vector form.
v(1) = A(2,1) - A(1,1);
v(2) = A(2,2) - A(1,2);
v = v(:);  % v is a col vector

% Unit vector along x-axis.
u = [1; 0];

% Find angle between major axis and x-axis. The angle is
% given by acos of the inner product of u and v divided by
% the product of their norms. Because the inputs are image
% points, they are in the first quadrant.
nv = norm(v);
nu = norm(u);
theta = acos(u'*v/nv*nu);
if theta > pi/2
  theta = -(theta - pi/2);
end
theta = theta*180/pi;  % Convert angle to degrees.

% Rotate by angle theta and crop the rotated image to original size.
C = imrotate(B,theta,'bilinear','crop');

% If the input was a boundary, re-extract it.
if  strcmp(type,'boundary')
  C = bwboundaries(C);
  C = C{1};
  % Shift so that centroid of the extracted boundary is
  % approx equal to the centroid of the original boundary:
  C(:,1) = C(:,1) - (min(C(:,1)) + max(C(:,1)))/2 + c(1);
  C(:,2) = C(:,2) - (min(C(:,2)) + max(C(:,2)))/2 + c(2);
end

author: dipum

repo: dipum-toolbox

save_path: github-repos/MATLAB/dipum-dipum-toolbox


text:
%  Figure 10.14      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%  fig10_14.m is a script to generate Fig. 10.14, the transient
%  response of the PD plus notch compensator of the
%  satellite position control, non-colocated case

% satellite system matrices
f =[0    1.0000         0         0;
   -0.9100   -0.0360    0.9100    0.0360;
         0         0         0    1.0000 ;
    0.0910    0.0036   -0.0910   -0.0036];
g =[0;
     0;
     0;
     1];

h =[1     0     0     0];

j =[0];
no3=conv(.25*[2 1],[1/.81 0 1]);
do3=conv([1/20 1],[1/625 2/25 1]);
sys3=tf(no3,do3);
sysG=ss(f,g,h,j);
[sysol]=series(sys3,sysG);
[Aol,Bol,Col,Dol]=ssdata(sysol);
Acl=Aol-Bol*Col;
hold off; clf
t=[0:0.1:40];
syscl=ss(Acl,Bol,Col,Dol)
[y,t]=step(syscl,t);
plot(t,y,'LineWidth',2);
grid on;
xlabel('Time (sec)');
ylabel('Amplitude');
title('Closed-loop step response of KD_3(s)G(s)')

%grid
nicegrid;


text:
% Function to derive audio features for silence and unvoiced regions
%
% Description
% Function to derive audio features for silence and unvoiced regions
%
%
% Inputs
%  x        : [samples] [Nx1]  Speech signal
%  fs       : [Hz]      [1x1]  Sampling frequency
%  frame_len_ms: [ms] [1x1]  Window length
%
% Outputs
%  Ts       : [samples] [Px1]  Locations of parameter values
%  Es       : [dB]      [Mx1]  Energy contour
%
% Example
%  Please see the HOWTO_glottalsource.m example file.
%
% References
%  [1] Drugman, T., Kane, J., Gobl, C., `Automatic Analysis of Creaky
%       Excitation Patterns', Submitted to Computer Speech and
%       Language.
%  [2] Kane, J., Drugman, T., Gobl, C., (2013) `Improved automatic
%       detection of creak', Computer Speech and Language 27(4), pp.
%       1028-1047.
%  [3] Drugman, T., Kane, J., Gobl, C., (2012) `Resonator-based creaky
%       voice detection', Interspeech 2012, Portland, Oregon, USA.
%
% Copyright (c) 2013 University of Mons, FNRS & 2013 Trinity College Dublin
%
% License
%  This code is a part of the GLOAT toolbox with the following
%  licence:
%  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 3 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Authors
%  Thomas Drugman <thomas.drugman@umons.ac.be> & John Kane <kanejo@tcd.ie>

function [Ts,Es,ZCs_ms,Xpos,pos] = sil_unv_features(x,fs,frame_len_ms)

if nargin < 3
    frame_len_ms=10;
end
frame_len=frame_len_ms/1000*fs;
frame_shift_ms=5;
Shift=frame_shift_ms/1000*fs;

Start=1; % intialise
Stop=Start+frame_len-1;
pos=mean([Start Stop]);

[ZCs,Xpos] = get_zero_x_rate(x,frame_len,Shift);
ZCs_ms=ZCs/1000*fs;

Ind=1;

Es=zeros(1,floor((length(x)-Stop+1)/Shift)+1);
Ts=zeros(size(Es));

win = hanning(Stop-Start+1);

while Stop<length(x)

    Mid=round((Start+Stop)/2);
    Ts(Ind)=Mid;

    Sig=x(Start:Stop); % Select frame segment
    Sig=Sig(:).*win; % Window segment

    Es(Ind)=mean(Sig.^2); % Get energy value

    Start=Start+Shift;
    Stop=Stop+Shift;
    pos(Ind)=mean([Start Stop]);
    Ind=Ind+1;
end

Es=log(Es);

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function f = fix(f)
%FIX   Round a CLASSICFUN pointwise toward zero.
%   G = FIX(F) returns the CLASSICFUN G such that G(x) = FIX(F(x)) for each x in
%   F.domain.
%
%   If F is complex, then the G = FIX(REAL(F)) + 1i*FIX(IMAG(F)).
%
%   Note that FIX() assumes the output G(X) is a constant. If it is not, then
%   garbage is returned with no warning.
%
% See also ROUND, CEIL, FLOOR.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

f.onefun = fix(f.onefun);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = psabr_2(a,b,r,n,f,k,t)
% SABR risk neutral density calculated from option prices
% calculated by density extrapolation technique sprice_2(...)
% as proposed by Kainth et al. using no specific cut off point

    eps = 1e-003;
    y1 = sprice_2(a,b,r,n,f,k+eps,t,-1);
    y2 = sprice_2(a,b,r,n,f,k,t,-1);
    y3 = sprice_2(a,b,r,n,f,k-eps,t,-1);
    y = (y1-2*y2+y3)/eps^2;
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MATLAB Code for                                                  %
%                                                                   %
%  Multi-Objective Particle Swarm Optimization (MOPSO)              %
%  Version 1.0 - Feb. 2011                                          %
%                                                                   %
%  According to:                                                    %
%  Carlos A. Coello Coello et al.,                                  %
%  "Handling Multiple Objectives with Particle Swarm Optimization," %
%  IEEE Transactions on Evolutionary Computation, Vol. 8, No. 3,    %
%  pp. 256-279, June 2004.                                          %
%                                                                   %
%  Developed Using MATLAB R2009b (Version 7.9)                      %
%                                                                   %
%  Programmed By: S. Mostapha Kalami Heris                          %
%                                                                   %
%         e-Mail: sm.kalami@gmail.com                               %
%                 kalami@ee.kntu.ac.ir                              %
%                                                                   %
%       Homepage: http://www.kalami.ir                              %
%                                                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Index SubIndex]=GetGridIndex(particle,G)

    c=particle.Cost;

    nobj=numel(c);
    ngrid=numel(G(1).Upper);

    str=['sub2ind(' mat2str(ones(1,nobj)*ngrid)];

    SubIndex=zeros(1,nobj);
    for j=1:nobj

        U=G(j).Upper;

        i=find(c(j)<U,1,'first');

        SubIndex(j)=i;

        str=[str ',' num2str(i)];
    end

    str=[str ');'];

    Index=eval(str);

end

author: sfvsfv

repo: Mathematical-modeling


text:
function degree = r8poly_degree ( na, a )

%*****************************************************************************80
%
%% R8POLY_DEGREE returns the degree of a polynomial.
%
%  Discussion:
%
%    The degree of a polynomial is the index of the highest power
%    of X with a nonzero coefficient.
%
%    The degree of a constant polynomial is 0.  The degree of the
%    zero polynomial is debatable, but this routine returns the
%    degree as 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NA, the dimension of A.
%
%    Input, real A(NA+1), the coefficients of the polynomials.
%
%    Output, integer DEGREE, the degree of A.
%
  degree = na;

  while ( 0 < degree )

    if ( a(degree+1) ~= 0.0 )
      return
    end

    degree = degree - 1;

  end

  return
end


text:
% Test file for trigtech/circconv

function pass = test_circconv(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

% Generate a random point to use as test values.
seedRNG(6178);
x = 2 * rand - 1;

%%
% Check operation in the face of empty arguments.

f = testclass.make();
g = testclass.make(@(x) sin(pi*x), [], pref);
pass(1) = (isempty(circconv(f,g)) && isempty(circconv(g,f)));

%%
% Simple checks

f_op = @(x) tanh(5*cos(pi*(x)));
g_op = @(x) ones(size(x));
f = testclass.make(f_op, [], pref);
g = testclass.make(g_op, [], pref);

hfg = circconv(f,g);
hgf = circconv(g,f);
% Answer should be zero since the functions are odd.
approx_fg = feval(hfg,x);
approx_gf = feval(hgf,x);
expected = 0;
err = abs(approx_fg - expected);
tol = 1e2*eps;
pass(2) = err < tol;

f_op = @(x) tanh(cos(pi*(x)));
a = x;
fa_op = @(x) tanh(cos(pi*(x-a)));

% Harder tests
f = testclass.make(f_op, [], pref);
g = circconv(f,f);

% Computed answer at zero:
approx = feval(g,0);
% Expected answer:
expected = sum(f.*f);
err = abs(approx - expected);
tol = 1e2*eps*vscale(g);
pass(3) = err < tol;

% Computed answer at x:
approx = feval(g,x);
% Expected answer:
h = testclass.make(fa_op, [], pref);
expected = sum(f.*h);
err = abs(approx - expected);
tol = 1e2*eps*vscale(g);
pass(4) = err < tol;

end


text:
function [V, converged, i] = newtonpf_S_cart(Ybus, Sbus, V0, ref, pv, pq, mpopt)
%NEWTONPF_S_CART  Solves power flow using full Newton's method (power/cartesian)
%   [V, CONVERGED, I] = NEWTONPF_S_CART(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
%
%   Solves for bus voltages using a full Newton-Raphson method, using nodal
%   power balance equations and cartesian coordinate representation of
%   voltages, given the following inputs:
%       YBUS  - full system admittance matrix (for all buses)
%       SBUS  - handle to function that returns the complex bus power
%               injection vector (for all buses), given the bus voltage
%               magnitude vector (for all buses)
%       V0    - initial vector of complex bus voltages
%       REF   - bus index of reference bus (voltage ang reference & gen slack)
%       PV    - vector of bus indices for PV buses
%       PQ    - vector of bus indices for PQ buses
%       MPOPT - (optional) MATPOWER option struct, used to set the
%               termination tolerance, maximum number of iterations, and
%               output options (see MPOPTION for details).
%
%   The bus voltage vector contains the set point for generator
%   (including ref bus) buses, and the reference angle of the swing
%   bus, as well as an initial guess for remaining magnitudes and
%   angles.
%
%   Returns the final complex voltages, a flag which indicates whether it
%   converged or not, and the number of iterations performed.
%
%   See also RUNPF, NEWTONPF, NEWTONPF_I_POLAR, NEWTONPF_I_CART.

%   MATPOWER
%   Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Baljinnyam Sereeter, Delft University of Technology
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% default arguments
if nargin < 7
    mpopt = mpoption;
end

%% options
tol         = mpopt.pf.tol;
max_it      = mpopt.pf.nr.max_it;
lin_solver  = mpopt.pf.nr.lin_solver;

%% initialize
converged = 0;
i = 0;
V = V0;
Vm = abs(V);
Vmpv = Vm(pv);

%% set up indexing for updating V
npv = length(pv);
npq = length(pq);
j1 = 1;         j2 = npq;           %% j1:j2 - Vr of pq buses
j3 = j2 + 1;    j4 = j2 + npv;      %% j3:j4 - Vr of pv buses
j5 = j4 + 1;    j6 = j4 + npq;      %% j5:j6 - Vi of pq buses
j7 = j6 + 1;    j8 = j6 + npv;      %% j7:j8 - Vi of pv buses

%% evaluate F(x0)
mis = V .* conj(Ybus * V) - Sbus(Vm);
F = [   real(mis([pq; pv]));
        imag(mis(pq));
        V(pv) .* conj(V(pv)) - Vmpv.^2  ];

%% check tolerance
normF = norm(F, inf);
if mpopt.verbose > 1
    fprintf('\n it    max P & Q mismatch (p.u.)');
    fprintf('\n----  ---------------------------');
    fprintf('\n%3d        %10.3e', i, normF);
end
if normF < tol
    converged = 1;
    if mpopt.verbose > 1
        fprintf('\nConverged!\n');
    end
end

%% attempt to pick fastest linear solver, if not specified
if isempty(lin_solver)
    nx = length(F);
    if nx <= 10 || have_feature('octave')
        lin_solver = '\';       %% default \ operator
    else    %% MATLAB and nx > 10 or Octave and nx > 2000
        lin_solver = 'LU3';     %% LU decomp with 3 output args, AMD ordering
    end
end

%% do Newton iterations
while (~converged && i < max_it)
    %% update iteration counter
    i = i + 1;

    %% evaluate Jacobian
    [dSbus_dVr, dSbus_dVi] = dSbus_dV(Ybus, V, 1);
    dV2_dVr = sparse(1:npv, npq+(1:npv), 2*real(V(pv)), npv, npv+npq);
    dV2_dVi = sparse(1:npv, npq+(1:npv), 2*imag(V(pv)), npv, npv+npq);

    %% handling of derivatives for voltage dependent loads
    %% (not yet implemented) goes here

    j11 = real(dSbus_dVr([pq; pv], [pq; pv]));
    j12 = real(dSbus_dVi([pq; pv], [pq; pv]));
    j21 = imag(dSbus_dVr(pq, [pq; pv]));
    j22 = imag(dSbus_dVi(pq, [pq; pv]));
    j31 = dV2_dVr;
    j32 = dV2_dVi;

    J = [   j11 j12;
            j21 j22;
            j31 j32;    ];

    %% compute update step
    dx = mplinsolve(J, -F, lin_solver);

    %% update voltage
    if npv
        V(pv) = V(pv) + dx(j3:j4) + 1j * dx(j7:j8);
    end
    if npq
        V(pq) = V(pq) + dx(j1:j2) + 1j * dx(j5:j6);
    end

    %% evalute F(x)
    mis = V .* conj(Ybus * V) - Sbus(Vm);
    F = [   real(mis([pq; pv]));
            imag(mis(pq));
            V(pv) .* conj(V(pv)) - Vmpv.^2  ];

    %% check for convergence
    normF = norm(F, inf);
    if mpopt.verbose > 1
        fprintf('\n%3d        %10.3e', i, normF);
    end
    if normF < tol
        converged = 1;
        if mpopt.verbose
            fprintf('\nNewton''s method power flow (power balance, cartesian) converged in %d iterations.\n', i);
        end
    end
end

if mpopt.verbose
    if ~converged
        fprintf('\nNewton''s method power flow (power balance, cartesian) did not converge in %d iterations.\n', i);
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
close all;
clearvars;
clc;
[LoD,HiD,LoR,HiR] = wfilters('db4');

load noisdopp;

% change the DWT mode to periodization.
old_dwt_mode = dwtmode('status','nodisp');
dwtmode('per');
% perform 4 level decomposition
[coefficients, levels] = wavedec(noisdopp,4, LoD,HiD);

% plot approximation and detail components
figure;
plot(coefficients); title('Coefficients');

export_fig images/noisdopp_db4_l4_decomposition_per.png -r120 -nocrop;

% reconstruct the signal
reconstructed = waverec(coefficients, levels, LoR, HiR);

% measure the maximum difference
max_abs_diff = max(abs(noisdopp-reconstructed))

figure;
for level=0:4
    level_app_coeffs = appcoef(coefficients, levels, LoR, HiR, level);
    subplot(511+level);
    plot(level_app_coeffs);
    fprintf('%d ', numel(level_app_coeffs));
    title(sprintf('Approximation coefficients @ level-%d', level));
end

% restore the old DWT mode
dwtmode(old_dwt_mode);

export_fig images/noisdopp_db4_l4_appcoeffs_per.png -r120 -nocrop;

author: indigits

repo: sparse-plex


text:
% solution of Poisson equation using sparse matrix
%
%  E.Holzbecher,    18.3.2011
%
%--------------------------------------------------
nx = 12; ny = 12;      % dimensions in x- and y-direction
h = 1/4;                % grid spacing

% boundary type indicators (1=Dirichlet, 0=Neumann no-flow)
ltop = logical(zeros(1,nx));         % top
lbottom = logical(zeros(1,nx));      % bottom
lleft = logical([ones(6,1) zeros(6,1)]);       % left
lright = logical([zeros(6,1) ones(6,1)]);      % right

% boundary values (Dirichlet only)
btop = ones(1,nx);                  % top
bbottom = zeros(1,nx);              % bottom
bleft = ones(ny,1);                 % left
bright = zeros(ny,1);               % right

q = 1*ones(nx,ny);                  % right hand side (source term)

N = nx*ny;
d = [-nx,-1,0,1,nx];
B = [ones(N,2) -4*ones(N,1) ones(N,2)];
q = reshape(q,N,1);
b = -h*h*q.*ones(N,1);
for i = 1:nx
    if ltop(i)
        b(i) = b(i)-btop(i);
    else
        B(i,3) = -3;
        %B(i-1,1) = 0;
    end
    if lbottom(i)
        b(N-nx+i) = b(N-nx+i)-bbottom(i);
    else
        B(N-nx+i,3) = -3;
        % B(N-nx+i,5) = 0;
    end
end
for i = 1:ny
    B(i*nx,2) = 0;
    if i<ny B(i*nx+1,4) = 0; end
    if lleft(i)
        b((i-1)*nx+1) = b((i-1)*nx+1)-bleft(i);
    else
        B((i-1)*nx+1,3) = B((i-1)*nx+1,3)+1;
    end
    if lright(i)
        b(i*nx) = b(i*nx)-bright(i);
    else
        B(i*nx,3) = B(i*nx,3)+1;
    end
end

A = spdiags(B,d,N,N);

% processing: solution
U = A\b;
U = reshape(U,nx,ny);

% check & visualize
4*del2(U)
surf (U)


text:
function [Hx,Hy,Ez,time] = MaxwellPNonCon2D(pinfo, Hx, Hy, Ez, FinalTime)

% function [Hx,Hy,Ez] = MaxwellPNonCon2D(pinfo, Hx, Hy, Ez, FinalTime)
% Purpose  : Integrate TM-mode Maxwell's until FinalTime starting with initial conditions Hx,Hy,Ez

Globals2D;
time = 0;

% Runge-Kutta residual storage
resHx = zeros(size(Hx)); resHy = resHx; resEz = resHx;

% compute time step size (taking into account variable polynomial order)
dt = 100;
for N=1:length(pinfo)
  dt = min(dt, 2./( N^2*max(pinfo(N).Fscale(:)/2)));
end

% outer time step loop
tstep = 1;
while (time<FinalTime)

   if(time+dt>FinalTime), dt = FinalTime-time; end

   for INTRK = 1:5
      % compute right hand side of TM-mode Maxwell's equations
      [rhsHx, rhsHy, rhsEz] = MaxwellPNonConRHS2D(pinfo, Hx,Hy,Ez);

      % initiate and increment Runge-Kutta residuals
      resHx = rk4a(INTRK)*resHx + dt*rhsHx;
      resHy = rk4a(INTRK)*resHy + dt*rhsHy;
      resEz = rk4a(INTRK)*resEz + dt*rhsEz;

      % update fields
      Hx = Hx+rk4b(INTRK)*resHx;
      Hy = Hy+rk4b(INTRK)*resHy;
      Ez = Ez+rk4b(INTRK)*resEz;
   end;

   % Increment time
   time = time+dt; tstep = tstep+1;
end
return

repo: Aero-matlab


text:
function test_failed=test_pbspline

Lr=[15,16,18,20];
ar=[ 3, 4, 6, 5];
or=[1, 1.5, 2,3];

%btypes={'ed','xd','stard','ec','xc','starc'};
btypes={'ed','xd','stard'};
centtypes={'wp','hp'};

test_failed=0;

disp(' ===============  TEST_PBSPLINE ============');

for ii=1:length(Lr)
  L=Lr(ii);
  a=ar(ii);
  N=L/a;

  for jj=1:length(or)
    order=or(jj);

    for kk=1:numel(btypes)
      btype=btypes{kk};

      for ll=1:2
        centstring=centtypes{ll};

        [g,nlen]=pbspline(L,order,a,btype,centstring);

        A=zeros(L,1);

        for n=0:N-1
          A=A+circshift(g,n*a);
        end;

        res=max(abs(A-1/sqrt(a)));
        [test_failed,fail]=ltfatdiditfail(res,test_failed);
        s=sprintf('PBSPLINE PU   %2s %s L:%3i a:%3i o:%3.5g %0.5g %s', ...
                  btype,centstring,L,a,order,res,fail);
        disp(s);

        gcutextend=middlepad(middlepad(g,nlen,centstring),L,centstring);

        res=norm(g-gcutextend);

        [test_failed,fail]=ltfatdiditfail(res,test_failed);
        s=sprintf('PBSPLINE NLEN %2s %s L:%3i a:%3i o:%3.5g %0.5g %s', ...
                  btype,centstring,L,a,order,res,fail);
        disp(s);

      end;

    end;
  end;

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [HR, HRS, HlR, HlRS, HiR, HiRS, ChiR, HshR, HshRS] = gaussian_method_v7_1_0(R, pars)
%GAUSSIAN_METHOD_V7_1_0 Estimate entropy values using the Gaussian Method.
%
% ---------
% ARGUMENTS
% ---------
% R    - response matrix
% pars - parameters structure

%   Copyright (C) 2009 Cesare Magri
%   Version 7.1.0

% -------
% LICENSE
% -------
% This software is distributed free under the condition that:
%
% 1. it shall not be incorporated in software that is subsequently sold;
%
% 2. the authorship of the software shall be acknowledged and the following
%    article shall be properly cited in any publication that uses results
%    generated by the software:
%
%      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
%      toolbox for the fast information analysis of multiple-site LFP, EEG
%      and spike train recordings. BMC Neuroscience 2009 10(1):81;
%
% 3.  this notice shall remain in place in each source file.

persistent previous_Nt previous_mask;

% HiR, HiRS and ChiR are not computed for the gaussian case, null values
% are returned instead:
HiR  = NaN;
HiRS = NaN;
ChiR = NaN;

% We have to recompute Nc, Nt, maxNt (in the case this routine is called by
% quadratic extrapolation)
Nt    = pars.Nt;
Nc    = size(R,1);
maxNt = size(R,2);
totNt = sum(Nt);
Ns    = pars.Ns;

logOf2times2 = 2*log(2);
logOf2piTimesExp1 = log(2*pi*exp(1));

if ~isequal(previous_Nt, Nt)
  mask = build_logical_mask(Nt, maxNt, Ns);
else
	mask = previous_mask;
end
R(:,~mask) = nan;

% H(R) and H_lin(R) -------------------------------------------------------
if pars.doHR || pars.doHlR

    % Removing mean:
    %meanRvec = mean(R(:,mask),2);
		meanRvec = (sum(sum(R,3),2)+maxNt*Ns-totNt)./totNt; %JM: this is faster, some juggling is needed to account for the -1's

	  %Rvec = R(:,mask) - meanRvec(:,ones(totNt,1));
    Rvec = R(:,mask) - meanRvec*ones(1,totNt);

    covPr = (Rvec * Rvec') ./ (totNt-1);
end

% H(R)
if pars.doHR
    HR = (Nc * logOf2piTimesExp1 + log(det(covPr))) / logOf2times2;

    % Bias correction:
    if pars.biasCorrNum==3
        HRbias = gaussian_bias(totNt, Nc);
        HR = HR - HRbias;
    end
else
    HR = 0;
end

% H_lin(R)
if pars.doHlR
    HlR = (Nc * logOf2piTimesExp1 + sum(log(diag(covPr)), 1)) / logOf2times2;

    % Bias correction:
    if pars.biasCorrNum==3
        HlRbias = Nc * gaussian_bias(totNt, 1);
        HlR = HlR - HlRbias;
    end
else
    HlR = 0;
end

% H(R|S) and H_lin(R|S) ---------------------------------------------------
if pars.doHRS || pars.doHlRS

    % Removing mean
    meanR = sum(R.*permute(mask(:,:,ones(Nc,1)), [3 1 2]), 2) ./ permute(Nt(:, ones(Nc,1)), [2 3 1]);
    Rzeromean = R - meanR(:, ones(maxNt,1), :);

    if Nc==1
        % HlRS is never computed for Nc==1, if we got here it's cause we
        % need to compute H(R|S):

				%covPrsMat = covPrsForSingletonNc(Rzeromean, Nt, mask, maxNt, Ns);
        covPrsMat = reshape(nansum(Rzeromean.*Rzeromean,2),1,[])./(Nt'-1); %JM:

        HRS  = (Nc * logOf2piTimesExp1 + log(covPrsMat(:))) / logOf2times2;

        HlRS = 0;

        % Bias correction:
        if pars.biasCorrNum==3
            HRSbias = gaussian_bias(Nt, 1);
            HRS = HRS - HRSbias;
        end
    else
        HRS  = zeros(Ns,1);
        HlRS = zeros(Ns,1);

        for s=1:Ns
            covPrs = Rzeromean(:, 1:Nt(s), s) * Rzeromean(:, 1:Nt(s), s).' ./ (Nt(s)-1);

            if pars.doHRS
                detCovPrs = det(covPrs);

                HRS(s) = (Nc * logOf2piTimesExp1 + log(detCovPrs)) / logOf2times2;

                % Bias correction:
                if pars.biasCorrNum==3
                    HRSbias = gaussian_bias(Nt(s), Nc);
                    HRS(s) = HRS(s) - HRSbias;
                end
            end;

            if pars.doHlRS
                diagCovPrs = diag(covPrs);

                HlRS(s) = (length(diagCovPrs) * logOf2piTimesExp1 + sum(log(diagCovPrs))) / logOf2times2;

                % Bias correction:
                if pars.biasCorrNum==3
                    HlRSbias = Nc * gaussian_bias(Nt(s), 1);
                    HlRS(s) = HlRS(s) - HlRSbias;
                end
            end
        end
    end

    if ~pars.doHRS
        HRS = 0;
    end

    if ~pars.doHlRS
        HlRS = 0;
    end

else
    HRS  = 0;
    HlRS = 0;
end

% H_sh(R) and H_sh(R|S) ---------------------------------------------------
if pars.doHshR || pars.doHshRS

    Rsh = shuffle_R_across_cells(R, Nt);

    if pars.doHshR
        % Removing mean:
        meanRshvec = mean(Rsh(:,mask),2);
        Rshvec = Rsh(:,mask) - meanRshvec(:,ones(totNt,1));

        covPshr = (Rshvec * Rshvec') ./ (totNt-1);

        HshR = (Nc * logOf2piTimesExp1 + log(det(covPshr))) / logOf2times2;

        % Bias correction:
        if pars.biasCorrNum==3
            HshRbias = gaussian_bias(totNt, Nc);
            HshR = HshR - HshRbias;
        end
    else
        HshR = 0;
    end

    % Remark: HshRS not computed for Nc==1
    if pars.doHshRS
        % Removing mean
        meanRsh = sum(Rsh.*permute(mask(:,:,ones(Nc,1)), [3 1 2]), 2) ./ permute(Nt(:, ones(Nc,1)), [2 3 1]);
        Rsh = Rsh - meanRsh(:,ones(maxNt,1),:);

        HshRS = zeros(Ns,1);
        for s=1:Ns
            covPshrs = Rsh(:, 1:Nt(s), s) * Rsh(:, 1:Nt(s), s).' ./ (Nt(s)-1);
            detCovPshrs = det(covPshrs);

            HshRS(s) = (Nc * logOf2piTimesExp1 + log(detCovPshrs)) / logOf2times2;
        end

        % Bias correction:
        if pars.biasCorrNum==3
            HshRSbias = gaussian_bias(Nt(:), Nc);
            HshRS = HshRS - HshRSbias;
        end
    else
        HshRS = 0;
    end

else
    HshR  = 0;
    HshRS = 0;
end

previous_Nt   = Nt;
previous_mask = mask;

function covPrs = covPrsForSingletonNc(R, Nt, mask, maxNt, Ns)

% Removing first singleton dimension:
R2D = zeros(maxNt, Ns);
R2D(mask) = R(1,mask);

% Computing the variance (R.*R faster than R.^2):
covPrs = (sum(R2D.*R2D,1)) ./ (Nt.' - 1);


text:
% IO utility functions to
% - save data to files.
% - read, edit mat file timestamps
%
% Binary Format for matrices and cell arrays (with matrix and cell array
% elements) inspired by Arvid Bottiger's write_matrix_bin,
% 1. 1*uint32               Indicator of cell array
% 2. 1*uint32               Dimensions of array
% 3. dimensions*uint32      Lengths of dimensions
% 4. a. Cell Array => recurse
%    b. Matrix
%       1*uint32            Indicator of character array
%       uint32 or float32   Dat%
%
%
%References:
%===============
%- write_matrix_bin by Arvid Bottiger
%   http://www.mathworks.com/matlabcentral/fileexchange/24483
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 7/14/2009
classdef IOUtil
    properties (Constant = true)
        memmapfileHeaderSize = 1024
        memmapfileDataField = 'data'
    end

    methods (Static)
        %Writes data to disk in binary format
        function writeBinary(filename, data)
            import edu.stanford.covert.util.IOUtil;

            %open file
            fid = fopen(filename, 'w');

            %write data
            IOUtil.writeDataBinary(fid, data);

            %close file
            fclose(fid);
        end

        %Reads data from disk in binary format
        function data = readBinary(filename)
            import edu.stanford.covert.util.IOUtil;

            %open file
            fid = fopen(filename, 'r');

            %read data
            data = IOUtil.readDataBinary(fid);

            %close file
            fclose(fid);
        end
    end

    methods (Static, Access = protected)
        %Writes data to disk in binary format
        function writeDataBinary(fid, data)
            import edu.stanford.covert.util.IOUtil;

            %cell array, number of dimensions, size
            dataIsCell = iscell(data);
            dataSize = size(data);
            dimensions = length(dataSize);
            fwrite(fid, [dataIsCell dimensions dataSize], 'uint32');

            %reshape data
            data = reshape(data,1,[]);

            %write data
            if dataIsCell
                for i = 1:length(data)
                    IOUtil.writeDataBinary(fid, data{i});
                end
            else
                switch class(data)
                    case 'char',    dataType = 1;
                    case 'uint8',   dataType = 2;
                    case 'uint16',  dataType = 3;
                    case 'uint32',  dataType = 4;
                    case 'uint64',  dataType = 5;
                    case 'int8',    dataType = 6;
                    case 'int16',   dataType = 7;
                    case 'int32',   dataType = 8;
                    case 'int64',   dataType = 9;
                    case 'single',  dataType = 10;
                    case 'double',  dataType = 11;
                    case 'logical', dataType = 12;
                    otherwise
                        throw(MException('IOUtil:unsupportedDatatype', 'Data type ''%s'' not supported', class(data)));
                end

                %character, write no. dimensions, size
                fwrite(fid,  dataType, 'uint32');

                %write data
                switch dataType
                    case 1,  fwrite(fid, data+0, 'uint16');
                    case 2,  fwrite(fid, data, 'uint16');  %Note not using uint8
                    case 3,  fwrite(fid, data, 'uint16');
                    case 4,  fwrite(fid, data, 'uint32');
                    case 5,  fwrite(fid, data, 'float64'); %Note not using uint64
                    case 6,  fwrite(fid, data, 'int16');   %Note not using int8
                    case 7,  fwrite(fid, data, 'int16');
                    case 8,  fwrite(fid, data, 'int32');
                    case 9,  fwrite(fid, data, 'float64'); %Note not using int64
                    case 10, fwrite(fid, data, 'float32');
                    case 11, fwrite(fid, data, 'float64');
                    case 12, fwrite(fid, data, 'uint16');  %Note not using ubit1
                    otherwise
                        throw(MException('IOUtil:unsupportedDatatype', 'Data type ''%s'' not supported', dataType));
                end
            end
        end

        %Reads data from disk in binary format
        function data = readDataBinary(fid)
            import edu.stanford.covert.util.IOUtil;

            %cell array, number of dimensions, size
            dataIsCell = fread(fid, 1, 'uint32');
            dimensions = fread(fid, 1, 'uint32');
            dataSize = fread(fid, dimensions, 'uint32')';
            count = prod(dataSize);

            %read data
            if dataIsCell
                data = cell(count, 1);
                for i = 1:count
                    data{i} = IOUtil.readDataBinary(fid);
                end
            else
                dataType = fread(fid, 1, 'uint32');
                switch dataType
                    case 1,  data = char(fread(fid, count, 'uint16')); %#ok<FREAD>
                    case 2,  data = uint8(fread(fid, count, 'uint16'));
                    case 3,  data = uint16(fread(fid, count, 'uint16'));
                    case 4,  data = uint32(fread(fid, count, 'uint32'));
                    case 5,  data = uint64(fread(fid, count, 'uint64'));
                    case 6,  data = int8(fread(fid, count, 'int16'));
                    case 7,  data = int16(fread(fid, count, 'int16'));
                    case 8,  data = int32(fread(fid, count, 'int32'));
                    case 9,  data = int64(fread(fid, count, 'int64'));
                    case 10, data = single(fread(fid, count, 'float32'));
                    case 11, data = double(fread(fid, count, 'float64'));
                    case 12, data = logical(fread(fid, count, 'uint16'));
                    otherwise
                        throw(MException('IOUtil:unsupportedDatatype', 'Data type ''%d'' not supported', dataType));
                end
            end

            if dimensions > 1
                data = reshape(data, dataSize);
            end
        end
    end

    methods (Static)
        function writeMemmapFile(data, fileName)
            import edu.stanford.covert.util.IOUtil;

            %write header
            format = {class(data) size(data) IOUtil.memmapfileDataField};
            IOUtil.writeBinary(fileName, format);

            %write spacing between header and data, and data
            fileInfo = dir(fileName);
            if fileInfo.bytes > IOUtil.memmapfileHeaderSize
                throw(MException('IOUtil:invalidHeader', 'Header size cannot be greater than %d bytes', IOUtil.memmapfileHeaderSize));
            end

            fid = fopen(fileName, 'a');
            fwrite(fid, zeros(IOUtil.memmapfileHeaderSize - fileInfo.bytes, 1, 'uint8'), 'uint8'); %spacing
            fwrite(fid, data, class(data)); %data
            fclose(fid);
        end

        function result = readMemmapFile(fileName)
            import edu.stanford.covert.util.IOUtil;

            format = IOUtil.readBinary(fileName);
            if all(format{2})
                result = memmapfile(fileName, ...
                    'Offset', IOUtil.memmapfileHeaderSize, ...
                    'Format', format, ...
                    'Repeat', 1, ...
                    'Writable', false);
            else
                result = struct('Data', struct(format{3}, []));
                result.Data.(format{3}) = zeros(format{2}, format{1});
            end
        end
    end

    methods (Static)
        function directories = getDirectoryNamesRecursively(directory)
            import edu.stanford.covert.util.IOUtil;
            if directory(end) == '/' || directory(end) == '\'
                directory = directory(1:end-1);
            end
            directories = {directory};
            files = dir(directory);
			files = files([files.isdir]);
            for i = 1:numel(files)
                if files(i).name(1) ~= '.'
                    directories = [directories;
                        IOUtil.getDirectoryNamesRecursively([directory filesep files(i).name])]; %#ok<AGROW>
                end
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function targNames = prtUtilGenerateDefaultTargetNames(indices2)
% targNames = prtUtilGenerateDefaultTargetNames(indices2)

% Previously static in prtDataSetBase
%     function targNames = generateDefaultTargetNames(indices2)
%         targNames = prtUtilCellPrintf('Target %d',num2cell(indices2));
%         targNames = targNames(:);
%     end
targNames = prtUtilCellPrintf('Target %d',num2cell(indices2));
targNames = targNames(:);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function val = get(w,prop_name)
   %GET Get waveform properties
   %   val = get(waveform,prop_name)
   %
   %   Valid property names:
   %       CHANNELTAG, NETWORK, STATION, LOCATION, CHANNEL,
   %       CHANNELINFO : (gets the info as a string 'Net.Sta.Loc.Cha')
   %       FREQ, START_STR,
   %       END_STR, DATA, NYQ, PERIOD
   %     also:
   %       START, END : return datenum format (MatLab's format)
   %       START_EPOCH, END_EPOCH : return epoch format.
   %
   %       DATA_LENGTH : return number of elements in data
   %
   %       The following return time duration of data
   %       DURATION_STR :  text format
   %       DURATION: matlab (datenum) format
   %       DURATION_EPOCH : antelope (epoch) format (just # of seconds)
   %
   %       TIMEVECTOR : return a vector of same length as data with
   %                    matlab-formatted times.
   %
   %       MISC_FIELDS: gets a list of fields that were added to this waveform
   %       UNITS : Find out what units the data is in (ex. counts, nm/s)
   %
   %   If waveform is N-dimensional, then VAL will be a cell of the same
   %   dimensions.  If GET would return single-values for the property, then
   %   VAL will be a matrix of type DOUBLE, arranged in the same dimensions.
   %
   %       If additional fields were added to waveform using ADDFIELD, then
   %       values from these can be retrieved using the fieldname
   %
   %       Example: Create a waveform, add a field, then get the field
   %           W = waveform;
   %           W = addfield(W,'P-pick', datenum('1/5/2007 02:12:15'));
   %           misc = get(W,'MISC_FIELDS'); %returns a cell with 'P-pick'
   %           PickTime = get(W,'P-pick'); % (same as) PickTime = get(W,misc);
   %
   %
   %   See also WAVEFORM/GETM, WAVEFORM/SET, WAVEFORM/ADDFIELD,
   %   WAVEFORM/DELFIELD

   %  11/28/2008 added WAVEFORM_OBJECT_VERSION, which returns a single version
   %  field from the waveform object.  This is set with the constructor, and
   %  is used in load_obj.

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   %global mep2dep

   %val_CELL = cell(size(w));
   %val = val_CELL;
   usedcell = false;
   singleWave = isscalar(w);
   switch upper(prop_name)

      % IDENTIFICATION PROPERTIES
      case {'STATION','CHANNEL','NETWORK','LOCATION'} %type:CELL
         val = get([w.cha_tag],prop_name);
         usedcell = true;

      case {'COMPONENT'} %type:CELL   GRANDFATHERED IN...
         val = get([w.cha_tag],'channel');
         usedcell = true;

         % DATA DESCRIBING PROPERTIES

      case {'FS', 'FREQ'} %type:DOUBLE
         val = [w.Fs];

      case 'DATA' %type:DOUBLE
         if singleWave
            val = w.data;
         else
            val = {w.data}; %cells of data values
            usedcell = true;
         end

      case 'NYQ' %type:DOUBLE
         val = [w.Fs] ./ 2;

      case 'PERIOD' %type:DOUBLE
         val = 1 ./ [w.Fs];

      case 'DATA_LENGTH' %type:DOUBLE
         val = zeros(size(w));
         for N = 1: numel(w)
            val(N) = numel(w(N).data);
         end

         % TIME PROPERTIES

      case {'START_STR'} %type:CELL
         val = stringDate([w.start]);
         usedcell = true;

      case {'START_EPOCH'} %type:DOUBLE
         val = mep2dep([w.start]);

      case {'START_MATLAB', 'START'} %type:DOUBLE
         val = [w.start];

      case {'END_STR'} %type:CELL
         val = stringDate(grabEndTime(w));
         usedcell = true;

      case {'END_EPOCH'} %type:DOUBLE
         val = mep2dep(grabEndTime(w));

      case {'END_MATLAB', 'END'} %type:DOUBLE
         val = grabEndTime(w);

      case {'DURATION_STR'} %type:CELL
         endvec = grabEndTime(w(:)) - [w(:).start];
         if numel(w) == 0
            val = durationDate(nan);
         elseif singleWave
            val = durationDate(endvec);
         else
            val = cell(size(w));
            for N = 1 : numel(w)
               val(N) = {durationDate(endvec(N))};
            end
            usedcell = true;
         end
      case {'DURATION_EPOCH'} %type:DOUBLE
         val = (grabEndTime(w) - [w.start]) * 86400;
         %epoch time is # of seconds.

      case {'DURATION_MATLAB', 'DURATION'} %type:DOUBLE
         val = grabEndTime(w) - [w.start];

      case {'TIMEVECTOR'} %type:CELL
         val = cell(size(w));
         for N = 1 : numel(w)
            Xvalues = linspace(w(N).start,get(w(N),'end'),length(w(N).data)+1);
            val{N} = Xvalues(1:end-1)';
         end
         usedcell = true;

      case {'MISC_FIELDS'} %type:CELL
         val = {w.misc_fields};
         usedcell = true;

      case {'UNITS'} %type : CELL
         val = {w.units};
         usedcell = true;

      case {'WAVEFORM_OBJECT_VERSION'}
         val = w(1).version;

      case {'SCNLOBJECT'} % legacy only
         % I'm making this whiny for now so it is obvious when things need to
         % be changed.  This warning should probably be pulled out once we
         % rerelease this.
         warning('waveform:get:scnlobjectsAreDepricated',...
            ['any usage of scnlobject should be replaced by ', ...
            '''channelinfo'' which uses the new ChannelTag class']);
         val = scnlobject([w.cha_tag]);
         %must add network & location, too.

      case {'CHANNELTAG'}
         val = [w.cha_tag];

      case {'CHANNELINFO'}
         val = string([w.cha_tag]);

      case {'HISTORY'}
         val = cell(size(w));
         for n=1:numel(w)
            val(n) = {w(n).history};
         end
         usedcell = true;

      otherwise
         %perhaps we're trying to get at one of the miscelleneous fields?
         val = cell(size(w));
         for n = 1 : numel(w)
            %loc is the position...
            %w(n).misc_fields should ALWAYS already be in uppercase
            mask = strcmp(upper(prop_name),w(n).misc_fields);
            %fieldwasfound = any(mask);
            %[fieldwasfound, loc] = ismember(prop_name, w(n).misc_fields);
            if any(mask)

               val{n} = w(n).misc_values{mask};
               %val{n} = w(n).misc_values{m};
            else
               error('Waveform:get:unrecognizedProperty',...
                  'Unrecognized property name : %s',  prop_name);
            end
         end
         %check to see if value can be returned as a numeric value instead
         %of cell.  Only if all values are numeric AND scalar
         numberize = true;
         for n=1:numel(val)
            if ~(isnumeric(val{n}) && isscalar(val{n}))
               numberize = false;
               usedcell = true;
               break
            end
         end
         if numberize,
            Z = val;
            val = nan(size(Z));
            for n=1:numel(Z)
               val(n) = Z{n};
            end
         end

   end;
   if (numel(val) == numel(w)) %THIS TEST CONDITIONALIZED 4/16/2008
      val = reshape(val,size(w)); %return values in proper shape
   end
   if usedcell && singleWave
      if isa(val,'cell')
         val = val{1}; % return the actual value, not a cell array
      end
   end;
end
%%%%%%%%%%%%%%%%
function val = grabEndTime(w)
   if isempty(w)
      %val = [];
      val = w.start;
      return
   end
   dlens = get(w,'data_length');
   dlens = dlens(:);

   myfrq = [w.Fs];
   myfrq = myfrq(:);

   seclen = dlens ./ myfrq;
   to_add = datenum([zeros(numel(w),5) seclen])';
   svals = [w(:).start];
   val = svals +  to_add;

   %endvec = datevec([w.start]) + [0 0 0 0 0 length(w.data)/w.Fs];
   %endvec = datevec([w.start]) + [0 0 0 0 0 get(w,'data_length') ./ get(w,'Freq')];
   %val = datenum(endvec);

end
%%%%%%%%%%%%%%%%
function val = stringDate(myDate)
   if isnan(myDate)
      warning('waveform:get:undefinedDate','Undefined date');
      val = 0;
      return
   end
   val = datestr(myDate,'yyyy-mm-dd HH:MM:SS.FFF');
end
%%%%%%%%%%%%%%%%
function val = durationDate(myDate)

   if isnan(myDate)
      %this likely happened because the frequency was undefined (NaN)
      val = '';
      return;
   end

   myDate = myDate(:);
   [yr, mo, da] = datevec(myDate);
   val = '';
   %if yr, val = [val sprintf('%2d years ',yr)]; end
   %if mo, val = [val sprintf('%2d months ',mo)]; end
   if (yr+mo+da)>1,
      val = [val sprintf('%4d days ',fix(myDate))];
   end
   val = [val datestr(myDate,'HH:MM:SS.FFF')];
end

author: geoscience-community-codes

repo: GISMO


text:
% $Id: WB08_Fig_2.m,v 1.1.1.1 2008/05/09 21:34:52 myself Exp $
%
% This prepares data for Figure 2 which is created with the corresponding
% shell script using GMT.  It also makes a plot in Matlab/Octave.
%
% Wessel, P. and J. M. Becker, 2008, Interpolation using a
%  generalized Green's function for a spherical surface spline
%  in tension, Geophys. J. Int., doi:10.1111/j.1365-246X.2008.03829.x
%
% Replicate Parker and then find the tension that minimizes the misfit at
% his 8 extra validation stations.

load mag_obs_1990.d
loni = mag_obs_1990(:,1);
lati = mag_obs_1990(:,2);
zi   = mag_obs_1990(:,3);
d=1;
% Set global 1x1 grid output coordinates

[X Y] = meshgrid (0:d:360, 0:d:90);

% First Parker's solution (p = 0)
Z = sphsplinet (loni, lati, zi, X, Y);

figure(1); clf
subplot (2,1,1)
contour (X, Y, Z)
drawnow
A = [X(:) Y(:) Z(:)];
save Fig_2_p0.d A -ascii -tabs

%Then used the wrong Oslo longitude to recreate Parker's figure
k = find (loni == 10.45)
loni(k) = 104.5;
Z = sphsplinet (loni, lati, zi, X, Y);
subplot (2,1,2)
contour (X, Y, Z)
A = [X(:) Y(:) Z(:)];
save Fig_2_orig.d A -ascii -tabs

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
clear,clc

addpath(pwd);
cd cvx;
addpath(genpath(pwd));
cd ..;

load('Ns=2.mat');
% load('Ns=8.mat');

Ns = 2;

NRF = [2,3,4,6,9,12,18];

SNR_dB = 0;
SNR = 10.^(SNR_dB./10);
realization = size(H,3);
smax = length(SNR);% enable the parallel

for r = 1:length(NRF)
    for reali = 1:realization
        [ FRF, FBB ] = SDR_AltMin( Fopt(:,:,reali), NRF(r) );
        [ WRF, WBB ] = Receiver( Wopt(:,:,reali), NRF(r) );
        R(r,reali) = log2(det(eye(Ns) + SNR/Ns * pinv(WRF * WBB) * H(:,:,reali) * FRF * FBB * FBB' * FRF' * H(:,:,reali)' * WRF * WBB));
    end
end
plot(NRF,sum(R,2)/realization,'Marker','diamond','LineWidth',1.5,'Color',[0.87058824300766 0.490196079015732 0]);
grid on
hold on

author: yuxianghao


text:
function [rpm] = MHz2rpm(MHz)
% Convert frequency from megahertz to revolutions per minute.
% Chad A. Greene 2012
rpm = MHz*60*1e+6;

author: Sable

repo: mcbench-benchmarks


text:
function [flag]=nonBlocking_isGoalReached(t_Kuka)
% this functions is used with nonBlocking motion functions to chekc whether
% the goal is reached.

%% Syntax:
% [flag]=nonBlocking_isGoalReached()

%% Arreguments
% t_Kuka: is the TCP/IP connection object

%% return value
% flag is a numerical value:
% flag=1: signifies the goal is reached
% flag=0: signifies the goal is not reached
% flag=-1: signifies an error

% Copy right: Mohammad SAFEEA 04th-April-2018

global paramName;
global paramVal;
x='eefpos';
%% if motion at EEF level
if(strcmp(x,paramName))
    % check size of array
    if(max(size(paramVal))==6)
    else
        flag=-1;
        return;
    end
% check if parameters are equal
flag=1;
    [eefpos,tempFlag]=GetActualEEFpos( t_Kuka );
    if tempFlag==false
        flag=1;
        return;
    end
    for i=1:6
        x1=paramVal{i};
        x2=eefpos{i};
        clearance=0.2; % 0.2 mm precission
        if((x1>(x2-clearance))&&(x1<(x2+clearance)))
        else
            flag=0;
            return;
        end
    end
    pause(0.1);
    return;
end
%% if motion at joint space level
x='jpos';
if(strcmp(x,paramName))
    % check size of array
    if(max(size(paramVal))==7)
    else
        flag=-1;
        return;
    end
    % check if parameters are equal
    flag=1;
    [eefpos,tempFlag]=GetActualJOINTSpos( t_Kuka );
    if tempFlag==false
        flag=1;
        return;
    end
    for i=1:6
        x1=paramVal{i};
        x2=eefpos{i};
        clearance=0.1*pi/180; % 0.1 degree precission
        if((x1>(x2-clearance))&&(x1<(x2+clearance)))
        else
            flag=0;
            return;
        end
    end
    pause(0.1);
    return;
end
end

function [ EEFpos,flag ] = GetActualEEFpos( t_Kuka )
EEEFpos={0,0,0,0,0,0};
theCommand='DcSeCarEEfP_';

for i=1:6
    x=EEEFpos{i};
    theCommand=[theCommand,sprintf('%0.5f',x),'_'];
end
% send the message through network
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);
[EEFpos,N]=getDoubleFromString(message);
if(N==0)
    flag=false;
    message=fgets(t_Kuka);
    [EEFpos,N]=getDoubleFromString(message);
else
    flag=true;
end
end

function [ EEFpos,flag ] = GetActualJOINTSpos( t_Kuka )
EEEFpos={0,0,0,0,0,0};
theCommand='DcSeCarJP_';

for i=1:6
    x=EEEFpos{i};
    theCommand=[theCommand,sprintf('%0.5f',x),'_'];
end
% send the message through network
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);
[EEFpos,N]=getDoubleFromString(message);
if(N==0)
    flag=false;
    message=fgets(t_Kuka);
    [EEFpos,N]=getDoubleFromString(message);
else
    flag=true;
end
end

function [jPos,j]=getDoubleFromString(message)
n=max(max(size(message)));
j=0;
numString=[];
for i=1:n
    if message(i)=='_'
        j=j+1;
        jPos{j}=str2num(numString);
        numString=[];
    else
        numString=[numString,message(i)];
    end
end
if j==0
    jPos={0};
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
% Copyright (C) 2012 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA
%
% You are free to use this software for academic purposes if you cite our paper:
% Quan Wang, Kim L. Boyer,
% The active geometric shape model: A new robust deformable shape model and its applications,
% Computer Vision and Image Understanding, Volume 116, Issue 12, December 2012, Pages 1178-1194,
% ISSN 1077-3142, 10.1016/j.cviu.2012.08.004.
%
% For commercial use, please contact the authors.

function GI = gaussianBlur(I,s)
%%  perform Gaussian blur
%   I: input image
%   s: standard deviation
%   GI: blurred image

M = gaussianMask(1,s);
if max(size(M))==0
    GI=I;
    return;
end
M = M/sum(sum(M));   % normalize the gaussian mask
GI=I;
for i=1:(2*s+1)
    GI=BoundMirrorExpand(GI);
end
GI = conv2(GI,M,'same');
for i=1:(2*s+1)
    GI=BoundMirrorShrink(GI);
end


text:
function msgString = getMessageString(messages)
% Returns a single string with the messages from cell array messages
if isempty(messages)
    msgString = '';
    return;
end;
msgString = num2str(messages{1});
for k = 2:length(messages)
    msgString = [msgString ' | ' num2str(messages{k})]; %#ok<AGROW>
end

author: VisLab

repo: EEG-Clean-Tools


text:
function [km2] = ha2km2(ha)
% Convert area from hectares to square kilometers.
% Chad A. Greene 2012
km2 = ha*0.01;