X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;ds=sidebyside;f=octave_packages%2Fsignal-1.1.3%2Fsftrans.m;fp=octave_packages%2Fsignal-1.1.3%2Fsftrans.m;h=0c6dc7d5e10b2a99b63616a2d1828e7de8a347d4;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/signal-1.1.3/sftrans.m b/octave_packages/signal-1.1.3/sftrans.m new file mode 100644 index 0000000..0c6dc7d --- /dev/null +++ b/octave_packages/signal-1.1.3/sftrans.m @@ -0,0 +1,183 @@ +## 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; if not, see . + +## usage: [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, 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: +## +## 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 +## ---------------- ------------------------- ------------------------ +## +## 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: +## +## H(S)=g prod(S-Xi)/prod(S-Xj) +## +## 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. + +function [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop) + + if (nargin != 5) + print_usage; + end + + 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"); + end + + 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))]; + end + end + 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(1,p); + 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(1, (p-z))]; + end + end + end + 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(1,p); + else + Sz = [C * Fc ./ Sz]; + if (p > z) + Sz = [Sz, zeros(1,p-z)]; + end + end + 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; + end + end +endfunction