眨眼信号分离及MATLAB实现

脑电中眨眼信号的分离原理：
通常采集到的脑电信号中会混有噪声和伪迹，而眨眼引起的肌肉电信号（以后简称眨眼信号）就是伪迹中的一种。眨眼电信号的频率分布较广泛，且幅值较大，所以采用了基于小波变换的多分辨率分析的方法进行眨眼信号的分离。
来说说我个人对小波变换多分辨率分析的理解。关于分辨率，最直观的的就是图像，所以先拿图片来举个例子：假设用100万像素的相机去拍照，得到一张图片A，再用50万像素的相机拍同一个地方，得到图片B，A与B内容相同，只是图片看上去的清晰程度不同。要是把B的细节在美化一下，就可以得到A。设A与B 之间的差异（细节）D=A-B，也就是说低分辨率下B加上细节部分D，就能得到高分辨率A，以此类推，B可以由更低分辨率的图片C加上C与B的细节部分D就能得到B，所以，最高分辨率A=D1+D2+D3+。。。+DN+A（n），其中A（n）为分辨率最低的图片，D1~DN为各图片与比他高级一的分辨率图片的细节部分。
同理，对于脑电信号分辨率越高相当于采样点越多，采样就越精确，信号越接近真实信号，分辨率越低采样点越少，图像越光滑。由于脑电信号是微弱的非线性的，所以包含在高分辨率下，而眨眼信号相对幅值较大，在低分辨率下也能分辨，所以只要提取出低分辨率下的信号，就能得到眨眼信号而去除脑电信号。具体代码如下：

function varargout = eeggui(varargin)
% EEGGUI MATLAB code for eeggui.fig
%      EEGGUI, by itself, creates a new EEGGUI or raises the existing
%      singleton*.
%
%      H = EEGGUI returns the handle to a new EEGGUI or the handle to
%      the existing singleton*.
%
%      EEGGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in EEGGUI.M with the given input arguments.
%
%      EEGGUI('Property','Value',...) creates a new EEGGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before eeggui_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to eeggui_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help eeggui

% Last Modified by GUIDE v2.5 29-Mar-2017 14:10:23

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @eeggui_OpeningFcn, ...
                   'gui_OutputFcn',  @eeggui_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before eeggui is made visible.
function eeggui_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to eeggui (see VARARGIN)

% Choose default command line output for eeggui
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes eeggui wait for user response (see UIRESUME)
% uiwait(handles.eeggui);


% --- Outputs from this function are returned to the command line.
function varargout = eeggui_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;



function ed_message_Callback(hObject, eventdata, handles)
% hObject    handle to ed_message (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ed_message as text
%        str2double(get(hObject,'String')) returns contents of ed_message as a double


% --- Executes during object creation, after setting all properties.
function ed_message_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ed_message (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in pb_readdata.
function pb_readdata_Callback(hObject, eventdata, handles)
% hObject    handle to pb_readdata (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global h10;
global l10;
global c10;
load('G:\桌面\脑电数据\c.mat');%脑电信号的校验和
load('G:\桌面\脑电数据\h.mat');%离散脑电信号高8位
load('G:\桌面\脑电数据\l.mat');%离散脑电信号低8位
for i=1:length(h)
     h10(i)=hex2dec(h(i));
end
for i=1:length(l)
    l10(i)=hex2dec(l(i));
end
for i=1:length(c)
    c10(i)=hex2dec(c(i));
end
str=['读取数据成功！']
set(handles.ed_message,'string',str);
% --- Executes on button press in pb_checksum.
function pb_checksum_Callback(hObject, eventdata, handles)
% hObject    handle to pb_checksum (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global h10;
global l10;
global c10;
global rawdata;
 for i=1:length(h10)
    sum=l10(i)+h10(i)+128+2;
    sum=bitxor(sum,4294967295);
    sumcheck=bitand(sum,255);
    if(sumcheck==c10(i))
        rawdata(i)=bitshift(h10(i),8)+l10(i);
    else
        rawdata(i)=0;
    end
 end
for i=1:length(rawdata)
      if(rawdata(i)>32768)
          rawdata(i)=rawdata(i)-65536;
      end
end
  str=['校验数据成功！']
set(handles.ed_message,'string',str);
plot(rawdata)
% --- Executes on button press in pb_wavelet.
function pb_wavelet_Callback(hObject, eventdata, handles)
% hObject    handle to pb_wavelet (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global aa4;
global rawdata;
global da;
global db;
global dc;
global dd;
global de;
global af;
global cc;
global ll;
[aa4,da,db,dc,dd,de,af,cc,ll]=waveletcontrol(rawdata);
str=['得到小波数据：各层的细节系数和平均系数']
set(handles.ed_message,'string',str);
%plot(eye);
% --- Executes on button press in pb_noise.
function pb_noise_Callback(hObject, eventdata, handles)
% hObject    handle to pb_noise (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global rawdata;
global da;
global db;
global dc;
global dd;
global de;
global af;
global cc;
global ll;
global aa;
dd1=zeros(size(da));
dd2=zeros(size(db));
dd3=zeros(size(dc));
dd4=zeros(size(dd));
dd5=zeros(size(de));
c=[af dd5 dd4 dd3 dd2 dd1];
aa=waverec(c,ll,'db4');%重构脑电信号
plot(aa);
str=['成功去除脑电噪声，但是包含眨眼信号！']
set(handles.ed_message,'string',str);
% --- Executes on button press in pb_eye.
function pb_eye_Callback(hObject, eventdata, handles)
% hObject    handle to pb_eye (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global aa4;
global y;
plot(aa4);
figure;
y=wthresh(aa4,'s',350);
plot(y);
% --- Executes on button press in pb_quchueye.
function pb_quchueye_Callback(hObject, eventdata, handles)
% hObject    handle to pb_quchueye (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global aa4;
global rawdata;
global y;
global aa;
for i=1:length(y)
   aaa4(i)=aa(i)-aa4(i);
end
plot(aaa4);
str=['去除眨眼信号！']
set(handles.ed_message,'string',str);

function varargout = waveletcontrol(varargin)
% WAVELETCONTROL MATLAB code for waveletcontrol.fig
%      WAVELETCONTROL, by itself, creates a new WAVELETCONTROL or raises the existing
%      singleton*.
%
%      H = WAVELETCONTROL returns the handle to a new WAVELETCONTROL or the handle to
%      the existing singleton*.
%
%      WAVELETCONTROL('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in WAVELETCONTROL.M with the given input arguments.
%
%      WAVELETCONTROL('Property','Value',...) creates a new WAVELETCONTROL or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before waveletcontrol_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to waveletcontrol_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help waveletcontrol

% Last Modified by GUIDE v2.5 29-Mar-2017 14:11:01

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @waveletcontrol_OpeningFcn, ...
                   'gui_OutputFcn',  @waveletcontrol_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before waveletcontrol is made visible.
function waveletcontrol_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to waveletcontrol (see VARARGIN)

% Choose default command line output for waveletcontrol
global rd;
handles.output = hObject;
rd=varargin{1};
% Update handles structure
guidata(hObject, handles);
uiwait(handles.waveletcontrol);
% UIWAIT makes waveletcontrol wait for user response (see UIRESUME)
% uiwait(handles.waveletcontrol);


% --- Outputs from this function are returned to the command line.
function varargout = waveletcontrol_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.aa4;
varargout{2} = handles.d1;
varargout{3} = handles.d2;
varargout{4} = handles.d3;
varargout{5} = handles.d4;
varargout{6} = handles.d5;
varargout{7} = handles.a5;
varargout{8} = handles.C;
varargout{9} = handles.L;
delete(handles.waveletcontrol);


function ed_mess_Callback(hObject, eventdata, handles)
% hObject    handle to ed_mess (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ed_mess as text
%        str2double(get(hObject,'String')) returns contents of ed_mess as a double


% --- Executes during object creation, after setting all properties.
function ed_mess_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ed_mess (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on selection change in pp_wavelet.
function pp_wavelet_Callback(hObject, eventdata, handles)
% hObject    handle to pp_wavelet (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns pp_wavelet contents as cell array
%        contents{get(hObject,'Value')} returns selected item from pp_wavelet


% --- Executes during object creation, after setting all properties.
function pp_wavelet_CreateFcn(hObject, eventdata, handles)
% hObject    handle to pp_wavelet (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on selection change in pp_N.
function pp_N_Callback(hObject, eventdata, handles)
% hObject    handle to pp_N (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns pp_N contents as cell array
%        contents{get(hObject,'Value')} returns selected item from pp_N


% --- Executes during object creation, after setting all properties.
function pp_N_CreateFcn(hObject, eventdata, handles)
% hObject    handle to pp_N (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in pb_start.
function pb_start_Callback(hObject, eventdata, handles)
% hObject    handle to pb_start (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%handles = guidata(hObject);
global rd;
global wave;
global N;
global aa4
% global a1;
% global d1;
[C,L]=wavedec(rd,5,'db4'); %db4小波5层分解
%重构各层的平均和细节
  a5 = appcoef(C,L,'db4',5);
  d5 = detcoef(C,L,5);
  a4 = appcoef(C,L,'db4',4);
  d4 = detcoef(C,L,4);
  a3 = appcoef(C,L,'db4',3);
  d3 = detcoef(C,L,3);
  a2 = appcoef(C,L,'db4',2);
  d2 = detcoef(C,L,2);
  a1 = appcoef(C,L,'db4',1);
  d1 = detcoef(C,L,1);

aa4 = wrcoef('a',C,L,'db4',4);%眨眼信号
handles.aa4=aa4;
handles.d1=d1;
handles.a1=a1;
handles.d2=d2;
handles.a2=a2;
handles.d3=d3;
handles.a3=a3;
handles.d4=d4;
handles.a4=a4;
handles.d5=d5;
handles.a5=a5;
handles.C=C;
handles.L=L;
guidata(hObject, handles);
axes(handles.axes1);
plot(a1);
axes(handles.axes2);
plot(a2);
axes(handles.axes3);
plot(a3);
axes(handles.axes4);
plot(a4)
axes(handles.axes5);
plot(a5)
str=['提取并重构各层的平均系数！']
set(handles.ed_mess,'string',str);
%plot(rd);
% --- Executes on button press in pb_end.
function pb_end_Callback(hObject, eventdata, handles)
% hObject    handle to pb_end (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
uiresume(handles.waveletcontrol);

这是使用了MATLAB的guide工具箱制作的一个带有软件界面的眨眼信号分离，结果如图


第一张图片为原始脑电信号，第二张左边是经过一定的去噪后的信号，右边是识别到的眨眼信号。