マルチスケールエントロピーモデル——MSE

MSE-matlabコード
シンボル説明:x入力信号ベクトル(例えば、脳電気信号または音声信号)nSfスケール係数の数mテンプレート長(暦長);Costaはm=2 r−マッチング閾値を用いる.通常、時間系列におけるサンプル偏差が選択される.xがz変換されると、公差は、rに標準偏差mse−マルチスケールエントロピーsf−を乗じたmseに対応するスケール係数解釈として定義される.Costaは、このスケール上のより多くの情報を反映するためにより高いエントロピー値を解釈する(ランダムであれば、あまり予測できない).天秤の上ではかなり安定している.
参考文献:Costaら.(2002)複雑な生理学的多尺度エントロピー分析時間系列.Costaら.(2005)生物信号のマルチスケールエントロピー解析.
MSE_Costa2005

function [mse sf] = MSE_Costa2005(x,nSf,m,r)

% pre-allocate mse vector
mse = zeros([1 nSf]);

% coarse-grain and calculate sample entropy for each scale factor
for ii = 1 : nSf
	% get filter weights
	f = ones([1 ii]);
	f = f/sum(f);

	% get coarse-grained time series (i.e., average data within non-overlapping time windows)
	y = filter(f,1,x);
	y = y(length(f):end);
	y = y(1:length(f):end);
	
	% calculate sample entropy
	mse(ii) = SampleEntropy(y,m,r,0);
end

% get sacle factors
sf = 1 : nSf;

colored_noise

function [y] = colored_noise(Sf,dur,b)

% [y] = colored_noise(Sf,dur,b)
%
% Sf  - sampling frequency
% dur - duration
% b   - exponent of the 1/f^b power law function
%       b = 0 (white), 1 (pink), 2 (brown), -1 (blue), -2 (violet)
%
% y - noise time domain signal
%     (y is normalized such that mean(y)==0 and std(y)==1)
%
% Description: The script calculates a colored noise.
% get number of samples
N = round(Sf * dur);
if mod(N,2)
    M = N + 1;
else
    M = N;
end

% get fourier spectrum and the corresponding frequency axis
F  = fft(rand([M 1]));
xf = (0:M/2-1)*(Sf/M);

% power law function
p = 1./xf(2:end).^b;

% get scaled magnitudes and phase
R = sqrt((abs(F(2:length(xf))).^2).*p');
P = angle(F(2:length(xf)));
	
% get symmetrical spectrum
F(2:length(xf))     = R .* (cos(P) + 1i.*sin(P));
F(length(xf)+2:end) = real(F(length(xf):-1:2)) - 1i*imag(F(length(xf):-1:2));

% IFFT
y = ifft(F);

% get correct number of samples
y = real(y(1:N));

% ensure unity standard deviation and zero mean value
y    = y - mean(y);
yrms = sqrt(mean(y.^2));
y    = y / yrms;

SampleEntropy

function SampEn = SampleEntropy(x,m,r,sflag)
% check inputs
SampEn = [];
if nargin < 1 || isempty(x), fprintf('Error: x needs to be defined!
'), return; end
if ~isvector(x), fprintf('Error: x needs to be a vector!
'), return; end
if nargin < 2 || isempty(m), m = 2; end
if nargin < 3 || isempty(r), r = 0.2; end
if nargin < 4 || isempty(sflag), sflag = 1; end

% force x to be column vector
x = x(:);
N = length(x);

 % normalize/standardize x vector
if sflag > 0, x = (x - mean(x)) / std(x); end  

% obtain subsequences of the signal
Xam = zeros(N-m,m+1); Xbm = zeros(N-m,m);
for ii = 1 : N-m                   % although for N-m+1 subsequences could be extracted for m,
	Xam(ii,:) = x(ii:ii+m);          % in the Richman approach only N-m are considered as this gives the same length for m and m+1
	Xbm(ii,:) = x(ii:ii+m-1);
end

% obtain number of matches
B = zeros(1,N-m); A = zeros(1,N-m);
for ii = 1 : N-m
	% number of matches for m
	d = abs(bsxfun(@minus,Xbm((1:N-m)~=ii,:),Xbm(ii,:)));
 	B(ii) = sum(max(d,[],2) <= r);
 	
	% number of matches for m+1
	d = abs(bsxfun(@minus,Xam((1:N-m)~=ii,:),Xam(ii,:)));
 	A(ii) = sum(max(d,[],2) <= r);
end

% get probablities for two sequences to match
B  = 1/(N-m-1) * B;                  % mean number of matches for each subsequence for m
Bm = 1/(N-m) * sum(B);               % probability that two sequences will match for m points (mean of matches across subsequences)
A  = 1/(N-m-1) * A;                  % mean number of matches for each subsequence for m+1
Am = 1/(N-m) * sum(A);               % probability that two sequences will match for m+1 points (mean of matches across subsequences)

cp = Am/Bm;                          % conditional probability
SampEn = -log(cp);                   % sample entropy

spectra

function [F xf] = spectra(X, Sf, nfft, win)

% check inputs
if nargin < 2, error('Error: Wrong number of inputs!
'); end;
if nargin < 3 || isempty(nfft), nfft = size(X,1); end;
if nargin < 4 || isempty(win),  win = @hann; end;
if nfft < size(X,1); fprintf('Info: nfft < size(X,1)! Using nfft = size(X,1).
'); nfft = size(X,1); end;

% reshape multidimensional data into 2-d matrix
siz = size(X);
n = siz(1); nrest = siz(2:end);
X = reshape(X,[n prod(nrest)]);

% get window function and multiply it with the data
xwin = window(win,n);
X = bsxfun(@times,X,xwin);

% compute FFT, F - complex numbers, fourier matrix
F = fft(X,nfft,1)*2/sum(xwin);

% get frequency axis
xf = (0:nfft/2)*(Sf/nfft);

% shorten matrix by half
F = F(1:length(xf),:);

% put data back into multi-d matrix
F = reshape(F,[size(F,1) nrest]);

main

clear all;
format compact

rand('seed',10)

% parameters
Sf  = 1000;
dur = 30;

% colors and b parameter for 1/f^b
ccol = [0 151 204; 0 0 0; 204 8 118; 204 122 0] / 255;
b    = [-1 0 1 2];
lab  = {'Blue' 'White' 'Pink' 'Brown'};

% get signals, mse, spectrum
nsamps = round(10.^linspace(log10(50),log10(2000),60));
[y mse mspe mswpe SampEn pe wpe alpha] = deal([]);
for ii = 1 : length(b)
	y(:,ii)        = colored_noise(Sf,dur,b(ii));
	[mse(:,ii) sf] = MSE_Costa2005(y(:,ii),20,2,std(y(:,ii))*0.15);
	SampEn(ii)     = SampleEntropy(y(:,ii),2,std(y(:,ii))*0.15,0);
end

% time vector
t = (0:length(y)-1)/Sf;

% spectrum
[F xf] = spectra(y, Sf, Sf*10, @hann);

% plot spectra
figure, set(gcf,'Color',[1 1 1]), hold on
for ii = 1 : length(b)
	subplot(1,length(b),ii)
	set(gca,'FontSize',12)
	plot(xf,abs(F(:,ii)),'-','Color',ccol(ii,:))
	xlim([0 400])
	ylim([0 0.2])
	xlabel('Frequency (Hz)')
	if ii == 1
		ylabel('Amplitude (a.u.)')
	end
	title([lab{ii} ' noise (1/f^{' num2str(b(ii)) '})'])
end


% plot time course
figure, set(gcf,'Color',[1 1 1]), hold on
for ii = 1 : length(b)
	subplot(length(b),1,ii)
	set(gca,'FontSize',12)
	plot(t,y(:,ii),'-','Color',ccol(ii,:))
	ylim([-3 3])
	if ii == length(b)
		xlabel('Time (s)')
	end
	ylabel('Amp. (a.u.)')
	title([lab{ii} ' noise (1/f^{' num2str(b(ii)) '})'])
end


% plot mse
figure, set(gcf,'Color',[1 1 1]), hold on, pp = []; labfull = {};
set(gca,'FontSize',12)
for ii = 1 : length(b)
	pp(ii) = plot(sf,mse(:,ii),'-','Color',ccol(ii,:),'LineWidth',2);
	labfull{ii} = [lab{ii} ' noise (1/f^{' num2str(b(ii)) '})'];
end
ylim([0 3])
xlabel('Scale')
ylabel('SampEn')
title('Multi-scale entropy')
legend(pp,labfull)
legend('boxoff')


% plot sample entropy
figure, set(gcf,'Color',[1 1 1]), hold on, pp = []; labfull = {};
set(gca,'FontSize',12)
for ii = 1 : length(b)
	rectangle('Position',[ii 0 1 SampEn(ii)],'FaceColor',ccol(ii,:))
end
set(gca,'XTick',1.5:length(b)+0.5,'XTickLabel',lab)
xlim([0.8 length(b)+1.2])
ylim([0 3])
xlabel('Noise type')
ylabel('SampEn')
title('Sample entropy')