Description of the included Matlab scripts

Compkautz.m

S=compkautz(a,N,x) produces tap-output signals of a complex Kautz filter, defined by a pole vector a, |a|<1, and excited by an input signal vector x. N is the desired length of the tap-output signals if N > length(x). Tap-output signals are the columns of the signal matrix S. Notice that S is real valued only in the case of real poles.

For example:

• G=compkautz(a,N,1) are the (length N) tap-output impulse responses of the Kautz filter defined by a. With a=.7*ones(M,1) we would have the Laguerre filter tap-output responses (Laguerre parameter 0.7 and filter order M).
• The orthonormatity of the impulse responses could be tested by forming the matrix product C=G'*G: for sufficient N, C is the (appropriate) unit matrix.
• More generally, C=S'*S is the correlation matrix which implements the (time domain) inner products of the tap-output signals.
• y=S*w, with S=compkautz(a,N,x), is the Kautz filter response to input x with Kautz filter coefficient vector w. The filter order is M=length(a).
• Identification: the least-square estimate corresponding to a particular Kautz filter structure (i.e. for a fixed a) and an given input-output-data pair (x,y) is formally w=inv(S'*S)*S'*y.
• Approximation: an infinite dimensional Kautz filter can be seen as an analysis-synthesis-structure implementing a particular Fourier transformation: any causal and finite-energy signal h has a representation h=G*c with Fourier coefficients c=G'*h. The partial sums of these representations are the corresponding least-square approximations.
• The synthesis filter coefficients c can be evaluated implementing convolutions by filtering: c=S(N,:), c=c(:), N=length(h), where S is the signal matrix produced by feeding the time-reversed signal, ih=h(N:-1:1) to S=compkautz(a,N,ih) and reading the tap-outputs at N. These are the lest-square optimal coefficients for approximating h with y=G*c, G=compkautz(a,N,1).

Mixkautz1.m

S=mixkautz1(a,N,x) produces tap-output signals of a real Kautz filter, defined by a pole vector a, excited by an input signal vector x, using a particular mixed structure for real and complex conjugate pole pairs. N is the desired length of the output, if N > length(X). Poles with imag(a)<0 are disregarded to avoid duplication produced by complex conjugate poles in a.

Mixkautz1.m works like compkautz.m (see the above examples), except for the filter order induced by the pole vector a. However, for a distribution of real and complex conjugate pair poles it is length(a).

For example:

• A complex pole a=r*exp(j*ang), 0<ang<pi, induces a Kautz structure corresponding to a complex conjugate pole pair (a,a’) (cc-Kautz). The Nth order tap-responses are produces by using pole vector a=a*ones(N/2,1).
• The cc-Kautz structure can be seen as a generalization of the Laguerre structure: the resolution of the modeling is focused to a frequency region defined by the pole angle ang with intensity r, 0<r<1.
• For a purely complex a, mixkautz1.m produces impulse responses whose magnitude response maximums correspond to the pole angles. Each resonance is represented by two responses with different shapes, due to the additional zeros at –1 and 1 of the second-order blocks, associated to the odd and even indexed signals, respectively.

Kautzapprx1.m

[c,E,h]=kautzapprx1(a,x) composes an approximation of a given response x using mixkautz1.m. The vector c contains the least-square optimal filter coefficients for approximating x with h=G*c, where G=mixkautz1(a,N,1) are the corresponding tap-output impulse responses, and N=length(x). E is the normalized root-mean-square error of the error signal x-h.

Coefficients c are the true orthonormal series-expansion coefficients (, the Fourier coefficient,) corresponding to the chosen Kautz structure. This fact is utilized in the evaluation of the coefficient: c’=S(N,:), S=mixkautz1(a,N,x(N:-1:1)), i.e., the convolution sums defining the coefficient are evaluated by filtering the time-inverted x and reading the tap-outputs at time N. Notice that the convolutions are in general infinite summations, but they can be implemented by finite filtering operations. This can be seen as a generalization of rectangular-window FIR filter design.

For example:

• We may use pole distributions reflecting the desired modeling resolution allocation, e.g. a=r.*exp(j*f) with logarithmically spaced f in [0,pi] for some choice of r.
• It is relatively easy to find usable pole sets by direct selection of prominent resonances and proper pole radius tuning (by inspecting the time and frequency responses).
• All-pole or pole-zero modeling can be applied in search for suitable poles. The Steiglitz-McBride method has proven to be the most useful of the standard (Matlab) methods.
• Easy and fast evaluation of the filter coefficients (and consequently the approximation error) enables many kinds of random and iterative search schemes.

Kautzapprx2.m

[c,e ,h]=kautzapprx2(a,M,x) approximates x with the impulse response h of a real Kautz filter with M identical blocks, each block defined by the pole vector a. It produces the same h as kautzapprx1.m with a=a*ones(1,M), a=a(:), but it preserves the structure of identical blocks , i.e., the order in the coefficient vector c. This is a way to essentially separate the pole optimization procedure and the model order selection. According to the orthogonality, we get the approximation error associated to the number of block used in the approximation, i=1,…,M, in the form e(i)=x’*x – c(1:p*i)’*c(1:p*i), where p is the order of one block. Observing e(i), i=1,…,M, gives a trivial way to perform model reduction with respect to a fixed substructure defined by a. 

For example:

• Laguerre and cc-Kautz filters are examples of structures with identical blocks. In these cases combined pole and filter order optimization methods based on random or iterative search are especially tempting.
• For higher block orders we probably first have to fix a substructure and then evaluate its contribution in the cascaded system.
• A typical approach is to choose a few prominent resonances from the magnitude response and represent them with appropriately radius tuned pole pairs. Effectively the repetitive use of the same block works two ways: it sharpens the "resonators" and the increase in filter order improves the off-resonance modeling.
• Any causal and stable all-pass filter induces a Kautz filter with identical blocks.

Braun1.m

An iterative method to optimize the pole positions defining a real Kautz filter. We have adopted this algorithm to the context of Kautz filters from a pure FIR-to-IIR filter conversion method proposed by Brandenstein and Unbehauen [1]. It resembles a lot the Steiglitz-McBride method of pole-zero modeling, but it genuinely optimizes the pole positions of the orthonormal representation, and not just the poles of a coupled rational overall structure. Furthermore, it produces always stable (and theoretically globally optimal) poles.

[e,Q]=braun1(N,I,x) uses I iterations to search for an Nth order Kautz structure to approximate x. Q(i,:) are the (Kautz filter common demonimator) polynomial coefficient and e(i) is the approximation (RMS) error as a function of the iterations, i=1,…,I. For sufficient I, we use roots(Q(k,:)), corresponding to minimum e from [m,k]=min(e), as the Kautz filter poles. Notice that this method doesn’t require intermediate evaluation of polynomial roots.

 For high orders N, and/or long signals x, braun1.m will cause warnings and take time. The error e may peak to NaN and Inf, but we should still get reasonable results. In our experience the algorithm works surprisingly well, even for really tough cases.

 For example:

• By choosing N=1,2, we optimize the Laguerre or cc-Kautz structures, respectively, but actually the optimization is done with respect to the one block structures, which is probably not what we want. It still gives a picture of the "center of gravity" of the poles.
• If we have observed e.g. five prominent resonances in the magnitude response we could choose N=10 in anticipation of five pole pairs corresponding to the resonances.
• However, a higher N is often required to capture the desired resonances. We have found that the poles corresponding to a chosen set of prominent resonances are always attainable by pruning a large set of poles produced by braun1.m.

Kautzdemo1.m

Kautzdemo1 is a demonstration of all the above Matlab functions. There is quite a lot of reading and at some points you may have to wait a while. One could argue that it would really need some chopping, but we have grabbed to this unique opportunity to awaken your interest in Kautz filters.

After an introduction, the user is requested to choose a pole distribution to work with. The consequent time and frequency responses produced by compkautz.m and mixkautz1.m are demonstrated. Kautzapprx1.m and kautzapprx2.m are presented using a measured acoustic guitar body response as the target response. The utilization of barun1.m is demonstrated in many ways.

Kautzdemo2.m

Kautzdemo2 is a case study type of Kautz filter design example concerning loudspeaker equalization. A pre-designed target response for the equalizer is used so the example deals with pure magnitude response equalization. The functions demonstrated are braun1.m and kautzapprx1.m, the latter implying the use of mixkautz1.m.

After presenting the measured, target equalized and equalizer responses, some straightforward FIR equalizers are presented for comparison. Results for Laguerre and two-pole Kautz equalizers are demonstrated for comparable filter orders. These are the most time-consuming calculations. Distributed pole set for the Kautz equalizers are composed using braun1.m. To really flatten a sharp and local deviation in the middle frequencies, 3 to 4 pole pairs are tuned manually and appended to the pole sets proposed by braun1.m. Finally the whole Kautz equalizer structure defined by 10 complex conjugate pole pairs is tuned manually. Equalization results accomplished during the process are compared those attained by Laguerre and FIR filtering.