From ab8ede08c7507ac439016e2028ec6228f12c675d Mon Sep 17 00:00:00 2001
From: marwan <>
Date: Thu, 14 Jun 2007 08:52:19 +0000
Subject: [PATCH] *** empty log message ***
---
corrgram.m | 147 ++++++++++++++++++++++++++++++++
migram.m | 239 +++++++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 386 insertions(+)
create mode 100644 corrgram.m
create mode 100644 migram.m
diff --git a/corrgram.m b/corrgram.m
new file mode 100644
index 0000000..a30c962
--- /dev/null
+++ b/corrgram.m
@@ -0,0 +1,147 @@
+function [c_out, l_out, t_out] = corrgram(varargin)
+%CORRGRAM Calculate windowed cross correlation between two signals.
+% C = CORRGRAM(A,B,MAXLAG,WINDOW,NOVERLAP) calculates the windowed cross
+% correlation between the signals in vector A and vector B. CORRGRAM splits
+% the signals into overlapping segments and forms the columns of C with
+% their cross correlation values up to maximum lag specified by scalar
+% MAXLAG. Each column of C contains the cross correlation function between
+% the short-term, time-localized signals A and B. Time increases linearly
+% across the columns of C, from left to right. Lag increases linearly down
+% the rows, starting at -MAXLAG. If lengths of A and B differ, the shorter
+% signal is filled with zeros. If N is the length of the signals, C is
+% a matrix with 2*MAXLAG+1 rows and
+% k = fix((N-NOVERLAP)/(WINDOW-NOVERLAP))
+% columns.
+%
+% [C,L,T] = CORRGRAM(...) returns a column of lag L and one of time T
+% at which the correlation coefficients are computed. L has length equal
+% to the number of rows of C, T has length k.
+%
+% C = CORRGRAM(A,B) calculates windowed cross correlation using defeault
+% settings; the defeaults are MAXLAG = floor(0.1*N), WINDOW = floor(0.1*N)
+% and NOVERLAP = 0. You can tell CORRGRAM to use the defeault for any
+% parameter by leaving it off or using [] for that parameter, e.g.
+% CORRGRAM(A,B,[],1000).
+%
+% CORRGRAM(A,B) with no output arguments plots the windowed cross
+% correlation using the current figure.
+%
+% EXAMPLE:
+% x = cos(0:.01:10*pi)';
+% y = sin(0:.01:10*pi)' + .5 * randn(length(x),1);
+% corrgram(x,y)
+%
+% See also CORRCOEF, CORR, XCORR, MIGRAM.
+
+% Copyright (c) 2007 by AMRON
+% Norbert Marwan, Potsdam University, Germany
+% http://www.agnld.uni-potsdam.de
+%
+% $Date$
+% $Revision$
+
+
+error(nargchk(2,5,nargin))
+verbose = 0;
+
+x = varargin{1}; y = varargin{2};
+x = x(:); y = y(:);
+
+
+% check input and inital setting of parameters
+
+nx = length(x); ny = length(y);
+if nx < ny % zero-pad x if it has length less than y
+ x(ny) = 0; nx = ny;
+end
+
+if ny < nx % zero-pad y if it has length less than x
+ y(nx) = 0;
+end
+
+maxlag = floor(nx/10);
+window = floor(nx/10);
+noverlap = 0;
+
+if length(varargin) > 2 & ~isempty(varargin{3})
+ maxlag = varargin{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 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 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
+
+
+% prepare time delayed signals
+X = buffer(x,maxlag+1,maxlag)';
+Y = fliplr(buffer(y,maxlag+1,maxlag)');
+
+% divide the delayed signals into overlapping windows
+% and compute the correlation coefficient
+cnt = 1;
+
+warning off
+C = zeros(2*maxlag+1, fix((nx-noverlap)/(window-noverlap)));
+if verbose, h = waitbar(0,'Compute cross correlation'); end
+
+% -MAXLAG:0
+[Yi dummy] = buffer(Y(:,1),window,noverlap,'nodelay');
+Yi = normalize(Yi);
+for i = 1:size(X,2), if verbose, waitbar(cnt/(2*size(X,2))), end
+ [Xi dummy] = buffer(X(:,i),window,noverlap,'nodelay');
+ Xi = normalize(Xi);
+ C(cnt,:) = mean(Xi .* Yi);
+ cnt = cnt + 1;
+end
+
+% 0:MAXLAG
+[Xi dummy] = buffer(X(:,end),window,noverlap,'nodelay');
+Xi = normalize(Xi);
+for i = 2:size(Y,2), if verbose, waitbar(cnt/(2*size(X,2))), end
+ [Yi dummy] = buffer(Y(:,i),window,noverlap,'nodelay');
+ Yi = normalize(Yi);
+ C(cnt,:) = mean(Xi .* Yi);
+ cnt = cnt + 1;
+end
+if verbose, delete(h), end
+
+warning on
+
+% create time scale for the windows
+t = (1:nx/size(Xi,2):nx)';
+l = (-maxlag:maxlag)';
+
+% display and output result
+if nargout == 0
+ newplot
+ imagesc(t, l, C)
+ xlabel('Time'), ylabel('Lag'), axis xy
+ title('Windowed cross correlation', 'fontweight', 'bold')
+ colorbar
+elseif nargout == 1,
+ c_out = C;
+elseif nargout == 2,
+ c_out = C;
+ l_out = l;
+elseif nargout == 3,
+ c_out = C;
+ l_out = l;
+ t_out = t;
+end
+
+
+function Y = normalize(X)
+ Y = X - repmat(mean(X), size(X,1), 1);
+ s = sqrt( sum(conj(Y) .* Y) / (size(Y,1) - 1) );
+ Y = Y ./ repmat(s, size(Y,1), 1);
diff --git a/migram.m b/migram.m
new file mode 100644
index 0000000..93d98a2
--- /dev/null
+++ b/migram.m
@@ -0,0 +1,239 @@
+function [c_out, l_out, t_out] = migram(varargin)
+%MIGRAM Calculate windowed mutual information between two signals.
+% I = MIGRAM(A,B,MAXLAG,WINDOW,NOVERLAP) calculates the windowed mutual
+% information between the signals in vector A and vector B. MIGRAM splits
+% the signals into overlapping segments and forms the columns of I with
+% their mutual information values up to maximum lag specified by scalar
+% MAXLAG. Each column of I contains the mutual information function
+% between the short-term, time-localized signals A and B. Time increases
+% linearly across the columns of I, from left to right. Lag increases
+% linearly down the rows, starting at -MAXLAG. If lengths of A and B
+% differ, the shorter signal is filled with zeros. If N is the length of
+% the signals, I is a matrix with 2*MAXLAG+1 rows and
+% k = fix((N-NOVERLAP)/(WINDOW-NOVERLAP))
+% columns.
+%
+% I = MIGRAM(A,B,MAXLAG,WINDOW,NOVERLAP,NBINS) calculates the mutual
+% information based on histograms with the number of bins NBINS.
+%
+% [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
+% for any parameter by leaving it off or using [] for that parameter, e.g.
+% MIGRAM(A,B,[],1000).
+%
+% MIGRAM(A,B) with no output arguments plots the mutual information
+% using the current figure.
+%
+% EXAMPLE:
+% x = cos(0:.01:10*pi)';
+% y = sin(0:.01:10*pi)' + .5 * randn(length(x),1);
+% migram(x,y)
+%
+% See also MI, CORRGRAM.
+
+% Copyright (c) 2007 by AMRON
+% Norbert Marwan, Potsdam University, Germany
+% http://www.agnld.uni-potsdam.de
+%
+% $Date$
+% $Revision$
+
+
+error(nargchk(2,6,nargin))
+verbose = 0;
+
+x = varargin{1}; y = varargin{2};
+x = x(:); y = y(:);
+
+
+% check input and inital setting of parameters
+
+nx = length(x); ny = length(y);
+if nx < ny % zero-pad x if it has length less than y
+ x(ny) = 0; nx = ny;
+end
+
+if ny < nx % zero-pad y if it has length less than x
+ y(nx) = 0;
+end
+
+maxlag = floor(nx/10);
+window = floor(nx/10);
+noverlap = 0;
+nbins = 10;
+
+if length(varargin) > 2 & ~isempty(varargin{3})
+ maxlag = varargin{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 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 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 nbins <= 0, error('Requires positive integer value for NBINS.'), end
+ if length(nbins) > 1, error('Requires NBINS to be a scalar.'), end
+end
+
+
+
+% prepare time delayed signals
+X = buffer(x,maxlag+1,maxlag)';
+Y = fliplr(buffer(y,maxlag+1,maxlag)');
+
+% divide the delayed signals into overlapping windows
+% and compute the correlation coefficient
+cnt = 1;
+
+warning off
+C = zeros(2*maxlag+1, fix((nx-noverlap)/(window-noverlap)));
+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');
+ C(cnt,:) = MI6(Xi, Yi, nbins);
+ cnt = cnt + 1;
+ end
+else
+ for i = 1:size(X,2), if verbose, waitbar(cnt/(2*size(X,2))), end
+ [Xi dummy] = buffer(X(:,i),window,noverlap,'nodelay');
+ C(cnt,:) = MI5(Xi, Yi, nbins);
+ cnt = cnt + 1;
+ end
+end
+
+% 0:MAXLAG
+[Xi dummy] = buffer(X(:,end),window,noverlap,'nodelay');
+Xi = normalize(Xi);
+if exist('accumarray','builtin') == 5
+ for i = 2:size(Y,2), if verbose, waitbar(cnt/(2*size(X,2))), end
+ [Yi dummy] = buffer(Y(:,i),window,noverlap,'nodelay');
+ C(cnt,:) = MI6(Xi, Yi, nbins);
+ cnt = cnt + 1;
+ end
+else
+ for i = 2:size(Y,2), if verbose, waitbar(cnt/(2*size(X,2))), end
+ [Yi dummy] = buffer(Y(:,i),window,noverlap,'nodelay');
+ C(cnt,:) = MI5(Xi, Yi, nbins);
+ cnt = cnt + 1;
+ end
+end
+
+if verbose, delete(h), end
+
+warning on
+
+% create time scale for the windows
+t = (1:nx/size(Xi,2):nx)';
+l = (-maxlag:maxlag)';
+
+% display and output result
+if nargout == 0
+ newplot
+ imagesc(t, l, C)
+ xlabel('Time'), ylabel('Lag'), axis xy
+ title('Windowed mutual information', 'fontweight', 'bold')
+ colorbar
+elseif nargout == 1,
+ c_out = C;
+elseif nargout == 2,
+ c_out = C;
+ l_out = l;
+elseif nargout == 3,
+ c_out = C;
+ t_out = t;
+ l_out = l;
+end
+
+
+% mutual information for Matlab version >= 6
+function Z = MI6(x, y, nbins)
+
+ % normalise the data and replace the values with integers
+ % in the range [1 nbins]
+ x = x - repmat(min(x), size(x,1), 1);
+ y = y - repmat(min(y), size(y,1), 1);
+ x = x ./ repmat(max(x) + eps, size(x,1), 1);
+ y = y ./ repmat(max(y) + eps, size(y,1), 1);
+
+ x = floor(x * nbins) + 1;
+ y = floor(y * nbins) + 1;
+
+ % compute probabilities
+ Z = zeros(1,size(x,2));
+ for i = 1:size(x,2)
+
+ Pxy = accumarray([x(:,i) y(:,i)] + 1, 1);
+ Px = sum(Pxy,1);
+ Py = sum(Pxy,2);
+
+ Pxy = Pxy / sum(Pxy(:));
+ Px = Px / sum(Px(:));
+ Py = Py / sum(Py(:));
+
+ % entropies
+ Ix = -sum((Px(Px ~= 0)) .* log(Px(Px ~= 0)));
+ Iy = -sum((Py(Py ~= 0)) .* log(Py(Py ~= 0)));
+ Ixy = -sum(Pxy(Pxy ~= 0) .* log(Pxy(Pxy ~= 0)));
+
+ % mutual information
+ Z(i) = Ix + Iy - Ixy;
+ end
+
+
+% mutual information for Matlab version < 6
+function Z = MI5(x, y, nbins)
+
+ % normalise the data and replace the values with integers
+ % in the range [1 nbins]
+ x = x - repmat(min(x), size(x,1), 1);
+ y = y - repmat(min(y), size(y,1), 1);
+ x = x ./ repmat(max(x) + eps, size(x,1), 1);
+ y = y ./ repmat(max(y) + eps, size(y,1), 1);
+
+ x = floor(x * nbins) + 1;
+ y = floor(y * nbins) + 1;
+
+ % compute probabilities
+ Z = zeros(1,size(x,2));
+ for i = 1:size(x,2)
+
+ Pxy = full(sparse(x(:,i) + 1, y(:,i) + 1, 1));
+ Px = sum(Pxy,1);
+ Py = sum(Pxy,2);
+
+ Pxy = Pxy / sum(Pxy(:));
+ Px = Px / sum(Px(:));
+ Py = Py / sum(Py(:));
+
+ % entropies
+ Ix = -sum((Px(Px ~= 0)) .* log(Px(Px ~= 0)));
+ Iy = -sum((Py(Py ~= 0)) .* log(Py(Py ~= 0)));
+ Ixy = -sum(Pxy(Pxy ~= 0) .* log(Pxy(Pxy ~= 0)));
+
+ % mutual information
+ Z(i) = Ix + Iy - Ixy;
+ end
+
+
--
GitLab