滤波器滤波过程用差分方程运算完成

clear all;
clc;
close all;
fs=1600;
f0=50;
N=400;
t=(0:N-1)/fs;
x=zeros(1,N);


for k=1:2:10
x=x+(10/pi/k)*sin(2*pi*k*f0*t);
end


fs2=fs/2;
wp=[40 60]/fs2;
ws=[30 80]/fs2;
Rp=1;
Rs=40;
[M,Wn]=ellipord(wp,ws,Rp,Rs);
[B,A]=ellip(M,Rp,Rs,Wn);
[H,f]=freqz(B,A,1000,fs);
gdy=grpdelay(B,A,1000,fs);


fprintf('B=%5.6e  %5.6e  %5.6e  %5.6e  \n',B);
fprintf('\n');
fprintf('A=%5.6e  %5.6e  %5.6e  %5.6e  \n',A);
fprintf('\n');


yy=filter(B,A,x);


for k=1:6
xx(k)=0;
y(k)=0;
end


for k=7:N+6
j=k-6;
xx(k)=x(j);
y(k)=B(1)*xx(k)+B(2)*xx(k-1)+B(3)*xx(k-2)+B(4)*xx(k-3)+B(5)*xx(k-4)...
    +B(6)*xx(k-5)+B(7)*xx(k-6)-A(2)*y(k-1)-A(3)*y(k-2)-A(4)*y(k-3)...
    -A(5)*y(k-4)-A(6)*y(k-5)-A(7)*y(k-6);
end
y=y(7:end);


figure(1);
subplot 211;
plot(f,20*log10(abs(H)));
title('椭圆形带通滤波器幅频响应曲线');
xlabel('频率/Hz');
ylabel('幅值/dB');


subplot 212;
plot(f,gdy);
title('椭圆形带通滤波器群延迟响应曲线');
xlabel('频率/Hz');
ylabel('群延迟');


figure(2);
subplot 211;
plot(t,x);
title('输入信号波形');
xlabel('时间/s');
ylabel('幅值');


subplot 212;
plot(t,yy,'color',[.6 .6 .6],'linewidth',3);
hold on;
plot(t,y,'k');
title('输出信号波形');
legend('1','2',2);
xlabel('时间/s');
ylabel('幅值');
box on;