matlab自己设计一个低通滤波器，要求滤出100Hz之外的频率，采样率为10000Hz

clc;clear all;
%归一化模拟切比雪夫I型低通滤波器的设计
Wp=2*pi*1000;Ws=2*pi*1500;rp=3;rs=30;%设计滤波器的参数
wp=1;ws=Ws/Wp; %频带变换得到归一化滤波器
[N,wc]=cheb1ord(wp,ws,rp,rs,'s'); %计算滤波器阶数和3dB截止频率
[z,p,k]=cheb1ap(N,rp); %计算归一化滤波器的零极点
[b,a]=zp2tf(z,p,k); %计算归一化滤波器系统函数的系数
w0=0:0.05*pi:2*pi;
[h0,w0]=freqs(b,a,w0); %求频率响应
figure(1)
plot(w0,20*log10(abs(h0)),'k');
title('归一化模拟且比雪夫I型低通滤波器');
xlabel('频率f/Hz');ylabel('幅度/dB');grid;
%一般低通切比雪夫低通滤波器的设计
[B,A]=lp2lp(b,a,Wp); %将归一化滤波器转换为一般模拟滤波器
w1=0:2*pi*100:2*pi*30000;
[h1,w1]=freqs(B,A,w1);
figure(2)
plot(w1/(2*pi),20*log10(abs(h1)),'k');
title('一般模拟且比雪夫I型低通滤波器');
xlabel('频率f/Hz');ylabel('幅度/dB');grid;
%冲激响应不变法设计低通巴特沃斯数字滤波器
Fs=10000; %采样频率
Wp1=Wp/Fs;Ws1=Ws/Fs;
rp1=3;rs1=30; %设计滤波器的参数
[N1,Wn]=cheb1ord(Wp,Ws,rp1,rs1,'s') %计算滤波器阶数
[b1,a1]=cheby1(N1,rp1,Wn,'s'); %样本AF的系数函数的分子分母系数
[bz,az]=impinvar(b1,a1,Fs); %冲激响应不变法从AF到DF变换
[C1,B1,A1]=dir2par(bz,az) %直接形式转换为并联型
w2=[Wp1,Ws1]; %数字临界频率
[H,f]=freqz(bz,az); %绘制数字滤波器的幅频特性和相频特性
[db,mag,pha,grd,f]=freqz_m(bz,az); %扩展函数检验滤波器的其他特性
figure(3)
plot(f/pi,db,'k');
title('冲激响应不变法设计低通切比雪夫I型数字滤波器');
xlabel('角频率w/pi');ylabel('振幅/dB');
axis([0,0.35,-30,5]);grid;
%用设计好的滤波器对信号进滤波处理
figure(4)
f1=500;f2=4000; %输入信号的频率
N=100; %数据长度
dt=1/Fs;n=0:N-1;t=n*dt; %采样间隔和时间序列
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号
subplot(2,1,1),plot(t,x),title('输入信号');
y=filtfilt(bz,az,x); %用滤波器进行滤波处理
y1=filter(bz,az,x); %进行滤波处理
subplot(2,1,2),plot(t,y,t,y1,':'),title('输出信号');xlabel('时间/s');
legend('filtfilt','filter') %加图例

freqz_m.m文件
function[db,mag,pha,grd,w]=freqz_m(b,a)
[H,w]=freqz(b,a,1000,'whole');
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(pha);

dir2par.m文件
function [C,B,A] = dir2par(b,a)
M=length(b);N=length(a);
[r1,p1,C]=residuez(b,a);
p=cplxpair(p1,10000000*eps);
I=cplxcomp(p1,p);
r=r1(I);
K=floor(N/2);B=zeros(K,2);A=zeros(K,3);
if K*2==N;
for i=1:2:N-2
Brow=r(i:1:i+1,:);
Arow=p(i:1:i+1,:);
[Brow,Arow]=residuez(Brow,Arow,[]);
B(fix((i+1)/2),:)=real(Brow);
A(fix((i+1)/2),:)=real(Arow);
end
[Brow,Arow]=residuez(r(N-1),p(N-1),[]);
B(K,:)=[real(Brow) 0];A(K,:)=[real(Arow) 0];
else
for i=1:2:N-1
Brow=r(i:1:i+1,:);
Arow=p(i:1:i+1,:);
[Brow,Arow]=residuez(Brow,Arow,[]);
B(fix((i+1)/2),:)=real(Brow);
A(fix((i+1)/2),:)=real(Arow);
end
end

cplxcomp.m文件
function I=cplxcomp(p1,p2)
I=[];
for j=1:1:length(p2)
for i=1:1:length(p1)
if(abs(p1(i)-p2(j))<0.0001)
I=[I,i];
end
end
end
I=I';