dsp.py

"""Digital Signal Processing capabilities for amodem."""
import functools

import numpy as np

import common
import exceptions
from kivy.logger import Logger as log


class FIR:
    def __init__(self, h):
        self.h = np.array(h)
        self.x_state = [0] * len(self.h)

    def __call__(self, x):
        x_ = self.x_state
        h = self.h
        for v in x:
            x_ = [v] + x_[:-1]
            yield np.dot(x_, h)
        self.x_state = x_

class Demux:
    def __init__(self, sampler, omegas, Nsym):
        self.Nsym = Nsym
        self.filters = [exp_iwt(-w, Nsym) / (0.5*self.Nsym) for w in omegas]
        self.filters = np.array(self.filters)
        self.sampler = sampler

    def __iter__(self):
        return self

    def next(self):
        frame = self.sampler.take(size=self.Nsym)
        if len(frame) == self.Nsym:
            symbol = np.dot(self.filters, frame)
            #log.info(symbol)
            magnitude = np.abs(symbol) #calculate distance between the symbol and the unit circle
            return symbol
        raise StopIteration

    __next__ = next

@functools.lru_cache(maxsize=64, typed=False) #improve the efficiency of this function using caching
def exp_iwt(omega, n):
    return np.exp(1j * omega * np.arange(n))

def norm(x):
    return np.sqrt(np.dot(x.conj(), x).real)

def rms(x):
    return np.mean(np.abs(x) ** 2, axis=0) ** 0.5

def coherence(x, omega):
    n = len(x)
    Hc = exp_iwt(-omega, n) / np.sqrt(0.5*n)
    norm_x = norm(x)
    
    if not norm_x:
        return 0.0
    return np.dot(Hc, x) / norm_x

def linear_regression(x, y):
    """ Find (a,b) such that y = a*x + b. """
    x = np.array(x)
    y = np.array(y)
    mean_x = np.mean(x)
    mean_y = np.mean(y)
    x_ = x - mean_x
    y_ = y - mean_y
    a = np.dot(y_, x_) / np.dot(x_, x_)
    b = mean_y - a * mean_x
    return a, b


class MODEM:

    def __init__(self, symbols):
        self.encode_map = {}
        symbols = np.array(list(symbols))
        bits_per_symbol = np.log2(len(symbols))
        bits_per_symbol = np.round(bits_per_symbol)
        N = (2 ** bits_per_symbol)
        assert N == len(symbols)
        bits_per_symbol = int(bits_per_symbol)

        #generate the bits that each symbol will represent
        for i, v in enumerate(symbols):
            bits = [int(i & (1 << j) != 0) for j in range(bits_per_symbol)]
            #print(bits)
            self.encode_map[tuple(bits)] = v

        self.symbols = symbols
        self.bits_per_symbol = bits_per_symbol
        #print('encode map')
        #print(self.encode_map)
        
        bits_map = dict(item[::-1] for item in self.encode_map.items())
        #print(bits_map)
        self.decode_list = [(s, bits_map[s]) for s in self.symbols]
        #print(self.symbols)
        #print(self.decode_list)
                
    def encode(self, bits):
        for bits_tuple in common.iterate(bits, self.bits_per_symbol, tuple):
            yield self.encode_map[bits_tuple]

    #take a stream of symbols and decode them into bits
    def decode(self, symbols, error_handler=None):
        """ Maximum-likelihood decoding, using naive nearest-neighbour. """
        symbols_vec = self.symbols
        _dec = self.decode_list
        for received in symbols:
            error = np.abs(symbols_vec - received)
            index = np.argmin(error)
            decoded, bits = _dec[index]
            if error_handler:
                error_handler(received=received, decoded=decoded)
            yield bits

##@brief convert an array of bits to a differentially encoded BPSK signal
##@param bits an array of 0 and 1's
##@param L the upsampling factor to be applied to the signal
##@param prev_bit the last differentially encoded bit from the previous call if the method is being called multiple times
def bits2baseband(bits,L,prev_bit=0):
    dif_encoded_bits = np.zeros(len(bits))
    for i in range(0,len(bits)):
        dif_encoded_bits[i] = (bits[i] + prev_bit)%2 
        prev_bit = dif_encoded_bits[i]


    toReturn = np.zeros(len(dif_encoded_bits)*L)
    for i in range(0,len(dif_encoded_bits)):
        #upsample the signal and apply NRZ encoding
        if (dif_encoded_bits[i] == 0):
            toReturn[i*L : (i*L)+L] = -1
        else:
            toReturn[i*L : (i*L)+L] = 1
    t = np.arange(start=0, stop=len(dif_encoded_bits)*L)
    return (toReturn,t,prev_bit)

def prbs(reg, poly, bits):
    """ Simple pseudo-random number generator. """
    mask = (1 << bits) - 1

    size = 0  # effective register size (in bits)
    while (poly >> size) > 1:
        size += 1

    while True:
        yield reg & mask
        reg = reg << 1
        if reg >> size:
            reg = reg ^ poly