chebyshev.m

chebyshev滤波器

function chebyshev ( mode , freq , alphap , alphas , fp , fs , fp2)
%用法：chebyshev ( mode , freq , alphap , alphas , fp , fs , fp2)
%mode用来选择滤波器模式，1为低通，2为高通，3为带通
%freq:采样频率 alphap:通带衰减 alphas:阻带衰减 
%fp:通带临界频率 fs:阻带临界频率 fp2:带通时上通带频率，低通和高通时无此输入
switch mode
    case 1
        T = 1/freq;
        wp = 2*pi*fp*T ; ws = 2*pi*fs*T;
        Op = tan(wp/2); Os = tan(ws/2);
        N = ceil(acosh(sqrt((10^(0.1*alphas)-1)/(10^(0.1*alphap)-1)))/acosh(Os/Op));
        fprintf ('N=%d\n' , N);
        e = sqrt(10^(0.1*alphap)-1);
        x(1:N) = -1/N*asinh(1/e);
        y(1:N) = 0;
        for i = 1:N
            y(i) = -(2*i-1)*pi/N/2;
        end;
        B(1:N) = 0;
        for k = 1:N
            B(k) = Op*(sinh(x(k))*sin(y(k))+j*cosh(x(k))*cos(y(k)));
        end;
        H(1:N+1) = 0;
        H(2) = 1; H(1) = B(1);
        for i = 2:N %求分母展开后系数H(0)~H(N)
            H(i+1) = 1;
            for k = i:-1:2
                H(k) = H(k)*B(i)+H(k-1);
            end;
            H(1) = H(1)*B(i);
        end;
        fprintf ('分子为%d\n分母为 \n' , Op^N/2^(N-1)/e);
        for i = 1:N+1
            fprintf ('H(%d) = %d\n' , i-1 , H(i));
        end;
        w = 0:pi/49:pi ;%w、s是向量
        s(1:50) = 0;
        for i=1:50
            s(i) = (1-exp(-j*w(i))) / (1+exp(-j*w(i)));
        end;
        fenmu(1:50) = 0;
        for i = 1:N+1
            for k=1:50
                fenmu(k) = fenmu(k)+H(i)*s(k)^(i-1);
            end;
        end;
        output(1:50)=0;
        output = Op^N/2^(N-1)/e./fenmu;
        stem ( w , abs(output) );
    case 2
        T = 1/freq;
        wp = 2*pi*fp*T ; ws = 2*pi*fs*T;
        Op = abs(cot(wp/2)); Os = abs(cot(ws/2));
        N = ceil(acosh(sqrt((10^(0.1*alphas)-1)/(10^(0.1*alphap)-1)))/acosh(Os/Op));
        fprintf ('N=%d\n' , N);
        e = sqrt(10^(0.1*alphap)-1);
        x(1:N) = -1/N*asinh(1/e);
        y(1:N) = 0;
        for i = 1:N
            y(i) = -(2*i-1)*pi/N/2;
        end;
        B(1:N) = 0;
        for k = 1:N
            B(k) = Op*(sinh(x(k))*sin(y(k))+j*cosh(x(k))*cos(y(k)));
        end;
        H(1:N+1) = 0;
        H(2) = 1; H(1) = B(1);
        for i = 2:N %求分母展开后系数H(0)~H(N)
            H(i+1) = 1;
            for k = i:-1:2
                H(k) = H(k)*B(i)+H(k-1);
            end;
            H(1) = H(1)*B(i);
        end;
        fprintf ('分子为%d\n分母为 \n' , Op^N/2^(N-1)/e);
        for i = 1:N+1
            fprintf ('H(%d) = %d\n' , i-1 , H(i));
        end;
        w = 0:pi/49:pi ;%w、s是向量
        s(1:50) = 0;
        for i=1:50
            s(i) = (1+exp(-j*w(i))) / (1-exp(-j*w(i)));
        end;
        fenmu(1:50) = 0;
        for i = 1:N+1
            for k=1:50
                fenmu(k) = fenmu(k)+H(i)*s(k)^(i-1);
            end;
        end;
        output(1:50)=0;
        output = Op^N/2^(N-1)/e./fenmu;
        stem ( w , abs(output) );
    case 3
        T = 1/freq;
        wp = 2*pi*fp*T ; ws = 2*pi*fs*T; wp2 = fp2*pi*2*T;
        cw0 = sin(wp+wp2)/(sin(wp)+sin(wp2));
        Op = abs((cw0-cos(wp))/sin(wp)); Os = abs((cw0-cos(ws))/sin(ws));
        N = ceil(acosh(sqrt((10^(0.1*alphas)-1)/(10^(0.1*alphap)-1)))/acosh(Os/Op));
        fprintf ('N=%d\n' , N);
        e = sqrt(10^(0.1*alphap)-1);
        x(1:N) = -1/N*asinh(1/e);
        y(1:N) = 0;
        for i = 1:N
            y(i) = -(2*i-1)*pi/N/2;
        end;
        B(1:N) = 0;
        for k = 1:N
            B(k) = Op*(sinh(x(k))*sin(y(k))+j*cosh(x(k))*cos(y(k)));
        end;
        H(1:N+1) = 0;
        H(2) = 1; H(1) = B(1);
        for i = 2:N %求分母展开后系数H(0)~H(N)
            H(i+1) = 1;
            for k = i:-1:2
                H(k) = H(k)*B(i)+H(k-1);
            end;
            H(1) = H(1)*B(i);
        end;
        fprintf ('分子为%d\n分母为 \n' , Op^N/2^(N-1)/e);
        for i = 1:N+1
            fprintf ('H(%d) = %d\n' , i-1 , H(i));
        end;
        w = 0:pi/49:pi ;%w、s是向量
        s(1:50) = 0;
        for i=1:50
            s(i) = (exp(j*w(i))^2-2*exp(j*w(i))*cw0+1)/(exp(j*w(i))^2-1);
        end;
        fenmu(1:50) = 0;
        for i = 1:N+1
            for k=1:50
                fenmu(k) = fenmu(k)+H(i)*s(k)^(i-1);
            end;
        end;
        output(1:50)=0;
        output = Op^N/2^(N-1)/e./fenmu;
        stem ( w , abs(output) );
end;