diff options
Diffstat (limited to 'macros/iirlp2mb.sci')
| -rw-r--r-- | macros/iirlp2mb.sci | 378 |
1 files changed, 341 insertions, 37 deletions
diff --git a/macros/iirlp2mb.sci b/macros/iirlp2mb.sci index 54f1a6b..09ec743 100644 --- a/macros/iirlp2mb.sci +++ b/macros/iirlp2mb.sci @@ -1,41 +1,345 @@ +// 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) -//This function does IIR Low Pass Filter to Multiband Filter Transformation. -//Calling Sequence -//[Num, Den, AllpassNum, AllpassDen] = iirlp2mb(B, A, Wo, Wt) -//[Num, Den, AllpassNum, AllpassDen] = iirlp2mb(B, A, Wo, Wt, Pass) -//Parameters -//B: real or complex value -//A: real or complex value -//Wo: scalar or vector -//Wt: scalar or vector, elements must be monotonically increasing and >= 0 and <= 1. -//Description -//This is an Octave function. -//This function does IIR Low Pass Filter to Multiband Filter Transformation. -//The first two parameters give the numerator and denominator of the prototype low pass filter. -//The third parameter is the normalized angular frequency/pi to be transformed. -//The fourth parameter is the normalized angular frequency/pi target vector. -//The first two output variables are the numerator and denominator of the transformed filter. -//The third and fourth output variables are the numerator and denominator of the allpass transform. -//The fifth parameter can have values pass or stop, default value is pass. -//Examples -//iirlp2mb(5,9,0.3,0.4) -//ans = 0.55556 - -B = varargin(1) -A = varargin(2) -Wo = varargin(3) -Wt = varargin(4) -rhs = argn(2) -lhs = argn(1) -if(rhs < 4 | rhs > 5) -error("Wrong number of input arguments.") -end -select(rhs) -case 4 then -[Num,Den,AllpassNum,AllpassDen] = callOctave("iirlp2mb",varargin(1),varargin(2),varargin(3),varargin(4)) -case 5 then -[Num,Den,AllpassNum,AllpassDen] = callOctave("iirlp2mb",varargin(1),varargin(2),varargin(3),varargin(4),varargin(5)) -end + + 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") +*/
\ No newline at end of file |