octave_packages/tsa-4.2.4/histo3.m

   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.
   4 %
   5 % R = HISTO3(Y)
   6 % R = HISTO3(Y, W)
   7 %       Y       data
   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
  17 %
  18 % Data compression can be performed in this way
  19 %       [R,tix] = histo3(Y)
  20 %               is the compression step
  21 %
  22 %       R.tix provides a compressed data representation.
  23 %       R.compressionratio estimates the compression ratio
  24 %
  25 %       R.X(tix) and R.X(R.tix)
  26 %               reconstruct the orginal signal (decompression)
  27 %
  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.
  30 %
  31 % see also: HISTO, HISTO2, HISTO3, HISTO4
  32 %
  33 % REFERENCE(S):
  34 %  C.E. Shannon and W. Weaver "The mathematical theory of communication" University of Illinois Press, Urbana 1949 (reprint 1963).
  35
  36
  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/
  43 %
  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.
  48 %
  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.
  53 %
  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/>.
  56
  57
  58 %%%%% check input arguments %%%%%
  59 [yr,yc] = size(Y);
  60 if nargin < 2,
  61         W = [];
  62 end;
  63 if ~isempty(W) && (yr ~= numel(W)),
  64         error('number of rows of Y does not match number of elements in W');
  65 end;
  66
  67 %%%%% identify all possible X's and generate overall Histogram %%%%%
  68 [sY, idx] = sort(Y(:),1);
  69 [tmp,idx1] = sort(idx);        % generate inverse index
  70
  71 ix  = diff(sY, [], 1) > 0;
  72 tmp = [find(ix); sum(~isnan(sY))];
  73
  74 R.datatype = 'HISTOGRAM';
  75 R.X = sY(tmp);
  76 R.N = sum(~isnan(Y), 1);
  77
  78 % generate inverse index
  79 if nargout>1,
  80         tix = cumsum([1; ix]);  % rank
  81         tix = reshape(tix(idx1), yr, yc);               % inverse sort rank
  82         cc  = 1;
  83         tmp = sum(ix) + 1;
  84         if exist('OCTAVE_VERSION') >= 5,
  85                 ; % NOP; no support for integer datatyp
  86         elseif tmp <= 2^8;
  87                 tix = uint8(tix);
  88                 cc = 8/1;
  89         elseif tmp <= 2^16;
  90                 tix = uint16(tix);
  91                 cc = 8/2;
  92         elseif tmp <= 2^32;
  93                 tix = uint32(tix);
  94                 cc = 8/4;
  95         end;
  96         R.compressionratio = (prod(size(R.X)) + (yr*yc)/cc) / (yr*yc);
  97         R.tix = tix;
  98 end;
  99
 100
 101 if yc==1,
 102         if isempty(W)
 103                 R.H = [tmp(1); diff(tmp)];
 104         else
 105                 C = cumsum(W(idx));     % cumulative weights
 106                 R.H = [C(tmp(1)); diff(C(tmp))];
 107         end;
 108         return;
 109
 110 elseif yc>1,
 111         % allocate memory
 112         H = zeros(size(R.X,1),yc);
 113
 114         % scan each channel
 115         for k = 1:yc,
 116                 if isempty(W)
 117                         sY = sort(Y(:,k));
 118                 else
 119                         [sY,ix] = sort(Y(:,k));
 120                         C = cumsum(W(ix));
 121                 end
 122                 ix = find(diff(sY,[],1) > 0);
 123                 if size(ix,1) > 0,
 124                         tmp = [ix; R.N(k)];
 125                 else
 126                         tmp = R.N(k);
 127                 end;
 128
 129                 t = 0;
 130                 j = 1;
 131                 if isempty(W)
 132                     for x = tmp',
 133                         acc = sY(x);
 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)
 137                         t = x;
 138                     end;
 139                 else
 140                     for x = tmp',
 141                         acc = sY(x);
 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)
 145                         t = C(x);
 146                     end;
 147                 end;
 148         end;
 149
 150         R.H = H;
 151 end;
 152
 153