butterworth.m

在matlab上实现通过巴特沃斯滤波器来设计数字滤波器。

function Butterworth(fp,fs,ap,as,Fs)
if nargin==5              %%对应数字滤波器设计
    wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
    Ts=1/Fs;
    Op=2*tan(wp/2)/Ts;Os=2*tan(ws/2)/Ts;
elseif nargin==4         %%对应模拟滤波器设计
    Op=fp;Os=fs;
end
rp=1;rs=Os/Op;           %归一化
C=sqrt(10^(ap/10)-1);
N=ceil(log10(sqrt(10^(as/10)-1)/C)/log10(rs));  %N向上取整
syms G p s z w g x;      %定义符号
G=1;                     
for i=1:N
    angle1=pi*(2*i+N-1)/(2*N);
    pk=exp(angle1*j);
    G=G*(p-pk);
end
G=1/G;                 %求得G(p)
if nargin==4
    g=s/Op;
    Gs=compose(G,g);    %用g=s/Op替代符号p
    g=j*s;
    Gr=compose(Gs,g);
    S=0:0.5:fs;
    S1=S;
    GS=subs(Gr,s,S1);
    figure(1);
    plot(S,20*log10(abs(GS)));grid on;title('对数幅频曲线');xlabel('f/Hz');ylabel('20lg|H(exp(jw))|');
    figure(2);
    plot(S,angle(GS));grid on;title('对数幅频曲线');xlabel('f/Hz');ylabel('angle(H(exp(jw))');
elseif nargin==5
    g=j*x;
    Gx=compose(G,g);
    g=s/Op;
    Gx=compose(Gx,g);%求得G(s)
    Gs=compose(G,g)
    g=2/Ts*(z-1)/(z+1);
    Hz=compose(Gs,g);%求得H(z)
    g=exp(j*w);
    Hw=compose(Hz,g)%求得H（exp(jw)）
    W=0:0.5:fs;
    W1=W*2*pi/Fs;
    W2=W*2*pi;
    GW=subs(Hw,w,W1);
    GS=subs(Gx,s,W2);
    figure(1);
    plot(W,20*log10(abs(GW)),'r-',W,20*log10(abs(GS)),'b-');grid on;title('对数幅频曲线');xlabel('f/Hz');ylabel('20lg|H(exp(jw))|');texlabel('r--DF，b--AF');
    figure(2);
    plot(W,angle(GW),'r-',W,angle(GS),'b-');grid on;title('对数相频曲线');xlabel('f/Hz');ylabel('angle(H(exp(jw))');
end