基于GRNN网络和小波变换的ECG信号睡眠监测matlab仿真

目录

一、理论基础

二、核心程序

三、仿真测试结果

---

一、理论基础

1.1原理概述

      在这个问题中，我们将会讲述使用GRNN（广义回归神经网络）网络和小波变换来进行ECG（心电图）信号睡眠监测的原理和相关的数学公式。

       首先，理解ECG和睡眠监测的关系是非常重要的。ECG是一种记录心脏电活动的技术，可以反映心脏的状态。在睡眠期间，心脏的活动模式会发生改变，因此，我们可以通过分析ECG信号来推断睡眠状态。

       小波变换是一种信号处理工具，能够将信号分解成不同的频率成分。在ECG分析中，小波变换常用于去除噪声，突出有用的信号特征。

       而GRNN是一种用于回归问题的神经网络。它的核心思想是通过学习输入与输出之间的映射关系来进行预测。在ECG睡眠监测中，GRNN可以被训练来根据ECG信号预测睡眠阶段。

下面是基于GRNN和小波变换的ECG信号睡眠监测的基本步骤：

1.预处理：首先，对原始的ECG信号进行小波变换，以去除噪声并提取有用的特征。

2. 特征提取：从小波变换后的ECG信号中提取特征。这可能包括心率、心率变异性、呼吸率等。这些特征可以通过计算小波系数的统计量（例如均值、标准差、峰值等）来获得。

3.训练GRNN：使用提取的特征和对应的睡眠阶段标签来训练GRNN。GRNN的核心是一个回归模型，它可以学习输入特征（即ECG信号的特征）和输出（即睡眠阶段）之间的映射关系。这个映射关系可以通过以下公式表示：

y = f(w,x)

其中y是预测的睡眠阶段，f是GRNN的回归函数，w是权重参数，x是输入的特征。
4. 预测：在训练完成后，我们可以使用GRNN来根据新的ECG信号预测睡眠阶段。这可以通过将新的ECG信号通过相同的预处理和特征提取步骤处理，然后将提取的特征输入到GRNN中进行预测来完成。

       这就是基于GRNN网络和小波变换的ECG信号睡眠监测的基本原理和数学公式。请注意，实际应用中可能还需要进一步的优化和调整，例如通过交叉验证选择最佳的小波函数和GRNN参数，以及通过其他技术（例如数据增强）来处理数据不足的问题等。

1.2具体实施过程

1， 不同睡眠阶段会发出不同频率的波，所以最好用频谱的方法来提取出各个睡眠阶段中的不同频率的波如：α波，β波，θ波等等。

2，算法不需要太复杂，简便易懂最好。

3， 有关XML文件的作用我在这解释下，在matlab中用xmltest这个function可以得到所有XML的值，化成图就可以得到这个

这个数据是8小时睡眠的数据，上面这个XML文件图是对这8小时睡眠数据分期的结论，这个是老师给我们的正确结果，可以看到，前6小时这个人都处于清醒状态，所以是阶段0，六小时后开始睡眠阶段的变化，Y轴代表睡眠阶段，X轴是小时。

4，

      这个是我们原始的EEG信号，这个项目就是把这段原始信号通过算法来判断出每个时刻都处于什么睡眠阶段，然后得出结果后跟上面的XML文件的标准分期做对比来评判这种方法的准确率。

5， 打开EDF文件后会有21个测量人体电信号的传感器所测得的数据，但是我们只用其中四个A4，A2，C3，C1。上面这个原始EEG信号就是用A2-C3的差值所画的图，另外再做一个A4-C1的图就行。

我们的算法方法如下所示：

首先对原始的信号通过小波变换提取

然后小波变换做8层小波变换，提取第二层，第三层，第7层计算对应的能量，并做归一化处理。得到S1，S2，S3两个特征值

然后将 S1，S2，S3这五个变量做为特征提取变量。

1.1信号的小波分解与重构原理

 

然后通过8层次的小波分解，具体代码如下所示：

然后小波分解之后每一层的信号对应着具体不同的脑电波，具体看我这里的程序注释。

1.2GRNN网络

        广义回归神经网络(Generalized regression neural network, GRNN)是一种建立在非参数核回归基础之上的神经网络，通过观测样本计算自变量和因变量之间的概率密度函数。GRNN结构如图1所示，整个网络包括四层神经元：输入层、模式层、求和层与输出层。

       GRNN神经网络的性能，主要通过对其隐回归单元的核函数的光滑因子来设置的，不同的光滑因子可获得不同的网络性能。 

二、核心程序

function [x_delta,x_theta,x_apha,x_beta,S2,S3,S5,S6,S7]=wavelete_energy_decompose(x)
%小波分解
[c,l]=wavedec(x,8,'db3');
%下面的程序是对信号用小波进行分解。
ca8=wrcoef('a',c,l,'db3',8);%心电、眼电干扰
cd8=wrcoef('d',c,l,'db3',8);%0.41~0.82   delta波
cd7=wrcoef('d',c,l,'db3',7);%0.82~1.64   delta波
cd6=wrcoef('d',c,l,'db3',6);%1.64~3.28   delta波
cd5=wrcoef('d',c,l,'db3',5);%3.28~6.56   theta波
cd4=wrcoef('d',c,l,'db3',4);%6.56~13.13  apha波
cd3=wrcoef('d',c,l,'db3',3);%13.13~26.25 beta波
cd2=wrcoef('d',c,l,'db3',2);%26.25~52.25 高频噪声
cd1=wrcoef('d',c,l,'db3',1);%52.25~105   高频噪声
%下面的程序时对小波分解的近似信号的系数提取
 
%%%对delta系数求绝对和
%x_delta=sum(abs(D8)^2)+sum(abs(D7)^2)+sum(abs(D6)^2);
x_delta=sum(abs((cd8)).^2)+sum(abs((cd7)).^2)+sum(abs((cd6)).^2);
%%%theta系数
%x_theta=sum(abs(D5)^2);
x_theta=sum(abs((cd5)).^2);
%%%apha系数
%x_apha=sum(abs(D4)^2);
x_apha=sum(abs((cd4)).^2);
%%%beta系数
%x_beta=sum(abs(D3)^2);
x_beta=sum(abs((cd3)).^2);
%%%能量归一化
energy_sum=x_delta+x_theta+x_apha+x_beta;
x_delta=x_delta/energy_sum;
x_theta=x_theta/energy_sum;
x_apha=x_apha/energy_sum;
x_beta=x_beta/energy_sum;



S2 = sum(abs((cd2)).^2);
S3 = sum(abs((cd3)).^2);
S5 = sum(abs((cd5)).^2);
S6 = sum(abs((cd6)).^2);
S7 = sum(abs((cd7)).^2);

S  = S2+S3+S5+S6+S7;

S2 = S2/S;
S3 = S3/S;
S5 = S5/S;
S6 = S6/S;
S7 = S7/S;

function theStruct = parseXML(filename)
% PARSEXML Convert XML file to a MATLAB structure.
try
   tree = xmlread(filename);
catch
   error('Failed to read XML file %s.',filename);
end

% Recurse over child nodes. This could run into problems 
% with very deeply nested trees.
try
   theStruct = parseChildNodes(tree);
catch
   error('Unable to parse XML file %s.',filename);
end


% ----- Local function PARSECHILDNODES -----
function children = parseChildNodes(theNode)
% Recurse over node children.
children = [];
if theNode.hasChildNodes
   childNodes = theNode.getChildNodes;
   numChildNodes = childNodes.getLength;
   allocCell = cell(1, numChildNodes);

   children = struct(             ...
      'Name', allocCell, 'Attributes', allocCell,    ...
      'Data', allocCell, 'Children', allocCell);

    for count = 1:numChildNodes
        theChild = childNodes.item(count-1);
        children(count) = makeStructFromNode(theChild);
    end
end

% ----- Local function MAKESTRUCTFROMNODE -----
function nodeStruct = makeStructFromNode(theNode)
% Create structure of node info.

nodeStruct = struct(                        ...
   'Name', char(theNode.getNodeName),       ...
   'Attributes', parseAttributes(theNode),  ...
   'Data', '',                              ...
   'Children', parseChildNodes(theNode));

if any(strcmp(methods(theNode), 'getData'))
   nodeStruct.Data = char(theNode.getData); 
else
   nodeStruct.Data = '';
end

% ----- Local function PARSEATTRIBUTES -----
function attributes = parseAttributes(theNode)
% Create attributes structure.

attributes = [];
if theNode.hasAttributes
   theAttributes = theNode.getAttributes;
   numAttributes = theAttributes.getLength;
   allocCell = cell(1, numAttributes);
   attributes = struct('Name', allocCell, 'Value', ...
                       allocCell);

   for count = 1:numAttributes
      attrib = theAttributes.item(count-1);
      attributes(count).Name = char(attrib.getName);
      attributes(count).Value = char(attrib.getValue);
   end
end

clc;
clear;
close all;
warning off;
addpath 'func\'

%读取数据选择
SEL = 2;%设置1，重新读取文件，设置2，调用已经读取的数据
 

if SEL == 1
   %读取XML文件
   File_Name = 'mros-visit1-aa0009-profusion.xml';
   value     = func_read_xml(File_Name);
   figure;
   subplot(211);
   plot(value);
   %读取edf数据
   File_Name           = 'mros-visit1-aa0009.edf';
   [sensorC3,sensorC4] = func_edf_read(File_Name);%3为已知分类的信号，用来进行数据的训练
   subplot(212);
   plot(sensorC3);
   save edf_data.mat value sensorC3 sensorC4
else
   load edf_data.mat 
   figure;
   subplot(211);
   plot(value);
   subplot(212);
   plot(sensorC3);
   axis([0,length(sensorC3),-600,600]);
end




hd        = fir1(63,0.3,'low');
sensorC3s = conv(sensorC3,hd);
sensorC4s = conv(sensorC4,hd);

%信号滤波 
sensorC3_filter = sensorC3s;     
        
%首先对已知的数据进行特征分析和提取，然后进行训练
NUM    = length(sensorC3_filter)/length(value);
for i = 1:length(value)-1
    T(i,1)    = value(i);
    data(i,:) = sensorC3_filter(NUM*(i-1)+1:NUM*i);
end
%信号特征提取
Feature = func_get_feature(data);



%睡眠大周期分割
%进行神经网络的训练
net_first    = newgrnn(Feature',T',0.1);
net_second   = newgrnn(Feature',T',0.016);
%然后对测试数据进行测试
%信号滤波 
sensorC4_filter = sensorC4s;    
 
NUM = length(sensorC4_filter)/length(value);
for i = 1:length(value)
    data2(i,:) = sensorC4_filter(NUM*(i-1)+1:NUM*i);
end
%信号特征提取
Feature2 = func_get_feature(data2);
T2       = sim(net_first,Feature2');
T2s      = func_smooth(T2,64);
T3       = round(T2s(1:length(sim(net_second,Feature2'))).*sim(net_second,Feature2'));

 
figure;
subplot(211);
plot(T3);title('测试信号识别结果');
subplot(212);
plot(sensorC4_filter);title('测试信号结果');
axis([0,length(sensorC4_filter),-600,600]);

三、仿真测试结果

A28-46