diff --git a/migram.m b/migram.m index 93d98a24d9f4f2723a4c76e1a0732ba1762a6bdf..f7f6897ab4760e8d4caf39f52dbefe6a61431918 100644 --- a/migram.m +++ b/migram.m @@ -16,13 +16,16 @@ function [c_out, l_out, t_out] = migram(varargin) % I = MIGRAM(A,B,MAXLAG,WINDOW,NOVERLAP,NBINS) calculates the mutual % information based on histograms with the number of bins NBINS. % +% I = MIGRAM(...,'norm') calculates the renormalised mutual +% information, which is I/log(NBINS) and ensures a value range [0 1]. +% % [I,L,T] = MIGRAM(...) returns a column of lag L and one of time T % at which the mutual information is computed. L has length equal % to the number of rows of I, T has length k. % % I = MIGRAM(A,B) calculates windowed mutual information using defeault % settings; the defeaults are MAXLAG = floor(0.1*N), WINDOW = floor(0.1*N), -% NOVERLAP = 0 and NBINS = 10. You can tell CORRGRAM to use the defeault +% NOVERLAP = 0 and NBINS = 10. You can tell MIGRAM to use the defeault % for any parameter by leaving it off or using [] for that parameter, e.g. % MIGRAM(A,B,[],1000). % @@ -44,7 +47,7 @@ function [c_out, l_out, t_out] = migram(varargin) % $Revision$ -error(nargchk(2,6,nargin)) +error(nargchk(2,7,nargin)) verbose = 0; x = varargin{1}; y = varargin{2}; @@ -66,32 +69,40 @@ maxlag = floor(nx/10); window = floor(nx/10); noverlap = 0; nbins = 10; +norm = 0; + +i_num = find(cellfun('isclass',varargin,'double')); +i_char = find(cellfun('isclass',varargin,'char')); -if length(varargin) > 2 & ~isempty(varargin{3}) - maxlag = varargin{3}; +if length(i_num) > 2 & ~isempty(varargin{i_num(3)}) + maxlag = varargin{i_num(3)}; if maxlag < 0, error('Requires positive integer value for maximum lag.'), end if length(maxlag) > 1, error('Requires MAXLAG to be a scalar.'), end end -if length(varargin) > 3 & ~isempty(varargin{4}) - window = varargin{4}; +if length(i_num) > 3 & ~isempty(varargin{i_num(4)}) + window = varargin{i_num(4)}; if window <= 0, error('Requires positive integer value for window length.'), end if length(window) > 1, error('Requires WINDOW to be a scalar.'), end end -if length(varargin) > 4 & ~isempty(varargin{5}) - noverlap = varargin{5}; +if length(i_num) > 4 & ~isempty(varargin{i_num(5)}) + noverlap = varargin{i_num(5)}; if noverlap < 0, error('Requires positive integer value for NOVERLAP.'), end if length(noverlap) > 1, error('Requires NOVERLAP to be a scalar.'), end if noverlap >= window, error('Requires NOVERLAP to be strictly less than the window length.'), end end -if length(varargin) > 5 & ~isempty(varargin{6}) - noverlap = varargin{6}; +if length(i_num) > 5 & ~isempty(varargin{i_num(6)}) + noverlap = varargin{i_num(6)}; if nbins <= 0, error('Requires positive integer value for NBINS.'), end if length(nbins) > 1, error('Requires NBINS to be a scalar.'), end end +% normalise the result +for i = 1:length(i_char) + if strcmpi(varargin(i_char(i)), 'norm'), norm = 1; end +end % prepare time delayed signals @@ -108,7 +119,6 @@ if verbose, h = waitbar(0,'Compute mutual information'); end % -MAXLAG:0 [Yi dummy] = buffer(Y(:,1),window,noverlap,'nodelay'); -Yi = normalize(Yi); if exist('accumarray','builtin') == 5 for i = 1:size(X,2), if verbose, waitbar(cnt/(2*size(X,2))), end [Xi dummy] = buffer(X(:,i),window,noverlap,'nodelay'); @@ -148,6 +158,12 @@ warning on t = (1:nx/size(Xi,2):nx)'; l = (-maxlag:maxlag)'; +% if result has to be normalised +if norm + C = C / log(nbins); +end + + % display and output result if nargout == 0 newplot @@ -236,4 +252,3 @@ function Z = MI5(x, y, nbins) Z(i) = Ix + Iy - Ixy; end -