gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\bilinear.m
function [zd, pd, kd, dd] = bilinear(z, p, k, fs, fp, fp1) %BILINEAR Bilinear transformation with optional frequency prewarping. % [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs) converts the s-domain transfer % function specified by Z, P, and K to a z-transform discrete % equivalent obtained from the bilinear transformation: % % H(z) = H(s) | % | s = 2*Fs*(z-1)/(z+1) % % where column vectors Z and P specify the zeros and poles, scalar % K specifies the gain, and Fs is the sample frequency in Hz. % [NUMd,DENd] = BILINEAR(NUM,DEN,Fs), where NUM and DEN are % row vectors containing numerator and denominator transfer % function coefficients, NUM(s)/DEN(s), in descending powers of % s, transforms to z-transform coefficients NUMd(z)/DENd(z). % [Ad,Bd,Cd,Dd] = BILINEAR(A,B,C,D,Fs) is a state-space version. % Each of the above three forms of BILINEAR accepts an optional % additional input argument that specifies prewarping. For example, % [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs,Fp) applies prewarping before % the bilinear transformation so that the frequency responses % before and after mapping match exactly at frequency point Fp % (match point Fp is specified in Hz). % % See also IMPINVAR. % Author(s): J.N. Little, 4-28-87 % J.N. Little, 5-5-87, revised % Copyright (c) 1988-98 by The MathWorks, Inc. % $Revision: 1.1 $ $Date: 1998/06/03 14:41:57 $ % Gene Franklin, Stanford Univ., motivated the state-space % approach to the bilinear transformation. [mn,nn] = size(z); [md,nd] = size(p); if (nd == 1 & nn < 2) & nargout ~= 4 % In zero-pole-gain form if mn > md error('Numerator cannot be higher order than denominator.') end if nargin == 5 % Prewarp fp = 2*pi*fp; fs = fp/tan(fp/fs/2); else fs = 2*fs; end z = z(finite(z)); % Strip infinities from zeros pd = (1+p/fs)./(1-p/fs); % Do bilinear transformation zd = (1+z/fs)./(1-z/fs); % real(kd) or just kd? kd = (k*prod(fs-z)./prod(fs-p)); zd = [zd;-ones(length(pd)-length(zd),1)]; % Add extra zeros at -1 elseif (md == 1 & mn == 1) | nargout == 4 % if nargout == 4 % State-space case a = z; b = p; c = k; d = fs; fs = fp; error(abcdchk(a,b,c,d)); if nargin == 6 % Prewarp fp = fp1; % Decode arguments fp = 2*pi*fp; fs = fp/tan(fp/fs/2)/2; end else % Transfer function case if nn > nd error('Numerator cannot be higher order than denominator.') end num = z; den = p; % Decode arguments if nargin == 4 % Prewarp fp = fs; fs = k; % Decode arguments fp = 2*pi*fp; fs = fp/tan(fp/fs/2)/2; else fs = k; % Decode arguments end % Put num(s)/den(s) in state-space canonical form. [a,b,c,d] = tf2ss(num,den); end % Now do state-space version of bilinear transformation: t = 1/fs; r = sqrt(t); t1 = eye(size(a)) + a*t/2; t2 = eye(size(a)) - a*t/2; ad = t2\t1; bd = t/r*(t2\b); cd = r*c/t2; dd = c/t2*b*t/2 + d; if nargout == 4 zd = ad; pd = bd; kd = cd; else % Convert back to transfer function form: p = poly(ad); zd = poly(ad-bd*cd)+(dd-1)*p; pd = p; end else error('First two arguments must have the same orientation.') end