060921bijaio3ceng.m

用于心电分析的小波分析程序

%-----------比较db(1-8)小波去噪优越性程序.---------------------------------
clc; clear all;
%----------从MIT-BIH库中读取心电异常第100组数据------------------------------
% This programm reads ECG data which are saved in format 212.

% (e.g., 100.dat from MIT-BIH-DB, cu01.dat from CU-DB,...)

% The data are displayed in a figure together with the
% annotations.(数据和说明在一起)

% The annotations are saved in the vector ANNOT, the corresponding

% times (in seconds) are saved in the vector ATRTIME.

% The annotations are saved as numbers, the meaning of the numbers can

% be found in the codetable "ecgcodes.h" available at www.physionet.org.

%

% ANNOT only contains the most important information, which is displayed

% with the program rdann (available on www.physionet.org) in the 3rd row.

% The 4th to 6th row are not saved in ANNOT.
%------ SPECIFY DATA ------------------------------------------------------

PATH= 'D:\matlab\matlab7.2\work'; % path, where data are saved

HEADERFILE= '100.hea';      % header-file in text format
                            %头文件.
ATRFILE= '100.atr';         % attributes-file in binary format
                            %说明文件.
DATAFILE='100.dat';         % data-file
                            %数据文件.
SAMPLES2READ=30000;         % number of samples to be read
%样本数.

                            % in case of more than one signal:

                            % 2*SAMPLES2READ samples are read



%------ LOAD HEADER DATA --------------------------------------------------

fprintf(1,'\\n$> WORKING ON %s ...\n', HEADERFILE);
%将格式化的文件写入字符串.

%signalh= fullfile(PATH, HEADERFILE);

%signalh= fullfile(HEADERFILE);

fid1=fopen('100.hea','r');
%'100.hea为打开的文件,'r'为打开的方式(为只读)
z= fgetl(fid1);
%从文件中读行,并丢掉换行符.

A= sscanf(z, '%*s %d %d %d',[1,3]);
%从格式化的字符串中读取(数据).(%*s 忽略,%d10进制)
nosig= A(1);  % number of signals(信号的数)

sfreq=A(2);   % sample rate of data(数的样品率).

clear A;

for k=1:nosig
%for循环(1=<k=<nosig).
    z= fgetl(fid1);
%从文件中读行,并丢掉换行符.
    A= sscanf(z, '%*s %d %d %d %d %d',[1,5]);
%从格式化的字符串中读取.
    dformat(k)= A(1);           % format; here only 212 is allowed
%[212 212]
    gain(k)= A(2);              % number of integers per mV
%[200 200]
    bitres(k)= A(3);            % bitresolution
%[11 11]
    zerovalue(k)= A(4);         % integer value of ECG zero point
%[1024 1024]
    firstvalue(k)= A(5);        % first integer value of signal (to test for errors)
%[995 1011]
end;
%结束循环.
fclose(fid1);
%关掉文件fid1.
clear A;
%清除A.

%------ LOAD BINARY DATA --------------------------------------------------

if dformat~= [212,212], error('this script does not apply binary formats different to 212.'); end;
%如果dformat~=[212 212]提示错误.
signald= fullfile(PATH, DATAFILE);            % data in format 212
%从部分文件中建立文件全名(或把DATAFILE文件放在PATH目录下).
fid2=fopen('100.dat','r');
%打开数据文件100.dat,'r'为只读.
A= fread(fid2, [3, SAMPLES2READ], 'uint8')';  % matrix with 3 rows, each 8 bits long, = 2*12bit
%读二进制文件fid2,生成3行,列为小于等于3000的矩阵.'uint8'无符号8字节型.(实际执行程序时生成为3000*3矩阵.)
fclose(fid2);
%关闭文件fid2.
M2H= bitshift(A(:,2), -4);
%把矩阵A第二列按位运算右移4位,即保留其高4位.
M1H= bitand(A(:,2), 15);
%把矩阵A第二列按位运算与15作与运算,即保留其低4位.
PRL=bitshift(bitand(A(:,2),8),9);     % sign-bit
%先做与运算保留了第4为,然后移9位,移到了第13位.
PRR=bitshift(bitand(A(:,2),128),5);   % sign-bit
%先做与运算保留了第8位,然后左移5位,移到了第13位.
M( : , 1)= bitshift(M1H,8)+ A(:,1)-PRL;
%把矩阵A第二列的低4位,左移8位,然后和其第1列求和,减去第13位(原第4位).
M( : , 2)= bitshift(M2H,8)+ A(:,3)-PRR;
%把矩阵A第二列的高4位,左移8位,然后和其第3列求和,减去第13位(原第8位).
if M(1,:) ~= firstvalue, error('inconsistency in the first bit values'); end;

switch nosig

case 2

    M( : , 1)= (M( : , 1)- zerovalue(1))/gain(1);

    M( : , 2)= (M( : , 2)- zerovalue(2))/gain(2);

    TIME=(0:(SAMPLES2READ-1))/sfreq;
%时间为采样点数除以样品率。
case 1

    M( : , 1)= (M( : , 1)- zerovalue(1));

    M( : , 2)= (M( : , 2)- zerovalue(1));

    M=M';

    M(1)=[];

    sM=size(M);

    sM=sM(2)+1;

    M(sM)=0;

    M=M';

    M=M/gain(1);

    TIME=(0:2*(SAMPLES2READ)-1)/sfreq;

otherwise  % this case did not appear up to now!

    % here M has to be sorted!!!

    disp('Sorting algorithm for more than 2 signals not programmed yet!');

end;

clear A M1H M2H PRR PRL;
%提取其中一个导联的数据进行去噪
%--------------------------小波分析心电信号---------------------------------
E=M(:,2);
E=E';
n=size(E);
s=E(1:2000);
%使用db(1-8)小波对E进行3层分解
for jj=1:8
    %jj表示小波的类型
    switch jj
        case 1 
            db='db1';
        case 2 
            db='db2';
        case 3 
            db='db3';
        case 4 
            db='db4';
        case 5 
            db='db5';
        case 6 
            db='db6';
        case 7 
            db='db7';
        case 8 
            db='db8';
        otherwise
    end;
[C L]=wavedec(E,3,db);

% 从c中提取尺度3下的近似小波系数
cA3=detcoef(C,L,3);
%从信号c中提取尺度1,2,3下的细节小波系数
[cD1,cD2,cD3]=detcoef(C,L,[1,2,3]);
%重构尺度3下的近似小波系数
A3=wrcoef('a',C,L,db,3);
% 重构尺度3,2,1下的细节小波系数. 

D1=wrcoef('d',C,L,db,1);
D2=wrcoef('d',C,L,db,2);
D3=wrcoef('d',C,L,db,3);

%-------------------使用阈值法去噪-----------------------------
%利用'ddencmp'得到除噪的默认参数
%thr,sorh,keepapp分别表示是阈值,软阈值(或硬阈值)和允许用户保存的低频系数.
%'den','wv',s分别为降噪('cmp'表示压缩),小波('wp'小波包),和降噪信号.
[thr,sorh,keepapp]=ddencmp('den','wv',E);
%去噪
clean=wdencmp('gbl',C,L,db,3,thr,sorh,keepapp);
%'gbl'为所有的级次都选择相同的阈值；'lvd'则对不同的级次选择不同的阈值。
%---------去噪效果衡量（SNR越大效果越好,MSE越小越好）--------------------------------
%选取信号的长度。
N=n(2);
x=E;%E为原始信号。
y=clean;%为去噪后的信号。
F=0;%SNR累加求和的初始值。
M=0;%MSE累加求和的初始值。
for ii=1:N
    %求原始信号与去噪后信号的差值平方。
   m(ii)=(x(ii)-y(ii))*(x(ii)-y(ii));
   %求去噪后信号的平方
   t(ii)=y(ii)*y(ii);
   f(ii)=t(ii)/m(ii);
   %累加求和
    F=F+f(ii);
    M=M+m(ii);
end;
SNR(jj)=10*log10(F);
MSE(jj)=M/N;
SM(jj)=SNR(jj)/MSE(jj);
end;
%------------------------比较小波函数db(1-8)优越性--------------------------
max_SNR = SNR(1);ts=1;
min_MSE=MSE(1);tm=1;
max_SM=SM(1);tsm=1;
%求SNR最大值
for ii=2:8
    if max_SNR==inf
        max_SNR=SNR(ii);
    end;
    if SNR==inf
        ii=ii+1;
    end;
    if max_SNR<=SNR(ii)
        max_SNR=SNR(ii);
        ts=ii;
    end;
end;
%求MSE最小值
for ii=2:8
    if min_MSE>=MSE(ii)
        min_MSE=MSE(ii);
        tm=ii;
    end;
end;
%求SM最大值
for ii=2:8
    if max_SM==inf
        max_SM=SM(ii);
    end;
    if SM(ii)==inf
        ii=ii+1;
    end;
    if max_SM<=SM(ii)
        max_SM=SM(ii);
        tsm=ii;
    end;
end;
disp('SNR的最大值为：');max_SNR
disp('最大SNR对应的db小波函数序号为：');ts
disp('MSE的最小值为：');min_MSE
disp('最小MSE对应的db小波函数序号为：');tm
disp('SM的最大值为：');max_SM
disp('最大SM对应的db小波函数序号为：');tsm
%以（即SM综合SNR,MSE的信息）为标准确定最佳的小波函数
switch tsm
        case 1 
            db='db1';
        case 2 
            db='db2';
        case 3 
            db='db3';
        case 4 
            db='db4';
        case 5 
            db='db5';
        case 6 
            db='db6';
        case 7 
            db='db7';
        case 8 
            db='db8';
        otherwise
    end;
%对比原始信号和除噪后的信号
subplot(2,1,1);
plot(s(1:1000));
title('原始信号')
%--------------------用运去噪效果最佳的小波函数除噪--------------------------
[C L]=wavedec(s,3,db);
% 从c中提取尺度3下的近似小波系数
cA3=detcoef(C,L,3);
%从信号c中提取尺度3，2，1下的细节小波系数
[cD1,cD2,cD3]=detcoef(C,L,[1,2,3]);
%重构尺度3下的近似小波系数
A3=wrcoef('a',C,L,db,3);
% 重构尺度3,2,1下的细节小波系数. 
D1=wrcoef('d',C,L,db,1);
D2=wrcoef('d',C,L,db,2);
D3=wrcoef('d',C,L,db,3);
[thr,sorh,keepapp]=ddencmp('den','wv',s);
clean=wdencmp('gbl',C,L,db,3,thr,sorh,keepapp);
subplot(2,1,2);
plot(clean(1:1000));
title('除噪后的信号')