Electronic filter design using RKFUN arithmetic
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 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:
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)])
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.
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:
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  (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.256702311079903e-01 - 3.812725100682988e-01i -4.256702311079903e-01 + 3.812725100682988e-01i -1.187966132528373e+00 - 8.330610885153928e-01i -1.187966132528373e+00 + 8.330610885153928e-01i
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 . 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 , has equiripple properties in both the stop- and passbands. Such filters are based on Zolotarev's equioscillating rational functions , which are also implemented in the gallery of the RKToolbox.
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 0 NaN 2.465190328815662e-31 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.
 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.
 M. Berljafa and S. Güttel. The RKFIT algorithm for nonlinear rational approximation, SIAM J. Sci. Comput., 39(5):A2049--A2071, 2017.
 H. Blinchikoff and A. Zverev. Filtering in the Time and Frequency Domains, John Wiley & Sons Inc., New York, 1976.
 W. Cauer. Ein Interpolationsproblem mit Funktionen mit positivem Realteil, Math. Z., 38(1):1-44, 1934.
 T. A. Driscoll, N. Hale, and L. N. Trefethen. Chebfun Guide, Pafnuty Publications, Oxford, 2014. http://www.chebfun.org
 M. Van Valkenburg. Analog Filter Design, Holt, Rinehart and Winston, 1982.
 Wikipedia, entry Electronic filter as of 17/09/2015. https://goo.gl/vkQwWG
 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.