www.gusucode.com > signal 工具箱matlab源码程序 > signal/+fspecs/@bscutoffwap/cheby1.m
function Hd = cheby1(this, varargin) %CHEBY1 Chebyshev Type I digital filter design. % Copyright 1988-2012 The MathWorks, Inc. N = this.FilterOrder; Ap = this.Apass; nfreq = this.NormalizedFrequency; normalizefreq(this, true); Fc1 = this.F3dB1; Fc2 = this.F3dB2; normalizefreq(this, nfreq); [isvalid, errmsg, errid] = checkincfreqs(this,{'F3dB1','F3dB2'}); if ~isvalid error(message(errid,errmsg)); end % Make up a 3dB point for digital lowpass (arbitrarily) Fc = (Fc1+Fc2)/2; % Compute theoretical lowpass passband edge Wc = tan(pi*Fc/2); % Find corresponding analog passband-edge frequency Wp = Wc/cosh(1/(N/2)*acosh(1/sqrt(10^(Ap/10)-1))); % Use half the order % Convert analog passband-edge frequency to digital Fp = 2*atan(Wp)/pi; % Compute alpha alpha = cos((Fc2+Fc1)*pi/2)/cos((Fc2-Fc1)*pi/2); % Compute k k = tan(Fc*pi/2)*tan((Fc2-Fc1)*pi/2); c1 = 2*alpha/(k+1); c2 = (1-k)/(1+k); % Solve LP to BS inverse Constantinides (Antoniou p.243) j = complex(0,1); z = exp(j*Fp*pi); b = c1*(z-1); ax2 = 0.5/(c2*z-1); b2_4ac = sqrt(c1^2*(1-z)^2-4*(c2*z-1)*(z-c2)); Fp1 = abs(real(log(ax2*(b - b2_4ac))/(j*pi))); Fp2 = abs(real(log(ax2*(b + b2_4ac))/(j*pi))); % Convert to bandstop with passband-edge specifications hs = fspecs.bspass(N,Fp1,Fp2,Ap); Hd = cheby1(hs,varargin{:});