text:
% [y,feature_names,cache] = ComputeChangeWindowFeatures(x,...)
%
% Computes the change window features y for the input one-dimensional
% per-frame time series data x. y(i,t) will correspond to change
% window feature i and frame t, and is the (transformation of) the
% difference between the average of the per-frame data at the end and start
% of the window defined by i and t.
%
% For window feature i, let rho_i be the change window radius, r_i be the
% window radius, and off_i be the window offset. If the transformation type
% for feature i is 'none', the window feature y(i,t) is the difference
% between the average per-frame data at the end of the window (from
% (t+r_i+off_i)-rho_i to (t+r_i+off_i)+rho_i) and at the start of the
% window (from (t-r_i+off_i)-rho_i to (t-r_i+off_i)+rho_i). If the
% transformation type for feature i is 'abs', then we compute the absolute
% value of this difference. If the transformation type for feature
% i is 'flip', then we compute the sign of x(t) times this difference.
%
% Input:
%
% x: 1 x NFRAMES array of per-frame data.
%
% Output:
%
% y: NWINDOWFEATURES x NFRAMES matrix of window data, where y(i,t)
% corresponds to window feature i and frame t. i indexes the radius and
% offset of the window as well as the transformation type.
% feature_names: 1 x NWINDOWFEATURES cell in which feature_names{i}
% describes the ith window feature computed. feature_names{i} is itself a
% cell that can be interpreted as pairs of a string description followed by
% a value, e.g.
% {'stat','change','trans','abs','radius',1,'offset',1,'change_window_radius',0}
% cache: if DOCACHE is set to true, then computations useful in other
% window feature calculations may be saved in the output cache.
%
% Optional inputs:
%
% 'change_window_radii': The 1 X NCHANGEWINDOWRADII change window radii to
% use (this defines the size of the sub-window averaged over at the end and
% start of the original window). default value: 0.
%
% DEFAULT WINDOW LOCATIONS:
% These window locations are used if window locations are not specified on
% a per-feature type basis. Default default values set by
% SetDefaultWindowParameters.
% Inputs interpreted by SetWindowParameters.
% The same window parameter interpretation is done in all
% Compute*WindowFeatures functions.
% 'windows': window offsets and radii to try ([radius1,offset1];...;[radiusn,offsetn])
% if empty, then window_radii and window_offsets are used to specify
% windows. if the window is [r,off], then the window feature at time t
% will be computed from the window at t-r+off to t+r+off. default value: [].
% 'window_radii': if windows is empty, then the cross-product of window_radii
% and window_offsets is used to set windows. if empty, then
% min_window_radius, max_window_radius, and nwindow_radii are used to set
% window_radii. default value: [].
% 'window_offsets': if windows is empty, then the cross-product of  window_radii
% and window_offsets is used to set windows. window_offsets are relative to
% radius, so the window corresponding to radius r_i and offset off_j_rel
% will be [r_i,off_i=r_i*off_j_rel]. default value: [-1,0,1].
% 'min_window_radius': if windows and window_radii are both empty, then
% window_radii is set to nwindow_radii evenly spaced radii between
% min_window_radii and max_window_radii. default value: 0
% 'max_window_radius': see 'min_window_radius'. default value: 20.
% 'nwindow_radii: see 'min_window_radius'. default value: 5.
%
% 'trans_types': default types of transformations to apply for each feature
% type, if not otherwise specified. Options include 'abs','flip', and
% 'none'. 'abs' corresponds to the absolute value, flip corresponds to
% flipping the sign of the feature if the per-frame feature at the frame t
% is negative, and 'none' corresponds to no transformation.
%
% 'sanitycheck': whether to compute all the features a second time in the
% obvious way to make sure that the optimized computations are correct.
% default value: false.
%
% 'docache': whether to cache computations that might be useful, e.g. the
% mean computations are useful when computing the change window features.
% default value: true.
%
% 'cache': input cached computations that might be useful in this
% computation.

function [y,feature_names,cache] = ComputeChangeWindowFeatures(x,varargin)

x = x(:)';
N = numel(x);
y = nan(0,N);
feature_names = {};

%% default parameters

[windows,window_offsets,...
  window_radii,min_window_radius,...
  max_window_radius,nwindow_radii] = ...
  SetDefaultWindowParameters();

% use all transformation types by default
%trans_types = 'all';
trans_types = uint8(15);

% for debugging purposes
SANITY_CHECK = true;

% whether to cache results
DOCACHE = true;

% initialize empty cache
cache = InitializeCache();

% initialize feature_types already computed to empty
%feature_types = {};

% change window radii (width = 2*radius + 1)
change_window_radii = 0;

relativeParams = [];

%% parse parameters

[...
  windows,...
  window_radii,window_offsets,...
  min_window_radius,max_window_radius,nwindow_radii,...
  trans_types,...
  SANITY_CHECK,...
  DOCACHE,...
  cache,...
  change_window_radii,...
  relativeParams,...
  ] = myparse(varargin,...
  'windows',windows,...
  'window_radii',window_radii,'window_offsets',window_offsets,...
  'min_window_radius',min_window_radius,'max_window_radius',max_window_radius,'nwindow_radii',nwindow_radii,...
  'trans_types',trans_types,...
  'sanitycheck',SANITY_CHECK,...
  'docache',DOCACHE,...
  'cache',cache,...
  'change_window_radii',change_window_radii,...
  'relativeParams',relativeParams);
%   'feature_types',feature_types,...

%% whether we've specified to use all trans types by default
%if ischar(trans_types) && strcmpi(trans_types,'all'),
%  trans_types = {'none','abs','flip','relative'};
%end

%% select default windows from various ways of specifying windows

[windows,window_radii,windowi2radiusi,nradii] = ...
  SetWindowParameters(...
  windows,window_offsets,...
  window_radii,min_window_radius,...
  max_window_radius,nwindow_radii);

%% main computation

%if ismember('relative',trans_types)
if bitand(8,trans_types)
  if DOCACHE && ~isempty(cache.relX)
    modX = cache.relX;
  else
    modX = convertToRelative(x,relativeParams);
    cache.relX = modX;
  end
end

for change_r_i = 1:numel(change_window_radii),

  % take the mean across all windows of radius change_r
  change_r = change_window_radii(change_r_i);
  change_w = 2*change_r+1;

  if DOCACHE && ismember(change_r,cache.mean.radii),
    cache_i = find(change_r == cache.mean.radii,1);
    res_mean = cache.mean.data{cache_i};
  else
    % average filter
    res_mean = MeanWindowCore(x,change_w);
    if DOCACHE,
      cache.mean.radii(end+1) = change_r;
      cache.mean.data{end+1} = res_mean;
    end
  end

%  if ismember('relative',trans_types)
  if bitand(8,trans_types)
    if DOCACHE && ismember(change_r,cache.meanRel.radii),
      cache_i = find(change_r == cache.meanRel.radii,1);
      resRel_mean = cache.meanRel.data{cache_i};
    else
      resRel_mean = MeanWindowCore(modX,change_w);
      % store for future computations
      if DOCACHE,
        cache.meanRel.radii(end+1) = change_r;
        cache.meanRel.data{end+1} = resRel_mean;
      end
    end
  end

  % loop over window radii
  for radiusi = 1:nradii,
    r = window_radii(radiusi);
    % don't do for r <= change_r -- the windows will overlap, cancel, ...
    if r <= change_r,
      continue;
    end
    w = 2*r+1;

    % take the difference between the end of the window and the start of
    % the window
    % res(t) corresponds to t+r-change_r
    % so res(t-r+change_r) corresponds to frame t
    res = res_mean(w:end) - res_mean(1:end-w+1);

%    if ismember('relative',trans_types),
    if bitand(8,trans_types),
      resRel = resRel_mean(w:end) - resRel_mean(1:end-w+1);
    end

    % offset
    windowis = find(windowi2radiusi == radiusi);
    for windowi = windowis',

      off = windows(windowi,2);
      % res(t-r+change_r) corresponds to frame t,
      % and we want to grab from [1+off,N+off] which is
      % [1+off-r+change_r,N+off-r+change_r], relative to res
      res1 = padgrab2(res,nan,1,1,1+off-r+change_r,N+off-r+change_r)/r;

%      if ismember('none',trans_types),
      if bitand(1,trans_types),
        y(end+1,:) = res1;
        feature_names{end+1} = {'stat','change','trans','none','radius',r,'offset',off,'change_window_radius',change_r}; %#ok<*AGROW>
      end

%      if ismember('abs',trans_types),
      if bitand(2,trans_types),
          y(end+1,:) = abs(res1);
          feature_names{end+1} = {'stat','change','trans','abs','radius',r,'offset',off,'change_window_radius',change_r};
      end

%      if ismember('flip',trans_types)
      if bitand(4,trans_types)
        res2 = res1;
        res2(x<0) = -res2(x<0);
        y(end+1,:) = res2;
        feature_names{end+1} = {'stat','change','trans','flip','radius',r,'offset',off,'change_window_radius',change_r};
      end

%      if ismember('relative',trans_types),
      if bitand(8,trans_types),
        resRel1 = padgrab2(resRel,nan,1,1,1+off-r+change_r,N+off-r+change_r)/r;
        y(end+1,:) = resRel1;
        feature_names{end+1} = {'stat','change','trans','relative','radius',r,'offset',off,'change_window_radius',change_r}; %#ok<*AGROW>
      end

      if SANITY_CHECK,
        extraStr = sprintf('change_r = %d',change_r);
%        if ismember('none',trans_types),
        if bitand(1,trans_types),
          fastY = res1;
        end
        res_dumb = nan(1,N);
        for n_dumb = 1:N,
          tmp1 = padgrab2(x,nan,1,1,n_dumb+off-r-change_r,n_dumb+off-r+change_r);
          tmp2 = padgrab2(x,nan,1,1,n_dumb+off+r-change_r,n_dumb+off+r+change_r);
          res_dumb(n_dumb) = (nanmean(tmp2) - nanmean(tmp1))/r;
        end
        checkSanity(fastY,res_dumb,r,off,'change','none',extraStr);

%        if ismember('abs',trans_types),
        if bitand(2,trans_types),
          fastY = abs(res1);
          res_dumb = nan(1,N);
          for n_dumb = 1:N,
            tmp1 = padgrab2(x,nan,1,1,n_dumb+off-r-change_r,n_dumb+off-r+change_r);
            tmp2 = padgrab2(x,nan,1,1,n_dumb+off+r-change_r,n_dumb+off+r+change_r);
            res_dumb(n_dumb) = abs(nanmean(tmp2) - nanmean(tmp1))/r;
          end
          checkSanity(fastY,res_dumb,r,off,'change','abs',extraStr);
        end

%        if ismember('flip',trans_types)
        if bitand(4,trans_types)
          res2 = res1;
          res2(x<0) = -res2(x<0);
          fastY = res2;
          res_dumb = nan(1,N);
          for n_dumb = 1:N,
            tmp1 = padgrab2(x,nan,1,1,n_dumb+off-r-change_r,n_dumb+off-r+change_r);
            tmp2 = padgrab2(x,nan,1,1,n_dumb+off+r-change_r,n_dumb+off+r+change_r);
            res_dumb(n_dumb) = (nanmean(tmp2) - nanmean(tmp1))/r;
            if x(n_dumb) < 0,
              res_dumb(n_dumb) = -res_dumb(n_dumb);
            end
          end
          checkSanity(fastY,res_dumb,r,off,'change','flip',extraStr);
        end

      end

    end
  end

end

author: kristinbranson

repo: JAABA


text:
function val = get(TC,varargin)

%GET    Get threecomp object properties.
%   V = GET(TC,'PropertyName') returns the value of the specified
%   property for the threecomp object specified by TC. Allowable properties
%   include:
%   WAVEFORM          waveform object containing the three component traces
%   BACKAZIMUTH       backazimuth to source
%   RECTILINEARITY    rectilinearity stored in a waveform object
%   PLANARITY         stored in a waveform object
%   ENERGY            energy stored in a waveform object
%   AZIMUTH           azimuth stored in a waveform object
%   INCLINATION       inclination stored in a waveform object
% See also threecomp object

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% Check inputs
if ~isa(TC,'threecomp')
    error('First argment must be a threecomp object');
end
if nargin~=2
    error('Incorrect number of arguments');
end
property = varargin{1};
if ~ischar(property)
    error('Second argument must be a character string');
end

[objSize1,objSize2] = size(TC);
numObj = numel(TC);

% Prepare property values
switch upper(property)

    case 'WAVEFORM'
        if numObj==1
            val = TC.traces;
        elseif objSize1>1 && objSize2==1
            for n = 1:numObj
                val(n,1:3) = TC(n).traces;
            end
        elseif objSize2>1
            for n1 = 1:objSize1
                for n2 = 1:objSize2
                    val(n1,n2,1:3) = TC(n).traces;
                end
            end
        else
            error('threecomp object is not 1x1, nx1, or nxm in size');
        end

    case 'TRIGGER'
        for n = 1:numObj
           val(n) = TC(n).trigger;
        end
        val = reshape(val,objSize1,objSize2);

    case 'BACKAZIMUTH'
        for n = 1:numObj
            if isempty(TC(n).backAzimuth)
                val(n) = NaN;
            else
                val(n) = TC(n).backAzimuth;
            end
        end
        val = reshape(val,objSize1,objSize2);
    case 'NSCL'
        for n = 1:numObj
            w = TC(n).waveform;
            net = get(w(1),'NETWORK');
            sta = get(w(1),'STATION');
            chan = get(w,'CHANNEL');
            chanString = [chan{1}(1:2) '[' chan{1}(3) chan{2}(3) chan{3}(3) ']' ];
            loc = get(w(1),'LOCATION');
            val(n) = {[net '_' sta '_' chanString '_' loc]};
        end
        val = reshape(val,objSize1,objSize2);
    case 'STATION'
        for n = 1:numObj
            w = TC(n).waveform;
            val(n) = {get(w(1),'STATION')};
        end
        val = reshape(val,objSize1,objSize2);
    case 'CHANNEL'
        for n = 1:numObj
            w = TC(n).waveform;
            val(n) = {get(w,'CHANNEL')};
        end
        val = reshape(val,objSize1,objSize2);
    case 'TYPE'
        for n = 1:numObj
            val(n) = gettype(TC(n));
        end
        val = reshape(val,objSize1,objSize2);

    case 'COMPLETENESS'
        for n = 1:numObj
            val(n) = getcompleteness(TC(n));
        end
        val = reshape(val,objSize1,objSize2);

    case 'ORIENTATION'
        for n = 1:numObj
            vals = TC(n).orientation;
            if isempty(vals)
               vals = [NaN NaN NaN NaN NaN NaN];
            end
            val(n,1:6) = vals;
        end
    case 'RECTILINEARITY'
        for n = 1:numObj
            if isempty(TC(n).rectilinearity)
                val(n) = waveform;
            else
                val(n) = TC(n).rectilinearity;
            end
        end
        val = reshape(val,objSize1,objSize2);

    case 'PLANARITY'
        for n = 1:numObj
            if isempty(TC(n).planarity)
                val(n) = waveform;
            else
                val(n) = TC(n).planarity;
            end
        end
        val = reshape(val,objSize1,objSize2);

    case 'ENERGY'
        for n = 1:numObj
            if isempty(TC(n).energy)
                val(n) = waveform;
            else
                val(n) = TC(n).energy;
            end
        end
        val = reshape(val,objSize1,objSize2);

    case 'AZIMUTH'
        for n = 1:numObj
            if isempty(TC(n).azimuth)
                val(n) = waveform;
            else
                val(n) = TC(n).azimuth;
            end
        end
        val = reshape(val,objSize1,objSize2);

    case 'INCLINATION'
        for n = 1:numObj
            if isempty(TC(n).inclination)
                val(n) = waveform;
            else
                val(n) = TC(n).inclination;
            end
        end
        val = reshape(val,objSize1,objSize2);

    otherwise
        error('Property not recognized.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DETERMINE TYPE PROPERTY
% where TC and val are scalars
%
function val = gettype(TC)

chan = get(TC,'CHANNEL');
chan = chan{1};
for n = 1:3
   val(1) = {[ chan{1}(3) chan{2}(3) chan{3}(3) ]};
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DETERMINE COMPLETENESS PROPERTY
% where TC and val are scalars
%
function val = getcompleteness(TC)

% 0 = no waveforms and no particle motions
% 1 = waveforms but no particle motions
% 2 = waves and all particle motion fields

val = 0;

if all(~isempty(TC.traces)) && ~isempty(TC.rectilinearity) && ~isempty(TC.planarity) && ~isempty(TC.energy) && ~isempty(TC.azimuth)  && ~isempty(TC.inclination)
    val = 2;
elseif all(~isempty(TC.traces))
   val = 1;
else
    val = 0;
end

author: geoscience-community-codes

repo: GISMO


text:
function sVF1 = rdivide(sVF1, sVF2)
%
% Syntax
%
%   sVF = sVF ./ a
%   sVF = sVF ./ sF
%
% Input
%  sVF - @S2VectorFieldHarmonic
%  sF  - @S2Fun
%  a   - double
%
% Output
%  sVF - @S2VectorFieldHarmonic
%

if isnumeric(sVF2)
  sVF1.sF = sVF1.sF .* (1./sVF2);
elseif isa(sVF2,'S2Fun')
  sVF1.sF = sVF1.sF ./ sVF2;
end

end

author: mtex-toolbox

repo: mtex

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


text:
% PHONE_SRC Creates a source for the phones in the TIMIT dataset
%
% Usage
%    src = PHONE_SRC(directory)
%
% Input
%    directory (char): The directory containing the TIMIT dataset.
%
% Output
%    src (struct): The TIMIT phone source.
%
% Description
%    The function searches through the directory structure, indexing each
%    audio file by adding it to the src structure and its constituent phones
%    as separate objects contained in that file. Phone types are reduced to
%    48 different phones [1]. Additionally, the clustering specified in [1] is
%    recorded in src.cluster. This is then used in when calculating the
%    classification error to allow confusions within the same cluster.
%
% References
%    [1] K.F. Lee and H-W. Hon, "Speaker-Independent Phone Recognition Using
%       Hidden Markov Models," IEEE TASSP 1989
%
% See also
%    CREATE_SRC, PHONE_PARTITION

function src = phone_src(directory)
	if nargin < 1
		directory = 'timit';
	end

	[phones,original,mapping] = phone_list(0);

	src = create_src(directory,@(x)(phone_objects_fun(x,phones,original,mapping)));

	[temp1,temp2,fold_mapping] = phone_list(1);

	[temp,reorder] = sort(phones);

	[temp,I] = ismember(reorder,mapping);
	src.cluster = fold_mapping(I);
end

function [segments,classes] = phone_objects_fun(file,phones,original,mapping)
	[path_str,name,ext] = fileparts(file);

	if strcmp(name,'SA1') || strcmp(name,'SA2')
		segments = [];
		classes = {};

		return;
	end

	set = phrase_set(file);

	if set == 3		% test but not in core, dev
		segments = [];
		classes = {};

		return;
	end

	phn = parse_phn([path_str filesep name '.PHN'],original,mapping);

	segments = struct();
	classes = {};

	for n = 1:length(phn)
		segments(n).u1 = phn(n).t1;
		segments(n).u2 = phn(n).t2;
		segments(n).subset = set;
		classes{n} = phones{phn(n).phone};
	end
end

function [phones,original,mapping] = phone_list(do_fold)
	if nargin < 1
		do_fold = 1;
	end

	stops = {'b','d','g','p','t','k','dx','q'};
	closures = {'cl','vcl','bcl','dcl','gcl','pcl','tcl','kcl','qcl'};
	affricates = {'jh','ch'};
	fricatives = {'s','sh','z','zh','f','th','v','dh'};
	nasals = {'m','n','ng','em','en','eng','nx'};
	semivowels_glides = {'l','r','w','y','hh','hv','el'};
	vowels = {'iy','ih','eh','ey','ae','aa','aw','ay','ah','ao','oy','ow','uh','uw','ux','er','ax','ix','axr','ax-h'};
	others = {'pau','epi','h#'};

	original = [stops closures affricates fricatives nasals semivowels_glides vowels others];

	% Kai-Fu Lee clustering, see K.F. Lee and H-W. Hon, "Speaker-Independent Phone Recognition Using Hidden Markov Models," IEEE TASSP 1989

	allophones = {{'uw','ux'},{'er','axr'},{'m','em'},{'n','nx'},{'ng','eng'},{'hh','hv'},{'cl','pcl','tcl','kcl','qcl'},{'vcl','bcl','dcl','gcl'},{'sil','h#','pau'}};

	folding = {{'sil','cl','vcl','epi'},{'el','l'},{'en','n'},{'sh','zh'},{'ao','aa'},{'ih','ix'},{'ah','ax'}};

	removed = {'q','ax-h'};

	phones = {};

	for k = 1:length(original)
		phone = original{k};
		for r = 1:length(allophones)
			if ~isempty(strmatch(phone,allophones{r},'exact'))
				phone = allophones{r}{1};
				break;
			end
		end
		if do_fold
			for r = 1:length(folding)
				if ~isempty(strmatch(phone,folding{r},'exact'))
					phone = folding{r}{1};
					break;
				end
			end
		end
		if ~isempty(strmatch(phone,removed,'exact'))
			phone = [];
		end
		if ~isempty(phone) && isempty(strmatch(phone,phones,'exact'))
			phones = [phones phone];
		end
		mapped{k} = phone;
	end

	for k = 1:length(mapped)
		if isempty(mapped{k})
			mapping(k) = 0;
		else
			ind = strmatch(mapped{k},phones,'exact');
			mapping(k) = ind;
		end
	end
end

function phn = parse_phn(filename,original,mapping)
	txt = fileread(filename);

	lines = regexp(txt,'\n','split');

	n = 1;
	phn = struct();

	for k = 1:length(lines)-1
		parts = regexp(lines{k},' ','split');
		phone_label = parts{3};
		phone = mapping(find(strcmp(phone_label,original),1));
		if phone == 0, continue; end
		if isempty(phone)
			error('empty phone');
		end

		phn(n).t1 = str2double(parts{1});
		phn(n).t2 = str2double(parts{2});
		phn(n).phone = phone;
		n = n+1;
	end
end

function set = phrase_set(file)
	core_set = {'FDHC0','FMGD0','FPAS0','MCMJ0','MJDH0','MKLT0','MNJM0','MTLS0', ...
				'FELC0','FMLD0','FPKT0','MDAB0','MJLN0','MLLL0','MPAM0','MWBT0', ...
				'FJLM0','FNLP0','MBPM0','MGRT0','MJMP0','MLNT0','MTAS1','MWEW0'};

	dev_set = {'FAKS0','MMDB1','MBDG0','FEDW0','MTDT0','FSEM0','MDVC0','MRJM4','MJSW0','MTEB0', ...
			   'FDAC1','MMDM2','MBWM0','MGJF0','MTHC0','MBNS0','MERS0','FCAL1','MREB0','MJFC0', ...
			   'FJEM0','MPDF0','MCSH0','MGLB0','MWJG0','MMJR0','FMAH0','MMWH0','FGJD0','MRJR0', ...
			   'MGWT0','FCMH0','FADG0','MRTK0','FNMR0','MDLS0','FDRW0','FJSJ0','FJMG0','FMML0', ...
			   'MJAR0','FKMS0','FDMS0','MTAA0','FREW0','MDLF0','MRCS0','MAJC0','MROA0','MRWS1'};

	path_str = fileparts(file);

	path_parts = regexp(path_str,filesep,'split');

	if strcmp(path_parts{end-2},'TRAIN')
		set = 0;
	elseif strcmp(path_parts{end-2},'TEST')
		if any(strcmp(path_parts{end},core_set))
			set = 1;
		elseif any(strcmp(path_parts{end},dev_set))
			set = 2;
		else
			set = 3;
		end
	else
		set = [];
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet