utilities.py

#! /usr/bin/env python3
###############################################################################
#
# Copyright 2024
# ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

###############################################################################

# Author: Eric Kooistra
# Purpose: Utilities and functions for DSP
# Description:
#
# References:
# [1] dsp_study_erko.txt

import numpy as np
from sys import exit

c_atol = 1e-15
c_rtol = 1e-8   # 1/2**32 = 2.3e-10


###############################################################################
# Utilities
###############################################################################

def verify_result(result):
    if result:
        print('PASSED')
    else:
        exit('FAILED')


def ceil_div(num, den):
    """ Return integer ceil value of num / den """
    return int(np.ceil(num / den))


def ceil_log2(num):
    """ Return integer ceil value of log2(num) """
    return int(np.ceil(np.log2(num)))


def ceil_pow2(num):
    """ Return power of 2 value that is equal or greater than num """
    return 2**ceil_log2(num)


def pow_db(volts):
    """Voltage to power in dB"""
    return 20 * np.log10(np.abs(volts) + c_atol)


def snr_db(signalData, noiseData):
    """Signal to noise ration in dB"""
    return pow_db(np.std(signalData) / (np.std(noiseData) + c_atol))


def is_integer_value(value):
    """Return true when value is sufficiently close to its integer value."""
    return np.isclose(value, np.round(value))


def is_even(n):
    """Return True if n is even, else False when odd.

    For all n in Z.
    """
    return n % 2 == 0


def is_odd(n):
    """Return True if n is odd, else False when even.

    For all n in Z, because result of n % +2 is 0 or +1.
    """
    return n % 2 == 1


def is_symmetrical(x, anti=False):
    """Return True when x[n] = +-x[N-1 - n], within tolerances, else False."""
    rtol = c_rtol
    atol = np.min(np.abs(x[np.nonzero(x)])) * rtol
    n = len(x)
    h = n // 2
    if is_even(n):
        if anti:
            return np.allclose(x[0:h], np.flip(-x[h:]), rtol=rtol, atol=atol)
        else:
            return np.allclose(x[0:h], np.flip(x[h:]), rtol=rtol, atol=atol)
    else:
        if anti:
            return np.allclose(x[0:h], np.flip(-x[h + 1:]), rtol=rtol, atol=atol) and np.abs(x[h]) < atol
        else:
            return np.allclose(x[0:h], np.flip(x[h + 1:]), rtol=rtol, atol=atol)


def read_coefficients_file(filepathname):
    coefs = []
    with open(filepathname, 'r') as fp:
        for line in fp:
            if line.strip():  # skip empty line
                s = float(line)   # one coef per line
                coefs.append(s)
    return np.array(coefs)


def one_bit_quantizer(x):
    """Quantize 0 and positive x to +1 and negative x to -1."""
    return np.signbit(x) * -1 + 2


def impulse_at_zero_crossing(x):
    """Create signed impulse at zero crossings of x."""
    diff = np.diff(one_bit_quantizer(x))
    return np.concatenate((np.array([0]), diff))