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
文档链接:公式及滤波器设计说明