clear all;
clc
wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;
wp=2*tan(wp/2);ws=2*tan(ws/2);
[N,wn]=buttord(wp,ws,Rp,As,'s');
[B,A]=butter(N,wn,'s'); %离散化处理
[num,den]=bilinear(B,A,1);%双线性
[h,w]=freqz(num,den);
figure(1)
plot(w/pi,20*log10(abs(h)/(max(h)))); %绘制滤波器幅频特性
grid on;
xlabel('w/Π');
ylabel('幅度(dB)');%显示网格线，给出横纵坐标标识
%原始心电图信号
xn1=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,...
-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,...
0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0,...
-2,-4,-2,0,-2,-4,-4,2,0,0,-2,-4,-2,0,0,-2,-4,-2,0,0,-4,-4,-2,-2,-4,-6,-6,-4,-4,8,-10,-8,...
-6,-6,-8,-12,-10,-8,-8,-10,-12,-10,-8,-8,-10,-10,-8,-6,-6,-8,-8,-4,-2,-4,-4,-4,...
0,0,-2,-4,-2,-2,0,-4];
n=0:length(xn1)-1;
figure(2);
plot(n,xn1);
title('原始心电信号时域波形'); grid on;
figure(3);
plot(n,abs(fft(xn1)));
title('原始心电信号频谱'); grid on;
y=filter(num,den,xn1);
figure(4)
plot(n,y);
title('IIR滤波器滤波后心电信号'); grid on;
figure(5);
plot(n,abs(fft(y)));
title('IIR滤波器滤波后信号心电频谱'); grid on;
%采用矩形窗设计
wp=0.2*pi;
ws=0.3*pi;
As=15;
wn=boxcar(18);%计算后N取18
N=18;
nn=0:N-1; alfa=(N-1)/2;
hd=sin(0.25*pi*(nn-alfa))./(pi*(nn-alfa));
h=hd.*wn';
[h1,w1]=freqz(h,1);
figure(6);
plot(w1/pi,20*log10(abs(h1)));
title('fir滤波器'); grid on;
xlabel('归一化频率/p');
ylabel('幅度/dB')
%
y1=filter(h,1,xn1);
figure(7);
plot(y1);
title('fir滤波器滤波后心电信号'); grid on;
figure(8);
plot(0:length(y1)-1,abs(fft(y1)));
title('fir滤波器滤波后信号心电频谱'); grid on;
