exa060603.m

心电信号加入0.5Hz的基线漂移

%----------------------------------------------------------------------------
% exa060603.m, for example 6.6.3;
% To design IIR Butteworth bandstop DF by analog-lowpass,
%通过模拟低通滤波器直接设计巴特沃斯阻带数字滤波器
% ---------------------------------------------------------------------------
clear all;
load data.txt;
data;
fp=[45 55];fs=[49 51];
%wp=[.19*pi 0.21*pi];ws=[.198*pi 0.202*pi];
Fs=1000;%采样频率
rp=3;rs=14;%通带的幅度误差和阻带内的最小衰减
wp=fp*2*pi/Fs;ws=fs*2*pi/Fs;
%
% Firstly to finish frequency prewarping;
wap=2*Fs*tan(wp./2);
was=2*Fs*tan(ws./2);
[n,wn]=buttord(wap,was,rp,rs,'s');
% Note: 's'!
[z,p,k]=buttap(n);% z p k分别表示滤波器的零点 极点 增益
[b,a]=zp2tf(z,p,k);
bw=wap(2)-wap(1);
w0=sqrt(wap(1)*wap(2));
[bt,at]=lp2bs(b,a,w0,bw);
%
% Note: z=(2/ts)(z-1)/(z+1);
[bz1,az1]=bilinear(bt,at,Fs);

%计算系统的幅频特性
[h,w]=freqz(bz1,az1,256,Fs);%[0 500]频率范围内选取256个频率点
y=filter(bz1,az1,data);
figure(1)
subplot(2,1,1)
stem(data);
subplot(2,1,2)
stem(y);
mag=abs(h);%计算幅度
ph=angle(h);%计算频率
ph=w*180/pi;
%figure(2)
%subplot(211)
%plot(w,mag);grid %幅度特性曲线
%xlabel('frequency(HZ)');
%ylabel('magnitude');
%subplot(212)
%plot(w,ph);grid %频率特性曲线
%xlabel('frequency(HZ)');
%ylabel('phaze');

% 利用MATLAB巴特沃斯函数直接设计数字滤波器
[n,wn]=buttord(wp/pi,ws/pi,rp,rs);
[b,a]=butter(n,wp/pi,'stop')%巴特沃斯阻带数字滤波器设计
[h1,w1]=freqz(b,a,256,Fs);
y=filter(b,a,data);
t=1:length(data);
data1=0.2*sin(25/pi*t)'+data;
subplot(311)
plot(data);grid on
axis([0 1000 -1 1 ])
title('干净的原信号')
subplot(312)
plot(data1);grid on 
axis([0 1000 -1 1 ])
title('加工频后的原信号')
subplot(313)
plot(y);grid on 
axis([0 1000 -1 1 ])
title('滤波后的信号')
mag1=abs(h1);%计算幅度
ph1=angle(h1);%计算频率
ph1=w1*180/pi;
figure(2)
subplot(211)
plot(w1,20*log10(abs(h1)))%幅度特性曲线
grid on
xlabel('frequency(HZ)');
ylabel('magnitude');
axis([ 0 500 -40 5])
title('系统幅度特性')
subplot(212)
plot(w1,ph1);grid %频率特性曲线
xlabel('frequency(HZ)');
ylabel('phaze');
title('系统频率特性')