## Copyright (C) 1999-2001 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{Sz}, @var{Sp}, @var{Sg}] =} sftrans (@var{Sz}, @var{Sp}, @var{Sg}, @var{W}, @var{stop}) ## ## Transform band edges of a generic lowpass filter (cutoff at W=1) ## represented in splane zero-pole-gain form. W is the edge of the ## target filter (or edges if band pass or band stop). Stop is true for ## high pass and band stop filters or false for low pass and band pass ## filters. Filter edges are specified in radians, from 0 to pi (the ## nyquist frequency). ## ## Theory: Given a low pass filter represented by poles and zeros in the ## splane, you can convert it to a low pass, high pass, band pass or ## band stop by transforming each of the poles and zeros individually. ## The following table summarizes the transformation: ## ## @example ## Transform Zero at x Pole at x ## ---------------- ------------------------- ------------------------ ## Low Pass zero: Fc x/C pole: Fc x/C ## S -> C S/Fc gain: C/Fc gain: Fc/C ## ---------------- ------------------------- ------------------------ ## High Pass zero: Fc C/x pole: Fc C/x ## S -> C Fc/S pole: 0 zero: 0 ## gain: -x gain: -1/x ## ---------------- ------------------------- ------------------------ ## Band Pass zero: b +- sqrt(b^2-FhFl) pole: b +- sqrt(b^2-FhFl) ## S^2+FhFl pole: 0 zero: 0 ## S -> C -------- gain: C/(Fh-Fl) gain: (Fh-Fl)/C ## S(Fh-Fl) b=x/C (Fh-Fl)/2 b=x/C (Fh-Fl)/2 ## ---------------- ------------------------- ------------------------ ## Band Stop zero: b +- sqrt(b^2-FhFl) pole: b +- sqrt(b^2-FhFl) ## S(Fh-Fl) pole: +-sqrt(-FhFl) zero: +-sqrt(-FhFl) ## S -> C -------- gain: -x gain: -1/x ## S^2+FhFl b=C/x (Fh-Fl)/2 b=C/x (Fh-Fl)/2 ## ---------------- ------------------------- ------------------------ ## Bilinear zero: (2+xT)/(2-xT) pole: (2+xT)/(2-xT) ## 2 z-1 pole: -1 zero: -1 ## S -> - --- gain: (2-xT)/T gain: (2-xT)/T ## T z+1 ## ---------------- ------------------------- ------------------------ ## @end example ## ## where C is the cutoff frequency of the initial lowpass filter, Fc is ## the edge of the target low/high pass filter and [Fl,Fh] are the edges ## of the target band pass/stop filter. With abundant tedious algebra, ## you can derive the above formulae yourself by substituting the ## transform for S into H(S)=S-x for a zero at x or H(S)=1/(S-x) for a ## pole at x, and converting the result into the form: ## ## @example ## H(S)=g prod(S-Xi)/prod(S-Xj) ## @end example ## ## The transforms are from the references. The actual pole-zero-gain ## changes I derived myself. ## ## Please note that a pole and a zero at the same place exactly cancel. ## This is significant for High Pass, Band Pass and Band Stop filters ## which create numerous extra poles and zeros, most of which cancel. ## Those which do not cancel have a "fill-in" effect, extending the ## shorter of the sets to have the same number of as the longer of the ## sets of poles and zeros (or at least split the difference in the case ## of the band pass filter). There may be other opportunistic ## cancellations but I will not check for them. ## ## Also note that any pole on the unit circle or beyond will result in ## an unstable filter. Because of cancellation, this will only happen ## if the number of poles is smaller than the number of zeros and the ## filter is high pass or band pass. The analytic design methods all ## yield more poles than zeros, so this will not be a problem. ## ## References: ## ## Proakis & Manolakis (1992). Digital Signal Processing. New York: ## Macmillan Publishing Company. ## @end deftypefn function [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop) if (nargin != 5) print_usage; endif Sz = Sz(:); Sp = Sp(:); C = 1; p = length(Sp); z = length(Sz); if z > p || p == 0 error("sftrans: must have at least as many poles as zeros in s-plane"); endif if length(W)==2 Fl = W(1); Fh = W(2); if stop ## ---------------- ------------------------- ------------------------ ## Band Stop zero: b ± sqrt(b^2-FhFl) pole: b ± sqrt(b^2-FhFl) ## S(Fh-Fl) pole: ±sqrt(-FhFl) zero: ±sqrt(-FhFl) ## S -> C -------- gain: -x gain: -1/x ## S^2+FhFl b=C/x (Fh-Fl)/2 b=C/x (Fh-Fl)/2 ## ---------------- ------------------------- ------------------------ if (isempty(Sz)) Sg = Sg * real (1./ prod(-Sp)); elseif (isempty(Sp)) Sg = Sg * real(prod(-Sz)); else Sg = Sg * real(prod(-Sz)/prod(-Sp)); endif b = (C*(Fh-Fl)/2)./Sp; Sp = [b+sqrt(b.^2-Fh*Fl); b-sqrt(b.^2-Fh*Fl)]; extend = [sqrt(-Fh*Fl); -sqrt(-Fh*Fl)]; if isempty(Sz) Sz = [extend(1+rem([1:2*p],2))]; else b = (C*(Fh-Fl)/2)./Sz; Sz = [b+sqrt(b.^2-Fh*Fl); b-sqrt(b.^2-Fh*Fl)]; if (p > z) Sz = [Sz; extend(1+rem([1:2*(p-z)],2))]; endif endif else ## ---------------- ------------------------- ------------------------ ## Band Pass zero: b ± sqrt(b^2-FhFl) pole: b ± sqrt(b^2-FhFl) ## S^2+FhFl pole: 0 zero: 0 ## S -> C -------- gain: C/(Fh-Fl) gain: (Fh-Fl)/C ## S(Fh-Fl) b=x/C (Fh-Fl)/2 b=x/C (Fh-Fl)/2 ## ---------------- ------------------------- ------------------------ Sg = Sg * (C/(Fh-Fl))^(z-p); b = Sp*((Fh-Fl)/(2*C)); Sp = [b+sqrt(b.^2-Fh*Fl); b-sqrt(b.^2-Fh*Fl)]; if isempty(Sz) Sz = zeros(p,1); else b = Sz*((Fh-Fl)/(2*C)); Sz = [b+sqrt(b.^2-Fh*Fl); b-sqrt(b.^2-Fh*Fl)]; if (p>z) Sz = [Sz; zeros((p-z),1)]; endif endif endif else Fc = W; if stop ## ---------------- ------------------------- ------------------------ ## High Pass zero: Fc C/x pole: Fc C/x ## S -> C Fc/S pole: 0 zero: 0 ## gain: -x gain: -1/x ## ---------------- ------------------------- ------------------------ if (isempty(Sz)) Sg = Sg * real (1./ prod(-Sp)); elseif (isempty(Sp)) Sg = Sg * real(prod(-Sz)); else Sg = Sg * real(prod(-Sz)/prod(-Sp)); endif Sp = C * Fc ./ Sp; if isempty(Sz) Sz = zeros(p,1); else Sz = [C * Fc ./ Sz]; if (p > z) Sz = [Sz; zeros(p-z, 1)]; endif endif else ## ---------------- ------------------------- ------------------------ ## Low Pass zero: Fc x/C pole: Fc x/C ## S -> C S/Fc gain: C/Fc gain: Fc/C ## ---------------- ------------------------- ------------------------ Sg = Sg * (C/Fc)^(z-p); Sp = Fc * Sp / C; Sz = Fc * Sz / C; endif endif endfunction