巴特沃斯滤波器的两种等价的表达公式

clc;
clear;
close all;

h=0.001;
gama=0.3;
T=1;
N=T/h;
s_f = 10;
n_f1 = 100;
n_f2 = 300;

M = 2;

%butterworth 2-order ,fs:1000Hz  wc:50Hz
filter_A = [1,   2,                 1               ] * 0.02008336556421;
filter_B = [1,   -1.561018075801,   0.6413515380576 ];

u1=0;u2=0;u3=0;
y1=0;y2=0;y3=0;

for i=1:N
   t(i) = i*h;
   
   s(i) = sin(2*pi*s_f*t(i));
   n1(i) = gama*sin(2*pi*n_f1*t(i));
   n2(i) = gama*sin(2*pi*n_f2*t(i));
   s_n(i) = s(i) + n1(i) + n2(i);
    
   u1 = s_n(i);
   if M==1
       y(i) = filter_A(1)*u1 + filter_A(2)*u2 + filter_A(3)*u3 - filter_B(2)*y2 - filter_B(3)*y3;

       u3 = u2;
       u2 = u1;
       y3 = y2;
       y2 = y(i);
   elseif M==2
       den = u1 - filter_B(2)*y2 - filter_B(3)*y3;
       y(i) = filter_A(1)*den + filter_A(2)*y2 + filter_A(3)*y3;
       
       y3 = y2;
       y2 = den;
   end
   
end

figure(1);
plot(t,s,t,s_n);legend('s','s\_n');

figure;
plot(t,s_n,t,y);legend('s\_n','y');

[abs_f,abs_mag] = abs_fft(s_n,h);
figure;
plot(abs_f,abs_mag);
[abs_f,abs_mag] = abs_fft(y,h);
hold on;
plot(abs_f,abs_mag);
legend("s_n","y");


function [abs_f,abs_mag] = abs_fft(x,dt)

N=length(x);
mag = fft(x,N);
f = (1:N)/N/dt;

abs_mag = abs(mag(1:(N/2+1)));
abs_f   = abs(f(1:(N/2+1)));

end

文档链接：公式及滤波器设计说明