1 简介
睡眠状况是评价人体健康状态的重要指标。本文提出一种基于枕下式的无扰睡眠监测系统，通过无扰获取的心率信号测算心率变异性（HRV），并结合隐马尔可夫模型（HMM），在对用户无扰无接触的环境下求解睡眠分期。针对现有 HMM 睡眠分期存在的问题，提出采用集合经验模态分解（EEMD）消除 HRV 个体差异导致的分期误差，再求解相应的睡眠分期。研究结果证明本文所提无扰式睡眠监测方案可实现 S1～S4 睡眠分期，正确率超过 60%，且性能优于现有 HMM 睡眠分期方案。
2 部分代码
function c=melcepst(s,fs,w,nc,p,n,inc,fl,fh)
%MELCEPST Calculate the mel cepstrum of a signal C=(S,FS,W,NC,P,N,INC,FL,FH)
%
%
% Simple use: c=melcepst(s,fs)% calculate mel cepstrum with 12 coefs, 256 sample frames
% c=melcepst(s,fs,'e0dD') % include log energy, 0th cepstral coef, delta and delta-delta coefs
%
% Inputs:
% s speech signal
% fs sample rate in Hz (default 11025)
% nc number of cepstral coefficients excluding 0'th coefficient (default 12)
% n length of frame (default power of 2 <30 ms))
% p number of filters in filterbank (default floor(3*log(fs)) )
% inc frame increment (default n/2)
% fl low end of the lowest filter as a fraction of fs (default = 0)
% fh high end of highest filter as a fraction of fs (default = 0.5)
%
%w any sensible combination of the following:
%
%'R' rectangular window in time domain
%'N'Hanning window in time domain
%'M'Hamming window in time domain (default)
%
% 't' triangular shaped filters in mel domain (default)
% 'n' hanning shaped filters in mel domain
% 'm' hamming shaped filters in mel domain
%
%'p'filters act in the power domain
%'a'filters act in the absolute magnitude domain (default)
%
% '0' include 0'th order cepstral coefficient
%'e' include log energy
%'d'include delta coefficients (dc/dt)
%'D'include delta-delta coefficients (d^2c/dt^2)
%
% 'z' highest and lowest filters taper down to zero (default)
% 'y' lowest filter remains at 1 down to 0 frequency and
% highest filter remains at 1 up to nyquist freqency
%
% If 'ty' or 'ny' is specified, the total power in the fft is preserved.
%
% Outputs:c mel cepstrum output: one frame per row
%
% Copyright (C) Mike Brookes 1997
%
% Last modified Thu Jun 15 09:14:48 2000
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 2 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 can obtain a copy of the GNU General Public License from
% ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to
% Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2 fs=11025; end
if nargin<3 w='M'; end
if nargin<4 nc=12; end
if nargin<5 p=floor(3*log(fs)); end
if nargin<6 n=pow2(floor(log2(0.03*fs))); end
if nargin<9
fh=0.5;
if nargin<8
fl=0;
if nargin<7
inc=floor(n/2);
end
end
end
if any(w=='R')
z=enframe(s,n,inc);
elseif any (w=='N')
z=enframe(s,hanning(n),inc);
else
z=enframe(s,hamming(n),inc);
end
f=rfft(z.');
[m,a,b]=melbankm(p,n,fs,fl,fh,w);
pw=f(a:b,:).*conj(f(a:b,:));
pth=max(pw(:))*1E-6;
if any(w=='p')
y=log(max(m*pw,pth));
else
ath=sqrt(pth);
y=log(max(m*abs(f(a:b,:)),ath));
end
c=rdct(y).';
nf=size(c,1);
nc=nc+1;
if p>nc
c(:,nc+1:end)=[];
elseif p<nc
c=[c zeros(nf,nc-p)];
end
if ~any(w=='0')
c(:,1)=[];
end
if any(w=='e')
c=[log(sum(pw)).' c];
end
% calculate derivative
if any(w=='D')
vf=(4:-1:-4)/60;
af=(1:-1:-1)/2;
ww=ones(5,1);
cx=[c(ww,:); c; c(nf*ww,:)];
vx=reshape(filter(vf,1,cx(:)),nf+10,nc);
vx(1:8,:)=[];
ax=reshape(filter(af,1,vx(:)),nf+2,nc);
ax(1:2,:)=[];
vx([1 nf+2],:)=[];
if any(w=='d')
c=[c vx ax];
else
c=[c ax];
end
elseif any(w=='d')
vf=(4:-1:-4)/60;
ww=ones(4,1);
cx=[c(ww,:); c; c(nf*ww,:)];
vx=reshape(filter(vf,1,cx(:)),nf+8,nc);
vx(1:8,:)=[];
c=[c vx];
end
if nargout<1
[nf,nc]=size(c);
t=((0:nf-1)*inc+(n-1)/2)/fs;
ci=(1:nc)-any(w=='0')-any(w=='e');
imh = imagesc(t,ci,c.');
axis('xy');
xlabel('Time (s)');
ylabel('Mel-cepstrum coefficient');
map = (0:63)'/63;
colormap([map map map]);
colorbar;
end
3 仿真结果
4 参考文献
[1]李翔等. "基于隐马尔可夫模型的枕下无扰式新型睡眠监测方案." 生物医学工程学杂志 35.2(2018):10.
部分理论引用网络文献，若有侵权联系博主删除。
5 MATLAB代码与数据下载地址
见博客主页
