path: root/macros/iirlp2mb.sci

// Copyright (C) 2018 - IIT Bombay - FOSSEE
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
// Original Source : https://octave.sourceforge.io/signal/
// Modifieded by: Abinash Singh Under FOSSEE Internship
// Last Modified on : 3 Feb 2024
// Organization: FOSSEE, IIT Bombay
// Email: toolbox@scilab.in
// IIR Low Pass Filter to Multiband Filter Transformation
//
// [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt)
// [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt,Pass)
//
// Num,Den:               numerator,denominator of the transformed filter
// AllpassNum,AllpassDen: numerator,denominator of allpass transform,
// B,A:                   numerator,denominator of prototype low pass filter
// Wo:                    normalized_angular_frequency/pi to be transformed
// Wt:                    [phi=normalized_angular_frequencies]/pi target vector
// Pass:                  This parameter may have values 'pass' or 'stop'.  If
//                        not given, it defaults to the value of 'pass'.
//
// With normalized ang. freq. targets 0 < phi(1) <  ... < phi(n) < pi radians
//
// for Pass == 'pass', the target multiband magnitude will be:
//       --------       ----------        -----------...
//      /        \     /          \      /            .
// 0   phi(1) phi(2)  phi(3)   phi(4)   phi(5)   (phi(6))    pi
//
// for Pass == 'stop', the target multiband magnitude will be:
// -------      ---------        ----------...
//        \    /         \      /           .
// 0   phi(1) phi(2)  phi(3)   phi(4)  (phi(5))              pi
//
function [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(varargin)

  usage = sprintf("iirlp2mb Usage: [Num,Den,AllpassNum,AllpassDen]=iirlp2mb(B,A,Wo,Wt[,Pass])\n");
  B = varargin(1);  // numerator polynomial of prototype low pass filter
  A = varargin(2);  // denominator polynomial of prototype low pass filter
  Wo = varargin(3); // (normalized angular frequency)/pi to be transformed
  Wt = varargin(4); // vector of (norm. angular frequency)/pi transform targets
                    // [phi(1) phi(2) ... ]/pi
  if(nargin < 4 || nargin > 5)
    error(usage)
  end

  // Checking proper input arguments
  if ~isvector(B) then 
    if ~isscalar(B) then // if numerator is having only one coeffcient then it can be passed directly
      error(sprintf("iirlp2mb : Argument B must be  a vector containing numerator polynomial  of the prototype low pass filter\n %s",usage));
    end
  end
  if ~isvector(A) then 
    if ~isscalar(A)then  // if Denominator is having only one coeffcient then it can be passed directly
      error(sprintf("iirlp2mb : A must be a  vector containing denominator polynomial of the prototype low pass filter\n %s",usage));
    end
  end
  if(nargin == 5)
    Pass = varargin(5);
    switch(Pass)
      case 'pass'
        pass_stop = -1;
      case 'stop'
        pass_stop = 1;
      otherwise
        error(sprintf("Pass must be pass or stop\n%s",usage))
    end
  else
    pass_stop = -1; // Pass == 'pass' is the default
  end
  if(abs(Wo) <= 0)
    error(sprintf("Wo is %f <= 0\n%s",abs(Wo),usage));
  end
  if(abs(Wo) >= 1)
    error(sprintf("Wo is %f >= 1\n%s",abs(Wo),usage));
  end
  oWt = 0;
  for i = 1 : length(Wt)
    if(abs(Wt(i)) <= 0)
      error(sprintf("Wt(%d) is %f <= 0\n%s",i,abs(Wt(i)),usage));
    end
    if(abs(Wt(i)) >= 1)
      error(sprintf("Wt(%d) is %f >= 1\n%s",i,abs(Wt(i)),usage));
    end
    if(abs(Wt(i)) <= oWt)
      error(sprintf("Wt(%d) = %f, not monotonically increasing\n%s",i,abs(Wt(i)),usage));
    else
      oWt = Wt(i);
    end
  end

  //                                                             B(z)
  // Inputs B,A specify the low pass IIR prototype filter G(z) = ---- .
  //                                                             A(z)
  // This module transforms G(z) into a multiband filter using the iterative
  // algorithm from:
  // [FFM] G. Feyh, J. Franchitti, and C. Mullis, "All-Pass Filter
  // Interpolation and Frequency Transformation Problem", Proceedings 20th
  // Asilomar Conference on Signals, Systems and Computers, Nov. 1986, pp.
  // 164-168, IEEE.
  // [FFM] moves the prototype filter position at normalized angular frequency
  // .5*pi to the places specified in the Wt vector times pi.  In this module,
  // a generalization allows the position to be moved on the prototype filter
  // to be specified as Wo*pi instead of being fixed at .5*pi.  This is
  // implemented using two successive allpass transformations.
  //                                         KK(z)
  // In the first stage, find allpass J(z) = ----  such that
  //                                         K(z)
  //    jWo*pi     -j.5*pi
  // J(e      ) = e                    (low pass to low pass transformation)
  //
  //                                          PP(z)
  // In the second stage, find allpass H(z) = ----  such that
  //                                          P(z)
  //    jWt(k)*pi     -j(2k - 1)*.5*pi
  // H(e         ) = e                 (low pass to multiband transformation)
  //
  //                                          ^
  // The variable PP used here corresponds to P in [FFM].
  // len = length(P(z)) == length(PP(z)), the number of polynomial coefficients
  //
  //        len      1-i           len       1-i
  // P(z) = SUM P(i)z   ;  PP(z) = SUM PP(i)z   ; PP(i) == P(len + 1 - i)
  //        i=1                    i=1              (allpass condition)
  // Note: (len - 1) == n in [FFM] eq. 3
  //
  // The first stage computes the denominator of an allpass for translating
  // from a prototype with position .5 to one with a position of Wo. It has the
  // form:
  //          -1
  // K(2)  - z
  // -----------
  //          -1
  // 1 - K(2)z
  //
  // From the low pass to low pass transformation in Table 7.1 p. 529 of A.
  // Oppenheim and R. Schafer, Discrete-Time Signal Processing 3rd edition,
  // Prentice Hall 2010, one can see that the denominator of an allpass for
  // going in the opposite direction can be obtained by a sign reversal of the
  // second coefficient, K(2), of the vector K (the index 2 not to be confused
  // with a value of z, which is implicit).

  // The first stage allpass denominator computation
  K = apd([%pi * Wo]);

  // The second stage allpass computation
  phi = %pi * Wt; // vector of normalized angular frequencies between 0 and pi
  P = apd(phi);  // calculate denominator of allpass for this target vector
  PP = flipdim(P,2); // numerator of allpass has reversed coefficients of P

  // The total allpass filter from the two consecutive stages can be written as
  //          PP
  // K(2) -  ---
  //          P         P
  // -----------   *   ---
  //          PP        P
  // 1 - K(2)---
  //          P
  AllpassDen = P - (K(2) * PP);
  AllpassDen = AllpassDen / AllpassDen(1); // normalize
  AllpassNum = pass_stop * flipdim(AllpassDen,2);
  [Num,Den] = transform(B,A,AllpassNum,AllpassDen,pass_stop);

endfunction

function [Num,Den] = transform(B,A,PP,P,pass_stop)

  // Given G(Z) = B(Z)/A(Z) and allpass H(z) = PP(z)/P(z), compute G(H(z))
  // For Pass = 'pass', transformed filter is:
  //                          2                   nb-1
  // B1 + B2(PP/P) + B3(PP/P)^  + ... + Bnb(PP/P)^
  // -------------------------------------------------
  //                          2                   na-1
  // A1 + A2(PP/P) + A3(PP/P)^  + ... + Ana(PP/P)^
  // For Pass = 'stop', use powers of (-PP/P)
  //
  na = length(A);  // the number of coefficients in A
  nb = length(B);  // the number of coefficients in B
  // common low pass iir filters have na == nb but in general might not
  n  = max(na,nb); // the greater of the number of coefficients
  //                              n-1
  // Multiply top and bottom by P^   yields:
  //
  //      n-1             n-2          2    n-3                 nb-1    n-nb
  // B1(P^   ) + B2(PP)(P^   ) + B3(PP^ )(P^   ) + ... + Bnb(PP^    )(P^    )
  // ---------------------------------------------------------------------
  //      n-1             n-2          2    n-3                 na-1    n-na
  // A1(P^   ) + A2(PP)(P^   ) + A3(PP^ )(P^   ) + ... + Ana(PP^    )(P^    )

  // Compute and store powers of P as a matrix of coefficients because we will
  // need to use them in descending power order
  global Ppower; // to hold coefficients of powers of P, access inside ppower()
  np = length(P);
  powcols = np + (np-1)*(n-2); // number of coefficients in P^(n-1)
  // initialize to "Not Available" with n-1 rows for powers 1 to (n-1) and
  // the number of columns needed to hold coefficients for P^(n-1)
  Ppower = %nan * ones(n-1,powcols);
  Ptemp = P;                   // start with P to the 1st power
  for i = 1 : n-1              // i is the power
    for j = 1 : length(Ptemp) // j is the coefficient index for this power
      Ppower(i,j)  = Ptemp(j);
    end
    Ptemp = conv(Ptemp,P);    // increase power of P by one
  end

  // Compute numerator and denominator of transformed filter
  Num = [];
  Den = [];
  for i = 1 : n
    //              n-i
    // Regenerate P^    (p_pownmi)
    if((n-i) == 0)
      p_pownmi = [1];
    else
      p_pownmi = ppower(n-i,powcols);
    end
    //               i-1
    // Regenerate PP^   (pp_powim1)
    if(i == 1)
      pp_powim1 = [1];
    else
      pp_powim1 = flipdim(ppower(i-1,powcols),2);
    end
    if(i <= nb)
      Bterm = (pass_stop^(i-1))*B(i)*conv(pp_powim1,p_pownmi);
      Num = polysum(Num,Bterm);
    end
    if(i <= na)
      Aterm = (pass_stop^(i-1))*A(i)*conv(pp_powim1,p_pownmi);
      Den = polysum(Den,Aterm);
    end
  end
  // Scale both numerator and denominator to have Den(1) = 1
  temp = Den(1);
  for i = 1 : length(Den)
    Den(i) = Den(i) / temp;
  end
  for i = 1 : length(Num)
    Num(i) = Num(i) / temp;
  end

endfunction

function P = apd(phi) // all pass denominator
  // Given phi, a vector of normalized angular frequency transformation targets,
  // return P, the denominator of an allpass H(z)
  lenphi = length(phi);
  Pkm1 = 1; // P0 initial condition from [FFM] eq. 22
  for k = 1:lenphi
    P = pk(Pkm1, k, phi(k)); // iterate
    Pkm1 = P;
  end

endfunction

function Pk = pk(Pkm1, k, phik) // kth iteration of P(z)

  // Given Pkminus1, k, and phi(k) in radians , return Pk
  //
  // From [FFM] eq. 19 :                     k
  // Pk =     (z+1  )sin(phi(k)/2)Pkm1 - (-1) (z-1  )cos(phi(k)/2)PPkm1
  // Factoring out z
  //              -1                         k    -1
  //    =   z((1+z  )sin(phi(k)/2)Pkm1 - (-1) (1-z  )cos(phi(k)/2)PPkm1)
  // PPk can also have z factored out.  In H=PP/P, z in PPk will cancel z in Pk,
  // so just leave out.  Use
  //              -1                         k    -1
  // PK =     (1+z  )sin(phi(k)/2)Pkm1 - (-1) (1-z  )cos(phi(k)/2)PPkm1
  // (expand)                                k
  //    =            sin(phi(k)/2)Pkm1 - (-1)        cos(phi(k)/2)PPkm1
  //
  //              -1                         k   -1
  //           + z   sin(phi(k)/2)Pkm1 + (-1)    z   cos(phi(k)/2)PPkm1
  Pk = zeros(1,k+1); // there are k+1 coefficients in Pk
  sin_k = sin(phik/2);
  cos_k = cos(phik/2);
  for i = 1 : k
    Pk(i)   = Pk(i)+ sin_k * Pkm1(i) - ((-1)^k * cos_k * Pkm1(k+1-i));
    //
    //                    -1
    // Multiplication by z   just shifts by one coefficient
    Pk(i+1) = Pk(i+1)+ sin_k * Pkm1(i) + ((-1)^k * cos_k * Pkm1(k+1-i));
  end
  // now normalize to Pk(1) = 1 (again will cancel with same factor in PPk)
  Pk1 = Pk(1);
  for i = 1 : k+1
    Pk(i) = Pk(i) / Pk1;
  end

endfunction
function p = ppower(i,powcols) // Regenerate ith power of P from stored PPower

  global Ppower
  if(i == 0)
    p  = 1;
  else
    p  = [];
    for j = 1 : powcols
      if(isnan(Ppower(i,j)))
        break;
      end
      p =  [p, Ppower(i,j)];
    end
  end

endfunction

function poly = polysum(p1,p2) // add polynomials of possibly different length

  n1 = length(p1);
  n2 = length(p2);
  if(n1 > n2)
    // pad p2
    p2 = [p2, zeros(1,n1-n2)];
  elseif(n2 > n1)
    // pad p1
    p1 = [p1, zeros(1,n2-n1)];
  end
  poly = p1 + p2;

endfunction
 
/*
test passed
// Butterworth filter of order 3 with cutoff frequency 0.5 
B = [0.1667   0.5000   0.5000   0.1667];
A = [1.0000e+00  -3.0531e-16   3.3333e-01  -1.8504e-17 ] ;
[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B, A, 0.5, [.2 .4 .6 .8]) // 5th argument default pass
[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B, A,[.2 .4 .6 .8],0.2)   // 5th argument default pass
[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B, A, 0.5, [.2 .4 .6 .8],"stop")
[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B, A,[.2 .4 .6 .8],0.5,"stop")

test passed
// scalar input for B and A
[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(0, 0, 0.5, [.2 .4 .6 .8])
[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(1, 2, [.2 .4 .6 .8],0.2,"stop")


test passed 
// complex B and A 
[Num,Den,AllpassNum,AllpassDen] = iirlp2mb([%i,2+%i],[1 2*%i 3], [.2 .4 .6 .8],0.2,"stop")

*/