butterworth.m

It is about Matlab using in DSP ,welcome to download.

function butterworth(wp,ws,ap,as,type,a,b)
N=0.5*log10((10.^(as/10)-1)/(10.^(ap/10)-1))/log10(ws/wp);
N=ceil(N)
Wc1=wp/(10.^(ap/10)-1).^(1/(2*N));
Wc2=ws/(10.^(as/10)-1).^(1/(2*N));
Wc=(Wc1+Wc2)/2
f=0:0.0001:0.5;
A=(1./(1+(2*pi*f/Wc).^(2*N)));
subplot(2,1,1);
plot(f,10*log10(A));
title('butterworth模拟低通滤波器');
for k=1:N
    s(k)=Wc*exp(j*(1/2+(2*k-1)/(2*N))*pi);
end
pp=poly(s);
k0=pp(length(pp));
fprintf('butterworth模拟低通滤波器H(s)：\n');
fprintf('                 %.4g\n',k0);
fprintf('-----------------------------------------------------------------------------------------------------------------\n');
pstr1(pp,'s');
fprintf('\n');

H=Wc.^N;
f=0.01:0.001:0.49;
if type=='l'
    c=1;
    for k=1:N
        H=H./(((1-exp(-j*2*pi*f))./(1+exp(-j*2*pi*f)))-s(k));
        c=c/(1-s(k));
    end
    v=(1+s)./(1-s);
    pz=poly(v);
    fprintf('数字低通滤波器H(z)：\n');
    fprintf('             %.4g(z+1).^%.0f\n',c,N);
    fprintf('----------------------------------------------------------------------------------------------------------------\n');
    pstr1(pz,'z');
    fprintf('\n');
    
end
if type=='h'
    c=1;
    for k=1:N
        H=H./(((1+exp(-j*2*pi*f))./(1-exp(-j*2*pi*f)))-s(k));
        c=c/(s(k)-1);
    end
    v=(s+1)./(s-1);
    pz=poly(v);
    fprintf('数字高通滤波器H(z)：\n');
    fprintf('             %.4g(z-1).^%.0f\n',c,N);
    fprintf('----------------------------------------------------------------------------------------------------------------\n');
    pstr1(pz,'z');
    fprintf('\n');
    
end
if type=='bp'
    c=1;
    for k=1:N
        H=H./(a*(exp(j*2*2*pi*f)-2*b*exp(j*2*pi*f)+1)./(exp(j*2*2*pi*f)-1)-s(k));
        c=c/(a-s(k));
    end
    p=[a-s -2*a*b a+s];
    v=roots(p);
    pz=poly(v);
    fprintf('数字带通滤波器H(z)：\n');
    fprintf('             %.4g(z.^2-1).^%.0f\n',c,N);
    fprintf('----------------------------------------------------------------------------------------------------------------\n');
    pstr1(pz,'z');
    fprintf('\n');
end
subplot(2,1,2);
plot(f,20*log10(abs(H)));
title('数字滤波器');