Author: Nalin Pithwa

Biomedical solutions via math

Nalin PithwaLeave a comment

Reproduced from newspaper, DNA, print edition, Mumbai, Sunday, Mar 5, 2017, (Section on Health): 

Biomedical solutions via math:

(A new model combines mathematics with biology to set the stage for cancer cure and other diseases):

Ann Arbor:

How do our genes give rise to proteins, proteins to cells, and cells to tissues and organs? What makes a cluster of cells become a liver or a muscle? The incredible complexity of these biological systems drives the work of biomedical scientists. But, two mathematicians have introduced a new way of thinking that may help set the stage for better understanding of our bodies and other living things.

The pair from University of Michigan Medical School and University of California, Berkeley talk of using math to understand how generic information and interactions between cells give rise to the actual function of a particular type of tissue. While the duo admit that it’s a highly idealized framework, which does not take into account every detail of this process, that’s what’s needed. By stepping back and making a simplified model based on mathematics, they hope to create a basis for scientists to understand the changes that happen over time within and between cells to make living tissues possible. It could also help with understanding of how diseases such as cancer can arise when things don’t go as planned.

Turning to Turing’s machine:

U-M Medical School Assistant Professor of Computational Medicine, Indika Rajapakse and Berkeley Professor Emeritus, Stephen Smale have worked on the concepts for several years. “All the time, this process is happening in our bodies, as cells are dying and arising, and yet, they keep the function of the tissue going,” says Rajapakse. “We need to use beautiful mathematics and beautiful biology together to understand the beauty of a tissue.”

For the new work, they even hearken back to the work of Alan Turing, the pioneering Btitish mathematician famous for his “Turing machine” computer that cracked the codes during World War II.

Toward the end of his life, Turing began looking at the mathematical underpinnings of morphogenesis — the process that allows natural patterns such as a zebra’s stripes to develop  as a living thing grows from an embryo to an adult.

“Our approach adapts Turing’s technique, combining genome dynamics within the cell and the diffusion dynamics between cells,” says Rajapakse, who leads the U-M 4D —- Genome Lab in the Department of Computational Medicine and Bio-Informatics.

His team of biologists and engineers conduct experiments that capture human genome dynamics to three dimensions using bio-chemical methods and high resolution imaging.

Bringing math and the genome together

Smale, who retired from Berkeley, but is still active in research, is considered a pioneer of modelling dynamic systems. Several years ago, Rajapakse approached him during a visit to U-M, where Smale earned his undergraduate and graduate degrees. They began exploring how to study the human genome — the set of genes in an organism’s DNA — as a dynamic system.

They based their work on the idea that while the genes of an organism remain the same throughout life, how cells use them does not.

Last spring, they published a paper that lays a mathematical foundation for gene regulation — the process that governs how often and when genes get “read” by cells in order to make proteins.

Instead of the nodes of those networks being static, as Turing assumed, the new work sees them as dynamic systems. The genes may be “hard-wired” into the cell, but how they are expressed depends on factors such as epigenetic tags added as a result of environmental factors, and more.

Next Step:

As a result of his work with Smale, Rajapakse now has funding from the Defense Advanced Research Projects Agency (DARPA), to keep exploring the issue of emergence of function — including what happens when the process changes.

Cancer, for instance, arises from a cell development and proliferation cycle gone awry. And the process by which induced pluripotent stem cells are made in a lab —- essentially turning back the clock on a cell type so that it regains the ability to become other cell types — is another example.

Rajapakse aims to use data from real world genome and cell biology experiments in his lab to inform future work, focused on cancer and cell reprogramming.

He’s also organizing a gathering of mathematicians from around the world to look at computational biology and the genome this summer in Barcelona.

**************************************************************************************

Thanks to DNA, Prof. Stephen Smale and Prof. Indika Rajapakse; this, according to me, is one of the several applications of math.

–Nalin Pithwa.

Parameters for selection of a processor — part 2

Nalin PithwaLeave a comment

Processors suitable for digital control range from standard microprocessors like 8051 to special purpose DSP processors, the primary difference being in the instruction sets and speed(s) of particular instruction(s), such as multiply. Standard microprocessors or general purpose processors are intended for laptops, workstations, and general digital data bookkeeping. Naturally, because digital control involves much numerical computation, the instruction sets for special-purpose DSP processors are rich in math capabilities and are better suited for control applications than the standard microprocessors or general purpose processors.

Processors, such as those found in microwave ovens have a broad range of requirements. For example, the speed, instruction set, memory, word length, and addressing mode requirements are all very minimal for the microwave oven. The consequences of a data error are minimal as well, especially, relative to a data error in a PC/laptop while it is calculating an income tax return. PC’s or laptops, on the other hand, require huge megabytes of memory, and they benefit from speed, error correction, larger word size, and sophisticated addressing modes.

DSP processors generally need speed, word length, and math instructions such as multiply, multiply-and-accumulate, and circular addressing. One typical feature of signal processors not found in general purpose processors is the use of a Harvard architecture, which consists of separate data and program memory. Although separate data and program memory offer significant speed advantages, the IC pin count is higher assuming external memory is allowed because instruction address, instruction data, data address, and data buses are separate. A modified Harvard architecture has been used which maintains some speed advantage, while eliminating the requirement for separate program and data buses, greatly reducing pin count in processors that have external memory capability (almost all have this feature).

While thinking of control versus signal processing applications, in the former, we often employ saturation and therefore absolutely require saturation arithmetic; whereas in the latter, to ensure signal fidelity, in most signal processing applications the algorithms must be designed to prevent saturation by scaling signals appropriately.

The consequences of numerical overflow in control computations can be serious, even destabilizing. In most forms of numerical computation, it is usually better to suffer the non-linearity of signal saturation than the effects of numerical overflow.

For most control applications, it is advantageous to select a processor that does not require much support hardware. One of the most commonly cited advantages of digital control is the freedom from noise in the control processor. Although it is true that controller noise is nominally limited to equalization noise, it is not true that the digital controller enjoys an infinite immunity from noise. Digital logic is designed with certain noise margins, which of course are finite. When electromagnetic radiation impinges on the digital control system, there ia finite probability of making an error. One of the consequences of digital control is that although it can have a very high threshold of immunity, without error detection and correction it is equally likely that the system will make a large error as a small one — the MSB and the LSB of a bus have equal margin against noise.

In addition to external sources of error-causing signals, the possibility for circuit failure exists. If a digital logic circuit threshold drifts outside the design range, the consequences are usually catastrophic.

For operational integrity, error detection is a very important feature.

I hope to compare, if possible, some families of Digital Control processors here, a bit later.

Regards,

Nalin Pithwa

Reference:

Digital Control of Dynamic Systems, Franklin, Powell and Workman.

Parameters for selection of a processor — part 1 : Some Notes

Nalin PithwaLeave a comment

Ref: (1) Digital Control of Dynamic Systems by Franklin, Powell and Workman(2) the internet wrb (3) Manuals of various DSP vendors like TI and Analog Devices.

Before examining some of the requirements placed on Digital Signal or Control Processors (DSP’s or DCP’s), let’s define some of the basic parametric considerations in choosing a processor:

  1. Speed: The rate at which a computer can execute instructions. Note that speed is important relative only to word length and instruction complexity. Although a bit of an oversimplification, the amount of work done per unit time can be loosely expressed as “power \approx Speed x Word length x Instruction power; where word length and instruction power are defined below. Basically, the aggregate measure of processor power must not only include the speed at which it can execute instructions but it must also to the complexity of the instructions and the sice of the word on which they operate. Quoting processor speed (often measured by MIPS) without word length and instruction complexity is like specifying a trenching machine by the length of the trench it can dig per unit time, without specifying how wide or how deep the resulting trench would be. Comparing versions of the same basic processor by comparing MIPS is valid, but most literature discusses microprocessors of the same type by referring to their clock rate.
  2. Arithmetic Logic Unit: (ALU):Usually referred to with the acronym ALU, the arithmetic logic unit is the hardware that performs the basic arithmetic operations in the processor, including add, subtract, and boolean logic operations.
  3. Multiplier Accumulator (MAC):Hardware dedicated to multiplying two numbers and adding into an accumulator register. In control and signal processors, the speed of the MAC is of paramount importance.
  4. Instruction Set: The basic operations that can be performed by the processor. These range from bit-flipping operations accomplished in the ALU to fixed- and floating-point arithmetic accomplished in the multiplier-accumulator. More complex instructions, such as multiply, can be measured in terms of the number of logical operations they require (such as shift and adds in the case of multiply without an MAC). Because many processors execute different instructions at different rates, processor power is best measured in terms of benchmark tasks. The key difference among benchmarks is the instruction mix they represent. Benchmarks can be misleading and must be examined carefully before they are used to make comparisons between machines.It is possible to construct a benchmark that favors one machine over another even though for the application being considered it may be inferior. Some common benchmarks include the Whetstone and the Dhrystone. (Refer to the web or the books Computer Organization and Computer Architecture by Hennessy and Patterson for more details of these benchmarks). The Whetstone benchmark is a measure of Fortran and floating point performance and records the time it takes to perform predefined floating point operations. The Dhrystone benchmark is a measure of integer processes per second. Integer processes include operations such as assignment, arithmetic and control. This benchmark uses integer, character, string, pointer and record operand types. Complex fast FFT’s are also a common benchmark. Again, care must be taken in comparing benchmarks because small details can be omitted which could change the outcome of the comparison.
  5. Word length:The width of the data stream measured in bits, usually an integral number of bytes. Although the power of a processor as defined above increases with word length,it should be clear that in most applications there is a point beyond which increasing word length has no benefit (just like the law of diminishing returns in economics/finance). For fixed point arithmetic, data rarely need more than 24 to 32 bits of precision (though, in my opinion, this could change with IOT processors).
  6. Memory Size:The amount of memory space or address space the processor supports.
    General purpose applications typically require far more memory than signal processing or control applications, and signal processing far more than control applications.
  7. Addressing:The sophistication or number of addressing modes available. They include direct, indirect, immediate, relative, and more. High level language compilers can be made more efficient with more sophisticated addressing modes, but DSP’s and DCP’s make little use of anything but the very basic addressing modes. Circular addressing is an auto-increment address register that has a selectable modulo base producing an address that circulates through a memory block of relatively small length, such as 32 words. FIR filters are efficiently implemented using circular addressing (not too common in control applications), because a finite block of data that contains the recent N samples of input can be addressed very easily.
  8. Digital I/O:Addressable logic ports for output or input of logical signals to peripheral devices, such as a switch or an LED indicator.
    Whereas general purpose processors usually employ an I/O space using a bus architecture, DSP’s and DSP’s usually employ employ a small number of digital I/O lines in addition to an I/O space (addresses reserved for input-output registers using the data bits).
  9. Error Detection and Correction: Some applications require no error detection or correction, others require both. Mobile devices invariably employ error correction. In control systems, parity is a common and simple error detection technique. If the number of logically high (or 1) bits in a byte is odd, the parity is said to be odd; else it is even. When a bit, byte or word of information is transmitted to a receiver along with the parity bit(s), it allows the receiver to check the parity of the received information against the parity generated by the sender. As long as only one bit was in error, parity generation and checking allows detection of the error. Parity can be incorporated in a bit-for-bit basis (100% overhead), a bit-for-byte basis, or a bit-for-word basis. A simple parity bit-per-byte of information transmitted is most common, and in such a case it is easy to picture that for each group of nine bits (eight data bits plus one parity), the action of the parity bit is to keep the parity of the nine bits constant (odd or even). It seems that there are some industrial processors (especially IOT processors) which contain internal error detection and correction, such as in the internal buses, ALU’s (arithmetic logic units), or multipliers. For real-time or critical applications, the error detection and correction feature is integrated in hardware.
  10. Interrupts:Processor interrupts provide a means of rapidly switching the program that is being executed from one task to another in response to an external (hardware) stimulus. Hardware interrupts are usually inputs that cause the processor to be vectored to a new location in memory. Some interrupts can be masked (ignored) by setting control registers inside the processor, and some purposely cannot be disabled. In control applications, the calculation of the next control often begins when an A/D converter signals that it has completed a conversion. Using the conversion complete signal as an interrupt input is a common practice. It should be noted, however, that some processors have variable response time to interrupts, often depending on the instruction being executed at the time of the interrupt request. Especially true in general purpose processors, the variance in response time can be large, and it is damaging to the controller in that it represents control delay jitter. One advantage of most DSP processors is that the interrupt response time is constant.
  11. ROM, RAM, PROM, EPROM, EEPROM:These are, of course, types of memory. Although non-volatile ROM is useful for storing fixed program code, RAM can store both variable and program data. RAM is volatile and must be loaded from somewhere upon power-up unless it contains only calculated and measured data. Although a PROM could be read and stored in RAM, PROM’s are typically very fast,, and there is little reason to stage them in staged memory structures. EPROMs are erasable PROMs; the erase cycle takes special hardware and is slow compared to the read access time. EPROMs are extremely useful to the development of digital control systems where the controller code is updated only infrequently. EEPROMs appeared after all the other memories listed above. EEPROMs are electrically erasable PROMs that have an in-circuit erase mode, which is v slow but requires no special hardware to accomplish. In a system requiring a small amount of non volatile storage and that needs to be updated infrequently from some external source, EEPROMs are inexpensive solid state candidates. Care must be taken in the design to prevent accidental overwrite of the program and/or data stored in EEPROMs, even under temporary fault conditions. Transient noise (glitches), addressing errors, and the like should be considered in designing hardware to protect the device from being put in the write mode inadvertently. In a nutshell, the three parameters are price, performance and power. For IoT processors, security too needs consideration.

We will examine the details soon.

Nalin Pithwa

Bilinear Transformation via Synthetic Division: The Matrix Method

Nalin PithwaLeave a comment

(Authors: Prof. Navdeep M. Singh, VJTI, University of Mumbai and Nalin Pithwa, 1992).

Abstract: The bilinear transformation can be achieved by using the method of synthetic division. A good deal of simplification is obtained when the process is implemented as a sequence of matrix operations. Besides, the matrices are found to have simple structures suitable for algorithmic implementation.

I) INTRODUCTION:

Davies [1] proposed a method for bilinear transformation using synthetic division. This method can be quite simplified when the synthetic division is carried out as a set of matrix operations because the operator matrices are found to have simple structures and so can be easily generated.

II) THE ALGORITHM:

Given a discrete time transfer function H(z), it is transformed to G(s) in the s-plane under the transformation:

z=\frac{s+1}{s-1}

This can be sequentially achieved as :-

H(z) \rightarrow H(s+1) \rightarrow H(\frac{1}{s}+1) \rightarrow H(\frac{2}{s}+1) \rightarrow H(\frac{2}{s-1}+1)=H(\frac{s+1}{s-1})=G(s)

The first step is to represent the given characteristic polynomial in the standard companion form. Since the companion form represents a monic polynomial appropriate scaling is required in the course of the algorithm to ensure that the polynomial generated is monic after each step of transformation.

The method is developed for a third degree polynomial and then generalized.

Step 1:

H(z) \rightarrow H_{1}(s+1) \rightarrow H(s+1) (decreasing of roots by one)

Given H(z)=z^{3}+a_{2}z^{2}+a_{1}z+a_{0} (a_{3}=1) (monic)

Then, H_{1}(s+1)=s^{3}+b_{2}s^{2}+b_{1}s+b_{0} where b_{2}=a_{2}+1, b_{1}=a_{1}+b_{2}, b_{0}=a_{0}+b_{1}

H(s+1)=s^{3}+(a_{2}+3)s^{2}+(a_{1}+2a_{2}+3)s+(a_{0}+a_{1}+a_{2}+1).

In the companion form, obviously the following transformation is sought:

A = \left[ \begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -a_{0} & -a_{1} & -a_{2} \end{array}\right] and B=\left[ \begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 &1 \\-b_{0} & -b_{1} & -b_{2} \end{array} \right] and C = \left[ \begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -a_{0}-a_{1}-a_{2}-1 & -a_{1} - 2a_{2}-3 & -a_{2}-3 \end{array} \right]

Performing elementary row and column transformations on A using matrix operators, the final row operator and column operator matrices, P_{1} and Q_{1}, respectively are found to be P_{1} = \left[ \begin{array}{ccc} (a_{0}+a_{1})/a_{0} & a_{2}/a_{0} & 1/a_{0} \\ -1 & 1 & 0 \\ 0 & -1 & 1\end{array} \right] and Q_{1} = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \end{array}\right].

Thus, P_{1}AQ_{1}=B.

In general, for a polynomial of degree n,

P_{1} = \left[ \begin{array}{ccccc} (a_{0}+a_{1})/a_{0} & a_{2}/a_{0} & a_{3}/a_{0} & \ldots & 1\\ -1 & 1 & 0 & \ldots & a_{0}\\ 0 & -1 & 1 & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \vdots & \vdots & \vdots & -1 & 1\end{array}\right] and Q_{1} = \left[ \begin{array}{ccccc} 1 & 0 & 0 \ldots & 0 \\ 1 & 1 & 0 & \ldots & \ldots \\ 1 & 1 & 1 & \ldots & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 1 & 1 & \ldots & 1\end{array}\right], where both the matrices are n \times n.

and B = P_{1}AQ_{1} = \left[ \begin{array}{ccccc} 0 & 1 & 0 & \ldots & 0\\ 0 & 0 & 1 & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ -b_{0} & -b_{1} & -b_{2} & \ldots & -b_{n-1} \end{array}\right], and B is also n \times n.

Where b_{j} = \sum_{i=j}^{n-1}a_{i} + 1, where i, j = 0, 1, 2, \ldots, (n-1).

Now, the transformation B \rightarrow C is sought respectively. The row and column operator matrices P_{c1} and Q_{c1} respectively are : P_{c1} = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ -1 & 1 & 0 \\ 1 & -2 & 1 \end{array}\right] and Q_{c1} = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 1 \end{array}\right].

P_{c1} and Q_{c1} are found to have the following general structures:

P_{c1} = \left[ \begin{array}{ccccc}1 & 0 & \ldots & \ldots & 0 \\ -1 & 1 & \ldots & \ldots & 0 \\ 1 & -2 & 1 & \ldots & 0 \\ -1 & 3 & -3 & 1 & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \ldots & \ldots & \ldots & \ldots & 1\end{array}\right] and Q_{c1} = \left[ \begin{array}{ccccc} 1 & 0 & \ldots & \ldots & 0 \\ 0 & 1 & 0 & \ldots & 0 \\ 0 & 1 & 1 & \ldots & 0 \\ 0 & 1 & 2 & 1 &\ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \ldots & \ldots & \ldots & \ldots & 1\end{array}\right], both general matrices P_{c1} and Q_{c1}, being of dimensions n \times n.

P_{c1} is lower triangular and can be generated as : P_{c1}(i,j) = [|P_{c1(i-j)}|+|P_{c1}(i-1,j-1)|] \times (-1)^{(i+j)}, where i=2, 3, 4, \ldots, n and so we get

P_{c1}(i,1)=(-1)^{(i-1)}, where i=1, 2, 3, \ldots, n, and P_{c1}(i,i)=1.

Similarly, Q_{c1} is lower triangular and can be generated as :

Q_{c1}(i,1)=0, where i=2,3,4, \ldots, n

Q_{c1}(i,j) = Q_{c1}(i-1,j) + Q_{c1}(i-1,j-1).

Thus, when A is the companion form of a polynomial of any degree n, then P_{c1}(P_{1}AQ_{1})Q_{c1} gives H(s+1) in the companion form.

Step 2:

H(s+1) \rightarrow (\frac{s^{3}}{b_{0}}).H(1/s + 1) (scaled inversion).

Let H(s+1) = s^{3} + b_{2}s^{2}+b_{1}s+b_{0} where b_{0}=a_{0} + a_{1} + a_{2} +1, b_{1} = a_{1} + 2a_{2}+3, b_{2}=a_{2}+3.

The scaling of the entire polynomial by b_{0} ensures that the polynomial generated is monic and hence, can be represented in the companion form.

The following transformation is sought: C = \left[ \begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -b_{0} & -b_{1} & -b_{2} \end{array}\right] \rightarrow D = \left[\begin{array}{ccc} 0 & 1 & 0\\ 0 & 0 & 1 \\ -1/b_{0} & (-b_{2}/b_{0}) & (-b_{1}/b_{0})\end{array}\right]

The row and column operator matrices P_{2} and Q_{2} respectively are:

P_{2}=\left[\begin{array}{ccc} b_{1}/b_{2} & 0 & 0 \\ 0 & b_{2}/b_{1} & 0 \\ 0 & 0 & 1/b_{0}\end{array}\right] and Q_{2} = \left[ \begin{array} {ccc}1/b_{0} & 0 & 0 \\ 0 & b_{2}/b_{1} & 0 \\ 0 & 0 & b_{1}/b_{2}\end{array}\right]

In general, P_{2} = \left[ \begin{array}{ccccc} b_{1}/b_{n-1} & 0 & 0 & 0 & \ldots \\ 0 & 0 & 0 & \ldots & \ldots \\ 0 & 0 & \ldots & \ldots & \ldots \\ 0 & \ldots & \ldots & b_{n-1}/b_{1} \\ \ldots & \ldots & \ldots & \ldots & 1/b_{0}\end{array}\right] and Q_{2} = \left[ \begin{array}{ccccc} 1/b_{0} & 0 & 0 & \ldots & 0\\ 0 & b_{n-1}/b_{0} & \ldots & \ldots & 0 \\ 0 & 0 & 0 & \ldots & 0 \\ 0 & \ldots & \ldots & \ldots & 0 \\ 0 & 0 & \ldots & \ldots & b_{1}/b_{n-1}\end{array}\right] , where both the general matrices P_{2} and Q_{2} are of dimensions n \times n.

So, we get D=P_{2}A_{1}Q_{2} = \left[ \begin{array}{ccccc} 0 & 1 & 0 & \ldots & 0 \\ 0 & 0 & 1 \ldots & 0 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1/b_{0} & (-b_{n-1}/b_{0}) & (-b_{n-2}/b_{0}) & \ldots & (-b_{1}/b_{0})\end{array}\right], which is also a matrix of dimension n \times n.

Step 3:

(s^{3}/b_{0})H(1/s + 1) \rightarrow (s^{3}/b_{0})H(1 + k/s), with k=2 (scaling of roots)

If F(x) = a_{n}x^{n} + a_{n-1}x^{n-1} + a_{n-2}x^{n-2} + \ldots + s_{0}, then F(x/k) = a_{n}x^{n} + ka_{n-1}x^{n-1} + k^{2}a_{n-2}x^{n-2}+ \ldots + a_{0}k^{n}.

The following transformation is sought:

D = \left[\begin{array}{ccc}0 & 1 & 0 \\ 0 & 0 & 1 \\ -1/b_{0} & (-b_{2}/b_{0}) & (-b_{1}/b_{0}) \end{array}\right] and E = \left[ \begin{array}{ccc}0 & 1 & 0 \\ 0 & 0 & 1 \\ -k^{3}/b_{0} & (-k^{2}b_{2}/b_{0}) & (-kb_{1}/b_{0})\end{array}\right].

The row and column operators, P_{3} and Q_{3}, respectively are:

P_{3}= \left[ \begin{array}{ccc} 1/k^{2} & 0 & 0 \\ 0 & 1/k & 0 \\ 0 & 0 &\  1\end{array}\right], and Q_{3}=\left[ \begin{array}{ccc} k^{3} & 0 & 0 \\ 0 & k^{2} & 0 \\ 0 & 0 & k \end{array}\right]

In general, P_{3} = \left[ \begin{array}{ccccc}1/k^{n-1} & 0 & 0 & \ldots & 0 \\ 0 & 1/k^{n-2} & 0 & \ldots & 0 \\ 0 & \ldots & \ldots & \ldots & 0 \\ 0 & \vdots & \vdots & \vdots & 1/k \\ 0 & \ldots & \ldots & \ldots & 1\end{array}\right], and Q_{3} = \left[ \begin{array}{ccccc}k^{n} & 0 & 0 & \ldots & 0 \\ 0 & k^{n-1} & 0 & \ldots & 0 \\ 0 & 0 & \ldots & \ldots & 0\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots &  k\end{array}\right], where the general matrices P_{3} and Q_{3} are both of dimensions n \times n

and E=P_{3}A_{2}Q_{3} = \left[ \begin{array}{ccccc}0 & 1 & 0 & \ldots & 0 \\ 0 & 0 & 1 & \ldots & 0 \\ 0 & \ldots & \ldots & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ (-k^{n}/b_{0}) & (-b_{n-1}k^{n-1}/b_{0}) & \ldots & \ldots & (-kb_{1}/b_{0})\end{array}\right]

Step 4:

(s^{3}/b_{0})H(1+k/s) \rightarrow (s^{3}/b_{0})H(1+k/(s-1)) (increasing of roots by one)

For the third degree case, the following transformation is sought:

E = \left[ \begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -A_{0} & -A_{1} & -A_{2}\end{array}\right] and F = \left[ \begin{array}{ccc}0 & 1 & 0 \\ 0 & 0 & 1 \\ -B_{0} & -B_{1} & -B_{2}\end{array}\right] and G = \left[ \begin{array}{ccc}0 & 1 & 0\\ 0 & 0 & 1 \\ (-A_{0}+A_{1}-A_{2}+1) & (-A_{1}+2A_{2}-3) & (3-A_{2}) \end{array}\right],

where B_{2}=A_{2}-1, B_{1}=A_{1}-B_{2}, B_{0}=A_{0}-B_{1}.

The row and column operators P_{4} and Q_{4} are :

P_{4}= \left[ \begin{array}{ccc} (A_{0}-A_{1})/A_{0} & (-A_{2}/A _{0}) & (-1/A_{0}) \\ 1 & 1 & 0 \\ 0 & 1 & 1\end{array}\right] and Q_{4}=\left[ \begin{array}{ccc} 1 & 0 & 0 \\ -1 & 1 & 0 \\ 1 & -1 & 1\end{array}\right]

In general, P_{4} = \left[ \begin{array}{ccccc} (C_{0}-C_{1})/C_{0} & (-C_{2}/C_{0}) & (-C_{3}/C_{0}) & \ldots & (1/C_{0}) \\ 1 & 1 & \ldots & \ldots & 0\\ 0 & 1 & 1 & \ldots & 0\\0 & \ldots & \ldots & 1 & 1 \\\vdots & \vdots & \vdots & \vdots & 1 \end{array}\right], and Q_{4}=\left[ \begin{array}{ccccc} 1 & 0 & \ldots & \ldots & 0 \\ -1 & 1 & 0 & \ldots & 0 \\ 1 & -1 & 1 & 0 & \ldots\\ \ldots & \ldots & \ldots & \ldots & 1 \\ \vdots & \vdots & \vdots & \vdots & 1\end{array}\right], both the general matrices P_{4} and Q_{4} being of dimensions n \times n.

Where C_{0}=k^{n}/b_{0} and C_{j}=k^{n-j}b_{n-j}/b_{0}, where j=1, 2, 3, \ldots, (n-1) and Q_{4} is a lower triangular matrix where q_{4}(i,j) = (-1)^{i+j}.

In general, we have F = P_{4}A_{3}Q_{4} = \left[ \begin{array}{ccccc} 0 & 1 & 0 & \ldots & 0 \\ 0 & 0 & 1 & \ldots & 0 \\ \ldots & \ldots & \ldots & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ -\hat{C_{0}} & -\hat{C_{1}} & -\hat{C_{2}} & \ldots & -\hat{C_{n-1}}\end{array}\right]

where -\hat{C_{0}} = -C_{0} + C_{1} - C_{2} + \ldots + (-1)^{n-1}

where -\hat{C_{1}} = -C_{1} + C_{2} - C_{3} + \ldots - (-1)^{n-1}

where -\hat{C_{2}} = -C_{2} + C_{3} - C_{4} + \ldots + (-1)^{n-1}

and so on \vdots

Now, the transformation F \rightarrow G is to be achieved. The row and column operators, P_{c4} and Q_{c4}, respectively are: P_{c4}=\left[ \begin{array}{ccc} 1 & 0 & 0 \\ 1 & 1 & 0\\ 1 & 2 & 1\end{array}\right] and Q_{c4}=\left[ \begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -1 & 1\end{array}\right]

In general, P_{c4}=\left[ \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0\\ 1 & 1 & 0 & \ldots & 0 \\ 1 & 2 & 1 & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & \ldots & \ldots & \ldots & 1\end{array}\right] and Q_{c4}= \left[ \begin{array} {ccccc} 1 & 0 & \ldots & \ldots & 0 \\ 0 & 1 & 0 & \ldots & 0 \\ 0 & -1  & 1 & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \ldots & \ldots & \ldots & \ldots & 1\end{array}\right], where both the general matrices P_{c4} and Q_{c4} are of dimensions n \times n.

P_{c4}, a lower triangular matrix can be easily generated as:

p_{c4}(i,i)=1; p_{c4}(i,j) = p_{c4}(i-1,j) + p_{c4}(i-1,j-1).

and Q_{c4}, also a lower triangular matrix can be easily generated as:

q_{c4}(i,i)=1; and q_{c4}(i,1)=0 when i \neq 1; and

q_{c4}(i,j) = [|| Q_{c4}(i-1,j)+ |Q_{c4}(i-1,j-1)|] \times (-1)^{i+j}.

Thus, P_{c4}(P_{4}EQ_{4})Q_{c4} finishes the process of bilinear transformation. Steps 2 and 3 can be combined so that the algorithm reduces to three steps only.

If the original polynomial is non-monic (that is, a_{n} \neq 1), then multiplying  the final tranformed polynomial by b_{0}a_{n} restores it to the standard form.

III. Stability considerations:

In the \omega plane, the Schwarz canonical approach can be applied as an algorithm directly to the canonical form of bilinear transformation of the polynomial obtained previously because a companion form is a non-derogatory matrix.

IV. An Example:

F_{0}(z) = 2z^{4} + 4z^{3} + 6z^{2} + 5z +1

F(z) = z^{4} + 2z^{3} + 3z^{2} + 2.5z + (0.5)z

b_{0} = a_{0} + a_{1} + a_{2} + 1 =9

b_{1} = a_{1} + a_{2} + a_{3} + 1 = 8.5

Similarly, b_{2}=6 and b_{3}=3

B=P_{1}AQ_{1}= \left[ \begin{array}{cccc} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\\ -9 & -8.5 & -6 & -3\end{array}\right]

C = \left[ \begin{array}{cccc}1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 1 & -2 & 1 & 0 \\ -1 & 3 & -3 & 1\end{array}\right] \times \left[ \begin{array}{cccc}0 & 1 & 0 & 0\\ 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 1\\ -9 & -8.5 & -6 & -3\end{array}\right] \times \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 1 & 2 & 1\end{array}\right]

Hence, C = \left[ \begin{array}{cccc} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ -9 & -18.5 & -15 & -6\end{array}\right]

H(s+1)=s^{4}+6s^{3} + 15s^{2}+18.5s+9

Steps 2 and 3:

E=\left[ \begin{array}{cccc} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\\ -2^{4}/9 & -2^{3}(6)/9 & -2^{2}(15)/9 & (-2)(18.5)/9\end{array}\right] = \left[ \begin{array}{cccc}0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\\ (-16/9) & (-48/9) & (-60/9) & (-37/9)\end{array}\right]

Step 4:

C_{0}=16/9, C_{1}=48/9, 60/9, 37/9.

-\hat{C_{0}}=0, -\hat{C_{1}}=-16/9, -\hat{C_{2}}=-32/9, -\hat{C_{3}}=-28/9

P_{c4}(P_{4}EQ_{4})Q_{c4}= \left[ \begin{array}{cccc}1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 2 & 1 & 0 \\ 1 & 3 & 3 & 1\end{array}\right] \times \left[ \begin{array}{cccc} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\\ 0 & -16/9 & -32/9 & -28/9\end{array}\right] \times \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 1 & -2 & 1\end{array}\right], so we get

P_{c4}(P_{4}EQ_{4})Q_{c4}=\left[ \begin{array}{cccc}0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\\ 0 & -3/9 & -3/9 & -1/9\end{array}\right]

The final monic polynomial is s^{4} + (1/9)s^{3}+(3/9)s^{2} + (3/9)s, and multiplying it by b_{0}a_{n}, that is, 9 \times 2 =18, restores it to the non-monic form:

18s^{4}+2s^{3}+6s^{2}+6s

V. Conclusion:

Since the operator matrices have lesser non-zero elements, storage requirements are lesser. The computational complexity should reduce for higher-order systems because the non-zero elements lesser manipulations are also lesser, besides lesser storage requirements. Additionally, the second and third steps can be combined giving a three step method only. Thus, the algorithm easily achieves bilinear transformation, especially, for higher systems compared to other available methods hitherto.

VI. References:

  1. Davies, A.C., “Bilinear Transformation of Polynomials”, IEEE Automatic Control, Nov. 1974.
  2. Barnett and Storey, “Matrix Methods in Stability Theory”, Thomas Nelson and Sons Ltd.
  3. Datta, B. N., “A Solution of the Unit Circle Problem via the Schwarz Canonical Form”, IEEE Automatic Control, Volume AC 27, No. 3, June 1982.
  4. Parthsarthy R., and Jaysimha, K. N., “Bilinear Transformation by Synthetic Division”, IEEE Automatic Control, Volume AC 29, No. 6, June 1986.
  5. Jury, E.I., “Theory and Applications of the z-Transform Method”, John Wiley and Sons Inc., 1984.