iir.m

来自「用matlab设计FIR和IIR滤波器的设计函数与例程」· M 代码 · 共 47 行

%实时检测QRS复波，采样率1000HZ
filename='D:\数据\蔡刿.mat';
load(filename);
ecgdata=data(:,2);
ecg=ecgdata(1:2000);

M=fs/50;
%------------lp,截至频率11HZ，延迟5个采样点，增益36
a=[1 zeros(1,M-1) -2 zeros(1,M-1) 1];
b=[1 -2 1];
ecgdata_denoise=filter(a,b,ecg);
 figure(1);subplot(2,1,1);plot(ecg(1:1000));title('orignal');
subplot(2,1,2);
plot(ecgdata_denoise(1:1000));title('低通滤波，去除60hz');
 figure(2);
 freqz(a,b,512,fs);
%------------hp,截至频率11HZ，延迟15.5个采样点，增益32
a=[ -1  zeros(1,15) 32 -32 zeros(1,15) 1];b=[32 -32 zeros(1,31)];
ecgdata_denoise1=filter(a,b,ecgdata_denoise);
 figure(3);
plot(ecgdata_denoise1(1:1000));title('高通滤波');grid on;
 figure(4);
 freqz(a,b,512,fs);
 %------------diff
a=[ 2 1 0 -1 -2];b=[10];
ecgdata_denoise2=filter(a,b,ecgdata_denoise1);
 figure(5);
 subplot(2,1,1);
 plot(ecgdata_denoise2(1:1000));title('diff');grid on;
 subplot(2,1,2);
 ecgdata_denoise3=ecgdata_denoise2.*ecgdata_denoise2;
plot(ecgdata_denoise3(1:1000));title('square');grid on;
 figure(6);
 freqz(a,b,512,fs);
 %-----------------------------sum 30p
 %-----------------------------sum 32p
 y=zeros(1,length(ecg)-30);
 for i=1:length(ecg)-30
     sum=0;
     for j=0:29
     sum=sum+ecgdata_denoise3(i+j);
     end
     y(i)=sum;
 end
 figure(7);
 plot(y(1:1000));grid on;