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;