Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

isolated matlab_htk from original monolithic ronw-matlab-tools reposi…

…tory
  • Loading branch information...
commit 090e264cf579e115db3be29715700dd4934feefa 1 parent e14a6df
Ron Weiss authored
Showing with 82 additions and 3,681 deletions.
  1. +82 −16 README
  2. 0  {matlab_htk → }/apply_mllr_transform.m
  3. +0 −6 celltools/@lazymap/display.m
  4. +0 −7 celltools/@lazymap/end.m
  5. +0 −3  celltools/@lazymap/iscell.m
  6. +0 −59 celltools/@lazymap/lazymap.m
  7. +0 −3  celltools/@lazymap/length.m
  8. +0 −3  celltools/@lazymap/numel.m
  9. +0 −15 celltools/@lazymap/set.m
  10. +0 −7 celltools/@lazymap/size.m
  11. +0 −17 celltools/@lazymap/subsref.m
  12. +0 −5 celltools/README
  13. +0 −14 celltools/cellcat.m
  14. +0 −35 celltools/cfilter.m
  15. +0 −25 celltools/crange.m
  16. +0 −31 celltools/czip.m
  17. +0 −39 celltools/flatten.m
  18. +0 −61 celltools/glob.m
  19. +0 −38 celltools/interleave.m
  20. +0 −45 celltools/map.m
  21. +0 −45 celltools/mapreduce.m
  22. +0 −51 celltools/reduce.m
  23. +0 −46 celltools/sfilter.m
  24. 0  {matlab_htk → }/compose_hmms.m
  25. 0  {matlab_htk → }/eval_htk_recognizer.m
  26. 0  {matlab_htk → }/get_htk_path.m
  27. 0  {matlab_htk → }/get_temp_filename.m
  28. +0 −60 hmm/README
  29. +0 −129 hmm/adapt_gmm.m
  30. +0 −67 hmm/convert_hmm_to_gaussian_emissions.m
  31. +0 −218 hmm/decode_hmm.m
  32. +0 −54 hmm/eval_gmm.m
  33. +0 −195 hmm/eval_hmm.m
  34. +0 −82 hmm/is_valid_gmm.m
  35. +0 −133 hmm/is_valid_hmm.m
  36. +0 −67 hmm/lmvnpdf.m
  37. +0 −15 hmm/logsum.m
  38. +0 −113 hmm/merge_states.m
  39. +0 −56 hmm/reorder_states.m
  40. +0 −24 hmm/sample_gaussian.m
  41. +0 −34 hmm/sample_gmm.m
  42. +0 −61 hmm/sample_hmm.m
  43. +0 −101 hmm/train_gmm.m
  44. 0  {matlab_htk → }/htkread.m
  45. 0  {matlab_htk → }/htkwrite.m
  46. 0  {matlab_htk → }/kmeans.m
  47. 0  {matlab_htk → }/learn_mllr_transform.m
  48. 0  {matlab_htk → }/logsum.m
  49. +0 −77 matlab_htk/README
  50. +0 −288 plottools/add_pan_and_zoom_controls_to_figure.m
  51. +0 −24 plottools/align_axes.m
  52. +0 −34 plottools/plot_on_same_axes.m
  53. +0 −39 plottools/plot_or_imagesc.m
  54. +0 −134 plottools/plot_pages.m
  55. +0 −387 plottools/plotall.m
  56. +0 −422 plottools/zooming_scrollbar.m
  57. +0 −132 process_options.m
  58. 0  {matlab_htk → }/read_htk_hmm.m
  59. 0  {matlab_htk → }/read_text_file.m
  60. +0 −164 source_sep/maxvq_bb.m
  61. 0  {matlab_htk → }/train_gmm_htk.m
  62. 0  {matlab_htk → }/train_hmm_htk.m
  63. 0  {matlab_htk → }/train_htk_recognizer.m
  64. 0  {matlab_htk → }/train_htk_recognizer.sh
  65. 0  {matlab_htk → }/write_htk_bnf.m
  66. 0  {matlab_htk → }/write_htk_hmm.m
  67. 0  {matlab_htk → }/write_htk_slf.m
  68. 0  {matlab_htk → }/write_text_file.m
98 README
View
@@ -1,17 +1,83 @@
-This directory contains a variety of Matlab functions that I've found
-to be useful for my work at LabROSA. The portions written by me are
-distributed under the terms of the GNU General Public License. See
-the file COPYING for details.
-
-Contents:
-hmm - functions for evaluating hidden Markov models and Gaussian
- mixture models
-matlab_htk - functions that interface with the HTK Speech Recognition
- Toolkit (http://htk.eng.cam.ac.uk/) for training
- HMMs, GMMs and simple speech recognizers
-source_sep - implementations of standard algorithms for source
- separation
-celltools - useful functional programming style functions for cell arrays
-plottools - Wrappers around basic plotting functions and code for some
- handy new GUI controls.
+This package contains a set of functions for calling interfacing with
+HTK from Matlab. Right now its mostly limited to training GMMs and
+HMMs. It converts your Matlab data into a format that HTK understands
+and calls HTK command line programs. The path to the HTK binaries is
+hardcoded in get_htk_path.m
+The portions written by me are distributed under the terms of the GNU
+General Public License. See the file COPYING for details.
+
+Unforuntately this has not been very thoroughly tested, but it has
+been working for me. If anyone wants to write some unit tests using
+http://mlunit.sf.net or something similar I won't say no.
+
+
+* Functions
+
+They are all reasonably commented...
+
+ - The important ones:
+ - train_hmm_htk - use HTK to train an HMM
+ - train_gmm_htk - use HTK to train a GMM
+ - train_htk_recognizer - train a simple HTK recognizer
+ - eval_htk_recognizer - recognize signal using HTK recognizer
+
+ - Utility functions:
+ - read_htk_hmm - read in an HTK HMM definition (text format only)
+ - write_htk_hmm - write an HMM structure as an HTK HMM definition
+ - htkread/htkwrite - read/write a Matlab matrix in HTK feature file format
+ (written by Mark Hasegawa-Johnson)
+ - compose_hmms - does FSM composition to form a big HMM from many
+ small HMMs (i.e. get an HMM to recognize a
+ sentence from a bunch of phone HMMs).
+ - kmeans - uses k-means to learn clusters from data. Used to
+ initialize GMMs and HMMs.
+ - logsum - takes the sum of a matrix of log likelihoods
+ - get_htk_path - centralized location to set the path to the HTK binaries
+
+* Data Structures
+
+The functions in this toolbox pass around the following structures:
+Note: all probabilities are stored as log probabilities
+
+** GMM
+ - gmm.nmix - number of components in the mixture
+ - gmm.priors - array of prior log probabilities over each state
+ - gmm.means - matrix of means (column x is mean of component x)
+ - gmm.covars - matrix of covariance (column x is the diagonal of the
+ covariance matrix of component x)
+
+** HMM with GMM observations
+ - hmm.name -
+ - hmm.nstates - number of states in the HMM
+ - hmm.emission_type - 'GMM'
+ - hmm.start_prob - array of log probs P(first observation is state x)
+ - hmm.end_prob - array of log probs P(last observation is state x)
+ - hmm.transmat - matrix of transition log probs (transmat(x,y)
+ = log(P(transition from state x to state y)))
+ - hmm.labels - optional cell array of labels for each state in the HMM
+ (for use in composing HMMs)
+ - hmm.gmms - array of GMM structures
+
+** HMM with Gaussian observations
+ - hmm.nstates - number of states in the HMM
+ - hmm.emission_type - 'gaussian'
+ - hmm.start_prob - array of log probs P(first observation is state x)
+ - hmm.end_prob - array of log probs P(last observation is state x)
+ - hmm.transmat - matrix of transition log probs (transmat(x,y)
+ = log(P(transition from state x to state y)))
+ - hmm.labels - optional cell array of labels for each state in the HMM
+ (for use in composing HMMs)
+ - hmm.means - matrix of means (column x is mean of state x)
+ - hmm.covars - matrix of means (column x is the diagonal of the
+ covariance matrix of component x)
+
+Note that each row of exp(hmm.transmat) does not necessarily sum to 1
+because for each state x there is some probability
+(exp(hmm.exit_prob(x))) that the next transition will be to a
+non-emitting exit state (i.e. the current observation is the last
+observation in the sequence). The correct invariant is:
+sum(exp(hmm.transmat, 2)) + exp(hmm.exit_prob) == ones(hmm.nstates, 1)
+
+
+2007-11-06 Ron Weiss <ronw@ee.columbia.edu>
0  matlab_htk/apply_mllr_transform.m → apply_mllr_transform.m
View
File renamed without changes
6 celltools/@lazymap/display.m
View
@@ -1,6 +0,0 @@
-function display(C)
-
-disp([inputname(1) ' = '])
-disp([' lazymap object:'])
-disp([' size: ' mat2str(size(C))])
-disp([' map function: ' func2str(C.func)])
7 celltools/@lazymap/end.m
View
@@ -1,7 +0,0 @@
-function e = end(C, k, n)
-
-if n == 1
- e = numel(C.cellarray);
-else
- e = size(C.cellarray, k);
-end
3  celltools/@lazymap/iscell.m
View
@@ -1,3 +0,0 @@
-function tf = iscell(C)
-
-tf = 1;
59 celltools/@lazymap/lazymap.m
View
@@ -1,59 +0,0 @@
-function Y = lazymap(fun, C)
-% Y = lazymap(fun, C)
-%
-% Behaves like map/cellfun, but with lazy evaluation. Returns an
-% object that behaves like a cell array with same dimensions as C.
-% Whenever an element in y is referenced (e.g. y(10) or y{1:end}) the
-% function fun is applied to the corresponding cells of C and the
-% result is returned.
-%
-% Example:
-% numbers = {0, 0.5, 1, pi/2, pi, 2*pi};
-% y = lazymap(@cos, numbers);
-% y(10) % returns cos(numbers(4))
-% y = set(y, 'Function', @sin);
-% y{10} % returns sin(numbers(4))
-% y{1:5} % returns {sin(numbers(1)), sin(numbers(2)), ..., sin(numbers(5))}
-%
-% This lazy behavior is useful when only some elements of fun(C) are
-% needed at any given time. For example, when extracting features
-% from a set of files, a call like:
-% features = lazymap(@(x) compute_features(x, param1, param2), filenames)
-% can be used to give convenient access to the features of a large
-% data set without having to store all of it in memory.
-%
-% Note that this class doesn't support memoization so subsequent
-% references to the same element, e.g. y{1}; y{1};, will result in
-% multiple calls to fun(C{1}).
-%
-% 2007-11-01 ronw@ee.columbia.edu
-
-% Copyright (C) 2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if ~isa(fun, 'function_handle')
- error('1st argument must be a function handle')
-end
-
-if ~iscell(C)
- error('2nd argument must be a cell array.')
-end
-
-Y.cellarray = C;
-Y.func = fun;
-
-% Can't inherit from built in types...
-% y = class(y, 'lazymap', cellarray);
-Y = class(Y, 'lazymap');
3  celltools/@lazymap/length.m
View
@@ -1,3 +0,0 @@
-function n = length(C)
-
-n = length(C.cellarray);
3  celltools/@lazymap/numel.m
View
@@ -1,3 +0,0 @@
-function n = numel(C, varargin);
-
-n = numel(C.cellarray, varargin{:});
15 celltools/@lazymap/set.m
View
@@ -1,15 +0,0 @@
-function C = set(C, varargin)
-
-propertyArgIn = varargin;
-while length(propertyArgIn) >= 2,
- prop = propertyArgIn{1};
- val = propertyArgIn{2};
- propertyArgIn = propertyArgIn(3:end);
-
- switch prop
- case 'Function'
- C.func = val;
- otherwise
- error('Invalid property')
- end
-end
7 celltools/@lazymap/size.m
View
@@ -1,7 +0,0 @@
-function varargout = size(C, dim)
-
-if nargin == 1
- [varargout{1:nargout}] = size(C.cellarray);
-else
- varargout{1} = size(C.cellarray, dim);
-end
17 celltools/@lazymap/subsref.m
View
@@ -1,17 +0,0 @@
-function varargout = subsref(C, S);
-
-if S.type == '.'
- error('Attempt to reference field of non-structure array.');
-end
-
-tmpS = S;
-tmpS.type = '()';
-B = subsref(C.cellarray, tmpS);
-B = cellfun(C.func, B, 'UniformOutput', 0);
-
-if S.type == '{}'
- varargout = B;
-else
- varargout{1} = B;
-end
-
5 celltools/README
View
@@ -1,5 +0,0 @@
-This package contains a set of standard functional programming tools
-that can be used on cell arrays. Contains implementations of map,
-reduce, filter (called cfilter so as not to conflict with the filter()
-function in the Signal Processing Toolbox), and lazymap which has the
-same semantics as map but uses lazy evaluation.
14 celltools/cellcat.m
View
@@ -1,14 +0,0 @@
-function c = cellcat(varargin)
-% c = cellcat(c1, c2, c3, ..., cn)
-%
-% Concatenates the concents of each cell array c1, ..., cn into one
-% big cell array.
-
-c = {};
-for n = 1:length(varargin)
- tmp = varargin{n};
- if ~iscell(tmp)
- tmp = {tmp};
- end
- c = {c{:} tmp{:}};
-end
35 celltools/cfilter.m
View
@@ -1,35 +0,0 @@
-function Y = cfilter(fun, C)
-% Y = cfilter(fun, C)
-%
-% Returns a cell array only containing elements of C for which fun(C{i}) == 1.
-%
-% 2007-11-06 ronw@ee.columbia.edu
-
-% Copyright (C) 2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if ~isa(fun, 'function_handle')
- error('1st argument must be a function handle')
-end
-
-if ~iscell(C)
- error('2nd argument must be a cell array.')
-end
-
-bool = zeros(size(C));
-for n = 1:numel(C)
- bool(n) = feval(fun, C{n});
-end
-Y = C(bool);
25 celltools/crange.m
View
@@ -1,25 +0,0 @@
-function c = crange(arg1, arg2, arg3)
-% c = crange(last)
-% c = crange(first, last)
-% c = crange(first, step, last)
-%
-% Creates a cell array containing a range of numbers. Analagous to
-% c = first:step:last
-% but creates a cell array instead of a normal matlab matrix. If
-% either first or step are not specified, they each default to 1.
-
-if nargin == 1
- first = 1;
- step = 1;
- last = arg1;
-elseif nargin == 2
- first = arg1;
- step = 1;
- last = arg2;
-elseif nargin == 3
- first = arg1;
- step = arg2;
- last = arg3;
-end
-
-c = num2cell(first:step:last);
31 celltools/czip.m
View
@@ -1,31 +0,0 @@
-function C = czip(varargin)
-% C = czip(cell1, cell2, cell3, ...)
-%
-% Like Python's zip function. Returns a cell array C containing a list of
-% cell arrays such that:
-% C{i} = {cell1{i}, cell2{i}, cell3{i}, ...}
-%
-% C is truncated to the length of the smallest of the input arrays.
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-len = mapreduce(@length, @(x,y) min(x,y), varargin);
-C = cell(len, 1);
-for l = 1:len
- for n = 1:nargin
- C{l}{n} = varargin{n}{l};
- end
-end
39 celltools/flatten.m
View
@@ -1,39 +0,0 @@
-function y = flatten(x)
-% y = flatten(x)
-%
-% Takes a cell array x that contains nested cell arrays and
-% flattens the contents into a single cell array.
-% E.g. flatten({1, {2, 3, {4}, 5}}) returns {1, 2, 3, 4, 5}
-
-% Copyright (C) 2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if ~iscell(x)
- error('flatten only works on cell arrays.');
-end
-
-y = inner_flatten(x);
-
-
-function y = inner_flatten(x)
-if ~iscell(x)
- y = {x};
-else
- y = {};
- for n = 1:length(x)
- tmp = inner_flatten(x{n});
- y = {y{:} tmp{:}};
- end
-end
61 celltools/glob.m
View
@@ -1,61 +0,0 @@
-function files = glob(pattern)
-% files = glob(pattern)
-%
-% Returns a cell array containing a list of files that match the given
-% pattern.
-%
-% 2008-04-30 ronw@ee.columbia.edu
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-files = glob_with_backoff(pattern);
-%files = dirglob(pattern);
-
-
-function files = dirglob(pattern)
-% This function doesn't produce the same output as lsglob:
-% 1. If pattern is 'some/path/' lsglob will return a single result while
-% this function will return all files in 'some/path/'.
-% I.e. dirglob('some/path/') is the same as lsglob('some/path/*')
-% 2. This doesn't work with nested globs (e.g. 'some/path/*/*.mat')
-path = fileparts(pattern);
-tmp = dir(pattern);
-files = cell(length(tmp), 1);
-for n = 1:length(tmp)
- files{n} = fullfile(path, tmp(n).name);
-end
-
-
-function files = lsglob(pattern)
-files = {};
-list_str = ls('-d1', pattern);
-idx = regexp(list_str, '\n');
-files = cell(length(idx), 1);
-last_idx = 1;
-for i = 1:length(idx)
- files{i} = list_str(last_idx:idx(i)-1);
- last_idx = idx(i) + 1;
-end
-
-
-function files = glob_with_backoff(pattern)
-try
- files = lsglob(pattern);
-catch
- err = lasterror();
- warning(sprintf('Backing off to dir because ls gave a %s error. Note that this has slightly different behavior so you might not get what you want.', err.identifier))
- files = dirglob(pattern);
-end
38 celltools/interleave.m
View
@@ -1,38 +0,0 @@
-function C = interleave(varargin)
-% C = interleave(cell1, cell2, cell3, ...)
-%
-% Interleaves the elements of cell1, cell2, cell3, ... in a lazy manner (i.e.
-% without making any copies of the elements of celli). This is similar to
-% Python's zip function (and roughly equivalent to
-% flatten(czip(cell1, cell2, cell3, ...)), except with lazy evaluation).
-%
-% Note that cell1, cell2, ... must be the same length.
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-lens = cell2mat(map(@length, varargin));
-if any(lens ~= lens(1))
- error('Inputs must all have the same length.');
-end
-
-len = sum(lens);
-C = lazymap(@(x) interleaved_index(varargin, x), crange(len));
-
-function y = interleaved_index(cell_arrays, x)
-narrays = length(cell_arrays);
-array = mod(x-1, narrays) + 1;
-offs = floor((x-1) / narrays) + 1;
-y = cell_arrays{array}{offs};
45 celltools/map.m
View
@@ -1,45 +0,0 @@
-function varargout = map(fun, varargin)
-% Y = map(fun, C)
-%
-% Wrapper around cellfun. Takes function handle fun and cell array
-% C and returns a new cell array that contains the result of
-% applying fun to each element of C.
-%
-% If fun takes N arguments, map must be passed N cell arrays
-% corresponding to those N arguments:
-% Y = map(fun, C1, C2, C3, ...);
-% Each input cell array must have the same size.
-%
-% Similarly, map can handle functions that have multiple outputs. If
-% called like this:
-% [Y1, Y2, Y3, ...] = map(fun, C)
-% Y1, Y2, ... will each be a cell array containing the analogous
-% outputs of fun.
-%
-% 2007-11-06 ronw@ee.columbia.edu
-
-% Copyright (C) 2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if ~isa(fun, 'function_handle')
- error('1st argument must be a function handle')
-end
-
-if ~iscell(varargin{1})
- error('2nd argument must be a cell array.')
-end
-
-
-[varargout{1:nargout}] = cellfun(fun, varargin{:}, 'UniformOutput', 0);
45 celltools/mapreduce.m
View
@@ -1,45 +0,0 @@
-function y = mapreduce(mapfun, reducefun, C, initial)
-% y = mapreduce(mapfun, reducefun, C, initial)
-%
-% Combines map and reduce into a single operation. First run mapfun
-% on each element of cell array C, then reduce the result to a single
-% element. This is essentially syntactic sugar for:
-% reduce(reducefun, map(mapfun, C), initial)
-%
-% 2007-11-13 ronw@ee.columbia.edu
-
-% Copyright (C) 2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if ~isa(reducefun, 'function_handle') || ~isa(mapfun, 'function_handle')
- error('1st two arguments must be function handles.')
-end
-
-if nargin(reducefun) > 0 && nargin(reducefun) ~= 2
- error('reduce function must take two arguments.');
-end
-
-if ~iscell(C)
- error('3rd argument must be a cell array.')
-end
-
-y = feval(mapfun, C{1});
-if nargin == 4
- y = reducefun(initial, y);
-end
-
-for n = 2:numel(C),
- y = feval(reducefun, y, feval(mapfun, C{n}));
-end
51 celltools/reduce.m
View
@@ -1,51 +0,0 @@
-function y = reduce(fun, C, initial)
-% y = reduce(fun, C, initial)
-%
-% Reduces the contents of cell array C to a single value by
-% accumulating the result of applying binary function fun to elements
-% from C.
-% e.g. reduce(@plus, {0, 1, 2, 3, 4}) will return 10.
-%
-% Optional argument initial can be used to set the initial value in
-% the accumulator, as if it was prepended to C.
-% e.g. reduce(@plus, {0, 1, 2, 3, 4}, 10) will return 20.
-%
-% 2007-11-06 ronw@ee.columbia.edu
-
-% Copyright (C) 2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if ~isa(fun, 'function_handle')
- error('1st argument must be a function handle')
-end
-
-if nargin(fun) > 0 && nargin(fun) ~= 2
- error('reduce function must take two arguments.');
-end
-
-if ~iscell(C)
- error('2nd argument must be a cell array.')
-end
-
-
-if nargin == 3
- y = feval(fun, initial, C{1});
-else
- y = C{1};
-end
-
-for n = 2:numel(C),
- y = feval(fun, y, C{n});
-end
46 celltools/sfilter.m
View
@@ -1,46 +0,0 @@
-function out_struct = sfilter(struct, varargin);
-% out_struct = sfilter(struct, 'field1', val1, 'field2', val2, ...)
-%
-% Filters the given struct array to only include those that match all
-% of the given field, value pairs.
-%
-% 2008-10-27 ronw@ee.columbia.edu
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-fields = varargin(1:2:end);
-values = varargin(2:2:end);
-nfields = length(fields);
-
-idx = 1;
-for x = 1:length(struct)
- include_this_struct = true;
- for y = 1:length(fields)
- tmp = getfield(struct(x), fields{y});
- if (isstr(tmp) && ~strcmp(tmp, values{y})) || any(tmp ~= values{y})
- include_this_struct = false;
- break;
- end
- end
- if include_this_struct
- out_struct(idx) = struct(x);
- idx = idx + 1;
- end
-end
-
-if ~exist('out_struct', 'var')
- out_struct = [];
-end
0  matlab_htk/compose_hmms.m → compose_hmms.m
View
File renamed without changes
0  matlab_htk/eval_htk_recognizer.m → eval_htk_recognizer.m
View
File renamed without changes
0  matlab_htk/get_htk_path.m → get_htk_path.m
View
File renamed without changes
0  matlab_htk/get_temp_filename.m → get_temp_filename.m
View
File renamed without changes
60 hmm/README
View
@@ -1,60 +0,0 @@
-This package contains a set of functions for evaluating HMMs and GMMs.
-
-* Functions
-
- - The important ones:
- - train_gmm - train a GMM from data
- - eval_gmm - compute the posterior probability of a GMM given data
- - eval_hmm - compute the posterior probabilities of all possible HMM
- state sequences given data
- - decode_hmm - find the most likely state sequence through the HMM
- given data
-
- - Utility functions:
- - logsum - takes the sum of a matrix of log likelihoods
- - lmvnpdf - compute the log probability of data under a
- multivariate Gaussian distribution
-
-* Data Structures
-
-The functions in this toolbox pass around the following structures:
-Note: all probabilities are stored as log probabilities
-
-** GMM
- - gmm.nmix - number of components in the mixture
- - gmm.priors - array of prior log probabilities over each state
- - gmm.means - matrix of means (column x is mean of component x)
- - gmm.covars - matrix of covariances (column x is the diagonal of the
- covariance matrix of component x)
-
-** HMM with GMM observations (does not work with all functions)
- - hmm.name -
- - hmm.nstates - number of states in the HMM
- - hmm.emission_type - 'GMM'
- - hmm.start_prob - array of log probs P(first observation is state x)
- - hmm.end_prob - array of log probs P(last observation is state x)
- - hmm.transmat - matrix of transition log probs (transmat(x,y)
- = log(P(transition from state x to state y)))
- - hmm.labels - optional cell array of labels for each state in the HMM
- (for use in composing HMMs)
- - hmm.gmms - array of GMM structures
-
-** HMM with Gaussian observations
- - hmm.nstates - number of states in the HMM
- - hmm.emission_type - 'gaussian'
- - hmm.start_prob - array of log probs P(first observation is state x)
- - hmm.end_prob - array of log probs P(last observation is state x)
- - hmm.transmat - matrix of transition log probs (transmat(x,y)
- = log(P(transition from state x to state y)))
- - hmm.labels - optional cell array of labels for each state in the HMM
- (for use in composing HMMs)
- - hmm.means - matrix of means (column x is mean of state x)
- - hmm.covars - matrix of means (column x is the diagonal of the
- covariance matrix of component x)
-
-Note that each row of exp(hmm.transmat) does not necessarily sum to 1
-because for each state x there is some probability
-(exp(hmm.end_prob(x))) that the next transition will be to a
-non-emitting exit state (i.e. the current observation is the last
-observation in the sequence). The correct invariant is:
-sum(exp(hmm.transmat), 2)' + exp(hmm.end_prob) == ones(hmm.nstates, 1)
129 hmm/adapt_gmm.m
View
@@ -1,129 +0,0 @@
-function gmm = adapt_gmm(pgmm, trdata, niter, adapt_weight, verb, MINCV)
-% gmmparams = adapt_gmm(prior_gmm, trdata, nmix, niter,
-% adapt_weight, verb, min_cv)
-%
-% Adapt the given prior_gmm to the given training data using MAP
-% adaptation.
-%
-% Inputs:
-% prior_gmm - initial GMM parameters to adapt.
-% trdata - training data (cell array of training sequences, each
-% column of the sequences arrays contains an observation)
-% niter - number of EM iterations to perform. Defaults to 10.
-% adapt_weight - weight of initial GMM in update equations (see HTK
-% book section 9.3). If adapt_weight is a vector, the
-% first element corresponds to the adaptation weight
-% for the Gaussian means (tau), and the second
-% correponse to the weight for the covariances (alpha).
-% Defaults to 10 (i.e. do not update covar parameters).
-% verb - set to 1 to output loglik at each iteration
-% min_cv - minimum covariance to avoid overfitting. Defaults to 1.
-%
-% Outputs:
-% gmmparams - structure containing hmm parameters learned from training
-% data (gmm.priors, gmm.means(:,1:nmix), gmm.covars(:,1:nmix))
-%
-% 2009-06-15 ronw@ee.columbia.edu
-
-DEBUG = false;
-
-if nargin < 3
- niter = 10;
-end
-if nargin < 4
- adapt_weight = 10;
-end
-if nargin < 5
- verb = 0;
-end
-if nargin < 6
- MINCV = 1;
-end
-
-if ~iscell(trdata)
- trdata = {trdata};
-end
-
-ndata = length(trdata);
-
-T = adapt_weight(1);
-adapt_covars = false;
-if length(adapt_weight) == 2
- adapt_covars = true;
- A = adapt_weight(2);
-end
-
-% Initialization
-gmm = pgmm;
-ppriors = exp(pgmm.priors);
-nmix = pgmm.nmix;
-
-if size(trdata{1}, 1) ~= size(gmm.means, 1)
- error(['Dimensionality of initial GMM not compatible with given ' ...
- 'training data']);
-end
-ndim = size(trdata{1}, 1);
-
-% sufficient statistics
-norm = zeros(size(gmm.priors));
-means = zeros(size(gmm.means));
-covars = zeros(size(gmm.covars));
-
-last_loglik = -Inf;
-for iter = 1:niter
- % E-step
- loglik = 0;
- norm(:) = 0;
- means(:) = 0;
- covars(:) = 0;
- for n = 1:ndata
- curr_data = trdata{n};
- [ll, posteriors] = eval_gmm(gmm, curr_data);
-
- loglik = loglik + sum(ll);
-
- norm = norm + sum(posteriors, 2)';
- means = means + curr_data * posteriors';
- covars = covars + curr_data.^2 * posteriors';
- end
-
- if verb,
- fprintf('Iteration %d: log likelihood = %f\n', iter, loglik)
-
- if DEBUG
- figure(1)
- plot_on_same_axes(pgmm.priors, gmm.priors)
- figure(2)
- cax = [-80 -10];
- subplot(211); imagesc(pgmm.means); axis xy; colorbar; caxis(cax);
- subplot(212); imagesc(gmm.means); axis xy; colorbar; caxis(cax);
- figure(3)
- cax = [0 100];
- subplot(211); imagesc(pgmm.covars); axis xy; colorbar; caxis(cax);
- subplot(212); imagesc(gmm.covars); axis xy; colorbar; caxis(cax);
- drawnow
- end
- end
-
- % Check for convergence
- if abs(loglik - last_loglik) < 1e-5
- fprintf('Converged at iteration %d\n', iter)
- break
- end
- last_loglik = loglik;
-
- % M-step
- % based on Huang, Acero, Hon, "Spoken Language Processing", p. 443 - 445
- npriors = (ppriors - 1 + norm) ./ sum(ppriors - 1 + norm);
- npriors(npriors < 0) = 0;
- gmm.priors = log(npriors);
- nrm = repmat(norm, [ndim 1]);
- gmm.means = (T * pgmm.means + means) ./ (T + nrm);
- if adapt_covars
- gmm.covars = ((A - 1) * pgmm.covars ...
- + T * (gmm.means - pgmm.means).^2 ...
- + (covars - 2*gmm.means.*means) .* nrm + gmm.means.^2) ...
- ./ (A - 1 + nrm);
- gmm.covars(gmm.covars < MINCV) = MINCV;
- end
-end
67 hmm/convert_hmm_to_gaussian_emissions.m
View
@@ -1,67 +0,0 @@
-function new_hmms = convert_hmm_to_gaussian_emissions(hmms)
-% new_hmms = convert_hmm_to_gaussian_emissions(hmms)
-%
-% Convert HMMs with GMM emissions in hmms to new_hmms with gaussian
-% emissions. hmms can be an array of hmm structures.
-%
-% 2007-01-18 ronw@ee.columbia.edu
-
-% Copyright (C) 2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-for n = 1:length(hmms)
- if strcmp(hmms(n).emission_type, 'gaussian')
- new_hmms(n) = hmms(n);
- continue;
- end
-
- new_hmms(n).name = hmms(n).name;
- new_hmms(n).emission_type = 'gaussian';
-
- nmix = [hmms(n).gmms(:).nmix];
- state_offset = cumsum([0, nmix(1:end-1)]);
-
- new_hmms(n).nstates = sum(nmix);
- new_hmms(n).transmat = repmat(-Inf, [new_hmms(n).nstates, new_hmms(n).nstates]);
-
- curr_state = 0;
- for s = 1:hmms(n).nstates
- ns = state_offset(s) + [1:nmix(s)];
-
- new_hmms(n).start_prob(ns) = hmms(n).start_prob(s);
- new_hmms(n).end_prob(ns) = hmms(n).end_prob(s);
-
- if(isfield(hmms(n), 'labels'))
- new_hmms(n).labels(ns) = hmms(n).labels(s);
- end
-
- new_hmms(n).means(:,ns) = hmms(n).gmms(s).means;
- new_hmms(n).covars(:,ns) = hmms(n).gmms(s).covars;
-
- for ss = 1:hmms(n).nstates
- nss = state_offset(ss) + [1:nmix(ss)];
- new_hmms(n).transmat(ns, nss) = hmms(n).transmat(s,ss);
- end
- end
-
- priors = [hmms(n).gmms(:).priors];
- % make sure priors is a column vector
- priors = priors(:);
- for r = 1:new_hmms(n).nstates
- new_hmms(n).transmat(r,:) = new_hmms(n).transmat(r,:) + priors';
- end
-
- new_hmms(n).start_prob = new_hmms(n).start_prob + priors';
-end
218 hmm/decode_hmm.m
View
@@ -1,218 +0,0 @@
-function [loglik, stateseq, recon, lattice, tb, mllattice] = decode_hmm(hmm, frameLogLike, maxRank, beamLogProb, normalize_lattice, verb)
-% [loglik, stateseq, recon, lattice, tb] = decode_hmm(hmm, seq, rank, beam)
-%
-% Performs Viterbi decode of seq. Does rank and beam pruning.
-% Assumes all hmm params are logprobs.
-%
-% 2007-02-26 ronw@ee.columbia.edu
-
-% Copyright (C) 2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-zeroLogProb = -1e200;
-
-% no rank pruning by default
-if nargin < 3
- maxRank = 0;
-end
-
-% no beam pruning by default
-if nargin < 4
- beamLogProb = -Inf;
-end
-
-if nargin < 5
- normalize_lattice = 0;
-end
-
-if nargin < 6
- verb = 0;
-end
-
-[nstates, nobs] = size(frameLogLike);
-if nstates ~= hmm.nstates && nstates == size(hmm.means, 1)
- seq = frameLogLike;
- ndim = nstates;
- nstates = hmm.nstates;
- if strcmp(hmm.emission_type, 'gaussian')
- frameLogLike = lmvnpdf(seq, hmm.means, hmm.covars);
- elseif strcmp(hmm.emission_type, 'GMM')
- for s = 1:hmm.nstates
- frameLogLike(s,:) = eval_gmm(hmm.gmms(s), seq);
- end
- else
- error('Unknown HMM emission distribution.');
- end
-end
-
-
-% how big should our rank pruning histogram be?
-histSize = 1000;
-
-stateseq = zeros(1, nobs);
-tb = zeros(hmm.nstates, nobs);
-
-if nargout > 3
- lattice = repmat(zeroLogProb, [hmm.nstates, nobs]);
-end
-% FIXME - there is a bug here in how the first frame is handled -
-% an extra transition probability from the non-existant 0th frame
-% to the 1st frame is added into lattice...
-prevLatticeFrame = hmm.start_prob(:);
-%prevLatticeFrame = hmm.start_prob(:) + frameLogLike(:,1);
-tb = zeros(hmm.nstates, nobs);
-
-% fill in the lattice...
-avg_nactive = 0;
-prevFrameMaxLogProb = 0;
-for obs = 1:nobs
- if verb >= 2
- tic
- end
-
- % beam pruning
- threshLogProb = prevFrameMaxLogProb + beamLogProb;
-
- % rank pruning
- if maxRank > 0
- tmp = prevLatticeFrame(:);
- min_tmp = 2*min(tmp(tmp > zeroLogProb));
- tmp(tmp < zeroLogProb) = min_tmp;
-
- [hst cdf] = hist(tmp, histSize);
-
- % want to look at the high ranks of the last frame
- hst = hst(end:-1:1);
- cdf = cdf(end:-1:1);
-
- hst = cumsum(hst);
- idx = min(find(hst >= maxRank));
- rankThresh = cdf(idx);
-
- % only change the threshold if it is stricter than the beam
- % threshold
- threshLogProb = max(threshLogProb, rankThresh);
-
- if verb >= 3
- %imgsc(prevLatticeFrame), colorbar, title(num2str(obs)), drawnow
-
- disp(['beam thresh = ' num2str(prevFrameMaxLogProb+beamLogProb) ...
- ', rank thresh = ' num2str(rankThresh) ...
- ', final thresh = ' num2str(threshLogProb)]);
- end
- end
-
- % which states are active?
- s_idx = find(prevLatticeFrame >= threshLogProb);
- nactive = numel(s_idx);
- avg_nactive = avg_nactive + nactive/nobs;
-
- vitPr = hmm.transmat(s_idx, :)' + repmat(prevLatticeFrame(s_idx), [1, hmm.nstates])';
-
- currllik = frameLogLike(:,obs);
-
-% v_idx = find(max(vitPr,[], 2) > zeroLogProb);
-% nv = length(v_idx);
-% currllik = repmat(zeroLogProb, [hmm.nstates, 1]);
-% if strcmp(hmm.emission_type, 'gaussian')
-% currllik(v_idx) = lmvnpdf(seq(:,obs), hmm.means(:, v_idx), ...
-% hmm.covars(:, v_idx));
-% else
-% for s = 1:nv
-% currllik(v_idx(s)) = eval_gmm(hmm.gmms(v_idx(s)), seq(:,obs));
-% end
-% end
-
- [prevLatticeFrame tb_tmp] = max(vitPr + repmat(currllik, [1, nactive]), [], 2);
- % This is equivalent to the above?
- %[prevLatticeFrame tb_tmp] = max(vitPr, [], 2);
- %prevLatticeFrame = prevLatticeFrame + currllik;
- tb(:,obs) = s_idx(tb_tmp);
-
- if nargout > 3
- lattice(:,obs) = prevLatticeFrame;
- end
-
- prevFrameMaxLogProb = max(prevLatticeFrame);
-
- if verb >= 2
- T = toc;
- disp(['frame ' num2str(obs), ' (' num2str(T) ' sec)' ...
- ': total active states: ' num2str(nactive)]);
- end
-end
-
-% include end_prob in lattice:
-ptmp = prevLatticeFrame;
-prevLatticeFrame = prevLatticeFrame + hmm.end_prob(:);
-
-%%%
-% do the traceback:
-[loglik s] = max(prevLatticeFrame(:));
-
-% we might have pruned too much, don't want to give up if end_prob
-% restrictions are too strong - so just ignore them in this case
-if loglik <= zeroLogProb
- warning(['decode_hmm: overpruned during decode,' ...
- ' ignoring probabilities that the hmms end in a' ...
- ' particular state']);
- prevLatticeFrame = ptmp;
- [loglik s] = max(prevLatticeFrame);
-end
-
-if nargout > 3
- lattice(:,end) = prevLatticeFrame;
-end
-
-for obs = nobs:-1:1
- stateseq(obs) = s;
- s = tb(s, obs);
-end
-
-
-% need to keep track of which gmm component was used in the
-% traceback for reconstruction - since this code doesn't do that,
-% can only guess as to the right reconstruction for GMM emissions -
-% should use convert_hmm_to_gaussian_emissions before calling this
-% function to get the full traceback
-if strcmp(hmm.emission_type, 'gaussian')
- recon = hmm.means(:,stateseq);
-else
- for s = 1:hmm.nstates
- means(:,s) = hmm.gmms(s).means * exp(hmm.gmms(s).priors)';
- end
-
- recon = means(:,stateseq);
-end
-
-if nargout > 3 & normalize_lattice
- nrm = logsum(lattice, 1);
- lattice = exp(lattice - repmat(nrm, [hmm.nstates, 1]));
- lattice(lattice < 1e-5) = 0;
-end
-
-if nargout > 4
- mllattice = sparse(hmm.nstates, nobs);
- for o = 1:nobs
- mllattice(stateseq(o),o) = 1;
- end
-end
-
-
-if verb
- disp(['decode_hmm: log likelihood: ' num2str(loglik) ...
- ', average number of active states per frame: ' ...
- num2str(avg_nactive)]);
-end
54 hmm/eval_gmm.m
View
@@ -1,54 +0,0 @@
-function [ll, post, mlg, mmserecon] = eval_gmm(gmm, data, norm)
-% [loglik, posteriors, mlseq, recon] = eval_gmm(gmm, data)
-%
-% Evaluate the log probability of each column of data given GMM gmm.
-%
-% Outputs:
-% loglik - log likelihood of each colimn of data
-% posteriors - posterior probability of each GMM component for each
-% column of data
-% mlseq - index of the most likely GMM component for each
-% column of data
-% recon - MMSE reconstruction of data given the GMM
-%
-% 2005-11-20 ronw@ee.columbia.edu
-
-% Copyright (C) 2005-2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if nargin < 3, norm = 1; end
-
-[ndim, ndat] = size(data);
-
-post = lmvnpdf(data, gmm.means, gmm.covars) + repmat(gmm.priors(:), [1, ndat]);
-ll = logsum(post, 1);
-
-if nargout > 1 && norm
- post = exp(post - repmat(logsum(post,1), gmm.nmix, 1));
-end
-
-if nargout > 2
- [mlg tmp] = ind2sub(size(post), find(post == repmat(max(post),gmm.nmix,1)));
-end
-
-if nargout > 3
- if norm
- mmserecon = gmm.means*post;
- else
- postnorm = exp(post - repmat(logsum(post,1), gmm.nmix, 1));
- mmserecon = gmm.means*postnorm;
- end
-end
-
195 hmm/eval_hmm.m
View
@@ -1,195 +0,0 @@
-function [loglik, lattice, alpha, beta, gamma] = eval_hmm(hmm, frameLogLike, maxRank, beamLogProb, do_backward, verb)
-% [loglik, lattice] = eval_hmm(hmm, seq, rank, beam)
-%
-% Performs forward-backward inference on seq. Does rank and beam
-% pruning. Assumes all hmm params are logprobs.
-%
-% 2008-08-11 ronw@ee.columbia.edu
-
-% Copyright (C) 2006-2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-
-% no rank pruning by default
-if nargin < 3
- maxRank = 0;
-end
-
-% No beam pruning by default.
-if nargin < 4
- beamLogProb = -Inf;
-end
-
-if nargin < 5
- do_backward = true;
-end
-
-if nargin < 6
- verb = 0;
-end
-
-% Don't bother doing backward calculation if all we want is log
-% likelihood.
-if nargout < 2
- do_backward = false;
-else
- do_backward = true;
-end
-
-zeroLogProb = -1e200;
-hmm.transmat(hmm.transmat < zeroLogProb) = zeroLogProb;
-
-% Verify type of observations. Can be observed sequence or
-% precomputed log likelihoods (i.e. for variational inference).
-[nstates, nobs] = size(frameLogLike);
-if nstates ~= hmm.nstates && nstates == size(hmm.means, 1)
- seq = frameLogLike;
- ndim = nstates;
- nstates = hmm.nstates;
- if strcmp(hmm.emission_type, 'gaussian')
- frameLogLike = lmvnpdf(seq, hmm.means, hmm.covars);
- elseif strcmp(hmm.emission_type, 'GMM')
- for s = 1:hmm.nstates
- frameLogLike(s,:) = eval_gmm(hmm.gmms(s), seq);
- end
- else
- error('Unknown HMM emission distribution.');
- end
-end
-
-
-%%%%%
-% Forward
-%%%%%
-alpha = zeros(nstates, nobs) - Inf;
-prevLatticeFrame = hmm.start_prob(:) + frameLogLike(:,1);
-alpha(:,1) = prevLatticeFrame;
-if verb >= 2
- fprintf('Starting forward pass...\n frame 1: ll = %f\n', ...
- logsum(prevLatticeFrame))
-end
-
-for obs = 2:nobs
- if verb >= 2; tic; end
-
- idx = prune_states(prevLatticeFrame, maxRank, beamLogProb, verb);
- pr = hmm.transmat(idx,:)' + repmat(prevLatticeFrame(idx), [1, hmm.nstates])';
- prevLatticeFrame = logsum(pr, 2) + frameLogLike(:, obs);
- alpha(:,obs) = prevLatticeFrame;
-
- if verb >= 2
- T = toc;
- fprintf(' frame %d: ll = %f (%f sec, %d active states)\n', obs, ...
- logsum(prevLatticeFrame), T, length(idx));
- end
-end
-alpha(alpha <= zeroLogProb) = -Inf;
-
-% Don't forget hmm.end_prob
-% This double counts frameLogLike(:,end)!!
-%nextLatticeFrame = hmm.end_prob(:) + frameLogLike(:,end);
-nextLatticeFrame = hmm.end_prob(:);
-loglik = logsum(prevLatticeFrame + nextLatticeFrame);
-if isinf(loglik) || isnan(loglik)
- nextLatticeFrame = frameLogLike(:,end);
- loglik = logsum(prevLatticeFrame + nextLatticeFrame);
-end
-
-if verb
- fprintf('eval_hmm: log likelihood = %f\n', loglik)
-end
-
-
-if ~do_backward
- return
-end
-
-%%%%%
-% Backward
-%%%%%
-beta = zeros(nstates, nobs) - Inf;
-beta(:,nobs) = nextLatticeFrame;
-if verb >= 2
- fprintf('Starting backward pass...\n frame %d: ll = %f\n', nobs, ...
- logsum(nextLatticeFrame));
-end
-
-for obs = nobs-1:-1:1
- if verb >= 2; tic; end
-
- % Do HTK style pruning (p. 137 of HTK Book version 3.4). Don't
- % bother computing backward probability if alpha*beta is more than a
- % certain distance from the total log likelihood.
- idx = prune_states(nextLatticeFrame + alpha(:,obs+1), 0, -20, verb);
- %idx = prune_states(nextLatticeFrame + alpha(:,obs+1), 10, -Inf, verb);
-
- pr = hmm.transmat(:,idx) + repmat(nextLatticeFrame(idx) ...
- + frameLogLike(idx,obs+1), [1, hmm.nstates])';
- nextLatticeFrame = logsum(pr, 2);
- beta(:,obs) = nextLatticeFrame;
-
- if verb >= 2
- T = toc;
- fprintf(' frame %d: ll = %f (%f sec, %d active states)\n', obs, ...
- logsum(nextLatticeFrame), T, length(idx));
- end
-end
-beta(beta <= zeroLogProb) = -Inf;
-
-gamma = alpha + beta;
-lattice = exp(gamma - repmat(logsum(gamma, 1), [hmm.nstates 1]));
-
-
-
-function [state_idx thresh] = prune_states(latticeFrame, ...
- maxRank, beamLogProb, verb)
-zeroLogProb = -1e200;
-frameLogProb = logsum(latticeFrame);
-
-% Beam pruning
-threshLogProb = frameLogProb + beamLogProb;
-
-% Rank pruning
-if maxRank > 0
- % How big should our rank pruning histogram be?
- histSize = 3*length(latticeFrame);
-
- tmp = latticeFrame(:);
- min_tmp = min(tmp(tmp > zeroLogProb)) - 1;
- tmp(tmp <= zeroLogProb) = min_tmp;
-
- [hst cdf] = hist(tmp, histSize);
-
- % Want to look at the high ranks of the last frame.
- hst = hst(end:-1:1);
- cdf = cdf(end:-1:1);
-
- hst = cumsum(hst);
- idx = min(find(hst >= maxRank));
- rankThresh = cdf(idx);
-
- % Only change the threshold if it is stricter than the beam
- % threshold.
- threshLogProb = max(threshLogProb, rankThresh);
-
- if verb >= 3
- fprintf('beam thresh = %f, rank thresh = %f, final thresh = %f\n', ...
- frameLogProb+beamLogProb, rankThresh, threshLogProb);
- end
-end
-
-% Which states are active?
-state_idx = find(latticeFrame >= threshLogProb);
-
82 hmm/is_valid_gmm.m
View
@@ -1,82 +0,0 @@
-function out = is_valid_gmm(gmm, verb)
-% y = is_valid_gmm(gmm)
-%
-% Returns 1 if and only if gmm is a valid GMM structure
-%
-% 2008-06-03 ronw@ee.columbia.edu
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if nargin < 2; verb = 0; end
-
-out = true;
-if ~isstruct(gmm)
- out = false;
- if verb; fprintf('is_valid_gmm: not a structure\n'); end
- return
-end
-
-if ~isfield(gmm, 'nmix')
- out = false;
- if verb; fprintf('is_valid_gmm: missing nmix field\n'); end
- return
-end
-
-if ~isfield(gmm, 'priors')
- out = false;
- if verb; fprintf('is_valid_gmm: missing priors field\n'); end
-else
- if length(gmm.priors) ~= gmm.nmix
- out = false;
- if verb; fprintf('is_valid_gmm: priors field is wrong length\n'); end
- end
- if abs(logsum(gmm.priors)) > 1e-3
- out = false;
- if verb; fprintf('is_valid_gmm: priors don''t sum to 1\n'); end
- end
-end
-
-if ~isfield(gmm, 'means')
- out = false;
- if verb; fprintf('is_valid_gmm: missing means field\n'); end
-else
- if size(gmm.means, 2) ~= gmm.nmix
- out = false;
- if verb; fprintf('is_valid_gmm: means field is wrong length\n'); end
- end
-end
-
-if ~isfield(gmm, 'covars')
- out = false;
- if verb; fprintf('is_valid_gmm: missing covars field\n'); end
-else
- if size(gmm.covars, 2) ~= gmm.nmix
- out = false;
- if verb; fprintf('is_valid_gmm: covars field is wrong length\n'); end
- end
- if isfield(gmm, 'means') && size(gmm.means, 1) ~= size(gmm.covars,1)
- out = false;
- if verb
- fprintf('is_valid_gmm: means and covars have inconsistent dimensions\n');
- end
- end
- if ~all(all(gmm.covars > 0))
- out = false;
- if verb; fprintf('is_valid_gmm: 0 or negative covars\n'); end
- end
-end
-
-
133 hmm/is_valid_hmm.m
View
@@ -1,133 +0,0 @@
-function out = is_valid_hmm(hmm, verb)
-% y = is_valid_hmm(hmm)
-%
-% Returns 1 if and only if hmm is a valid HMM structure
-%
-% 2008-06-03 ronw@ee.columbia.edu
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if nargin < 2; verb = 0; end
-
-out = true;
-if ~isstruct(hmm)
- out = false;
- if verb; fprintf('is_valid_hmm: not a structure\n'); end
- return
-end
-
-if ~isfield(hmm, 'nstates')
- out = false;
- if verb; fprintf('is_valid_hmm: missing nstates field\n'); end
- return
-end
-if ~isfield(hmm, 'emission_type')
- out = false;
- if verb; fprintf('is_valid_hmm: missing emissionn type field\n'); end
-end
-
-
-if ~isfield(hmm, 'start_prob')
- out = false;
- if verb; fprintf('is_valid_hmm: missing start_prob field\n'); end
-else
- if length(hmm.start_prob) ~= hmm.nstates
- out = false;
- if verb; fprintf('is_valid_hmm: start_prob field is wrong length\n'); end
- end
- if abs(logsum(hmm.start_prob)) > 1e-3
- out = false;
- if verb; fprintf('is_valid_hmm: start_prob doesn''t sum to 1\n'); end
- end
-end
-
-if ~isfield(hmm, 'end_prob')
- out = false;
- if verb; fprintf('is_valid_hmm: missing end_prob field\n'); end
-else
- if length(hmm.end_prob) ~= hmm.nstates
- out = false;
- if verb; fprintf('is_valid_hmm: end_prob field is wrong length\n'); end
- end
- if ~isfield(hmm, 'transmat')
- out = false;
- if verb; fprintf('is_valid_hmm: missing transmat field\n'); end
- else
- if ~all(size(hmm.transmat) == [hmm.nstates hmm.nstates])
- out = false;
- if verb; fprintf('is_valid_hmm: transmat is wrong size\n'); end
- end
- if ~all(abs(sum(exp(hmm.transmat), 2) + exp(hmm.end_prob(:)) - 1) < 1e-3)
- out = false;
- if verb;
- fprintf('is_valid_hmm: transmat is not normalized properly\n');
- end
- end
- end
-end
-
-
-if strcmp(hmm.emission_type, 'GMM')
- if ~isfield(hmm, 'gmms')
- out = false;
- if verb; fprintf('is_valid_hmm: missing gmms field\n'); end
- else
- for s = 1:length(hmm.gmms)
- tmp = is_valid_gmm(hmm.gmms(s));
- if ~tmp
- out = false;
- if verb
- fprintf('is_valid_hmm: Error in state %d:\n ', s);
- is_valid_gmm(hmm.gmms(s), verb);
- end
- break
- end
- end
- end
-end
-if strcmp(hmm.emission_type, 'gaussian')
- if ~isfield(hmm, 'means')
- out = false;
- if verb; fprintf('is_valid_hmm: missing means field\n'); end
- else
- if size(hmm.means, 2) ~= hmm.nstates
- out = false;
- if verb; fprintf('is_valid_hmm: means field is wrong length\n'); end
- end
- end
-
- if ~isfield(hmm, 'covars')
- out = false;
- if verb; fprintf('is_valid_hmm: missing covars field\n'); end
- else
- if size(hmm.covars, 2) ~= hmm.nstates
- out = false;
- if verb; fprintf('is_valid_hmm: covars field is wrong length\n'); end
- end
- if isfield(hmm, 'means') && size(hmm.means, 1) ~= size(hmm.covars,1)
- out = false;
- if verb
- fprintf(['is_valid_hmm: means and covars have inconsistent ' ...
- 'dimensions\n']);
- end
- end
- if ~all(all(hmm.covars > 0))
- out = false;
- if verb; fprintf('is_valid_hmm: 0 or negative covars\n'); end
- end
- end
-end
-
67 hmm/lmvnpdf.m
View
@@ -1,67 +0,0 @@
-function lpr = lmvnpdf(obs, mu, cv);
-% lpr = lmvnpdf(obs, mu, cv)
-%
-% Return the log probability of obs under the Gaussian distribution
-% parameterized by mu and cv.
-%
-% obs is an array of column vectors (DxO). mu and cv are also arrays
-% of column vectors (this only supports diagonal covariance matrices,
-% so mu and cv must both be DxC where C is the number of Gaussians).
-% lpr will be a CxO matrix where each row contains the log probability
-% of each observation given one of the C Gaussians.
-%
-% 2006-06-19 ronw@ee.columbia.edu
-
-% Copyright (C) 2006-2007 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if nargin < 3
- cv = 1;
-end
-
-[ndim, nobs] = size(obs);
-[ndim_mu, nmu] = size(mu);
-[ndim_cv, ncv] = size(cv);
-
-% make sure all the arguments are consistent
-if ndim ~= ndim_mu
- error('lmvnpdf: obs and mu must have the same number of dimensions.');
-end
-if nmu ~= ncv
- if ncv == 1
- % use the same diagonal covariance for each distribution
- cv = repmat(cv, 1, nmu);
- else
- error('lmvnpdf: mu and cv must have the same number of components.');
- end
-end
-ngauss = nmu;
-
-% are covariances scalar?
-if ndim_cv == 1
- cv = repmat(cv, ndim, 1);
-end
-
-
-% vectorized like there is no tomorrow:
-% ||x-y|| = x'x - 2*x'y + y'y
-% x'x = repmat(sum(x.^2),xc,1);
-% y'y = repmat(sum(y.^2),yc,1);
-%
-% but here, its ||(x-y)/cv||:
-% where cv has the same size as x (mu), but not the same as y (obs)...
-
-lpr = -0.5*(repmat(sum((mu.^2)./cv, 1)' + sum(log(cv))', [1 nobs]) ...
- - 2*(mu./cv)'*obs + (1./cv)'*(obs.^2) + ndim*log(2*pi));
15 hmm/logsum.m
View
@@ -1,15 +0,0 @@
-function s = logsum(array, dim)
-% s = logsum(array, dim)
-%
-% returns log(sum(exp(array), dim)) minimizing possibility of over/underflow
-
-
-if nargin < 2
- amax = max(array(:));
- s = log(sum(exp(array - amax))) + amax;
-else
- amax = max(array,[],dim);
- rep = ones(1, length(size(array)));
- rep(dim) = size(array,dim);
- s = log(sum(exp(array - repmat(amax, rep)), dim)) + amax;
-end
113 hmm/merge_states.m
View
@@ -1,113 +0,0 @@
-function new_hmm = merge_states(hmm, idx)
-% new_hmm = merge_states(hmm, idx)
-%
-% Merge states of the given HMM. idx is a list of state indices to
-% merge. If idx is a matrix, the indices in each column will be
-% merged into a single state. Also works on GMM structures.
-%
-% Examples:
-% - Merge states 1, 3, and 5: merge_states(hmm, [1 3 5])
-% - Merge states 1:5 and 10:20: merge_states(hmm, [1:5; 10:20]')
-% - Merge succesive pairs of states:
-% merge_states(hmm, reshape(1:hmm.nstates, [2, hmm.nstates/2]))
-%
-% Note that the merging isn't at all correct - it just dumbly takes
-% the weighted average of the given states.
-%
-% 2008-06-03 ronw@ee.columbia.edu
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-[nr nc] = size(idx);
-if nr == 1 || nc == 1
- % Turn into a column vector
- idx = idx(:);
-end
-
-if is_valid_gmm(hmm)
- new_hmm = merge_states_gmm(hmm, idx);
-else
- new_hmm = merge_states_hmm(hmm, idx);
-end
-
-
-
-function new_gmm = merge_states_gmm(gmm, idx)
-[nr nc] = size(idx);
-new_gmm = gmm;
-states_to_delete = zeros(gmm.nmix, 1);
-for c = 1:nc
- i = idx(:,c);
- states_to_delete(i(2:end)) = 1;
-
- lp = gmm.priors(i);
- new_gmm.priors(i(1)) = logsum(lp);
-
- p = exp(lp(:) - logsum(lp));
- new_gmm.means(:,i(1)) = gmm.means(:,i) * p;
- new_gmm.covars(:,i(1)) = gmm.covars(:,i) * p;
-end
-
-i = find(~states_to_delete);
-new_gmm.nmix = length(i);
-new_gmm.priors = new_gmm.priors(i);
-new_gmm.means = new_gmm.means(:,i);
-new_gmm.covars = new_gmm.covars(:,i);
-
-
-
-function new_hmm = merge_states_hmm(hmm, idx)
-if strcmp(hmm.emission_type, 'GMM')
- error('HMMs with GMM emissions are not supported.');
-end
-[nr nc] = size(idx);
-new_hmm = hmm;
-states_to_delete = zeros(hmm.nstates, 1);
-for c = 1:nc
- i = idx(:,c);
- states_to_delete(i(2:end)) = 1;
-
- lp = logsum(hmm.transmat(:,i), 1);
- p = exp(lp(:) - logsum(lp));
-
- new_hmm.start_prob(i(1)) = logsum(hmm.start_prob(i));
- new_hmm.transmat(i(1),i(1)) = logsum(logsum(hmm.transmat(i,i), 2) + p);
- for s = 1:hmm.nstates
- new_hmm.transmat(s,i(1)) = logsum(hmm.transmat(s,i));
- end
- new_hmm.end_prob(i(1)) = logsum(hmm.end_prob(i));
-
- new_hmm.means(:,i(1)) = hmm.means(:,i) * p;
- new_hmm.covars(:,i(1)) = hmm.covars(:,i) * p;
-end
-
-i = find(~states_to_delete);
-new_hmm.nstates = length(i);
-new_hmm.start_prob = new_hmm.start_prob(i);
-new_hmm.transmat = new_hmm.transmat(i,i);
-new_hmm.end_prob = new_hmm.end_prob(i);
-new_hmm.means = new_hmm.means(:,i);
-new_hmm.covars = new_hmm.covars(:,i);
-
-% Get rid of NaNs introduced by logsum(-Inf)
-new_hmm.start_prob(isnan(new_hmm.start_prob)) = -Inf;
-new_hmm.end_prob(isnan(new_hmm.end_prob)) = -Inf;
-new_hmm.transmat(isnan(new_hmm.transmat)) = -Inf;
-
-% make sure transmat and end_prob are normalized properly
-norm = logsum(cat(2, logsum(new_hmm.transmat, 2), new_hmm.end_prob'), 2);
-new_hmm.transmat = new_hmm.transmat - repmat(norm, [1 new_hmm.nstates]);
-new_hmm.end_prob = new_hmm.end_prob - norm';
56 hmm/reorder_states.m
View
@@ -1,56 +0,0 @@
-function new_hmm = reorder_states(hmm, idx)
-% new_hmm = reorder_states(hmm, idx)
-%
-% Rearrange the states of hmm using the given indices.
-%
-% 2008-09-12 ronw@ee.columbia.edu
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if is_valid_gmm(hmm)
- new_hmm = reorder_states_gmm(hmm, idx);
-else
- new_hmm = reorder_states_hmm(hmm, idx);
-end
-
-function new_gmm = reorder_states_gmm(gmm, idx)
-new_gmm = gmm;
-new_gmm.nmix = length(idx);
-new_gmm.means = gmm.means(:,idx);
-new_gmm.covars = gmm.covars(:,idx);
-% Ensure priors are normalized (if components are deleted and not
-% just rearranged).
-new_gmm.priors = gmm.priors(idx) - logsum(gmm.priors(idx));
-
-
-function new_hmm = reorder_states_hmm(hmm, idx)
-new_hmm = hmm;
-new_hmm.nstates = length(idx);
-if strcmp(hmm.emission_type, 'GMM')
- new_hmm.gmms = hmm.gmms(idx);
-else
- new_hmm.means = hmm.means(:,idx);
- new_hmm.covars = hmm.covars(:,idx);
-end
-new_hmm.start_prob = hmm.start_prob(idx);
-new_hmm.transmat = hmm.transmat(idx,idx);
-new_hmm.end_prob = hmm.end_prob(idx);
-
-% Make sure everything is normalized properly.
-new_hmm.start_prob = new_hmm.start_prob - logsum(new_hmm.start_prob);
-norm = logsum(cat(2, logsum(new_hmm.transmat, 2), new_hmm.end_prob'), 2);
-new_hmm.transmat = new_hmm.transmat - repmat(norm, [1 new_hmm.nstates]);
-new_hmm.end_prob = new_hmm.end_prob - norm';
24 hmm/sample_gaussian.m
View
@@ -1,24 +0,0 @@
-function y = sample_gaussian(mu, cv)
-% sample = sample_gaussian(mu, cv)
-%
-% Generate a random sample from a Gaussian distribution with the given
-% parameters.
-%
-% 2008-06-04 ronw@ee.columbia.edu
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-y = randn(size(mu)).*sqrt(cv) + mu;
34 hmm/sample_gmm.m
View
@@ -1,34 +0,0 @@
-function samples = sample_gmm(gmm, nsamp)
-% samples = sample_gmm(gmm, N)
-%
-% Generate N random samples from the given GMM.
-%
-% 2008-06-04 ronw@ee.columbia.edu
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-prior_pdf = exp(gmm.priors);
-prior_cdf = cumsum(prior_pdf);
-
-ndim = size(gmm.means, 1);
-
-samples = zeros(ndim, nsamp);
-for n = 1:nsamp
- p = rand(1);
- c = min(find(prior_cdf >= p));
- samples(:,n) = sample_gaussian(gmm.means(:,c), gmm.covars(:,c));
-end
-
61 hmm/sample_hmm.m
View
@@ -1,61 +0,0 @@
-function [samples state_seq] = sample_hmm(hmm)
-% [samples state_seq] = sample_hmm(hmm)
-%
-% Generate a random sample from the given HMM.
-%
-% 2008-06-04 ronw@ee.columbia.edu
-
-% Copyright (C) 2008 Ron J. Weiss
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-sp_pdf = exp(hmm.start_prob);
-sp_cdf = cumsum(sp_pdf);
-trans_pdf = exp([hmm.transmat hmm.end_prob']);
-trans_cdf = cumsum(trans_pdf, 2);
-
-% Initial state
-p = rand(1);
-s = min(find(sp_cdf >= p));
-state_seq = s;
-samples = sample_from_state(hmm, s);
-
-i = 1;
-while true
- % Select next component or exit.
- p = rand(1);
- os = s;
- s = min(find(trans_cdf(os,:) >= p));
- if isempty(s)
- s = length(trans_cdf(os,:));
- end
-
- if s > hmm.nstates
- break
- else
- i = i + 1;
- samples(:,i) = sample_from_state(hmm, s);
- state_seq(i) = s;
- end
-end
-
-
-function y = sample_from_state(hmm, s)
-if strcmp(hmm.emission_type, 'GMM')
- y = sample_gmm(hmm.gmms(s), 1);
-elseif strcmp(hmm.emission_type, 'gaussian')
- y = sample_gaussian(hmm.means(:,s), hmm.covars(:,s));
-else
- error(['Invalid HMM emission type: ' hmm.emission_type]);
-end
101 hmm/train_gmm.m
View
@@ -1,101 +0,0 @@
-function gmm = train_gmm(trdata, nmix, niter, verb, CVPRIOR, mu0);
-% gmmparams = train_gmm(trdata, nmix, niter, verb, cvprior, mu0);
-%
-% Train a GMM with diagonal covariance.
-%
-% Inputs:
-% trdata - training data (cell array of training sequences, each
-% column of the sequences arrays contains an
-% observation)
-% nmix - number of mixture components. Defaults to 3.
-% niter - number of EM iterations to perform. Defaults to 10.
-% verb - set to 1 to output loglik at each iteration
-% cvprior -
-%
-% Outputs:
-% gmmparams - structure containing hmm parameters learned from training
-% data (gmm.priors, gmm.means(:,1:nmix), gmm.covars(:,1:nmix))
-%
-% 2007-11-06 ronw@ee.columbia.edu
-
-if nargin < 2
- nmix = 3;
-end
-if nargin < 3
- niter = 10;
-end
-
-if nargin < 4
- verb = 0;
-end
-
-% prior on observation covariances to avoid overfitting:
-if nargin < 5
- CVPRIOR = 1;
-end
-
-if ~iscell(trdata)
- trdata = {trdata};
-end
-
-ndata = length(trdata);
-
-
-% Initialization
-gmm.priors = log(ones(1, nmix)/nmix);
-gmm.nmix = nmix;
-
-if nargin < 6 | numel(mu0) == 1 & mu0 == 1
- gmm.means = kmeans(cat(2, trdata{:}), nmix, niter/2);
-else
- if size(mu0, 2) == nmix
- gmm.means = mu0;
- end
-end
-
-ndim = size(trdata{1}, 1);
-%gmm.covars = ones(ndim, nmix);
-gmm.covars(:,1:nmix) = repmat(var(trdata{1}')', [1 nmix]);
-
-
-% sufficient statistics
-norm = zeros(size(gmm.priors));
-means = zeros(size(gmm.means));
-covars = zeros(size(gmm.covars));
-
-last_loglik = 0;
-for iter = 1:niter
- % E-step
- loglik = 0;
- norm(:) = 0;
- means(:) = 0;
- covars(:) = 0;
- for n = 1:ndata
- curr_data = trdata{n};
- [ll, posteriors] = eval_gmm(gmm, curr_data);
-
- loglik = loglik + sum(ll);
-
- norm = norm + sum(posteriors, 2)';
- means = means + curr_data * posteriors';
- covars = covars + curr_data.^2 * posteriors';
- end
-
- if verb,
- fprintf('Iteration %d: log likelihood = %f\n', iter, loglik);
- end
-
- % Check for convergence
- if abs(loglik - last_loglik) < 1e-5
- break
- end
- last_loglik = loglik;
-
- % M-step
- gmm.priors = log(norm/sum(norm));
-
- nrm = repmat(1./norm, [ndim 1]);
- gmm.means = means .* nrm;
- gmm.covars = (covars - 2*gmm.means.*means) .* nrm + gmm.means.^2;
- gmm.covars(gmm.covars < CVPRIOR) = CVPRIOR;
-end
0  matlab_htk/htkread.m → htkread.m
View
File renamed without changes
0  matlab_htk/htkwrite.m → htkwrite.m
View
File renamed without changes
0  matlab_htk/kmeans.m → kmeans.m
View
File renamed without changes
0  matlab_htk/learn_mllr_transform.m → learn_mllr_transform.m
View
File renamed without changes
0  matlab_htk/logsum.m → logsum.m
View
File renamed without changes
77 matlab_htk/README
View
@@ -1,77 +0,0 @@
-This package contains a set of functions for calling interfacing with
-HTK from Matlab. Right now its mostly limited to training GMMs and
-HMMs. It converts your Matlab data into a format that HTK understands
-and calls HTK command line programs. The path to the HTK binaries is
-hardcoded in get_htk_path.m
-
-Unforuntately this has not been very thoroughly tested, but it has
-been working for me. If anyone wants to write some unit tests using
-http://mlunit.sf.net or something similar I won't say no.
-
-
-* Functions
-
-They are all reasonably commented...
-
- - The important ones:
- - train_hmm_htk - use HTK to train an HMM
- - train_gmm_htk - use HTK to train a GMM
- - train_htk_recognizer - train a simple HTK recognizer
- - eval_htk_recognizer - recognize signal using HTK recognizer
-
- - Utility functions:
- - read_htk_hmm - read in an HTK HMM definition (text format only)
- - write_htk_hmm - write an HMM structure as an HTK HMM definition
- - htkread/htkwrite - read/write a Matlab matrix in HTK feature file format
- (written by Mark Hasegawa-Johnson)
- - compose_hmms - does FSM composition to form a big HMM from many
- small HMMs (i.e. get an HMM to recognize a
- sentence from a bunch of phone HMMs).
- - kmeans - uses k-means to learn clusters from data. Used to
- initialize GMMs and HMMs.
- - logsum - takes the sum of a matrix of log likelihoods
- - get_htk_path - centralized location to set the path to the HTK binaries
-
-* Data Structures
-
-The functions in this toolbox pass around the following structures:
-Note: all probabilities are stored as log probabilities
-
-** GMM
- - gmm.nmix - number of components in the mixture
- - gmm.priors - array of prior log probabilities over each state
- - gmm.means - matrix of means (column x is mean of component x)
- - gmm.covars - matrix of covariance (column x is the diagonal of the
- covariance matrix of component x)
-
-** HMM with GMM observations
- - hmm.name -
- - hmm.nstates - number of states in the HMM
- - hmm.emission_type - 'GMM'
- - hmm.start_prob - array of log probs P(first observation is state x)
- - hmm.end_prob - array of log probs P(last observation is state x)
- - hmm.transmat - matrix of transition log probs (transmat(x,y)
- = log(P(transition from state x to state y)))
- - hmm.labels - optional cell array of labels for each state in the HMM
- (for use in composing HMMs)
- - hmm.gmms - array of GMM structures
-
-** HMM with Gaussian observations
- - hmm.nstates - number of states in the HMM
- - hmm.emission_type - 'gaussian'
- - hmm.start_prob - array of log probs P(first observation is state x)
- - hmm.end_prob - array of log probs P(last observation is state x)
- - hmm.transmat - matrix of transition log probs (transmat(x,y)
- = log(P(transition from state x to state y)))
- - hmm.labels - optional cell array of labels for each state in the HMM
- (for use in composing HMMs)
- - hmm.means - matrix of means (column x is mean of state x)
- - hmm.covars - matrix of means (column x is the diagonal of the
- covariance matrix of component x)
-
-Note that each row of exp(hmm.transmat) does not necessarily sum to 1
-because for each state x there is some probability
-(exp(hmm.exit_prob(x))) that the next transition will be to a
-non-emitting exit state (i.e. the current observation is the last
-observation in the sequence). The correct invariant is:
-sum(exp(hmm.transmat, 2)) + exp(hmm.exit_prob) == ones(hmm.nstates, 1)
288 plottools/add_pan_and_zoom_controls_to_figure.m
View
@@ -1,288 +0,0 @@
-function add_pan_and_zoom_controls_to_figure(h_fig, all_axes)
-% add_pan_and_zoom_controls_to_figure(fig_handle, axes)
-%
-% Adds controls to the given figure window that let you pan and zoom
-% the vertical and horizontal axes and change the color limits of all
-% of the given axes at once. If no axes are specified, the controls
-% apply to all axes in the given figure. Separate controls are added
-% for each of the x, y, and color axes.
-%
-% If the limits are consistent across all of the specified axes
-% scrollbar is added at the edge of the figure. Otherwise simple
-% pan and zoom buttons are used.
-
-% Copyright (C) 2008 Ron J. Weiss (ronw@ee.columbia.edu)
-%
-% 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.
-%
-% You should have received a copy of the GNU General Public License
-% along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-if nargin < 1; h_fig = gcf; end
-
-% Axes are specified with functions to allow controls to behave
-% dynamically. I.e. if a new subplot is added after this function is
-% called, the controls can operate on the new subplot as well. Note
-% that this will lead to inconsistent scrollbar behavior if the new
-% subplot is not aligned with the subplots present at the time this
-% function was called.
-if nargin < 2
- xaxes_fun = @() get_axes(h_fig, 'axes');
- all_axes = xaxes_fun();
-else
- xaxes_fun = @() all_axes;
-end
-yaxes_fun = @() get_axes(feval(xaxes_fun), 'image');
-
-all_image_axes = get_axes(all_axes, 'image');
-all_images = findobj(all_axes, 'Type', 'image', 'Tag', '');
-if isempty(all_image_axes)
- yaxes_fun = xaxes_fun;
- all_image_axes = all_axes;
-end
-
-% Check to see if xlim, ylim, and clim are consistent across all axes.
-align_x_axes = are_limits_consistent(get(all_axes, 'XLim'));
-align_y_axes = are_limits_consistent(get(all_image_axes, 'YLim'));
-align_c_axes = are_limits_consistent(get(all_image_axes, 'CLim'));
-
-% Find current settings and absolute limits for xlim, ylim, and clim.
-% The parsing is kind of nasty because images and 1D plots need to be
-% handled differently and because the output type of get(axes) is
-% different for scalar axes and array axes.
-if ~isempty(all_images)
- current_xlim = get(all_axes, 'XLim');
- if iscell(current_xlim)
- current_xlim = cat(1, current_xlim{:});
- current_xlim = [min(current_xlim(1,:)) max(current_xlim(2,:))];
- end
- current_ylim = get(all_image_axes, 'YLim');
- if iscell(current_ylim)
- current_ylim = cat(1, current_ylim{:});
- current_ylim = [min(current_ylim(1,:)) max(current_ylim(2,:))];
- end
- current_clim = get(all_image_axes, 'CLim');
- if iscell(current_clim)
- current_clim = cat(1, current_clim{:});
- current_clim = [min(current_clim(1,:)) max(current_clim(2,:))];
- end
-
- cd = get(all_images, 'CData');
- if ~iscell(cd)
- cd = {cd};
- end
- sz_x = max(cellfun(@(x) size(x, 2), cd));
- sz_y = max(cellfun(@(x) size(x, 1), cd));
-
- for n = 1:length(cd)
- cd{n} = cd{n}(~isinf(cd{n}) & ~isnan(cd{n}));
- if isempty(cd{n})
- cd{n} = 0;
- end
- end
- clim_min = min(cellfun(@(x) min(x(:)), cd));
- clim_max = max(cellfun(@(x) max(x(:)), cd));
-
- original_xlim = [min(1, current_xlim(1)) max(sz_x, current_xlim(2))];
- original_ylim = [min(1, current_ylim(1)) max(sz_y, current_ylim(2))];
- original_clim = [clim_min clim_max];
-
- % Don't bother adding colorbar controls if the image is all
- % zeros. (Matlab sets the CLim of such image [-1 1] for some
- % reason, but there is no information there).
- if all(original_clim == 0)
- align_c_axes = false;
- end
-else
- align_c_axes = false;
- original_clim = [0 1];
- if length(all_axes) > 1
- min_x = min(cellfun(@(x) x(1), get(all_axes, 'XLim')));
- max_x = max(cellfun(@(x) x(2), get(all_axes, 'XLim')));
- original_xlim = [min_x max_x];
- min_y = min(cellfun(@(x) x(1), get(all_axes, 'YLim')));
- max_y = max(cellfun(@(x) x(2), get(all_axes, 'YLim')));
- original_ylim = [min_y max_y];
- else
- original_xlim = get(all_axes, 'XLim');
- original_ylim = get(all_axes, 'YLim');
- end
- current_xlim = original_xlim;
- current_ylim = original_ylim;
- current_clim = original_clim;
-end
-
-% Make sure the reset button resets the limits to the extremes
-% (i.e. make sure it works as expected if current_lim falls outside
-% of original_lim).
-original_xlim(1) = min(original_xlim(1), current_xlim(1));
-original_xlim(2) = max(original_xlim(2), current_xlim(2));
-original_ylim(1) = min(original_ylim(1), current_ylim(1));
-original_ylim(2) = max(original_ylim(2), current_ylim(2));
-original_clim(1) = min(original_clim(1), current_clim(1));
-original_clim(2) = max(original_clim(2), current_clim(2));
-
-setup_axis_controls(h_fig, 'XLim', original_xlim, current_xlim, ...
- xaxes_fun, align_x_axes);
-setup_axis_controls(h_fig, 'YLim', original_ylim, current_ylim, ...
- yaxes_fun, align_y_axes);
-setup_axis_controls(h_fig, 'CLim', original_clim, current_clim, ...
- yaxes_fun, align_c_axes);
-
-% Link axes so that standard zoom/pan controls apply to all axes (but
-% not colorbars).
-%linkaxes(all_axes, 'xy')
-
-% Make sure toolbar is still displayed.
-set(h_fig, 'Toolbar', 'figure')
-
-
-function y = are_limits_consistent(lim)
-y = true;
-if iscell(lim) && any(cat(2, lim{:}) ~= repmat(lim{1}, [1 length(lim)]))
- y = false;
-end
-
-
-function setup_axis_controls(h_fig, property, orig_lim, curr_lim, axes_fun, align_axes)
-
-if align_axes
- % Add a scrollbar control for the given property of the given axes.
- switch lower(property(1))
- case 'x'
- hidden_panel_pos = [0.05 -1.0 0.9 0.9];
- control_location = 'Bottom';
- case 'y'
- hidden_panel_pos = [-1.0 0.05 0.9 0.9];
- control_location = 'Left';
- case 'c'
- hidden_panel_pos = [0 0.05 1.05 0.9];
- control_location = 'Right';
- otherwise
- return
- end
-
- % Create hidden panel to control scrollbar placement.
- h_hidden_panel = uipanel('Parent', h_fig, 'Tag', mfilename, ...
- 'Units', 'normalized', 'Position', hidden_panel_pos, ...
- 'Visible', 'off', 'HitTest', 'off');
- normalized_lim = (curr_lim - orig_lim(1)) / (orig_lim(2) - orig_lim(1));
- h = zooming_scrollbar(h_hidden_panel, normalized_lim, ...
- @(a,b,c,d) scrollbar_update_axes(a, b, c, d, property, orig_lim, ...
- axes_fun, h_fig), control_location);
- set(h, 'Tag', mfilename);
- set(feval(axes_fun), property, curr_lim)
-else
- % Since axes have inconsistent axes, we cannot create a single
- % scrollbar control for all of them. Instead add simple button
- % controls that zoom and pan all axes equally.
-
- switch lower(property(1))
- case 'x'
- pos = [40 0 20 20];
- dim = 1;
- scroll_strings = '<>';
- case 'y'
- pos = [ 0 40 20 20];
- dim = 2;
- scroll_strings = 'v^';
- otherwise
- return
- end
-
- ZOOM = 0.9;
- PAN = 0.05;
- args = {property, axes_fun};
- uicontrol(h_fig, 'Style', 'pushbutton', 'String', 'R', 'Position', pos, ...
- 'Tag', mfilename, ...
- 'Callback', @(a,b) simple_reset_callback(args{:}));
- pos(dim) = pos(dim) + pos(dim+2);
- uicontrol(h_fig, 'Style', 'pushbutton', 'String', '+', 'Position', pos, ...
- 'Tag', mfilename, ...
- 'Callback', @(a,b) simple_zoom_callback(args{:}, ZOOM));
- pos(dim) = pos(dim) + pos(dim+2);
- uicontrol(h_fig, 'Style', 'pushbutton', 'String', '-', 'Position', pos, ...
- 'Tag', mfilename, ...
- 'Callback', @(a,b) simple_zoom_callback(args{:}, 1/ZOOM));
- pos(dim) = pos(dim) + pos(dim+2);
- uicontrol(h_fig, 'Style', 'pushbutton', 'String', scroll_strings(1), ...
- 'Tag', mfilename, ...
- 'Position', pos, ...
- 'Callback', @(a,b) simple_scroll_callback(args{:}, -PAN));
- pos(dim) = pos(dim) + pos(dim+2);
- uicontrol(h_fig, 'Style', 'pushbutton', 'String', scroll_strings(2), ...
- 'Tag', mfilename, ...
- 'Position', pos, ...
- 'Callback', @(a,b) simple_scroll_callback(args{:}, PAN));
-end
-
-
-function scrollbar_update_axes(h_scrollbar, h_hidden_panel, min_val, ...
- max_val, prop, original_lim, axes_fun, h_fig)
-minimum = original_lim(1);
-range = original_lim(2) - original_lim(1);
-new_lim = [min_val max_val]*range + minimum;
-new_range = new_lim(2) - new_lim(1);
-
-all_axes = feval(axes_fun);
-set(all_axes, prop, new_lim);
-
-if strcmpi(prop, 'CLim')
- for cb_axis = findobj(h_fig, 'Tag', 'Colorbar')'
- yt = get(cb_axis, 'YTick');
- ntick = length(yt);
- yl = get(cb_axis, 'YLim');
- cb_range = yl(2) - yl(1);
-
- ytl = num2cell((yt - yl(1)) * new_range/cb_range + new_lim(1));
- for n = 1:length(ytl)
- ytl{n} = str2num(sprintf('%.2e ', ytl{n}));
- end
-
- set(cb_axis, 'YTickLabel', ytl);
- end
-end
-
-</