Evaluation of loss functions for discrete time systems with a pulse transfer function B(z)/A(z).

signal_variance.py

import numpy as np
from numpy.polynomial.polynomial import Polynomial
import copy

# Define the coefficients of the polynomials A(z) and B(z)

# Introduction to Stochastic Control theory. Karl J. Astrom example
# A_coeffs = [1, 0.7, 0.5, -0.3]
# B_coeffs = [1, 0.3, 0.2, 0.1]

# 
A_coeffs = [1, 0.4, 0.1]
B_coeffs = [1, 0.9, 0.8]

A_coeffs_reverse = copy.deepcopy(A_coeffs)
B_coeffs_reverse = copy.deepcopy(B_coeffs)
A_coeffs_reverse.reverse()
B_coeffs_reverse.reverse()

# Create polynomial objects
A = Polynomial(A_coeffs_reverse)
B = Polynomial(B_coeffs_reverse)

print(f"A(x) = {A}")
print(f"B(x) = {B}")
# Check the roots of A(z) to ensure all are inside the unit circle
roots = A.roots()
if np.any(np.abs(roots) >= 1):
    raise ValueError("A(z) has roots on or outside the unit circle")

# Check the leadi2ng coefficient of A(z)
if A_coeffs_reverse[-1] <= 0:
    raise ValueError("The leading coefficient of A(z) is non-positive")

# Function to evaluate the integral I based on recursive formulas
def evaluate_integral(A_coeffs, B_coeffs):
    I = 0
    order = len(A_coeffs)
    # print(f"order {order}")
    A_table = np.zeros(((2*order)-1 , order+1))
    B_table = np.zeros(((2*order)-1 , order+1))
    A_table[0,0:order] = A_coeffs
    A_table[1,0:order] = A_coeffs_reverse
    B_table[0,0:order] = B_coeffs
    B_table[1,0:order] = A_coeffs_reverse
    # print(f"table shape {A_table.shape}")
    I += B_table[0,-2]**2 / A_table[0, 0]
    for k in range(order):
        A_table[2*k+1, -1] = A_table[2*k+1, 0] / A_table[2*k, 0]
        B_table[2*k+1, -1] = B_table[2*k, order-k-1] / B_table[2*k+1, order-k-1]

        for i in range(order-k-1):
            A_table[(2*k)+2, i] = A_table[2*k, i] - (A_table[2*k+1, -1] * A_table[2*k+1, i])
            B_table[(2*k)+2, i] = B_table[2*k, i] - (B_table[2*k+1, -1] * B_table[2*k+1, i])

        I += B_table[(2*k)+2, 0:(order-k-1)][-1]**2 / A_table[(2*k)+2, 0]

        is_last_k = (2*k)+3 >= A_table.shape[0]
        if is_last_k:
            break
        A_table[(2*k)+3, 0:(order-k-1)] = np.flip(A_table[(2*k)+2, 0:(order-k-1)])
        B_table[(2*k)+3, 0:(order-k-1)] = np.flip(A_table[(2*k)+2, 0:(order-k-1)])
        
            
        
    print("A-table :")
    print(A_table)
    print("\n")
    print("B-table :")
    print(B_table)

    I *= 1/A_table[0, 0]
    return I


# Evaluate the integral
I = evaluate_integral(A_coeffs, B_coeffs)
print("Evaluated integral:", I)