1 function [R, tix] = histo3(Y, W)
2 % HISTO3 calculates histogram for multiple columns with common bin values
3 % among all data columns, and can be useful for data compression.
8 % W weight vector containing weights of each sample,
9 % number of rows of Y and W must match.
10 % default W=[] indicates that each sample is weighted with 1.
11 % R struct with these fields
12 % R.X the bin-values, bin-values are equal for each channel
13 % thus R.X is a column vector. If bin values should
14 % be computed separately for each data column, use HISTO2
15 % R.H is the frequency of occurence of value X
16 % R.N are the number of valid (not NaN) samples
18 % Data compression can be performed in this way
20 % is the compression step
22 % R.tix provides a compressed data representation.
23 % R.compressionratio estimates the compression ratio
25 % R.X(tix) and R.X(R.tix)
26 % reconstruct the orginal signal (decompression)
28 % The effort (in memory and speed) for compression is O(n*log(n)).
29 % The effort (in memory and speed) for decompression is O(n) only.
31 % see also: HISTO, HISTO2, HISTO3, HISTO4
34 % C.E. Shannon and W. Weaver "The mathematical theory of communication" University of Illinois Press, Urbana 1949 (reprint 1963).
37 % $Id: histo3.m 8383 2011-07-16 20:06:59Z schloegl $
38 % Copyright (C) 1996-2002,2008,2011 by Alois Schloegl <a.schloegl@ieee.org>
39 % This is part of the TSA-toolbox. See also
40 % http://hci.tugraz.at/schloegl/matlab/tsa/
41 % http://octave.sourceforge.net/
42 % http://biosig.sourceforge.net/
44 % This program is free software: you can redistribute it and/or modify
45 % it under the terms of the GNU General Public License as published by
46 % the Free Software Foundation, either version 3 of the License, or
47 % (at your option) any later version.
49 % This program is distributed in the hope that it will be useful,
50 % but WITHOUT ANY WARRANTY; without even the implied warranty of
51 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
52 % GNU General Public License for more details.
54 % You should have received a copy of the GNU General Public License
55 % along with this program. If not, see <http://www.gnu.org/licenses/>.
58 %%%%% check input arguments %%%%%
63 if ~isempty(W) && (yr ~= numel(W)),
64 error('number of rows of Y does not match number of elements in W');
67 %%%%% identify all possible X's and generate overall Histogram %%%%%
68 [sY, idx] = sort(Y(:),1);
69 [tmp,idx1] = sort(idx); % generate inverse index
71 ix = diff(sY, [], 1) > 0;
72 tmp = [find(ix); sum(~isnan(sY))];
74 R.datatype = 'HISTOGRAM';
76 R.N = sum(~isnan(Y), 1);
78 % generate inverse index
80 tix = cumsum([1; ix]); % rank
81 tix = reshape(tix(idx1), yr, yc); % inverse sort rank
84 if exist('OCTAVE_VERSION') >= 5,
85 ; % NOP; no support for integer datatyp
96 R.compressionratio = (prod(size(R.X)) + (yr*yc)/cc) / (yr*yc);
103 R.H = [tmp(1); diff(tmp)];
105 C = cumsum(W(idx)); % cumulative weights
106 R.H = [C(tmp(1)); diff(C(tmp))];
112 H = zeros(size(R.X,1),yc);
119 [sY,ix] = sort(Y(:,k));
122 ix = find(diff(sY,[],1) > 0);
134 while R.X(j)~=acc, j=j+1; end;
135 %j = find(sY(x)==R.X); % identify position on X
136 H(j,k) = H(j,k) + (x-t); % add diff(tmp)
142 while R.X(j)~=acc, j=j+1; end;
143 %j = find(sY(x)==R.X); % identify position on X
144 H(j,k) = H(j,k) + C(x)-t; % add diff(tmp)