www.gusucode.com > signal 工具箱matlab源码程序 > signal/ellip.m

    function [num, den, z, p] = ellip(n, Rp, Rs, Wn, varargin)
%ELLIP Elliptic or Cauer digital and analog filter design.
%   [B,A] = ELLIP(N,Rp,Rs,Wp) designs an Nth order lowpass digital elliptic
%   filter with Rp decibels of peak-to-peak ripple and a minimum stopband
%   attenuation of Rs decibels. ELLIP returns the filter coefficients in
%   length N+1 vectors B (numerator) and A (denominator). The passband-edge
%   frequency Wp must be 0.0 < Wp < 1.0, with 1.0 corresponding to half the
%   sample rate.  Use Rp = 0.5 and Rs = 20 as starting points, if you are
%   unsure about choosing them.
%
%   If Wp is a two-element vector, Wp = [W1 W2], ELLIP returns an order 2N
%   bandpass filter with passband  W1 < W < W2. [B,A] =
%   ELLIP(N,Rp,Rs,Wp,'high') designs a highpass filter. [B,A] =
%   ELLIP(N,Rp,Rs,Wp,'low') designs a lowpass filter. [B,A] =
%   ELLIP(N,Rp,Rs,Wp,'stop') is a bandstop filter if Wp = [W1 W2].
%
%   When used with three left-hand arguments, as in [Z,P,K] = ELLIP(...),
%   the zeros and poles are returned in length N column vectors Z and P,
%   and the gain in scalar K. 
%
%   When used with four left-hand arguments, as in [A,B,C,D] = ELLIP(...),
%   state-space matrices are returned. 
%
%   ELLIP(N,Rp,Rs,Wp,'s'), ELLIP(N,Rp,Rs,Wp,'high','s') and 
%   ELLIP(N,Rp,Rs,Wp,'stop','s') design analog elliptic filters. In this
%   case, Wp is in [rad/s] and it can be greater than 1.0.
%
%   % Example 1:
%   %   For data sampled at 1000 Hz, design a sixth-order lowpass  
%   %   elliptic filter with a passband edge frequency of 300Hz, 3 dB of 
%   %   ripple in the passband, and 50 dB of attenuation in the stopband.
%
%   Wn = 300/500;               % Normalized passband edge frequency
%   [z,p,k] = ellip(6,3,50,Wn);
%   [sos] = zp2sos(z,p,k);      % Convert to SOS form
%   h = fvtool(sos)             % Plot magnitude response
%
%   % Example 2:
%   %   Design a 6th-order Elliptic band-pass filter which passes 
%   %   frequencies between 0.2 and 0.5, and with 5 dB of ripple in the
%   %   passband, and 80 dB of attenuation in the stopband
%
%   [b,a]=ellip(6,5,80,[.2,.5]);    % Bandpass digital filter design         
%   h = fvtool(b,a);                % Visualize filter
%
%   See also ELLIPORD, CHEBY1, CHEBY2, BUTTER, BESSELF, FREQZ, FILTER,
%   DESIGNFILT.

%   Author(s): L. Shure, 4-29-88
%   	   T. Krauss, 3-25-93, revised
%   Copyright 1988-2013 The MathWorks, Inc.

%   References:
%     [1] T. W. Parks and C. S. Burrus, Digital Filter Design,
%         John Wiley & Sons, 1987, chapter 7, section 7.3.3.


[btype,analog,errStr,msgobj] = iirchk(Wn,varargin{:});
if ~isempty(errStr) 
  error(msgobj);
end

% Cast to enforce precision rules. 
n = signal.internal.sigcasttofloat(n,'double','ellip','N','allownumeric');
Rp = signal.internal.sigcasttofloat(Rp,'double','ellip','Rp','allownumeric');
Rs = signal.internal.sigcasttofloat(Rs,'double','ellip','Rs','allownumeric');
Wn = signal.internal.sigcasttofloat(Wn,'double','ellip','Wp','allownumeric');

if n>500
	error(message('signal:ellip:InvalidRange'))
end

% step 1: get analog, pre-warped frequencies
if ~analog,
	fs = 2;
	u = 2*fs*tan(pi*Wn/fs);
else
	u = Wn;
end

% step 2: convert to low-pass prototype estimate
if btype == 1	% lowpass
	Wn = u;
elseif btype == 2	% bandpass
	Bw = u(2) - u(1);
	Wn = sqrt(u(1)*u(2));	% center frequency
elseif btype == 3	% highpass
	Wn = u;
elseif btype == 4	% bandstop
	Bw = u(2) - u(1);
	Wn = sqrt(u(1)*u(2));	% center frequency
end

% step 3: Get N-th order elliptic lowpass analog prototype
[z,p,k] = ellipap(n, Rp, Rs);

% Transform to state-space
[a,b,c,d] = zp2ss(z,p,k);

% step 4: Transform to lowpass, bandpass, highpass, or bandstop of desired Wn
if btype == 1		% Lowpass
	[a,b,c,d] = lp2lp(a,b,c,d,Wn);

elseif btype == 2	% Bandpass
	[a,b,c,d] = lp2bp(a,b,c,d,Wn,Bw);

elseif btype == 3	% Highpass
	[a,b,c,d] = lp2hp(a,b,c,d,Wn);

elseif btype == 4	% Bandstop
	[a,b,c,d] = lp2bs(a,b,c,d,Wn,Bw);
end

% step 5: Use Bilinear transformation to find discrete equivalent:
if ~analog,
	[a,b,c,d] = bilinear(a,b,c,d,fs);
end

if nargout == 4
	num = a;
	den = b;
	z = c;
	p = d;
else	% nargout <= 3
% Transform to zero-pole-gain and polynomial forms:
	if nargout == 3
		[z,p,k] = ss2zp(a,b,c,d,1);
		num = z;
		den = p;
		z = k;
	else % nargout <= 2
		den = poly(a);
    if analog
      zeroLimitFlag = true;
    else
      zeroLimitFlag = false;
    end
    zinf = ltipack.getTolerance('infzero',zeroLimitFlag);
    [z,k] = ltipack.sszero(a,b,c,d,[],zinf);
    num = k * poly(z);
    num = [zeros(1,length(den)-length(num))  num];
	end
end