# Electronic filter design using RKFUN arithmetic

### Mario Berljafa and Stefan Güttel, September 2015Download PDF or m-file

## Contents

## Introduction

With Version 2.2 of the RKToolbox the `rkfun` class and its methods have been substantially extended. It is now possible to perform basic arithmetic operations on rational functions that go beyond scalar multiplication and addition of constants. For example, one can now add, multiply, divide, or exponentiate `rkfun` objects. Composition of an `rkfun` with a Moebius transform is also possible. Furthermore, the `rkfun` constructor now comes with the option to construct an object from a symbolic expression, or a string specifying a rational function from the newly introduced `gallery`. We now demonstrate some of these features. An application to electronic filter design is given below.

## The RKFUN gallery

The `rkfun` class provides a gallery function which allows for the quick construction of some useful rational functions. A list of the functions currently implemented can be obtained by typing:

```
help rkfun.gallery
```

GALLERY Collection of rational functions. obj = rkfun.gallery(funname, param1, param2, ...) takes funname, a case-insensitive string that is the name of a rational function family, and the family's input parameters. See the listing below for available function families. constant Constant function of value param1. cheby Chebyshev polynomial (first kind) of degree param1. cayley Cayley transformation (1-z)/(1+z). moebius Moebius transformation (az+b)/(cz+d) with param1 = [a,b,c,d]. sqrt Zolotarev sqrt approximation of degree param1 on the positive interval [1,param2]. invsqrt Zolotarev invsqrt approximation of degree param1 on the positive interval [1,param2]. sqrt0h balanced Remez approximation to sqrt(x+(h*x/2)^2) of degree param3 on [param1,param2], where param1 <= 0 <= param2 and h = param4. sqrt2h balanced Zolotarev approximation to sqrt(x+(hx/2)^2) of degree param5 on [param1,param2]U[param3,param4], param1 < param2 < 0 < param3 < param4, h = param6. invsqrt2h balanced Zolotarev approximation to 1/sqrt(x+(hx/2)^2) of degree param5 on [param1,param2]U[param3,param4], param1 < param2 < 0 < param3 < param4, h = param6. sign Zolotarev sign approximation of degree 2*param1 on the union of [1,param2] and [-param2,-1]. step Unit step function approximation for [-1,1] of degree 2*param1 with steepness param2.

For example, we can construct an `rkfun` corresponding to a Moebius transform as follows:

```
r = rkfun.gallery('moebius', [4, 3, 2, -1])
```

r = RKFUN object of type (1, 1). Real-valued Hessenberg pencil (H, K) of size 2-by-1. coeffs = [0.000, 1.000]

As always with `rkfun` objects, we can perform several computations on `r`, such as computing its roots and poles:

format long e disp([roots(r) poles(r)])

-7.500000000000000e-01 5.000000000000000e-01

## The symbolic constructor

Symbolic strings can also be used to construct `rkfun` objects, provided that the required MATLAB symbolic toolboxes are installed. Here we construct an `rkfun` corresponding to the rational function :

```
s = rkfun('3*(z-3)*(z^2-1)/(z^5+1)')
```

s = RKFUN object of type (2, 4). Real-valued Hessenberg pencil (H, K) of size 5-by-4. coeffs = [0.000, 0.000, -0.000, -0.000, 23.570]

This works by first finding the roots and poles of the input function symbolically and then constructing the `rkfun` with these roots and poles numerically. Note that `s` is only of type as one of the five denominator roots has been cancelled out by symbolic simplifications prior to the numerical construction. The constructor should issue an error if the provided string fails to be parsed by the symbolic engine, or if it does not represent a rational function.

## Arithmetic operations with rational functions

Since Version 2.2 of the RKToolbox one can add, multiply, divide, and exponentiate `rkfun` objects. All these operations are implemented via transformations on generalized rational Krylov decompositions [1, 2]. For example, here is the product of the two rational functions and from above:

disp(r.*s)

RKFUN object of type (3, 5). Real-valued Hessenberg pencil (H, K) of size 6-by-5. coeffs = [0.000, 0.000, 0.000, -0.000, -0.000, ...]

Note that it is necessary to use point-wise multiplication `.*`, not the matrix-multiplication operator `*`. This is to be consistent with the MATLAB notation, and also with the notation used in the Chebfun system [5] (where the `*` operator returns the inner product of two polynomials; something we have not yet implemented for `rkfun` objects).

It is also possible to compose `rkfun` objects as long as the inner function is a Moebius transform, i.e., a rational function of type at most . The rational function from above *is* a Moebius transform, hence we can form the function :

f = 1./s(r)

f = RKFUN object of type (4, 4). Real-valued Hessenberg pencil (H, K) of size 5-by-4. coeffs = [-0.000, -0.002, -0.000, -0.005, -0.042]

The roots of should correspond to the poles of , which are the poles of mapped under the inverse function . The latter inverse function is indeed well defined in the whole complex plane as is an invertible Moebius transform. We can compute it via the `inv` command. Let's verify that the roots/poles agree numerically:

rts = sort(roots(f)); pls = sort(feval(inv(r), poles(s))); disp([rts, pls])

Column 1 -4.256702311079903e-01 - 3.812725100682988e-01i -4.256702311079903e-01 + 3.812725100682988e-01i -1.187966132528373e+00 - 8.330610885153919e-01i -1.187966132528373e+00 + 8.330610885153919e-01i Column 2 -4.256702311079904e-01 - 3.812725100682987e-01i -4.256702311079904e-01 + 3.812725100682987e-01i -1.187966132528373e+00 - 8.330610885153925e-01i -1.187966132528373e+00 + 8.330610885153925e-01i

## Constructing rational filters

Rational filter functions are ubiquitous in scientific computing and engineering. For example, in signal processing [3, 6] one is interested in deriving rational functions that act as filters on selected frequency bands. Using the `gallery` of the RKToolbox and some `rkfun` transformations, we can construct meaningful rational filters in just a few lines of MATLAB code. Below is a plot of four popular filter types, which are obtained by rational transforms of Chebyshev polynomials or the `step` function from `rkfun`'s `gallery`.

x = rkfun; butterw = 1./(1 + x.^16); cheby = rkfun('cheby',8); cheby1 = 1./(1 + 0.1*cheby.^2); cheby2 = 1./(1 + 1./(0.1*cheby(1./x).^2)); ellip = rkfun('step'); figure subplot(221),ezplot(butterw,[0,2]); title('Butterworth') subplot(222),ezplot(cheby1, [0,2],'r');title('Chebyshev type 1') subplot(223),ezplot(cheby2, [0,2],'g');title('Chebyshev type 2') subplot(224),ezplot(ellip, [0,2],'m');title('Elliptic')

The reader may compare this plot with that at the bottom of the Wikipedia page on *electronic filters* [7]. For example, the filter `cheby2` involves multiple inverses of a Chebyshev polynomial in the transformed variable . It has the so-called *equiripple property* in the stopband, which is the region where the filter value is close to zero. The elliptic filter `ellip`, also known as Cauer filter [4], has equiripple properties in both the stop- and passbands. Such filters are based on Zolotarev's equioscillating rational functions [8], which are also implemented in the `gallery` of the RKToolbox.

## Limitations

Although we hope that the new `rkfun` capabilities demonstrated above are already useful for many practical purposes, there are still some short-comings one has to be aware of. The main problem is that combinations of `rkfun` objects may have degrees higher than theoretically necessary, which may lead to an unnecessarily fast growth of parameters. For example, when subtracting an `rkfun` from itself,

```
r = rkfun('(1-x)/(1+x)');
d = r - r
```

d = RKFUN object of type (2, 2). Real-valued Hessenberg pencil (H, K) of size 3-by-2. coeffs = [0.000, -1.000, 1.000]

we currently obtain a type rational function, instead of the expected type . This is because the sum (or difference) of two type rational functions is of type *in the worst case*, and we currently do not perform any degree reduction on the sum (or difference). (A similar problem is encountered with multiplication or division.) The roots and poles of the function `d` are meaningless, but the evaluation works fine, except when we evaluate at (or nearby) a "legacy" pole:

d([-1, -1+1e-16, 0, 1, 2])

ans = Columns 1 through 3 NaN 0 0 Columns 4 through 5 0 -1.850371707708596e-17

A numerical degree reduction would probably require the concept of a "domain of evaluation." To illustrate, consider the rational function . The exact type of this function is , but the residue associated with is tiny, so one may conclude that this pole could be removed. However, when evaluating for , the removal of this pole would lead to an inaccurate result (try ). Hence, reliable degree reduction may only be possible when a relevant "domain of evaluation" is specified.

## References

[1] M. Berljafa and S. Güttel. *Generalized rational Krylov decompositions with an application to rational approximation,* SIAM J. Matrix Anal. Appl., 36(2):894--916, 2015.

[2] M. Berljafa and S. Güttel. *The RKFIT algorithm for nonlinear rational approximation,* SIAM J. Sci. Comput., 39(5):A2049--A2071, 2017.

[3] H. Blinchikoff and A. Zverev. *Filtering in the Time and Frequency Domains,* John Wiley & Sons Inc., New York, 1976.

[4] W. Cauer. *Ein Interpolationsproblem mit Funktionen mit positivem Realteil,* Math. Z., 38(1):1-44, 1934.

[5] T. A. Driscoll, N. Hale, and L. N. Trefethen. *Chebfun Guide,* Pafnuty Publications, Oxford, 2014. http://www.chebfun.org

[6] M. Van Valkenburg. *Analog Filter Design,* Holt, Rinehart and Winston, 1982.

[7] Wikipedia, entry *Electronic filter* as of 17/09/2015. https://goo.gl/vkQwWG

[8] E. I. Zolotarev. *Application of elliptic functions to questions of functions deviating least and most from zero,* Zap. Imp. Akad. Nauk St. Petersburg, 30:1-59, 1877. In Russian.