This file is indexed.

/usr/share/octave/packages/signal-1.2.2/iirlp2mb.m is in octave-signal 1.2.2-1build1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
## Copyright (C) 2011 Alan J. Greenberger <alanjg@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)
## [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
##
## Example of use:
## [B, A] = butter(6, 0.5);
## [Num, Den] = iirlp2mb(B, A, 0.5, [.2 .4 .6 .8]);

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
   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; # Pass == 'pass' is 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
   oWt = 0;
   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
      if(Wt(i) <= oWt)
         error("Wt(%d) = %f, not monotonically increasing\n%s",i,Wt(i),usage);
      else
         oWt = Wt(i);
      endif
   endfor

   #                                                             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 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 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 = revco(P); # 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(1); # normalize
   AllpassNum = pass_stop * revco(AllpassDen);
   [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 = 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 = 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;
   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