Matrix Functions for Digital Signal Processing. C Code and Octave Script.
1. Abbreviation
2. Introduction
A mathematical model of an object (plant) can be described by a system of differential equations. Let the model be described by a system of linear differential equations:
Equations (16) can be written in a simpler matrix form:
Now let’s consider a onedimensional (nonmatrix) version of equation (19):
where x(t) state function, u(t) control function
The solution of equation (25) has:
In this case the first term of the solution (26):
The second term of the solution:
The solution of the matrix equation (19) has a similar form as the onedimensional equation (25) and is written as:
Unfortunately it is impossible to write a simple formula for Φ(t,τ) as it was in the onedimensional case (27). Here are a few properties for this matrix:
If the matrix F in (19) does not depend on time then the matrix equation will take the form:
And solution of equation (38):
Note.

Thus in order to find a solution to equation (38) we must be able to calculate the matrix exponential function (40)!
So far we have been considering a continuous model. In the case of DSP it becomes necessary to write down the solution (33) for discrete moments of time:
Then using (33a):
We assume that the control u(t) remains constant for a small Ts time:
And introducing notation:
Then:
(47) is used to describe discrete linear dynamical systems in DSP applications!
For the case F independent of time, formulas (46) and (47) will be rewritten:
Conclusions

It becomes necessary to calculate matrix functions in DSP tasks for describing dynamic systems

The matrix exponential function is an important task for DSP

The calculation of the inverse matrix and matrix functions implies the use of square nxn matrices
3. Example
Consider as an example L (induction), R (resistor), C (capacitor) a circuit whose input voltage is u(t) and it is required to calculate the output voltage u_out(t). See the Figure 31.
Figure 31: LRC circuit
This circuit is described by a second order linear differential equation:
Further introducing the usual notation:
Rewrite (50) as:
Now let’s introduce the variables of the state vector:
And rewrite (52) in matrix form (38):
The solution (41) will take the form:
Take the following values for our example:
Then
Below I will show that:
Then:
Now we solve this task using a discrete form (48). Let’s choose the sampling time:
Equations (59) and (61) show the reaction of our dynamical system on the step function (56). See the Figure 31:
Figure 31: Digital State Model: x1(t) = u_out(t). Step Response
4. Inverse Matrix
There is a need to calculate the inverse matrix in DSP tasks. See the (58). The inverse matrix was also in the article “Kalman Filter with an example. C code and Octave script” for calculating the Kalman gain:
The definition of the inverse matrix is given by the formula (7):
The inverse matrix can be calculated by the formula:
Note: The inverse matrix exists for nondegenerate (nonsingular) matrices: det(A) ≠ 0
Using (62)(66) we can calculate the inverse matrix (58) from our example:
Then
Now let’s check the correctness of the calculations:
4.1. Inverse Matrix. Cayley–Hamilton Method. Faddeev Algorithm.
The characteristic polynomial of the nxn square matrix A is called a polynomial:
Using the Cayley–Hamilton theorem which states:
It can be proved that the inverse matrix is calculated by the formula:
The formula (71) takes the simple form for a 2×2 matrix:
Using (72) we can calculate (58) again:
The formula (71) takes a complex form for large matrices. It is convenient to use D. K. Faddeev’s recurrent algorithm in this case:
Notes:

At the last step n of the Faddeev’s algorithm the matrix Bn is always zero
This result can check of the calculations correctness 
The Faddeev’s algorithm also calculate the determinant of the matrix (75) and its characteristic polynomial (76)
Now we calculate the inverse matrix, determinant and characteristic polynomial by the Faddeev method for the matrix A:
Then:
Note: see the [1] for more information
4.2. Inverse Matrix. Gauss–Jordan Elimination Method
Now introduce elementary operations on the rows of the matrix:
Now let there be some square matrix A for which it is necessary to find the inverse matrix. For example:
For it we compose an auxiliary rectangular matrix:
And then performing elementary operations on the Help matrix we bring it to the form:
The matrix B will be the desired inverse matrix for A:
Let’s continue the example (78):
As a result we get the inverse matrix:
Let’s check the calculations correctness:
5. Matrix Functions
5.1 Matrix Function: exp(At) using Laplace Transformation
We had to calculate the matrix exponential in the above example (57) in which there is a time variable:
To solve this task we will use the Laplace transform:
Here are some images of Laplace functions:
Table 51: Table of Originals and Images by Laplace
The Laplace image of the function (83) corresponds to the second line of Table 51. Writing it in matrix form:
Then:
Using the formula (72) we can calculate the inverse matrix (86):
Obviously the inverse Laplace transform from (87) gives us the value of the desired function:
Now we can calculate the inverse Laplace transform for each element of the matrix (88) using lines 5 and 6 from Table 51:
Finally we get:
Now we prove the formula (60) using (89)
5.2 Eigenvalues and Eigenvectors of the Matrix. Transformation of Matrices to Diagonal or Jordan Form
Any matrix A:
can be interpreted as an operator that maps vector ‘x’ to vector ‘y’:
If the matrix (operator) A maps the ‘x’ vector into itself then such a vector is called the eigenvector of the matrix (operator) A and the proportionality coefficient is called the eigenvalue of the matrix A:
where x – the eigenvector of the matrix A
λ – the eigenvalue of the matrix A
To find the eigenvalues of the matrix it is necessary to solve the equation:
It means to find the roots of the characteristic polynomial.
Further solving the matrix equation:
for ‘x’ we find the eigenvectors of the matrix.
Example
Let there be a matrix
Then such a matrix can be converted to a diagonal form while the eigenvalues of the matrix will be located on the diagonal. An auxiliary matrix V is constructed for such transformation the columns of which are the eigenvectors of this matrix:
Using the matrix V we obtain a diagonal matrix:
Or in the back direction:
Let’s make a matrix V for our example from eigenvectors:
By the formula (72) we calculate the inverse matrix:
In the back direction:
Now the eigenvalues and vectors are calculated for the matrix (60) from the DSP example.
Then
In the back direction:
It means that the characteristic polynomial has multiple (repeated) roots. Then such a matrix can be transformed to Jordan (quasidiagonal) form. For example, let the matrix have elementary divisors:
Then the Jordan form will take the form:
Elementary blocks of the Jordan form can be distinguished:
To convert the matrix to Jordan form an auxiliary matrix T is constructed:
The transformation is done according to the same formulas as for the diagonal matrix case:
The algorithm for constructing the matrix T is quite hard and I will not give it here. See the [1] for more information.
5.3 Matrix Functions. Taylor and Maclaurin Series
Matrix functions can be determined using the Taylor and Maclaurin series.
Taylor ‘s series:
Let ‘s write the Maclaurin series in matrix form for few functions:
Take for example a diagonal matrix:
and substitute it into the formula (110):
Similarly substituting (114) into (112) and (113):
In general for an arbitrary function f and a diagonal matrix D:
Example
We can write for the zero matrix:
Example
For an identity matrix:
Now consider Jordan block (matrix) of the form:
Using (109) it can be proved that the matrix function from (121):
Example
Substituting (98) or (108) into (110) we can get a general result:
Example
Using the example for the matrix (95) and the result (100), as well as the formula (124):
5.4 Matrix Function: exp(A). Pade Aproximation
Notes:

Using eigenvalues and eigenvectors of the matrix to calculate matrix functions with the (124), (125) is a timeconsuming task.

Using a Maclaurin (or Taylor) series and limiting it to a finite terms number also requires a lot of computational costs.
The use of Pade Approximation is more effective. For the exponential function this approximation has the form:
Now let’s rewrite (126) in matrix form:
This approximation works well in the case when the matrix norm (see the (13)(15)) is small. In the future, we will require that the norm of the matrix is less or equal than to 0.5:
Note: next I will use the matrix norm 1, see the (13)
If (130) is fulfilled, then the relative approximation error (see the [3], [4]) can be estimated by the formula:
The matrix formula (129) for q=4:
If the matrix has a norm greater than 0.5, then the norm of the matrix decreases as follows:
Finally the algorithm for matrix exponential function calculating will take the form:
Example
Let’s take as an example the matrix from formula (60):
Notes:

The result is obtained with good accuracy

More information about it you can see in the [3], [4]
6. Conclusions

The task of calculating a matrix function arises when solving a system of differential equations (16)

Matrix functions are calculated only for square matrices

The classical method of calculating a matrix function is performed using formulas (124), (125), which is based on the matrix transform to a diagonal or Jordan form. This method is laborious and difficult to apply in practice

If the matrix function depends on time, then the Laplace transform can be applied to calculate the matrix function. See the example using formulas (85)(89)

The simplest way to calculate a matrix function is Pade approximation. See the example for calculating the matrix exponential function with formulas (135) – (137)

In the article I did not describe two more ways to determine and calculate a matrix function

using the minimal annulling polynomial of the matrix, see the Appendix 1

using the residues of a meromorphic function (complex analysis), see the Appendix 2


Frequent task in matrix arithmetic is the inverse matrix calculation. In this case the Faddeev algorithm is effective, see the (73)(76). You can use the formula (72) for a simple variant of 2×2 matrices
7. Octave GNU MatrixFunctionsSupport.m
The Octave script supports all plots of the article, generates the tests for the C code and calculates the Matrix Functions.
Script Tests:
% Test 1
% expm(0), exp(I)
% Test 2
% Pade Approximation Example: exp(C)=Rq(C) with q=5
% Without Matrix Norm Check
% Test 3
% The matrices exp(D) and D are commuting
% Test 4
% Faddeev Algorithm Example
% Test 5
% Calculate the inverse matrix
% Support the Gauss Algorithm from the article
% Test 6
% Model of the L,C,R circuit. See the article.
% d^2(u_out(t))/dt^2 + 10 d(u_out(t))/dt + 100 = 100 u(t)
% Variant 1 (test 6)
% Continuous Model
% for the equation: d^2(u_out(t))/dt^2 + 10 d(u_out(t))/dt + 100 = 100 u(t)
% In Laplace form: 100/(s^2+10s+100)
% Variant 2 (test 6)
% Continuous State Model
% for the equation: d^2(u_out(t))/dt^2 + 10 d(u_out(t))/dt + 100 = 100 u(t)
% Variant 3 (test 6)
% Discrete State Model
% for the equation: d^2(u_out(t))/dt^2 + 10 d(u_out(t))/dt + 100 = 100 u(t)
% Test 7
% [0 0.01;1 0.1]
% Calculate the matrix eigenvalues
% Transform and Diagonal matrices
% Test 8
% expm([1 2;3 2])
% article example
% Test 9
% Pade approximation exp(x)=Rq(x) for q_val=4
% Scalar Case
% Test 10
% Pade approximation for q_val=4
% Matrix Case
% Article example
8. C Code file MatrixFunctions.c/h
The file is supported the matrix with the 2×2 size.
Functions:
/*
Matrix Inverse
Input
float* matrix_in – pointer on the input float matrix_in[4] with the 2×2 size
float* matrix_out – pointer on the output float matrix_out[4] with the 2×2 size
Note: A11 = matrix[0], A12 = matrix[1], A21 = matrix[2], A22 = matrix[3]
Return void
matrix_out – inverse of the matrix_in: matrix_out = (matrix_in)^1
Note: if the matrix_in is not invertible (determinant=0)
then the zero matrix_out: [0.0 0.0; 0.0 0.0] will be returned
*/
void MatrixInverseFunc(float* matrix_in, float* matrix_out);
/*
Matrix Add
Input
float* matrix_in1 – pointer on the input float matrix_in1[4] with the 2×2 size
float* matrix_in2 – pointer on the input float matrix_in2[4] with the 2×2 size
float* matrix_out – pointer on the output float matrix_out[4] with the 2×2 size
Note: A11 = matrix[0], A12 = matrix[1], A21 = matrix[2], A22 = matrix[3]
Return void
matrix_out = matrix_in1 + matrix_in2
*/
void MatrixAddFunc(float* matrix_in1, float* matrix_in2, float* matrix_out);
/*
Matrix Subtraction
Input
float* matrix_in1 – pointer on the input float matrix_in1[4] with the 2×2 size
float* matrix_in2 – pointer on the input float matrix_in2[4] with the 2×2 size
float* matrix_out – pointer on the output float matrix_out[4] with the 2×2 size
Note: A11 = matrix[0], A12 = matrix[1], A21 = matrix[2], A22 = matrix[3]
Return void
matrix_out = matrix_in1 – matrix_in2
*/
void MatrixSubFunc(float* matrix_in1, float* matrix_in2, float* matrix_out);
/*
Matrix Multiply
Input
float* matrix_in1 – pointer on the input float matrix_in1[4] with the 2×2 size
float* matrix_in2 – pointer on the input float matrix_in2[4] with the 2×2 size
float* matrix_out – pointer on the output float matrix_out[4] with the 2×2 size
Note: A11 = matrix[0], A12 = matrix[1], A21 = matrix[2], A22 = matrix[3]
Return void
matrix_out = matrix_in1 x matrix_in2
*/
void MatrixMultiplyFunc(float* matrix_in1, float* matrix_in2, float* matrix_out);
/*
Number and Matrix Multiply
Input
float number – number
float* matrix_in – pointer on the input float matrix_in[4] with the 2×2 size
float* matrix_out – pointer on the output float matrix_out[4] with the 2×2 size
Note: A11 = matrix[0], A12 = matrix[1], A21 = matrix[2], A22 = matrix[3]
Return void
matrix_out = number x matrix_in2
*/
void NumberAndMatrixMultiplyFunc(float number, float* matrix_in, float* matrix_out);
/*
Matrix Norm I
Input
float* matrix_in – pointer on the input float matrix_in[4] with the 2×2 size
Note: A11 = matrix[0], A12 = matrix[1], A21 = matrix[2], A22 = matrix[3]
Return float – Matrix Norm I = max(ABS(A11)+ABS(A21), ABS(A12)+ABS(A22))
*/
float MatrixNormI_Func(float* matrix_in);
/*
Matrix Exponential Function
Input
float* matrix_in – pointer on the input float matrix_in1[4] with the 2×2 size
float* matrix_out – pointer on the output float matrix_out[4] with the 2×2 size
Note: A11 = matrix[0], A12 = matrix[1], A21 = matrix[2], A22 = matrix[3]
Return void
matrix_out = exp(matrix_in)
*/
void MatrixExpFunc(float* matrix_in, float* matrix_out);