matlab低通滤波器

clc; % 清除命令窗口
clear; % 清除所有变量
close all; % 关闭所有的图形窗口

N=2^10; % 定义一个变量N，值为2的10次方
n=0:N-1; % 生成一个从0到N-1的序列
Fs=2000; % 采样频率
tn=n*(1/Fs); % 时间序列
Fn=Fs/N*n; % 频率序列
Wk=2*pi*Fn; % 角频率序列
Wn=Wk/Fs; % 标准化角频率序列

f=100; % 信号频率
f_noise=250; % 噪声频率
f_noise2=20; % 第二个噪声频率
xn=sin(2*pi*f*tn); % 生成正弦信号

xn_noise = xn+0.6*sin(2*pi*f_noise*tn)+0.2*sin(2*pi*f_noise2*tn); % 添加噪声到原始信号


Xk_noise=fft(xn_noise,N); % FFT变换
Xk_noise=abs(Xk_noise); % 获取FFT结果的幅度

subplot(2,2,4); % 切换到第四个子图
plot(Fn,Xk_noise); % 绘制频域的噪声信号
xlim([0,Fs/2]);

Wp=200;Ws=300; % 设定通带截止频率和阻带截止频率
Wp=2*pi*Wp;Ws=2*pi*Ws; % 将截止频率转化为角频率
Rp=1;Rs=40; % 设定通带最大衰减和阻带最小衰减

[M,Wc]=buttord(Wp,Ws,Rp,Rs,"s"); % 计算滤波器的阶数和截止角频率
[B,A]=butter(M,Wc,"low","s"); % 设计一个模拟Butterworth低通滤波器
fprintf('The order of the filter is %d\n', M);%显示滤波器的阶数。
Hs=freqs(B,A,Wk); % 计算滤波器的频率响应
Hs=abs(Hs);
Hs_copy=Hs;
Hs=20*log10(Hs);

[Bz,Az] = impinvar(B,A,Fs);
Hz=freqz(Bz,Az,Wn);
Hz=abs(Hz);
Hz_copy=Hz;
Hz=20*log10(Hz);

subplot(2,2,3); % 切换到第三个子图
plot(Fn,Hs,"r",Fn,Hz,"b"); % 绘制滤波器的频率响应
legend("模拟滤波器","数字滤波器");
xlim([0,Fs/2]);

Yk=Hs_copy.*Xk_noise;
subplot(2,2,2); % 切换到第三个子图
plot(Fn,Yk); % 绘制滤波器的频率响应
xlim([0,Fs/2]);

% yn=ifft(Yk,N);
yn=filtfilt(Bz,Az,xn_noise);
subplot(2,2,1); % 切换到第二个子图
plot(tn,xn,"r",tn,xn_noise,"b",tn,yn,"g"); % 绘制添加了噪声的信号
legend("原信号","加噪信号","滤波器输出信号");
xlim([0,20/f]);