program.m

用于心电信号的处理

%212 Format:360Fre -> 540(1 Channal)|1080(2 Channal).
%%%%------------->心电信号的读取!
%readheader('E:\MATLAB\develop\100.hea')
clear all;
NN={100:350 200:250 100:300 100:150};           %各层选取的只含噪声点
fid1 = fopen('100.dat','rb');
[A] = fread(fid1,1080*4,'uchar');
j=1;
for i=1:360*4
   %212->12bits.
   B(i)=bitand(A(j+1),15)*2^8+A(j);             %取低4位，左移8位，成12位
   C(i)=bitand(A(j+1),240)*2^4+A(j+2);          %取高4位，左移4位，成12位
   j=j+3;
   %Signed.
   if B(i)>=2048
      B(i)=(4096-B(i))*(-1);
   end   
   
   if C(i)>=2048
      C(i)=(4096-C(i))*(-1);
   end      
end   
figure;
subplot(211);plot(B);
subplot(212);plot(C);

fid2 = fopen('100A.txt','wt');                  %文件写入
fprintf(fid2,'%d\n',B);
fid3 = fopen('100B.txt','wt');
fprintf(fid3,'%d\n',C);

fclose(fid1);
fclose(fid2);
fclose(fid3);

%%
signal=B;
level=5;wf='bior1.5';                           %滤波选用的小波函数

%%%%%%%%%%%%%%%%%%%%%各种滤波方法
%-------->1.软阈值滤波
S_thr1=wden(signal,'rigrsure','s','mln',5,'sym8'); %规则:rigrsure,尺度改变比例:mln
figure;
subplot(211);
plot(B);axis tight;grid on;axis([0 1500 800 1300]);
title('原始信号');
subplot(212);
plot(S_thr1);axis tight;grid on;axis([0 1500 800 1300]);
title('软阈值滤波后波形');

%-------->2.硬阈值滤波
S_thr2=wden(signal,'rigrsure','h','mln',5,'sym8'); %规则:rigrsure,尺度改变比例:mln
figure;
subplot(211);
plot(B);axis tight;grid on;axis([0 1500 800 1300]);
title('原始信号');
subplot(212);
plot(S_thr2);axis tight;grid on;axis([0 1500 800 1300]);
title('硬阈值滤波后波形');

%-------->3.FIR滤波器进行滤波(切比雪夫Ⅱ型 低通滤波器)
in=B;
fs=360;                                         %采样频率
fn=fs/2;                                        %奈奎斯特采样率
Wp=10/fn;Ws=100/fn;                             %通带截止频率及阻带截止频率
Rp=3;As=60;                                     %通带最大衰减Rp=3 dB-----阻带最小衰减As=60
[n,Wn]=cheb2ord(Wp,Ws, Rp,As)   
[b, a]=cheby2(n,Rp,Wn,'low');
%%绘制信号
figure;
subplot(3, 2, [1 2]);
plot(in); 
grid on;
ylabel('幅度');
title('输入信号');
out1=filter(b, a, in);
subplot(3, 2, [5 6]);
plot(out1); 
grid on;axis([0 1500 800 1300]);
ylabel('幅度');
title('FIR低通滤波后输出信号');

%-------->4.IIR滤波器进行滤波(巴特沃斯 低通数字滤波器)
in=B;
fp=40;fs=310;Fs=1000;Rp=3;Rs=60;T=1/Fs;         %指标设计
Wlp=2*tan(2*pi*fp*T/2)/pi;Wls=2*tan(2*pi*fs*T/2)/pi;    %求归一化频率
[N,Wn]=buttord(Wlp,Wls,Rp,Rs);                  %确定最小阶数N和频率参数Wn
[bp,ap]=butter(N,1,'s');                        %得归一化低通原型
[bs,as]=lp2lp(bp,ap,Wn*pi*Fs);                  %转换为模拟低通
[b,a]=bilinear(bs,as,Fs);                       %双线性法进行模数转换
out2=filter(b,a,in);                             %滤波
                %%绘制图形
figure;
subplot(3, 2, [1 2]);
plot(in); 
grid on;
ylabel('幅度');
title('输入信号');
subplot(3, 2, [5 6]);
plot(out2); 
grid on;axis([0 1500 800 1300]);
ylabel('幅度');
title('IIR低通滤波后输出信号');

%-------->5.空域相关法滤波

[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);             %设计指定滤波器
[swa,swd] = swt(signal,level,Lo_D,Hi_D);        %离散平稳小波变换
Swd_n=swd;                                      %小波系数的处理
swd_org=swd;
mask_n=zeros(size(Swd_n));                      %先把系数处理矩阵设置为全0。

for j=1:(level-1) 
                                                %在1:(level-1)分解层次上对高频系数处理，最后一层无法求相关系数，所以不作处理。
    Noise_d1=swd_org(j,:);
    Noise_d1=Noise_d1(NN{j});
    Noise_var=var(Noise_d1);                    %以信号的前只含有噪声的点估计噪声在各层的方差，个数为NN数组中定义
    Pw_var=var(swd_org(j,:));
    Corr=swd_org(j,:).*swd_org(j+1,:);          %定义相关系数为相邻两层的乘积。
      
    cc=1.2;                                     %===》用以设定停止迭代的噪声能量阈值，需要根据情况调节
    while Pw_var>cc*Noise_var
    Pw=sum(abs(swd(j,:)).^2);                   %计算小波能量
    Pcorr=sum(abs(Corr).^2);                    %计算相关系数能量
    Corr_new=Corr.*((Pw/Pcorr)^0.5);            %归一化
     
    corr_mod=abs(Corr_new);                     %取模
    w_mod=abs(swd(j,:));
    swd_n=swd(j,:).*(corr_mod>w_mod);           %corr_mod>w_mod成立，则认为相应的点为边缘
    swd_n1=(swd_n~=0);
    mask_n(j,:)=mask_n(j,:)+swd_n1;             %将选出的点赋给系数处理矩阵相应位置。
    swd_n0=ones(size(swd_n1));
    swd_n0=swd_n0-swd_n1;
    
    swd(j,:)=swd(j,:).*swd_n0;                  %将高频系数选出大值后的地方置0。
    Pw_var=var(swd(j,:));
    Corr_new=Corr_new.*swd_n0;                  %将相关系数选出大值后的地方置0。
    Corr=Corr_new;
    end
end  
    
mask_max=ones(1,length(mask_n));
mask_n=[mask_n((1:(level-1)),:);mask_max];      %最后一层系数处理矩阵全置1。
Swd_reg=swd_org.*mask_n;

signal_n=iswt(swa,Swd_reg,wf);
xcrr=signal_n-B;                                %求滤波误差信号。


figure;                                         %信号滤波前后比较。
subplot(3,1,1);
plot(signal); axis tight;grid on;
title('原始信号');
subplot(3,1,2);
plot(signal_n); axis tight;grid on;
title('空域法滤波后信号');
subplot(3,1,3);
plot(xcrr); axis tight;grid on;
title('滤波误差信号');

figure;                                         %空域法滤波与其他滤波方法的比较。
subplot(6,1,1);
plot(signal); axis tight;grid on;axis([0 1500 900 1230]);
title('原始信号');
subplot(612);
plot(S_thr1);axis tight;grid on;axis([0 1500 900 1230]);
title('软阈值滤波');
subplot(6,1,3);
plot(S_thr2); axis tight;grid on;axis([0 1500 900 1230]);
title('硬阈值滤波');
subplot(614);
plot(out1);axis tight;grid on;axis([0 1500 900 1230]);
title('FIR滤波');
subplot(615);
plot(out2);axis tight;grid on;axis([0 1500 900 1230]);
title('IIR滤波');
subplot(616);
plot(signal_n);axis tight;grid on;axis([0 1500 900 1230]);
title('空域相关滤波');

xzb1=10*log10(var(S_thr1)/var(signal-S_thr1));
xzb2=10*log10(var(S_thr2)/var(signal-S_thr2));
xzb3=10*log10(var(out1)/var(signal-out1));
xzb4=10*log10(var(out2)/var(signal-out2));
xzb5=10*log10(var(signal_n)/var(signal-signal_n));
fprintf('软阈值滤波信噪比为：%f',xzb1)
fprintf('硬阈值滤波信噪比为：%f',xzb2)
fprintf('FIR滤波信噪比为：%f',xzb3)
fprintf('IIR滤波信噪比为：%f',xzb4)
fprintf('空域相关滤波信噪比为：%f',xzb5)