%% Copyright (C) 1999 Paul Kienzle %% %% This program is free software: you can redistribute it and/or modify %% it under the terms of the GNU General Public License as published by %% the Free Software Foundation, either version 3 of the License, or %% (at your option) any later version. %% %% This program is distributed in the hope that it will be useful, %% but WITHOUT ANY WARRANTY; without even the implied warranty of %% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the %% GNU General Public License for more details. %% %% You should have received a copy of the GNU General Public License %% along with this program; see the file COPYING. If not, see %% . %% -*- texinfo -*- %% @deftypefn {Function File} {[@var{numd}, @var{dend}] =} bilinear2 (@var{num}, @var{den}, @var{fs}) %% @deftypefnx {Function File} {[@var{zd}, @var{pd}, @var{kd}] =} bilinear2 (@var{z}, @var{p}, @var{k}, @var{fs}) %% @deftypefnx {Function File} {[@var{Ad}, @var{Bd}, @var{Cd}, @var{Dd}] =} bilinear2 (@var{A}, @var{B}, @var{C}, ,@var{D}, @var{fs}) %% @deftypefnx {Function File} {[__] =} bilinear2 (__, @var{fp}) %% %% Transform a s-plane filter specification into a z-plane specification. Filters %% can be specified in either transfer function, zero-pole-gain or %% state-space form. fs is the sampling frequency. The %% fp parameter does prewarping to match fp in the s and z planes %% %% blinear2 differs from the bilinear function in the signal processing %% toolbox, which uses T as an input argument rather than fs (1/T). bilinear2 %% is compatible with the MATLAB function. The bilinear function does the %% transformation on the poles and zeros while bilinear2 does the transformation %% on the coefficients in the transfer function representation. %% %% Demos are included, try "demo bilinear2" %% %% @end deftypefn %% Author: Tony Richardson %% Created: 2021-04-11 function [Ad, Bd, Cd, Dd] = bilinear2(Ass, Bss, Css, Dss, xarg1, xarg2) % Transform all representations to transfer function (coefficient form) prewarp = false; if nargout<=2 if nargin==3 B = Ass; A = Bss; fs = Css; elseif nargin==4 B = Ass; A = Bss; fs = Css; fp = Dss; prewarp = true; else print_usage; endif elseif nargout==3 if nargin==4 [B, A] = zp2tf(Ass, Bss, Css); fs = Dss; elseif nargin==5 [B, A] = zp2tf(Ass, Bss, Css); fs = Dss; fp = xarg1; prewarp = true; else print_usage; endif elseif nargout==4 if nargin==5 [B, A] = ss2tf(Ass, Bss, Css, Dss); fs = xarg1; elseif nargin==6 [B, A] = ss2tf(Ass, Bss, Css, Dss); fs = xarg1; fp = xarg2; prewarp = true; else print_usage; endif else print_usage; endif % Initialization B = B(:)'; A = A(:)'; if prewarp C = 2*pi*fp/tan(pi*fp/fs); else C = 2*fs; endif M = length(B); N = length(A); % Numerator and denominator of BLT biltn = [1 -1]; biltd = [1 1]; % Transform numerator (B) coefficients if M > 1 Bt = zeros(1, M); for n = 1:M zp = 1; for k = 1:(n-1) zp = conv(zp, biltd); end for k = 1:(M-n) zp = conv(zp, C*biltn); end Bt = Bt + B(n)*zp; end else Bt = B(1); end for k = 1:(N-M) Bt = conv(Bt, biltd); end % Transform denominator (A) coefficients if N > 1 At = zeros(1, N); for n = 1:N zp = 1; for k = 1:(n-1) zp = conv(zp, biltd); end for k = 1:(N-n) zp = conv(zp, C*biltn); end At = At + A(n)*zp; end else At = A(1); end for k = 1:(M-N) At = conv(At, biltd); end %Normalize the coefficient vectors Bt = Bt/At(1); At = At/At(1); % Convert transfer function to desired form if nargout==3 % Zero-pole [Ad, Bd, Cd] = tf2zp(Bt, At); elseif nargout==4 % State-space [Ad, Bd, Cd, Dd] = tf2ss(Bt, At); else % Transfer function Ad = Bt; Bd = At; end end %!demo %! %! [z,p,k] = cheb1ap(10,3); %! [A,B,C,D] = zp2ss(z,p,real(k)); %! %! Fs = 2e3; %! u1 = 2*Fs*tan(100*(2*pi/Fs)/2); %! u2 = 2*Fs*tan(500*(2*pi/Fs)/2); %! Bw = u2 - u1; %! Wo = sqrt(u1*u2); %! [At,Bt,Ct,Dt] = lp2bp(A,B,C,D,Wo,Bw); %! [b,a] = ss2tf(At,Bt,Ct,Dt); %! %! figure(1) %! clf %! w = 0:50:10000; %! [h] = freqs(b,a,w); %! plot(w,mag2db(abs(h))) %! hold on %! ylim([-165 5]) %! [U1,U2] = meshgrid([u1 u2],ylim); %! plot(U1,U2) %! legend('Magnitude response','Lower Passband Edge','Upper Passband Edge') %! hold off %! xlabel('Angular Frequency (rad/s)') %! ylabel('Magnitude (dB)') %! grid %! %! [Ad,Bd,Cd,Dd] = bilinear2(At,Bt,Ct,Dt,Fs); %! [bz,az] = ss2tf(Ad,Bd,Cd,Dd); %! [hz, fz] = freqz(bz, az, 512, Fs); %! %! figure(2) %! clf %! plot(fz, mag2db(abs(hz))) %! ylim([-260 10]) %! grid on %! xlabel('Frequency (Hz)') %! ylabel('Magnitude (dB)') %! title('Magnitude Response (dB)') %!demo %! %! Fc = 20; %! Fs = 200; %! [z,p,k] = ellip(6,3,90,2*pi*Fc,'s'); %! [num,den] = zp2tf(z,p,k); %! %! f = 0:0.5:50; %! [h] = freqs(num,den,2*pi*f); %! figure(1) %! clf %! plot(f,mag2db(abs(h))) %! hold on %! xlim([0 50]) %! [l1,l2] = meshgrid(Fc,[-120 0]); %! plot(l1,l2) %! grid %! legend('Magnitude response','Passband Edge') %! xlabel('Frequency (Hz)') %! ylabel('Magnitude (dB)') %! %! [numd,dend] = bilinear2(num,den,Fs,20); %! [hz, fz] = freqz(numd, dend, 512, Fs); %! %! figure(2) %! clf %! plot(fz, mag2db(abs(hz))) %! ylim([-110 5]) %! grid on %! xlabel('Frequency (Hz)') %! ylabel('Magnitude (dB)') %! title('Magnitude Response (dB)')