octave_packages/signal-1.1.3/iirlp2mb.m

   1 ## Copyright (C) 2011 Alan J. Greenberger <alanjg@ptd.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## IIR Low Pass Filter to Multiband Filter Transformation
  17 ##
  18 ## [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt)
  19 ## [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt,Pass)
  20 ##
  21 ## Num,Den:               numerator,denominator of the transformed filter
  22 ## AllpassNum,AllpassDen: numerator,denominator of allpass transform,
  23 ## B,A:                   numerator,denominator of prototype low pass filter
  24 ## Wo:                    normalized_angular_frequency/pi to be transformed
  25 ## Wt:                    [phi=normalized_angular_frequencies]/pi target vector
  26 ## Pass:                  This parameter may have values 'pass' or 'stop'.  If
  27 ##                        not given, it defaults to the value of 'pass'.
  28 ##
  29 ## With normalized ang. freq. targets 0 < phi(1) <  ... < phi(n) < pi radians
  30 ##
  31 ## for Pass == 'pass', the target multiband magnitude will be:
  32 ##       --------       ----------        -----------...
  33 ##      /        \     /          \      /            .
  34 ## 0   phi(1) phi(2)  phi(3)   phi(4)   phi(5)   (phi(6))    pi
  35 ##
  36 ## for Pass == 'stop', the target multiband magnitude will be:
  37 ## -------      ---------        ----------...
  38 ##        \    /         \      /           .
  39 ## 0   phi(1) phi(2)  phi(3)   phi(4)  (phi(5))              pi
  40 ##
  41 ## Example of use:
  42 ## [B, A] = butter(6, 0.5);
  43 ## [Num, Den] = iirlp2mb(B, A, 0.5, [.2 .4 .6 .8]);
  44
  45 function [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(varargin)
  46    usage = sprintf(
  47     "%s: Usage: [Num,Den,AllpassNum,AllpassDen]=iirlp2mb(B,A,Wo,Wt[,Pass])\n"
  48     ,mfilename());
  49    B = varargin{1};  # numerator polynomial of prototype low pass filter
  50    A = varargin{2};  # denominator polynomial of prototype low pass filter
  51    Wo = varargin{3}; # (normalized angular frequency)/pi to be transformed
  52    Wt = varargin{4}; # vector of (norm. angular frequency)/pi transform targets
  53                      # [phi(1) phi(2) ... ]/pi
  54    if(nargin < 4 || nargin > 5)
  55       error("%s",usage)
  56    endif
  57    if(nargin == 5)
  58       Pass = varargin{5};
  59       switch(Pass)
  60          case 'pass'
  61             pass_stop = -1;
  62          case 'stop'
  63             pass_stop = 1;
  64          otherwise
  65             error("Pass must be 'pass' or 'stop'\n%s",usage)
  66       endswitch
  67    else
  68       pass_stop = -1; # Pass == 'pass' is the default
  69    endif
  70    if(Wo <= 0)
  71       error("Wo is %f <= 0\n%s",Wo,usage);
  72    endif
  73    if(Wo >= 1)
  74       error("Wo is %f >= 1\n%s",Wo,usage);
  75    endif
  76    oWt = 0;
  77    for i = 1 : length(Wt)
  78       if(Wt(i) <= 0)
  79          error("Wt(%d) is %f <= 0\n%s",i,Wt(i),usage);
  80       endif
  81       if(Wt(i) >= 1)
  82          error("Wt(%d) is %f >= 1\n%s",i,Wt(i),usage);
  83       endif
  84       if(Wt(i) <= oWt)
  85          error("Wt(%d) = %f, not monotonically increasing\n%s",i,Wt(i),usage);
  86       else
  87          oWt = Wt(i);
  88       endif
  89    endfor
  90
  91    #                                                             B(z)
  92    # Inputs B,A specify the low pass IIR prototype filter G(z) = ---- .
  93    #                                                             A(z)
  94    # This module transforms G(z) into a multiband filter using the iterative
  95    # algorithm from:
  96    # [FFM] G. Feyh, J. Franchitti, and C. Mullis, "All-Pass Filter
  97    # Interpolation and Frequency Transformation Problem", Proceedings 20th
  98    # Asilomar Conference on Signals, Systems and Computers, Nov. 1986, pp.
  99    # 164-168, IEEE.
 100    # [FFM] moves the prototype filter position at normalized angular frequency
 101    # .5*pi to the places specified in the Wt vector times pi.  In this module,
 102    # a generalization allows the position to be moved on the prototype filter
 103    # to be specified as Wo*pi instead of being fixed at .5*pi.  This is
 104    # implemented using two successive allpass transformations.
 105    #                                         KK(z)
 106    # In the first stage, find allpass J(z) = ----  such that
 107    #                                         K(z)
 108    #    jWo*pi     -j.5*pi
 109    # J(e      ) = e                    (low pass to low pass transformation)
 110    #
 111    #                                          PP(z)
 112    # In the second stage, find allpass H(z) = ----  such that
 113    #                                          P(z)
 114    #    jWt(k)*pi     -j(2k - 1)*.5*pi
 115    # H(e         ) = e                 (low pass to multiband transformation)
 116    #
 117    #                                          ^
 118    # The variable PP used here corresponds to P in [FFM].
 119    # len = length(P(z)) == length(PP(z)), the number of polynomial coefficients
 120    #
 121    #        len      1-i           len       1-i
 122    # P(z) = SUM P(i)z   ;  PP(z) = SUM PP(i)z   ; PP(i) == P(len + 1 - i)
 123    #        i=1                    i=1              (allpass condition)
 124    # Note: (len - 1) == n in [FFM] eq. 3
 125    #
 126    # The first stage computes the denominator of an allpass for translating
 127    # from a prototype with position .5 to one with a position of Wo. It has the
 128    # form:
 129    #          -1
 130    # K(2)  - z
 131    # -----------
 132    #          -1
 133    # 1 - K(2)z
 134    #
 135    # From the low pass to low pass tranformation in Table 7.1 p. 529 of A.
 136    # Oppenheim and R. Schafer, Discrete-Time Signal Processing 3rd edition,
 137    # Prentice Hall 2010, one can see that the denominator of an allpass for
 138    # going in the opposite direction can be obtained by a sign reversal of the
 139    # second coefficient, K(2), of the vector K (the index 2 not to be confused
 140    # with a value of z, which is implicit).
 141
 142    # The first stage allpass denominator computation
 143    K = apd([pi * Wo]);
 144
 145    # The second stage allpass computation
 146    phi = pi * Wt; # vector of normalized angular frequencies between 0 and pi
 147    P = apd(phi);  # calculate denominator of allpass for this target vector
 148    PP = revco(P); # numerator of allpass has reversed coefficients of P
 149
 150    # The total allpass filter from the two consecutive stages can be written as
 151    #          PP
 152    # K(2) -  ---
 153    #          P         P
 154    # -----------   *   ---
 155    #          PP        P
 156    # 1 - K(2)---
 157    #          P
 158    AllpassDen = P - (K(2) * PP);
 159    AllpassDen /= AllpassDen(1); # normalize
 160    AllpassNum = pass_stop * revco(AllpassDen);
 161    [Num,Den] = transform(B,A,AllpassNum,AllpassDen,pass_stop);
 162 endfunction
 163
 164 function [Num,Den] = transform(B,A,PP,P,pass_stop)
 165    # Given G(Z) = B(Z)/A(Z) and allpass H(z) = PP(z)/P(z), compute G(H(z))
 166    # For Pass = 'pass', transformed filter is:
 167    #                          2                   nb-1
 168    # B1 + B2(PP/P) + B3(PP/P)^  + ... + Bnb(PP/P)^
 169    # -------------------------------------------------
 170    #                          2                   na-1
 171    # A1 + A2(PP/P) + A3(PP/P)^  + ... + Ana(PP/P)^
 172    # For Pass = 'stop', use powers of (-PP/P)
 173    #
 174    na = length(A);  # the number of coefficients in A
 175    nb = length(B);  # the number of coefficients in B
 176    # common low pass iir filters have na == nb but in general might not
 177    n  = max(na,nb); # the greater of the number of coefficients
 178    #                              n-1
 179    # Multiply top and bottom by P^   yields:
 180    #
 181    #      n-1             n-2          2    n-3                 nb-1    n-nb
 182    # B1(P^   ) + B2(PP)(P^   ) + B3(PP^ )(P^   ) + ... + Bnb(PP^    )(P^    )
 183    # ---------------------------------------------------------------------
 184    #      n-1             n-2          2    n-3                 na-1    n-na
 185    # A1(P^   ) + A2(PP)(P^   ) + A3(PP^ )(P^   ) + ... + Ana(PP^    )(P^    )
 186
 187    # Compute and store powers of P as a matrix of coefficients because we will
 188    # need to use them in descending power order
 189    global Ppower; # to hold coefficients of powers of P, access inside ppower()
 190    np = length(P);
 191    powcols = np + (np-1)*(n-2); # number of coefficients in P^(n-1)
 192    # initialize to "Not Available" with n-1 rows for powers 1 to (n-1) and
 193    # the number of columns needed to hold coefficients for P^(n-1)
 194    Ppower = NA(n-1,powcols);
 195    Ptemp = P;                   # start with P to the 1st power
 196    for i = 1 : n-1              # i is the power
 197       for j = 1 : length(Ptemp) # j is the coefficient index for this power
 198          Ppower(i,j)  = Ptemp(j);
 199       endfor
 200       Ptemp = conv(Ptemp,P);    # increase power of P by one
 201    endfor
 202
 203    # Compute numerator and denominator of transformed filter
 204    Num = [];
 205    Den = [];
 206    for i = 1 : n
 207       #              n-i
 208       # Regenerate P^    (p_pownmi)
 209       if((n-i) == 0)
 210          p_pownmi = [1];
 211       else
 212          p_pownmi = ppower(n-i,powcols);
 213       endif
 214       #               i-1
 215       # Regenerate PP^   (pp_powim1)
 216       if(i == 1)
 217          pp_powim1 = [1];
 218       else
 219          pp_powim1 = revco(ppower(i-1,powcols));
 220       endif
 221       if(i <= nb)
 222          Bterm = (pass_stop^(i-1))*B(i)*conv(pp_powim1,p_pownmi);
 223          Num = polysum(Num,Bterm);
 224       endif
 225       if(i <= na)
 226          Aterm = (pass_stop^(i-1))*A(i)*conv(pp_powim1,p_pownmi);
 227          Den = polysum(Den,Aterm);
 228       endif
 229    endfor
 230    # Scale both numerator and denominator to have Den(1) = 1
 231    temp = Den(1);
 232    for i = 1 : length(Den)
 233       Den(i) = Den(i) / temp;
 234    endfor
 235    for i = 1 : length(Num)
 236       Num(i) = Num(i) / temp;
 237    endfor
 238 endfunction
 239
 240 function P = apd(phi) # all pass denominator
 241    # Given phi, a vector of normalized angular frequency transformation targets,
 242    # return P, the denominator of an allpass H(z)
 243    lenphi = length(phi);
 244    Pkm1 = 1; # P0 initial condition from [FFM] eq. 22
 245    for k = 1 : lenphi
 246       P = pk(Pkm1, k, phi(k)); # iterate
 247       Pkm1 = P;
 248    endfor
 249 endfunction
 250
 251 function Pk = pk(Pkm1, k, phik) # kth iteration of P(z)
 252    # Given Pkminus1, k, and phi(k) in radians , return Pk
 253    #
 254    # From [FFM] eq. 19 :                     k
 255    # Pk =     (z+1  )sin(phi(k)/2)Pkm1 - (-1) (z-1  )cos(phi(k)/2)PPkm1
 256    # Factoring out z
 257    #              -1                         k    -1
 258    #    =   z((1+z  )sin(phi(k)/2)Pkm1 - (-1) (1-z  )cos(phi(k)/2)PPkm1)
 259    # PPk can also have z factored out.  In H=PP/P, z in PPk will cancel z in Pk,
 260    # so just leave out.  Use
 261    #              -1                         k    -1
 262    # PK =     (1+z  )sin(phi(k)/2)Pkm1 - (-1) (1-z  )cos(phi(k)/2)PPkm1
 263    # (expand)                                k
 264    #    =            sin(phi(k)/2)Pkm1 - (-1)        cos(phi(k)/2)PPkm1
 265    #
 266    #              -1                         k   -1
 267    #           + z   sin(phi(k)/2)Pkm1 + (-1)    z   cos(phi(k)/2)PPkm1
 268    Pk = zeros(1,k+1); # there are k+1 coefficients in Pk
 269    sin_k = sin(phik/2);
 270    cos_k = cos(phik/2);
 271    for i = 1 : k
 272       Pk(i)   += sin_k * Pkm1(i) - ((-1)^k * cos_k * Pkm1(k+1-i));
 273       #
 274       #                    -1
 275       # Multiplication by z   just shifts by one coefficient
 276       Pk(i+1) += sin_k * Pkm1(i) + ((-1)^k * cos_k * Pkm1(k+1-i));
 277    endfor
 278    # now normalize to Pk(1) = 1 (again will cancel with same factor in PPk)
 279    Pk1 = Pk(1);
 280    for i = 1 : k+1
 281       Pk(i) = Pk(i) / Pk1;
 282    endfor
 283 endfunction
 284
 285 function PP = revco(p) # reverse components of vector
 286    l = length(p);
 287    for i = 1 : l
 288       PP(l + 1 - i) = p(i);
 289    endfor
 290 endfunction
 291
 292 function p = ppower(i,powcols) # Regenerate ith power of P from stored PPower
 293    global Ppower
 294    if(i == 0)
 295       p  = 1;
 296    else
 297       p  = [];
 298       for j = 1 : powcols
 299          if(isna(Ppower(i,j)))
 300             break;
 301          endif
 302          p =  horzcat(p, Ppower(i,j));
 303       endfor
 304    endif
 305 endfunction
 306
 307 function poly = polysum(p1,p2) # add polynomials of possibly different length
 308    n1 = length(p1);
 309    n2 = length(p2);
 310    if(n1 > n2)
 311       # pad p2
 312       p2 = horzcat(p2, zeros(1,n1-n2));
 313    elseif(n2 > n1)
 314       # pad p1
 315       p1 = horzcat(p1, zeros(1,n2-n1));
 316    endif
 317    poly = p1 + p2;
 318 endfunction