matlab测试IIR低通和陷波滤波器

下面的代码直接可以用于仿真测试

clear all;
clc;

f0=50;
Ts=0.0025;%采样率1khz
len=400;

n=0:len-1;
x=sin(2*pi*30*n*Ts)+sin(2*pi*150*n*Ts)+sin(2*pi*100*n*Ts) + sin(2*pi*70*n*Ts);
y=zeros(1,len);

%B=[0.93721 -1.32542 0.93721];%50Hz陷波系数
%A=[1 -1.3254 0.87443];
%Gain=1;

B=[1,2,1];%2阶低通巴特沃斯，FS=1000,FC=60
A=[1, -1.47548044, 0.586919508];
Gain = 0.027859;

w0=[0 0 0];
w1=[0 0 0];
 
for i=1:len
    w0(1)=x(i);
    c=B(1)*w0(1)+B(2)*w0(2)+B(3)*w0(3);
    y(i)=c*Gain-A(2)*w1(2)-A(3)*w1(3);
    w1(1)=y(i);
    w0(3)=w0(2);w0(2)=w0(1);
    w1(3)=w1(2);w1(2)=w1(1);
end
%或者下面的代码
% for i=1:len
%     w0(1)=A(1)*x(i)-A(2)*w0(2)-A(3)*w0(3);
%     y(i)=(B(1)*w0(1)+B(2)*w0(2)+B(3)*w0(3))*Gain;
%     w0(3)=w0(2);w0(2)=w0(1);
% end
figure
subplot(221);
plot(x);
subplot(222);
plot(y);
subplot(223);
a = fft(x,len);
plot(n,abs(a));
subplot(224);
b = fft(y,len);
plot(n,abs(b));