n=2300;
k=0:1/n:1;
x=randn(1,n+1)+5*sin(2*pi*50*k);
x1=randn(1,n+1);
w=linspace(0,pi,512);
X=abs(freqz(x,1,w));
figure:subplot(3,1,1);
plot(k,x);
title('滤波前的有干扰的肌电信号波形');
subplot(3,1,2);
plot(k,x1);
title('无干扰的原始肌电信号波形');
subplot(3,1,3);
plot(w/pi,X);
title('滤波前的幅频特性');
%滤出50HZ噪声的滤波器，采用IIR带阻滤波器
Fn=50;
Fs=2300;
W3=0.4;
W0=2*pi*Fn/Fs;
beta=cos(W0);
alpha=min(roots([1,-2/cos(W3),1]));
a=[1,-beta*(1+alpha),alpha];
b=[1,-2*beta,1]*(1+alpha)/2;
y1=filter(b,a,x);
figure:subplot(3,1,1);
plot(k,y1);
title('经带阻滤波后的波形图');
X=abs(freqz(y1,1,w));
subplot(3,1,2);
plot(w/pi,X);
title('经带阻滤波后的幅频特性');
subplot(3,1,3);
H1=abs(freqz(b,a));
plot(w/pi,H1);
title('带阻滤波器的幅频特性');
%采用10-1000HZ的带通滤波器
Wp1=10;
Wp2=1000;
Ws1=5;
Ws2=1100;
Rp=3;
Rs=30;
Fs=2300;
Wp=[Wp1 Wp2];
Ws=[Ws1 Ws2];
[n1,Wn1]=buttord(Wp/(Fs/2),Ws/(Fs/2),Rp,Rs);
[b1,a1]=butter(n1,Wn1);
y2=filter(b1,a1,y1);
figure:subplot(3,1,1);
plot(k,y2);
title('再经带通滤波后的波形图');
Y=abs(freqz(y2,1,w));
subplot(3,1,2);
plot(w/pi,Y);
title('再经带通滤波后的幅频特性');
subplot(3,1,3);
H2=abs(freqz(b1,a1));
plot(w/pi,H2);
title('带通滤波器的幅频特性');
