I am requesting registration as an octave developer. My sourceforge username
is ajgreenberger.
I needed to design a multiband filter, so I looked at the literature and
programmed the Feyh, Franchitti, Mullis [FFM] algorithm in Octave. I wish to
submit this but have some issues:
1 There is a matlab function iirlp2mb for this described in
http://www.mathworks.com/help/toolbox/filterdesign/ref/iirlp2mb.html [matl].
It has four forms of invocation:
[Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt)
[Num,Den,AllpassNum,AllpassDen]=iirlp2mb(B,A,Wo,Wt,Pass)
[G,AllpassNum,AllpassDen] = iirlp2mb(Hd,Wo,Wt)
[G,AllpassNum,AllpassDen] = iirlp2mb(...,Pass)
I have implemented the first two. I don't understand the last two, which are
related to dfilt. I don't think Octave supports dfilt. Should I call my
implementation iirlp2mb or name it something else?
2 [matl] says, "Wo | Frequency value to be transformed from the prototype filter". I
interpreted this to be the cutoff frequency of the prototype low pass filter defined by B,A. In
the [FFM] algorithm, the low pass cutoff must be .5 . So if Wo is not "close" to .5, I
perform an initial transformation of the low pass to a low pass with cutoff of .5 . If I run either
[b, a] = ellip(3, 0.1, 30, 0.409);
[num2, den2] = iirlp2mb(b, a, 0.409, [2 4 6 8]/10, 'pass');
[num2, den2] = iirlp2mb(b, a, 0.409, [2 4 6 8]/10, 'stop');
or
[b, a] = ellip(3, 0.1, 30, 0.5);
[num2, den2] = iirlp2mb(b, a, 0.5, [2 4 6 8]/10, 'pass');
[num2, den2] = iirlp2mb(b, a, 0.5, [2 4 6 8]/10, 'stop');
I get figures similar to that shown in [matl]. But if I run their example
[b, a] = ellip(3, 0.1, 30, 0.409);
[num2, den2] = iirlp2mb(b, a, 0.5, [2 4 6 8]/10, 'pass');
[num2, den2] = iirlp2mb(b, a, 0.5, [2 4 6 8]/10, 'stop');
I get nonsense. Is there someone in this group that has access to a matlab
license and is willing to run some tests of mine against theirs? Is there
someone who can explain the real meaning of Wo?
I have attached an svn diff of the code. I am a novice at this and am
expecting to have to do modifications until you think it is ready.
Index: main/signal/inst/iirlp2mb.m
===================================================================
--- main/signal/inst/iirlp2mb.m (revision 0)
+++ main/signal/inst/iirlp2mb.m (revision 0)
@@ -0,0 +1,283 @@
+## Copyright (C) 2011 Alan J. Greenberger <ala...@ptd.net>
+##
+## 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 <http://www.gnu.org/licenses/>.
+
+## IIR Low Pass Filter to Multiband Filter Transformation
+##
+## [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt) uses Pass='pass'
+## [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt,Pass) 'stop' or 'pass'
+##
+## Num,Den numerator,denominator of the transformed filter
+## AllpassNum,AllpassDen numerator,denominator of transforming allpass filter
+## B,A numerator,denominator of prototype low pass filter
+## Wo (nomalized angular frequency)/pi of prototype cutoff (preferably .5)
+## Wt vector of (norm. angular frequency)/pi transformed edge targets
+
+function [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(varargin)
+ usage = sprintf(
+ "%s: Usage: [Num,Den,AllpassNum,AllpassDen]=iirlp2mb(B,A,Wo,Wt[,Pass])\n"
+ ,mfilename());
+ 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
+ # Algorithm used assumes one starts with a low pass filter with a cutoff
+ # of .5, e.g. [b,a] = butter(4, .5) . This implementation tests whether
+ # Wo is "close" to .5 . If it is not, it adds an initial step of
+ # transforming the given low pass to a filter with a Wo of .5 and issues
+ # a warning.
+ Wt = varargin{4}; # vector of (norm. angular frequency)/pi transform targets
+ # [phi(1) phi(2) ... ]/pi
+ if(nargin < 4 || nargin > 5)
+ error("%s",usage)
+ endif
+ if(nargin == 5)
+ Pass = varargin{5};
+ switch(Pass)
+ case 'pass'
+ pass_stop = 1;
+ case 'stop'
+ pass_stop = -1;
+ otherwise
+ error("Pass must be 'pass' or 'stop'\n%s",usage)
+ endswitch
+ else
+ pass_stop = 1; # the default
+ endif
+ if(Wo <= 0)
+ error("Wo is %f <= 0\n%s",Wo,usage);
+ endif
+ if(Wo >= 1)
+ error("Wo is %f >= 1\n%s",Wo,usage);
+ endif
+ for i = 1 : length(Wt)
+ if(Wt(i) <= 0)
+ error("Wt(%d) is %f <= 0\n%s",i,Wt(i),usage);
+ endif
+ if(Wt(i) >= 1)
+ error("Wt(%d) is %f >= 1\n%s",i,Wt(i),usage);
+ endif
+ endfor
+ if(abs(Wo - .5) > .005) # Wo is not within 1% of .5
+ # Need to do an extra translation from a given low pass prototype filter
+ # with cutoff Wo to a filter with cutoff .5. First compute the
+ # denominator of an allpass for translating from a prototype with cutoff
+ # .5 to one with a cutoff of Wo.
+ P = p([Wo*pi]);
+ # From the low pass to low pass tranformation 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 from Wo to .5 can be
+ # obtained merely by a sign reversal.
+ P(2) = -P(2);
+ PP = revco(P);
+ [B,A] = transform(B,A,PP,P,1); # move low pass filter cutoff to .5
+ s1 = "A first step was added to transform the prototype to Wo = .5\n";
+ s2 = "Better to start with a prototype low pass filter with cutoff .5\n";
+ warning("%s: Wo = %f\n %s %s",mfilename(),Wo,s1,s2);
+ endif
+
+ phi = pi * Wt; # vector of normalized angular frequencies between 0 and pi
+ # Specified multiband filter edges 0 < phi(1) < phi(2) ... < phi(n) < pi
rads
+ #
+ # For Pass == 'pass', the target multiband will be:
+ # ------- --------- ----------
+ # | | | | | ...
+ # 0 phi(1) phi(2) phi(3) phi(4) phi(5) pi
+ #
+ # For Pass == 'stop', the target multiband will be:
+ # -------- ---------- -----------
+ # | | | | | | ...
+ # 0 phi(1) phi(2) phi(3) phi(4) phi(5) phi(6) pi
+
+ # This module takes a low pass IIR filter and transforms it 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
+
+ # PP(z)
+ # Find allpass H(z) = ---- ([FFM] eq. 2) to transform Z of low pass filter
+ # P(z)
+ # G(Z) to G(H(z)).
+ #
+ # ^
+ # Variable PP used here for 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
+ # note: (len - 1) == n in [FFM] eq. 3
+ P = p(phi); # calculate denominator of allpass for this vector of edges
+ PP = revco(P); # numerator of allpass has reversed coefficients of P
+ [Num,Den] = transform(B,A,PP,P,pass_stop);
+ AllpassNum = PP;
+ AllpassDen = P;
+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 = NA(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);
+ endfor
+ Ptemp = conv(Ptemp,P); # increase power of P by one
+ endfor
+
+ # 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);
+ endif
+ # i-1
+ # Regenerate PP^ (pp_powim1)
+ if(i == 1)
+ pp_powim1 = [1];
+ else
+ pp_powim1 = revco(ppower(i-1,powcols));
+ endif
+ if(i <= nb)
+ Bterm = (pass_stop^(i-1))*B(i)*conv(pp_powim1,p_pownmi);
+ Num = polysum(Num,Bterm);
+ endif
+ if(i <= na)
+ Aterm = (pass_stop^(i-1))*A(i)*conv(pp_powim1,p_pownmi);
+ Den = polysum(Den,Aterm);
+ endif
+ endfor
+ # Scale both numerator and denominator to have Den(1) = 1
+ temp = Den(1);
+ for i = 1 : length(Den)
+ Den(i) = Den(i) / temp;
+ endfor
+ for i = 1 : length(Num)
+ Num(i) = Num(i) / temp;
+ endfor
+endfunction
+
+function P = p(phi)
+ # Given phi, a vector of normalized angular frequency transformation
targets,
+ # return P, the denominator of an allpass H(z)
+ len = length(phi);
+ Pkm1 = 1; # P0 initial condition from [FFM] eq. 22
+ for k = 1 : len
+ P = pk(Pkm1, k, phi(k)); # iterate
+ Pkm1 = P;
+ endfor
+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) += sin_k * Pkm1(i) - ((-1)^k * cos_k * Pkm1(k+1-i));
+ #
+ # -1
+ # Multiplication by z just shifts by one coefficient
+ Pk(i+1) += sin_k * Pkm1(i) + ((-1)^k * cos_k * Pkm1(k+1-i));
+ endfor
+ # 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;
+ endfor
+endfunction
+
+function PP = revco(p) # reverse components of vector
+ l = length(p);
+ for i = 1 : l
+ PP(l + 1 - i) = p(i);
+ endfor
+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(isna(Ppower(i,j)))
+ break;
+ endif
+ p = horzcat(p, Ppower(i,j));
+ endfor
+ endif
+endfunction
+
+function poly = polysum(p1,p2) # add polynomials of possibly different length
+ n1 = length(p1);
+ n2 = length(p2);
+ if(n1 > n2)
+ # pad p2
+ p2 = horzcat(p2, zeros(1,n1-n2));
+ elseif(n2 > n1)
+ # pad p1
+ p1 = horzcat(p1, zeros(1,n2-n1));
+ endif
+ poly = p1 + p2;
+endfunction
------------------------------------------------------------------------------
Special Offer-- Download ArcSight Logger for FREE (a $49 USD value)!
Finally, a world-class log management solution at an even better price-free!
Download using promo code Free_Logger_4_Dev2Dev. Offer expires
February 28th, so secure your free ArcSight Logger TODAY!
http://p.sf.net/sfu/arcsight-sfd2d
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
