X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Ftsa-4.2.4%2Fhisto3.m;fp=octave_packages%2Ftsa-4.2.4%2Fhisto3.m;h=edf515080cda244e359689911ab987c2dcd19082;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/tsa-4.2.4/histo3.m b/octave_packages/tsa-4.2.4/histo3.m
new file mode 100644
index 0000000..edf5150
--- /dev/null
+++ b/octave_packages/tsa-4.2.4/histo3.m
@@ -0,0 +1,153 @@
+function [R, tix] = histo3(Y, W)
+% HISTO3 calculates histogram for multiple columns with common bin values 
+%    among all data columns, and can be useful for data compression. 
+%
+% R = HISTO3(Y)
+% R = HISTO3(Y, W)
+%	Y	data
+%	W	weight vector containing weights of each sample, 
+%		number of rows of Y and W must match.
+%		default W=[] indicates that each sample is weighted with 1. 
+% 	R 	struct with these fields 
+%       R.X  	the bin-values, bin-values are equal for each channel
+%		thus R.X is a column vector. If bin values should 
+%		be computed separately for each data column, use HISTO2
+%       R.H  	is the frequency of occurence of value X 
+%  	R.N  	are the number of valid (not NaN) samples 
+%
+% Data compression can be performed in this way
+%   	[R,tix] = histo3(Y) 
+%      		is the compression step
+%
+%	R.tix provides a compressed data representation. 
+%	R.compressionratio estimates the compression ratio
+%
+% 	R.X(tix) and R.X(R.tix) 
+%		reconstruct the orginal signal (decompression) 
+%
+% The effort (in memory and speed) for compression is O(n*log(n)).
+% The effort (in memory and speed) for decompression is O(n) only. 
+%
+% see also: HISTO, HISTO2, HISTO3, HISTO4
+%
+% REFERENCE(S):
+%  C.E. Shannon and W. Weaver "The mathematical theory of communication" University of Illinois Press, Urbana 1949 (reprint 1963).
+
+
+%	$Id: histo3.m 8383 2011-07-16 20:06:59Z schloegl $
+%	Copyright (C) 1996-2002,2008,2011 by Alois Schloegl <a.schloegl@ieee.org>	
+%       This is part of the TSA-toolbox. See also 
+%       http://hci.tugraz.at/schloegl/matlab/tsa/
+%       http://octave.sourceforge.net/
+%       http://biosig.sourceforge.net/
+%
+%    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/>.
+
+
+%%%%% check input arguments %%%%%
+[yr,yc] = size(Y);
+if nargin < 2, 
+	W = []; 
+end; 
+if ~isempty(W) && (yr ~= numel(W)),
+	error('number of rows of Y does not match number of elements in W');
+end;
+
+%%%%% identify all possible X's and generate overall Histogram %%%%%
+[sY, idx] = sort(Y(:),1);
+[tmp,idx1] = sort(idx);        % generate inverse index
+
+ix  = diff(sY, [], 1) > 0;
+tmp = [find(ix); sum(~isnan(sY))];
+
+R.datatype = 'HISTOGRAM';
+R.X = sY(tmp);
+R.N = sum(~isnan(Y), 1);
+
+% generate inverse index
+if nargout>1,
+        tix = cumsum([1; ix]);	% rank 
+        tix = reshape(tix(idx1), yr, yc);		% inverse sort rank
+        cc  = 1;
+        tmp = sum(ix) + 1;
+	if exist('OCTAVE_VERSION') >= 5,
+		; % NOP; no support for integer datatyp 
+        elseif tmp <= 2^8;
+                tix = uint8(tix);
+                cc = 8/1;
+        elseif tmp <= 2^16;
+                tix = uint16(tix);
+                cc = 8/2;
+        elseif tmp <= 2^32;
+                tix = uint32(tix);
+                cc = 8/4;
+        end;
+        R.compressionratio = (prod(size(R.X)) + (yr*yc)/cc) / (yr*yc);
+	R.tix = tix;        
+end;
+
+
+if yc==1, 
+	if isempty(W)
+		R.H = [tmp(1); diff(tmp)];
+	else
+		C = cumsum(W(idx));  	% cumulative weights  
+		R.H = [C(tmp(1)); diff(C(tmp))];
+	end;
+	return;
+
+elseif yc>1,
+        % allocate memory
+        H = zeros(size(R.X,1),yc);
+        
+        % scan each channel
+        for k = 1:yc,
+		if isempty(W)
+			sY = sort(Y(:,k));
+		else
+			[sY,ix] = sort(Y(:,k));
+			C = cumsum(W(ix));
+		end
+		ix = find(diff(sY,[],1) > 0);
+                if size(ix,1) > 0,
+                        tmp = [ix; R.N(k)];
+                else
+                        tmp = R.N(k);
+                end;
+
+                t = 0;
+                j = 1;
+		if isempty(W)
+                    for x = tmp',
+                        acc = sY(x);
+                        while R.X(j)~=acc, j=j+1; end;
+                        %j = find(sY(x)==R.X);   % identify position on X 
+                        H(j,k) = H(j,k) + (x-t);  % add diff(tmp)
+                        t = x;
+                    end;
+		else
+                    for x = tmp',
+                        acc = sY(x);
+                        while R.X(j)~=acc, j=j+1; end;
+                        %j = find(sY(x)==R.X);   % identify position on X 
+                        H(j,k) = H(j,k) + C(x)-t;  % add diff(tmp)
+                        t = C(x);
+                    end;
+		end; 
+        end;
+        
+	R.H = H;
+end;
+
+