bandstop.m

自己做的一个滤波程序

clear;
w1=2*50*pi;
w2=2*49*pi;
w3=2*51*pi;
k=1.039;

%产生随机信号
n=1281;
a=rand(n,1);
m=0:1:n-1;
figure(1)
stem(m,a);
xlabel('时间序列n');
ylabel('振幅');
title('噪声');

%采样间隔时间,产生各频率信号
dt=1/6400;
t=0:dt:0.2;
c=100*sin(w1*t);
d=5*sin(k*w2*t);
e=2*sin(k*w3*t);
f=5*sin(5*w1*t);
%绘制50HZ 49*K 51*K 250HZ信号
figure(2) 
  subplot(4,1,1);
plot(t,c,'b');
xlabel('时间');
ylabel('振幅');
title('50HZ输入信号');
grid;
axis([0 0.2,-110 110]);
  subplot(4,1,2);
plot(t,d,'b');
xlabel('时间');
ylabel('振幅');
title('49*k HZ输入信号');
grid;
axis([0 0.2,-6 6]);
subplot(4,1,3);
plot(t,e,'b');
xlabel('时间');
ylabel('振幅');
title('51*k HZ输入信号');
grid;
axis([0 0.2,-2.5 2.5]);
  subplot(4,1,4);
plot(t,f,'b');
xlabel('时间');
ylabel('振幅');
title('250HZ输入信号');
grid;
axis([0 0.2,-6 6]);
%计算总输入U
u=10*a'+f+c+d+e;
%计算N,和Wn
[N,Wn]=buttord(1/3200*[49*k*0.8 51*k*1.2],1/3200*[49*k 51*k],1,40);
%采用 IIR 滤波(Third-orders Elliptic Filter),
[b,a]=butter(N,Wn,'stop');%巴特沃斯滤波器,带阻虑波
i=filter(b,a,u);
%绘制总输入,和滤波输出信号
figure(3) 
  subplot(2,1,1);
plot(t,u,'b');
xlabel('时间');
ylabel('振幅');
title('总输入信号');
grid;
axis([0 0.2,-150 150]);
  subplot(2,1,2);
plot(t,i,'b');
xlabel('时间');
ylabel('振幅');
title('滤波输出信号');
grid;
axis([0 0.2,-50 95]);
%滤波特性
[h,omega]=freqz(b,a,128);
figure(4)
plot(omega/pi,20*log10(abs(h)));
xlabel('OMEGA');
ylabel('衰减')
title('滤波特性');
grid;