www.gusucode.com > signal 工具箱matlab源码程序 > signal/ellipord.m
function [order,wn] = ellipord(wp,ws,rp,rs,opt) %ELLIPORD Elliptic filter order selection. % [N, Wp] = ELLIPORD(Wp, Ws, Rp, Rs) returns the order N of the lowest % order digital elliptic filter which has a passband ripple of no more % than Rp dB and a stopband attenuation of at least Rs dB. Wp and Ws are % the passband and stopband edge frequencies, normalized from 0 to 1 % (where 1 corresponds to pi radians/sample). For example, % Lowpass: Wp = .1, Ws = .2 % Highpass: Wp = .2, Ws = .1 % Bandpass: Wp = [.2 .7], Ws = [.1 .8] % Bandstop: Wp = [.1 .8], Ws = [.2 .7] % ELLIPORD also returns Wp, the elliptic natural frequency to use with % ELLIP to achieve the specifications. % % [N, Wp] = ELLIPORD(Wp, Ws, Rp, Rs, 's') does the computation for an % analog filter, in which case Wp and Ws are in radians/second. % % NOTE: If Rs is much much greater than Rp, or Wp and Ws are very close, % the estimated order can be infinite due to limitations of numerical % precision. % % % Example 1: % % For 1000 Hz data, design a lowpass filter with less than 3 dB of % % ripple in the passband defined from 0 to 40 Hz and at least 60 dB % % of ripple in the stopband defined from 150 Hz to the Nyquist % % frequency (500 Hz): % % Wp = 40/500; Ws = 150/500; % Rp = 3; Rs = 60; % [n,Wp] = ellipord(Wp,Ws,Rp,Rs) % Gives mimimum order of filter % [b,a] = ellip(n,Rp,Rs,Wp); % Elliptic filter design % freqz(b,a,512,1000); % Plots the frequency response % title('n=4 Elliptic Lowpass Filter') % % % Example 2: % % Design a bandpass filter with a passband of 60 Hz to 200 Hz, with % % less than 3 dB of ripple in the passband, and 40 dB attenuation in % % the stopbands that are 50 Hz wide on both sides of the passband. % % Wp = [60 200]/500; Ws = [50 250]/500; % Rp = 3; Rs = 40; % [n,Wp] = ellipord(Wp,Ws,Rp,Rs) % Gives mimimum order of filter % [b,a] = ellip(n,Rp,Rs,Wp); % Elliptic filter design % freqz(b,a,512,1000) % Plots the frequency response % % See also ELLIP, BUTTORD, CHEB1ORD, CHEB2ORD, DESIGNFILT. % Author(s): L. Shure, 6-9-88 % T. Krauss, 11-18-92, updated % Copyright 1988-2013 The MathWorks, Inc. % Reference(s): % [1] Rabiner and Gold, p 241. narginchk(4,5); nargoutchk(0,2); % Cast to enforce precision rules wp = signal.internal.sigcasttofloat(wp,'double','ellipord','Wp','allownumeric'); ws = signal.internal.sigcasttofloat(ws,'double','ellipord','Ws','allownumeric'); rp = signal.internal.sigcasttofloat(rp,'double','ellipord','Rp','allownumeric'); rs = signal.internal.sigcasttofloat(rs,'double','ellipord','Rs','allownumeric'); if nargin == 4 opt = 'z'; elseif nargin == 5 if ~strcmp(opt,'z') && ~strcmp(opt,'s') error(message('signal:ellipord:InvalidParam')); end end [msg,msgobj]=freqchk(wp,ws,opt); if ~isempty(msg), error(msgobj); end ftype = 2*(length(wp) - 1); if wp(1) < ws(1) ftype = ftype + 1; % low (1) or reject (3) else ftype = ftype + 2; % high (2) or pass (4) end % first, prewarp frequencies from digital (unit circle) to analog (imag. axis): if strcmp(opt,'z') % digital WP=tan(pi*wp/2); WS=tan(pi*ws/2); else % don't have to if analog already WP=wp; WS=ws; end % next, transform to low pass prototype with passband edge of 1 and stopband % edges determined by the following: (see Rabiner and Gold, p.258) if ftype == 1 % low WA=WS/WP; order = findelliporder(WA,rp,rs); elseif ftype == 2 % high WA=WP/WS; order = findelliporder(WA,rp,rs); elseif ftype == 3 % stop if strcmp(opt,'s'), % For analog bandstop, convert back to digital wp = 2*atan(wp)/pi; ws = 2*atan(ws)/pi; end Fpass1 = wp(1); Fpass2 = wp(2); Fstop1 = ws(1); Fstop2 = ws(2); c = sin(pi*(Fpass1+Fpass2))/(sin(pi*Fpass1)+sin(pi*Fpass2)); wpa = abs(sin(pi*Fpass2)/(cos(pi*Fpass2)-c)); ws1 = sin(pi*Fstop1)/(cos(pi*Fstop1)-c); ws2 = sin(pi*Fstop2)/(cos(pi*Fstop2)-c); wsa = min(abs([ws1,ws2])); order = ellipord(wpa,wsa,rp,rs,'s'); elseif ftype == 4 % pass WA=(WS.^2 - WP(1)*WP(2))./(WS*(WP(1)-WP(2))); order = findelliporder(WA,rp,rs); end % natural frequencies are simply the passband edges (WP). % finally, transform frequencies from analog to digital if necessary: if strcmp(opt,'z') % digital wn = wp; else wn = WP; end %-------------------------------------------------------------------------- function order = findelliporder(WA,rp,rs) % find the minimum order elliptic filter to meet the more demanding spec: WA = min(abs(WA)); epsilon = sqrt(10^(0.1*rp)-1); k1 = epsilon/sqrt(10^(0.1*rs)-1); k = 1/WA; capk = ellipke([k^2 1-k^2]); capk1 = ellipke([(k1^2) 1-(k1^2)]); order = ceil(capk(1)*capk1(2)/(capk(2)*capk1(1))); % if both warnings are in effect, only print the first one if (1-k1^2) == 1 warning(message('signal:ellipord:MustBeFiniteAttn', 'ellipord')) elseif k^2 == 1 warning(message('signal:ellipord:MustBeFiniteBandEdges', 'ellipord')) end