figure

fs=500;%采样率
f1=30;%信号频率
f2=70;%信号频率
f3 =100;%信号频率
f4 =150;%信号频率
T=2;%时宽2s
n=round(T*fs);%采样点个数
t=linspace(0,T,n);%时域横坐标
x = sin(2*pi*f1*t) + sin(2*pi*f2*t) + sin(2*pi*f3*t) + 2*cos(2*pi*f4*t);
subplot(221);
plot(t,x);%画时域图
xlabel("t/s")
grid on

X = fftshift(fft(x./(n))); %用fft得出离散傅里叶变换
f=linspace(-fs/2,fs/2-1,n);%频域横坐标，注意奈奎斯特采样定理，最大原信号最大频率不超过采样频率的一半
subplot(222);
plot(f,abs(X));%画双侧频谱幅度图
xlabel("f/Hz")
ylabel("幅度")
grid on

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:n
    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

subplot(223);
plot(t,y);%画双侧频谱幅度图

subplot(224);
Y = fftshift(fft(y./(n)));
plot(f,abs(Y))