%----------------------------------------------------------------- % *** 2.161 Signal Processing - Continuous and Discrete *** % Tutorial MATLAB Function % % lpcheby2 - Chebyshev Type 2 continuous lowpass filter design. % % Source: Class Handout: Introduction to Continuous Time Filter Design % % Usage: [z,p,k,n] = lpcheby2(wc, Rc, ws, Rs) or % [z,p,k,n] = lpcheby2(wc, Rc, ws, Rs,'zpk') returns the zeros, % poles, gain and order for a Chebyshev Type 1 % continuous lowpass filter. % [b,a,n] = lpcheby2(wc, Rc, ws, Rs,'tf') returns the transfer % function numerator and denominator coefficients, % and order for a Chebyshev Type 1 continuous lowpass % filter. % where wc - passband cut-off frequency (rad/s) % Rc - attenuation at passband cut-off (dB) (Rc > 0) % ws - stopband frequency (rad/s) % Rs - attenuation at stopband (dB) (Rs > 0). % % Author: D. Rowell % Revision: 2.0 9-23-2007 % %--------------------------------------------------------------- % function [varargout] = lpcheby2(wc,Rc, ws, Rs,varargin) % % Compute the required order n % epsilon = sqrt(10^(Rc/10)-1); lambda = sqrt(10^(Rs/10)-1); n = ceil(acosh(lambda/epsilon)/acosh(ws/wc)); wprod = wc*ws; epsilonhat= 1/(epsilon*cos(n*acos(ws/wc))); % % Create the finite (imaginary) zeros % Handle even and odd orders separately to eliminate the zero at infinity % for n odd if (rem(n,2)) m = n - 1; z = cos([1:2:n-2 n+2:2:2*n-1]*pi/(2*n))'; else m = n; z = cos((1:2:2*n-1)*pi/(2*n))'; end z = (z - flipud(z))./2; % Transform back to the s-plane z = i*ws./z; % Organize zeros in complex pairs: j = [1:m/2; m:-1:m/2+1]; z = z(j(:)); % % Determine the poles - choosing first only complex conjugates in the % left-half plane % alpha = asinh(1/epsilonhat)/n; p = []; for j=1:floor(n/2) gamma = (2*j - 1)*pi/(2*n); newp = -wc*(sinh(alpha)*sin(gamma)+ i*cosh(alpha)*cos(gamma)); p = [p; newp; conj(newp)]; end % If n is odd - add a real pole if rem(n,2)==1 % n is odd gamma = (2*(floor(n/2)+1) - 1)*pi/(2*n); p = [p; -wc*(sinh(alpha)*sin(gamma))]; end p=wprod./p; % % compute the gain % k = real(prod(-p)/prod(-z)); % % Return the appropriate system representation % if (nargin == 5 ) if (varargin{1}(1:2) == 'tf') sys = zpk(z,p,k); [b,a] = tfdata(sys); varargout(1) = b; varargout(2) = a; varargout(3) = {n}; elseif (varargin{1}(1:3) == 'zpk') varargout(1)={z}; varargout(2)={p}; varargout(3)={k}; varargout(4)={n}; end else varargout(1)={z}; varargout(2)={p}; varargout(3)={k}; varargout(4)={n}; end