veeresht/commpy

Digital Communication with Python

CommPy

CommPy is an open source toolkit implementing digital communications algorithms in Python using NumPy and SciPy.

Objectives

• To provide readable and useable implementations of algorithms used in the research, design and implementation of digital communication systems.

Available Features

Channel Coding

• Encoder for Convolutional Codes (Polynomial, Recursive Systematic). Supports all rates and puncture matrices.
• Viterbi Decoder for Convolutional Codes (Hard Decision Output).
• MAP Decoder for Convolutional Codes (Based on the BCJR algorithm).
• Encoder for a rate-1/3 systematic parallel concatenated Turbo Code.
• Turbo Decoder for a rate-1/3 systematic parallel concatenated turbo code (Based on the MAP decoder/BCJR algorithm).
• Binary Galois Field GF(2^m) with minimal polynomials and cyclotomic cosets.
• Create all possible generator polynomials for a (n,k) cyclic code.
• Random Interleavers and De-interleavers.
• Belief Propagation (BP) Decoder and triangular systematic encoder for LDPC Codes.

Channel Models

• SISO Channel with Rayleigh or Rician fading.
• MIMO Channel with Rayleigh or Rician fading.
• Binary Erasure Channel (BEC)
• Binary Symmetric Channel (BSC)
• Binary AWGN Channel (BAWGNC)

Wifi 802.11 Simulation Class

• A class to simulate the transmissions and receiving parameters of physical layer 802.11 (currently till VHT (ac)).

Filters

• Rectangular
• Raised Cosine (RC), Root Raised Cosine (RRC)
• Gaussian

Impairments

• Carrier Frequency Offset (CFO)

Modulation/Demodulation

• Phase Shift Keying (PSK)
• Quadrature Amplitude Modulation (QAM)
• OFDM Tx/Rx signal processing
• MIMO Maximum Likelihood (ML) Detection.
• MIMO K-best Schnorr-Euchner Detection.
• MIMO Best-First Detection.
• Convert channel matrix to Bit-level representation.
• Computation of LogLikelihood ratio using max-log approximation.

Sequences

• PN Sequence
• Zadoff-Chu (ZC) Sequence

Utilities

• Decimal to bit-array, bit-array to decimal.
• Hamming distance, Euclidean distance.
• Upsample
• Power of a discrete-time signal

Links

• Estimate the BER performance of a link model with Monte Carlo simulation.
• Link model object.
• Helper function for MIMO Iteration Detection and Decoding scheme.

FAQs

Why are you developing this?

During my coursework in communication theory and systems at UCSD, I realized that the best way to actually learn and understand the theory is to try and implement ''the Math'' in practice :). Having used Scipy before, I thought there should be a similar package for Digital Communications in Python. This is a start!

What programming languages do you use?

CommPy uses Python as its base programming language and python packages like NumPy, SciPy and Matplotlib.

How can I contribute?

Implement any feature you want and send me a pull request :). If you want to suggest new features or discuss anything related to CommPy, please get in touch with me ([email protected]).

How do I use CommPy?

Requirements/Dependencies

• python 3.2 or above
• numpy 1.10 or above
• scipy 0.15 or above
• matplotlib 1.4 or above
• nose 1.3 or above
• sympy 1.7 or above

Installation

• To use the released version on PyPi, use pip to install as follows::

$ pip install scikit-commpy

• To work with the development branch, clone from github and install as follows::

$ git clone https://github.com/veeresht/CommPy.git
$ cd CommPy
$ python setup.py install

• conda version is curently outdated but v0.3 is still available using::

$ conda install -c https://conda.binstar.org/veeresht scikit-commpy

Citing CommPy

If you use CommPy for a publication, presentation or a demo, a citation would be greatly appreciated. A citation example is presented here and we suggest to had the revision or version number and the date:

V. Taranalli, B. Trotobas, and contributors, "CommPy: Digital Communication with Python". [Online]. Available: github.com/veeresht/CommPy

I would also greatly appreciate your feedback if you have found CommPy useful. Just send me a mail: [email protected]

For more details on CommPy, please visit https://veeresht.info/CommPy/

commpy's Issues

rrcosfilter with Fs=3

I've noticed a bizzar frequency response from the rrcosfilter. I suspect this is an edge case due to the if/else in the module. The time waveform doesn't look unusual.

import matplotlib.pyplot as plt
import scipy.fftpack as fftpack
from commpy import rrcosfilter
import numpy as np
import scipy

def padded_centered_fft(taps,mintaps=1024,db=True,db_min=-300,fs=1):
	""" returns the fft of the zero padded taps, and centers it """
	if len(taps) < mintaps:
		taps=np.concatenate((taps,np.zeros(mintaps-len(taps))))
	dfreq, fft = fftpack.fftshift(fftpack.fftfreq(len(taps))),fftpack.fftshift(scipy.fft(taps))
	fft_abs=np.abs(fft)
	fft_db = 20*np.log10(fft_abs+max(np.abs(fft))*(10**((db_min-50)/20)) )
	fft_db[np.where(fft_db < db_min)]=db_min
	dfreq*=fs
	return (dfreq, fft_db) if db else (dfreq,fft)

plt.plot(*padded_centered_fft(rrcosfilter(N=400,alpha=0.15,Ts=1,Fs=3.01)[1]),label='3.01')
plt.plot(*padded_centered_fft(rrcosfilter(N=400,alpha=0.15,Ts=1,Fs=3)[1]),label='3')
plt.plot(*padded_centered_fft(rrcosfilter(N=400,alpha=0.15,Ts=1,Fs=2.99)[1]),':',label='2.99')
plt.legend(loc='best')
plt.show()

convcode.py throws a TypeError when running if installed numpy version is 1.12 or greater

Appears that as of 1.12, numpy float64s can no longer be used as indices. This was encountered on a system with Python 2.7.14 and numpy 1.13.3:

  File "/usr/lib/python2.7/site-packages/commpy/channelcoding/convcode.py", line 337, in conv_encode
    puncture_matrix[0:].sum()/np.size(puncture_matrix, 1), 'int')
TypeError: 'numpy.float64' object cannot be interpreted as an index

Rolling back to 1.11.0 shows a warning:

/usr/lib/python2.7/site-packages/commpy/channelcoding/convcode.py:337: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  puncture_matrix[0:].sum()/np.size(puncture_matrix, 1), 'int')
/home/threexc/Projects/Recorder/T9_FrameGenerator/make_symbol.py:1047: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  (_,curr_data_time,_,_) = make_symbol(scr_coded[i*ncbps:(i+1)*ncbps],rate_dict[rate][1],i+1)

NumPy log2 bug in QAMModem class

There's a sort-of bug in the QAMModem class (at least on some python implementations). From your example above:
self.num_bits_symbol returns 5 instead of 6 (for QAMModem(64))
log2(64) = 5.99999...
then int(log2(64)) = 5
Note: log(64)/log(2) = 6.0.
I'm using python-2.7.5.amd64 on Win10.

• reported by Keith Douglas (http://veeresht.info/blog/cfo.html)

Using conv encoder for rate other than 1/2

I am having trouble while using conv_encoder for rate other than 1/2. Also the only example that was given was for rate 1/2 encoder. It would really help a lot if you can also add examples related to using conv_encoder for rate other than 1/2 (ex. 2/3, 3/4 etc.)

Travis CI build failing for Python 3.3

• File "/home/travis/build/veeresht/CommPy/commpy/init.py", line 17, in
  from filters import *
  ImportError: No module named 'filters'
• commpy/channelcoding/algcode NotPython: Couldn't parse '/home/travis/build/veeresht/CommPy/commpy/channelcoding/algcode.py' as Python source: 'invalid syntax' at line 35
• commpy/channelcoding/gfields NotPython: Couldn't parse '/home/travis/build/veeresht/CommPy/commpy/channelcoding/gfields.py' as Python source: 'invalid syntax' at line 63
• commpy/channelcoding/tests/test_convcode NotPython: Couldn't parse '/home/travis/build/veeresht/CommPy/commpy/channelcoding/tests/test_convcode.py' as Python source: 'invalid token' at line 17

conv_encode exceeds the expected code rate

Please, could you resolve my misunderstanding:
when I try to use the method conv_encode() for the binary stream array([1., 0., 1., 0., 0., 1., 1.]) the code rate exceed the 1/2:

from numpy import *

memory = array([6])
g_matrix = array([[0o171, 0o133]])
trellis = Trellis(memory, g_matrix)
print(conv_encode([1., 0., 1., 0., 0., 1., 1.], trellis))

output: [1 1 1 0 0 0 0 1 1 1 0 1 1 0 0 0 1 1 1 1 0 1 1 0 1 1]

So, code rate k/n = 7 / 26

Do I use the considered method incorrectly?

Constellation mapping implemented for MQAM is not the Gray mapping

I have noted that the constellation mapping implemented is not actually the Gray mapping. Because of that, e.g. for high order QAM, there is a performance mismatch with respect to the theoretical BER vs SNR curves for the AWGN channel, if one tries to use the QAMModem methods to modulate and demodulate a sequence of bits and compute the BER afterward.

Wrong direct push & Python 2 in Travis

Veeresh,

I mistakenly made a direct push rather than a PR (bad configuration of my new damn IDE...). The new version is no longer compatible with python 2.

Could you please check that the push suits you (if not you can do a revert), I prefer not to revert myself in order not to add more mess in the tree. Indeed, as you will see the addition is quite basic.

You should also remove Python 2 from Travis to avoid compatibility errors that no longer make sense with the end of Python 2.

ImportError: No module named 'filters'

I just tried installing Commpy on my Macbook Pro (2015) using pip.

When I type import commpy I get the following error message (Python 3.5) :

/Library/Frameworks/Python.framework/Versions/3.5/bin/python3.5 "/Users/Eric/Dropbox/Python TIPE/Bruit/test pillow"
Traceback (most recent call last):
  File "/Users/User_name/path/file", line 1, in <module>
    import commpy.channels
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/commpy/__init__.py", line 17, in <module>
    from filters import *
ImportError: No module named 'filters'

Do anybody know how I can solve this issue ?

Thank you,

Eric

Warning messages in tests

Unittests features few warning messages, it would be good to fix them.

<me>/CommPy$ python commpy/tests/test_*
...<me>/CommPy/commpy/links.py:324: ComplexWarning: Casting complex values to real discards the imaginary part
  received_msg[receive_size * i:receive_size * (i + 1)] = \
.<me>/CommPy/commpy/modulation.py:633: RuntimeWarning: divide by zero encountered in true_divide
  return -LLR / (2 * noise_var)
.................................................<me>/CommPy/commpy/channels.py:239: RuntimeWarning: divide by zero encountered in long_scalars
  return absolute(self.fading_param[0]) ** 2 / absolute(self.fading_param[1])
<me>/CommPy/commpy/channels.py:239: RuntimeWarning: divide by zero encountered in double_scalars
  return absolute(self.fading_param[0]) ** 2 / absolute(self.fading_param[1])
...
----------------------------------------------------------------------
Ran 56 tests in 257.543s

OK

QAM Error

Hi. I'm using Anaconda (Python 3.5.0 - x64) on Linux.
When I try this example code:

from numpy import real, imag, arange
from numpy.fft import fft, ifft
from numpy.random import randint
from commpy.modulation import QAMModem
from commpy.impairments import add_frequency_offset
import matplotlib.pyplot as plt

# System Parameters 
# Sampling Frequency (Hz)
Fs = 15.36e6 
# Frequency Offset (Hz)
delta_f = 500
# FFT Size for OFDM 
N = 1024

# Generate message bits
msg_bits = randint(0, 2, 6144)
# 64-QAM Modulation of the message bits
qam64 = QAMModem(64)
tx_sc_symbols = qam64.modulate(msg_bits)
# Add Frequency offset to the modulated symbols
rx_sc_symbols = add_frequency_offset(tx_sc_symbols, Fs, delta_f)

# Ideal 64-QAM constellation
plt.figure(1)
plt.plot(real(tx_sc_symbols), imag(tx_sc_symbols), '.b')
plt.xlabel('I-component')
plt.ylabel('Q-component')
plt.xticks(arange(-10, 12, 2))
plt.yticks(arange(-10, 12, 2))
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.grid(True)

# 64-QAM constellation for Single Carrier with Frequency offset
plt.figure(2)
plt.plot(real(rx_sc_symbols), imag(rx_sc_symbols), '.b', alpha=0.5)
plt.xlabel('I-component')
plt.ylabel('Q-component')
plt.xticks(arange(-10, 12, 2))
plt.yticks(arange(-10, 12, 2))
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.grid(True)

# Generate OFDM signal (single symbol) using modulated symbols
tx_ofdm_waveform = ifft(tx_sc_symbols)
# Add frequency offset to the OFDM signal
rx_ofdm_waveform = add_frequency_offset(tx_ofdm_waveform, Fs, delta_f)
# OFDM Demodulation to extract the 64-QAM symbols
rx_ofdm_symbols = fft(rx_ofdm_waveform)

# 64-QAM Constellation for OFDM with Frequency offset
plt.figure(3)
plt.plot(real(rx_ofdm_symbols), imag(rx_ofdm_symbols), '.b', alpha=0.6)
plt.xlabel('I-component')
plt.ylabel('Q-component')
plt.xticks(arange(-10, 12, 2))
plt.yticks(arange(-10, 12, 2))
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.grid(True)

plt.show()

I receive this error:

File "...../modulation.py", line 45, in modulate, 
baseband_symbols = self.constellation[index_list]

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

Wrong documentation for conv_encode

The docstring says it expects a generator_matrix, however, a trellis is expected.

Also, "M" is not an input, while code_type and puncture_matrix is not documented.

PR #74 introduce a bug in viterbi_decode

I just checked the merged PR #74 and it comes out that it introduces a bug. You can reproduce the issue by running test_conv_encode_viterbi_decode in test_convcode.py.

This bug is undetected by Travis due to the flag -a '!slow' in the configuration file. I'm about to open a PR to remove this flag and we could discuss about this change there (see #77).

CommPy/.travis.yml

Line 34 in b104c99

- nosetests -a '!slow' --with-coverage --cover-package=commpy --logging-level=INFO

Let's now discuss on the bug introduced in #74. The new _compute_branch_metric with LRU caching is called as follow:

CommPy/commpy/channelcoding/convcode.py

Line 623 in b104c99

branch_metric = _compute_branch_metrics(decoding_type, tuple(r_codeword), tuple(i_codeword_array))

However, with a hard decoding type, the function calls .astype on the arguments. This is obviously not possible on tuples.

CommPy/commpy/channelcoding/convcode.py

Lines 575 to 578 in b104c99

	@functools.lru_cache(maxsize=128, typed=False)
	def _compute_branch_metrics(decoding_type, r_codeword, i_codeword_array):
	if decoding_type == 'hard':
	return hamming_dist(r_codeword.astype(int), i_codeword_array.astype(int))

It is not possible to just keep the tuples provided since hamming_dist relies on numpy capabilities.

CommPy/commpy/utilities.py

Lines 112 to 132 in b104c99

	def hamming_dist(in_bitarray_1, in_bitarray_2):
	"""
	Computes the Hamming distance between two NumPy arrays of bits (0 and 1).
	Parameters
	----------
	in_bit_array_1 : 1D ndarray of ints
	NumPy array of bits.
	in_bit_array_2 : 1D ndarray of ints
	NumPy array of bits.
	Returns
	-------
	distance : int
	Hamming distance between input bit arrays.
	"""
	distance = np.bitwise_xor(in_bitarray_1, in_bitarray_2).sum()
	return distance

I don't have time now to explore if this change should be reverted or if there is a better solution. Indeed, converting back and forth from numpy arrays to tuples may lead to a overhead greater than the gain obtained with the LRU cache. I'll try to have a look in the next days but I can't be sure on when it will be possible. @eSoares do you have any suggestions ?

Specification of decimal equivalent of polynomial (g_matrix entry) for trellis

Hey,

I want to know how to specify the g_matrix entries as decimals. The problem is, I have the generator polynomials, but don't know how to convert them to decimal numbers.

Say one of the polynomials is 1+D^2+D^3, for a memory=3 encoder. Does this correspond to a binary number of 1011 (LSB first format) or 1101 (MSB first format)? In both the examples given in the documentation, the binary number is symmetric.

Trying to Demodulate in soft mode fails

Trying to change the "hard" demodulate in the example in commpy/examples/conv_encode_decode.py to "soft" fails with AttributeError: 'QAMModem' object has no attribute 'symbol_mapping'.

Wrong output table for trellis in convcode.py (Probably in other parts too)

When creating the output table for Convolutional Codes, you have to notice that the generator matrix polynomials and shift registers don't match up in their direction (polynomials are expressed left-to-right, while the corresponding shift register is right-to-left. To calculate the output, you need to reverse either one first before multiplication.

A quick fix:
In lines 157-158 of convcode.py replace:
generator_array = dec2bitarray(g_matrix[l][r],
memory[l] + 1)
with:
generator_array = dec2bitarray(g_matrix[l][r],
memory[l] + 1)[::-1]
to reverse the generator's corresponding bit array to match with the shift register's order.

Note: I didn't check the feedback loop. It might have the same problem.

Gray Coding

I checked out 16-QAM modulation and saw that the constellation is not using gray coding. My question is that will be available in next versions?

filters.py - fix variable name in gaussian filter function

Fix variable name 'time_index' to 'time_idx'.

Reported by: Benjamin Sackenreuter

How to test pnsequence in sequences.py file

How to test commpy.sequences.py file

commpy.modulation.mimo_ml don't work with more than 2 receive antennas

The following line generates obviously an array with shape (2,...).

x_ideal = array([tile(constellation, m), repeat(constellation, m)])

README tutorial is broken for channel coding

README tutorial is broken

Please replace trellis = cc.Trellis(M, generator_matrix) # Trellis structure by trellis = cc.Trellis(np.array([M]), generator_matrix) # Trellis structure, or so it seems.

AttributeError                            Traceback (most recent call last)
<ipython-input-5-d8ad681885c7> in <module>
      3 m = np.array([L-1]) # number of delay elements
      4 generator_matrix = np.array([[0o171, 0o133]]) # generator branches
----> 5 trellis = cc.Trellis([M], generator_matrix) # Trellis structure

<ME>/CommPy/commpy/channelcoding/convcode.py in __init__(self, memory, g_matrix, feedback, code_type, polynomial_format)
    120         self.code_type = code_type
    121 
--> 122         self.total_memory = memory.sum()
    123         self.number_states = pow(2, self.total_memory)
    124         self.number_inputs = pow(2, self.k)

AttributeError: 'list' object has no attribute 'sum'

import numpy as np
import commpy.channelcoding.convcode as cc
import commpy.modulation as modulation

def BER_calc(a, b):
    num_ber = np.sum(np.abs(a - b))
    ber = np.mean(np.abs(a - b))
    return int(num_ber), ber

N = 100000 #number of symbols per the frame
message_bits = np.random.randint(0, 2, N) # message

M = 4 # modulation order (QPSK)
k = np.log2(M) #number of bit per modulation symbol
modem = modulation.PSKModem(M) # M-PSK modem initialization 
trellis = cc.Trellis(M, generator_matrix) # Trellis structure

rate = 1/2 # code rate
L = 7 # constraint length
m = np.array([L-1]) # number of delay elements
generator_matrix = np.array([[0o171, 0o133]]) # generator branches
trellis = cc.Trellis(M, generator_matrix) # Trellis structure

awgn noise is still not correct when input signals are complex values

I found 2 bugs from awgn(), using complex input signals.

https://github.com/veeresht/CommPy/blob/master/commpy/channels.py#L97
I think 'is' should'nt be used for comparing dtype, instead of '=='.

For example...

>>> a
array([ 0.+1.j,  2.+0.j,  3.+0.j])
>>> a.dtype
dtype('complex128')
>>> a.dtype is complex
False
>>> a.dtype == complex
True
>>> a.dtype is np.dtype('complex128')
True

In python reference , 'x is y is true if and only if x and y are the same object.'
so You want to use 'is', compare with np.dtype('complex128').

https://github.com/veeresht/CommPy/blob/master/commpy/channels.py#L93
The Average energy of input signals is calculated from this equation.

avg_energy = sum(input_signal * input_signal)/len(input_signal)

But in the case input_signal is complex, you must calculate signal's absolute value before you square it.

avg_energy = sum(abs(input_signal) * abs(input_signal))/len(input_signal)

help to improve library

Hello,
In case of a desire to adding new modulation and so on, I will be glad to help you.
Ashkan.

Wrong BER results?

Using the example "conv_encode_decode" and lowering the SNR to even negative values shows that BER have a limit of 0.5 (50%). Shouldn't this value by 1 (100%)?

Or I'm doing something wrong?

Thanks for the help,

AWGN noise_variance dependence on FEC code rate

Hey,
This isn't an issue with the code, but something I want to understand the theory behind, and I couldn't find a better place to ask.

CommPy/commpy/channels.py

Line 699 in 1d26b07

noise_variance = avg_energy / (2 * rate * snr_linear)

In an AWGN channel, why does the noise variance depend on the FEC code rate? If you can share some references as to where you found this, it will be very helpful. (I am trying to do a self-study course on Coding theory where I am trying to code and implement the concepts I learn... and I'm trying to implement different channel models to understand better.)

MAP decoder for convolutional codes

In the README it is stated, that there is a MAP decoder for convolutional codes, based on BCJR. But I can't find it. Only the one for turbo codes. Did I miss it?

Modification of documentation for convolutional encoding

The Trellis function has the argument for the generator polynomial marked as

g_matrix : 2D ndarray of ints (octal representation)

i.e. base 8 polynomial. This is extremely misleading as the actual argument expected is that of a standard base 10 integer (the function calls the utility function dec2bin). This doesn't cause problems in the examples as the polynomials used there are less than 8, but upon implementing anything with longer memory the output was all wrong and caused a big headache.

Suggestion: amend the documentation to denote g_matrix as base 10.

awgn noise is not correct

https://github.com/veeresht/CommPy/blob/master/commpy/channels.py#L98

The complex angle of the noise is always the same (45 degree angle), which is totally wrong. one way to do it is to take the code for the real-only noise, and then apply a random complex rotation. something like this (for a single sample):

angle = np.random.random_sample(1)
noise = noise * np.exp(1j_angle_2.0*np.pi)

CAZAC sequence is incorrect

The definition of the CAZAC sequence is imprecise leading to failure modes when seq_length is an even number. Following https://en.wikipedia.org/wiki/Zadoff%E2%80%93Chu_sequence and tests to verify the value of the auto-correlation the following code

CommPy/commpy/sequences.py

Line 92 in c7f202a

zcseq = exp((-1j * pi * u * arange(seq_length) * (arange(seq_length)+1)) / seq_length)

should be replaced by:

def zcsequence(u, seq_length, q):
    """
    Generate a Zadoff-Chu (ZC) sequence.
    Parameters
    ----------
    u : int
        Root index of the the ZC sequence.
    seq_length : int
        Length of the sequence to be generated. Usually a prime number.
    q : int
        Cyclic shift of the sequence
    Returns
    -------
    zcseq : 1D ndarray of complex floats
        ZC sequence generated.
    """
    assert u>0
    assert u<seq_length
    assert np.gcd(u,seq_length)==1
    for el in [u,seq_length,q]:
        assert float(el).is_integer()

    cf = seq_length%2
    n = arange(seq_length)
    zcseq = exp( -1j * pi * u * n * (n+cf+2.*q) / seq_length)
    return zcseq

Here is an example:

import numpy as np

from numpy import empty, exp, pi, arange, int8, fromiter, sum

def zcsequence(u, seq_length):
    """
    Generate a Zadoff-Chu (ZC) sequence.
    Parameters
    ----------
    u : int
        Root index of the the ZC sequence.
    seq_length : int
        Length of the sequence to be generated. Usually a prime number.
    Returns
    -------
    zcseq : 1D ndarray of complex floats
        ZC sequence generated.
    """
    zcseq = exp((-1j * pi * u * arange(seq_length) * (arange(seq_length)+1)) / seq_length)

    return zcseq
def zcsequence_corrected(u, seq_length, q=0):
    """
    Generate a Zadoff-Chu (ZC) sequence.
    Parameters
    ----------
    u : int
        Root index of the the ZC sequence.
    seq_length : int
        Length of the sequence to be generated. Usually a prime nu![example-failure](https://user-images.githubusercontent.com/3166436/84265615-62527e80-ab23-11ea-829d-e37df16d2527.png)mber.
    Returns
    -------
    zcseq : 1D ndarray of complex floats
        ZC sequence generated.
    """
    assert u>0
    assert u<seq_length
    assert np.gcd(u,seq_length)==1
    for el in [u,seq_length,q]:
        assert float(el).is_integer()

    cf = seq_length%2
    n = arange(seq_length)
    zcseq = exp( -1j * pi * u * n * (n+cf+2.*q) / seq_length)

    return zcseq

if __name__ == '__main__':
    
    import matplotlib.pyplot as plt
    
    u = 3
    seq_length = 20

    seq = zcsequence(u,seq_length)
    x_F = np.fft.fft(seq) / np.sqrt(seq.size)
    h_chan_est_f = (np.fft.ifft(np.conj(x_F) * x_F)*np.sqrt(seq.size)).T
    plt.plot(np.absolute(h_chan_est_f)**2,marker='o',label='before')
    
    seq = zcsequence_corrected(u,seq_length)
    x_F = np.fft.fft(seq) / np.sqrt(seq.size)
    h_chan_est_f = (np.fft.ifft(np.conj(x_F) * x_F)*np.sqrt(seq.size)).T
    plt.plot(np.absolute(h_chan_est_f)**2,marker='o',label='after')
    
    plt.xlabel('Delta symbol')
    plt.ylabel('auto-correlation')
    plt.legend()
    plt.grid()
    plt.show()

Plotting modulated signals?

Hello everyone,

I'm working on an encoding/modulating setup much like the example in https://github.com/veeresht/CommPy/blob/master/commpy/examples/conv_encode_decode.py. I think I have a good start on understanding the process, but my background in signals processing is a bit lacking. I'm hoping someone here can help me with the next steps.

For context, my goal is to form a relatively reliable communication scheme over a noisy (and potentially transcoded/compressed) audio channel. I've been referencing this library and another called amodem, but it seems like this library is missing the component for actually generating audio and the other library does not incorporate any error correcting codes and does not survive my tests through the audio channel.

I'm at a point where, with the code here, I can take binary data and generate complex QAM symbols. My understanding is that each QAM symbol represents the amplitude and and phase of one baud's signal (at the carrier's frequency). I'd like to graph this signal as a raw waveform, but I'm unsure how to do that. This is as far as I got, but it seems cumbersome to be calculating the amplitude and phase explicitly for each symbol:

f_w = 2*np.pi*carrier_hz
t_series = np.arange(0, 1/baud, 1/sample_hz)
wave_samples = np.abs(curr_qam) * np.cos(f_w * np.angle(curr_qam) * t_series)

plt.plot(t_series, wave_samples)
plt.show()

If I do this for every QAM symbol and concatenate the results, I think that would be exactly the signal transmitted over some analog medium (after passed through some DAC), but I wonder if there is a better/simpler or more elegant way of doing this with the CommPy library.

modulaiton.py PSKModem constellation not using numpy array

self.constellation = map(self._constellation_symbol,
self.symbol_mapping)

in the constructor, should be

self.constellation = array(map(self._constellation_symbol,
self.symbol_mapping))

similar to what QAMModem is doing.

Otherwise, this type error occurs.
TypeError: list indices must be integers, not list

Trying to run turbo_encode

import numpy as np
import commpy.channelcoding as cx
import commpy.channelcoding.convcode as cc
memory = np.array([2])
g_matrix = np.array([[5, 7]])
trellis = cc.Trellis(memory, g_matrix)
trellis.k =1 
trellis.n = 2
trellis.number_states=0
trellis.number_inputs = 0
trellis.total_memory = 0
x = np.array([1,0,1,0,1,1])
x=np.array(x);
print(cx.turbo_encode(x,trellis,trellis,cx.RandInterlv))

I get this error:

File "C:\Users\sankalp.1\AppData\Local\Programs\Python\Python37\lib\site-packages\commpy\utilities.py", line 50, in dec2bitarray
   bitarray[bit_width - i - 1] = int(binary_string[length - i - 1])
IndexError: index -1 is out of bounds for axis 0 with size 0

Can anyone help me figure out how to eliminate this error? It seems to be working fine in conv_encode.

Thank you.

Example plotConsModem.py isn't working...

FYI, I went ahead and tried your example.

Get this error.

'PSKModem' object has no attribute 'plot_constellation'

I'm not the greatest at programming so I could've missed something.

RootRaisedCosine impulse response seem incorrect

When comparing a RRC impulse response generated by CommPy to those of other reference implementations, the results seem incorrect -- e.g., not even symmetric.
I don't have screenshots at the moment, since I changed implementations, but may try to show evidence here.

Just a warning to those using it for pulse shaping.. it's not perfect.

interleaver python3 compatibility

Using an interleaver in a turbo encoder, the map operator was causing python3 compatibility problems. The following seems to fix it, though it may not be optimal:

    out_array = array(list(map(lambda x: in_array[x], self.p_array)))

size mismatch in channel coding README

size mismatch in channel coding README

ValueError                                Traceback (most recent call last)
<ipython-input-8-e1f7b72f48a6> in <module>
     36 
     37 
---> 38     NumErr, BER_soft[cntr] = BER_calc(message_bits, decoded_soft[:-(L-1)]) # bit-error ratio (soft decision)
     39     NumErr, BER_hard[cntr] = BER_calc(message_bits, decoded_hard[:-(L-1)]) # bit-error ratio (hard decision)
     40     NumErr, BER_uncoded[cntr] = BER_calc(message_bits, demodulated_uncoded) # bit-error ratio (uncoded case)

<ipython-input-1-a3bc5a317da7> in BER_calc(a, b)
      4 
      5 def BER_calc(a, b):
----> 6     num_ber = np.sum(np.abs(a - b))
      7     ber = np.mean(np.abs(a - b))
      8     return int(num_ber), ber

ValueError: operands could not be broadcast together with shapes (100000,) (99998,)

EbNo = 5 # energy per bit to noise power spectral density ratio (in dB)
snrdB = EbNo + 10*np.log10(k*rate) # Signal-to-Noise ratio (in dB)
noiseVar = 10**(-snrdB/10) # noise variance (power)

N_c = 10 # number of trials

BER_soft = np.empty((N_c,))
BER_hard = np.empty((N_c,))
BER_uncoded = np.empty((N_c,))

for cntr in range(N_c):
    
    message_bits = np.random.randint(0, 2, N) # message
    coded_bits = cc.conv_encode(message_bits, trellis) # encoding
    
    modulated = modem.modulate(coded_bits) # modulation
    modulated_uncoded = modem.modulate(message_bits) # modulation (uncoded case)

    Es = np.mean(np.abs(modulated)**2) # symbol energy
    No = Es/((10**(EbNo/10))*np.log2(M)) # noise spectrum density

    noisy = modulated + np.sqrt(No/2)*\
        (np.random.randn(modulated.shape[0])+\
         1j*np.random.randn(modulated.shape[0])) # AWGN
    
    noisy_uncoded = modulated_uncoded + np.sqrt(No/2)*\
        (np.random.randn(modulated_uncoded.shape[0])+\
         1j*np.random.randn(modulated_uncoded.shape[0])) # AWGN (uncoded case)

    demodulated_soft = modem.demodulate(noisy, demod_type='soft', noise_var=noiseVar) # demodulation (soft output)
    demodulated_hard = modem.demodulate(noisy, demod_type='hard') # demodulation (hard output)
    demodulated_uncoded = modem.demodulate(noisy_uncoded, demod_type='hard') # demodulation (uncoded case)

    decoded_soft = cc.viterbi_decode(demodulated_soft, trellis, tb_depth, decoding_type='unquantized') # decoding (soft decision)
    decoded_hard = cc.viterbi_decode(demodulated_hard, trellis, tb_depth, decoding_type='hard') # decoding (hard decision)


    NumErr, BER_soft[cntr] = BER_calc(message_bits, decoded_soft[:-(L-1)]) # bit-error ratio (soft decision)
    NumErr, BER_hard[cntr] = BER_calc(message_bits, decoded_hard[:-(L-1)]) # bit-error ratio (hard decision)
    NumErr, BER_uncoded[cntr] = BER_calc(message_bits, demodulated_uncoded) # bit-error ratio (uncoded case)

mean_BER_soft = np.mean(BER_soft) # averaged bit-error ratio (soft decision)
mean_BER_hard = np.mean(BER_hard) # averaged bit-error ratio (hard decision)
mean_BER_uncoded = np.mean(BER_uncoded) # averaged bit-error ratio (uncoded case)

print("Soft decision:\n"+str(mean_BER_soft)+"\n")
print("Hard decision:\n"+str(mean_BER_hard)+"\n")
print("Uncoded message:\n"+str(mean_BER_uncoded)+"\n")

Strange results in wifi 802.11

Hello,

I know I was the one who submitted the code in first place, and as far as I know, I implemented it according to the IEEE 802.11 standard.
According to the IEEE 802.11 section 17.3.5.6 the data is coded with a convolutional encoded of rate 1/2. Higher rates are achieved by then puncturing. In the receiver side, there is depuncturing by adding a dummy "zero" on the place of the omitted zeros.

This works good in some MCS (0,1,2,3,4) where the parameters have coding of 1/2 and 3/4 and modulation BPSK, QPSK and 16-QAM.
Graph with the BER vs SNR bellow (message size was 600bits, 2000 runs each SNR):

When higher MCS are tried (5, 6,7, 8, 9) the behaviour becomes strange as can be seen in a similar plot for those MCS bellow (same conditions of experiment):

• The first and more visible problem is with MCS 6, after a particular SNR (about 24) the BER stops decreasing at the same rate as others MCS (it still decreases but only slightly, maybe rounding differences).
• Second the MCS 7 achieves better BER than MCS 5 at most SNRs, notice both are 64-QAM, the only difference are the punturing performed. Similar behaviour happens with MCS 9 in relation with MCS 8.

Does someone have any idea why this simulation results are happening?
Some error in code or some error in my interpretation of the results?

If it helps debug, I notice some warnings in my console:

/commpy/commpy/modulation.py:131: RuntimeWarning: divide by zero encountered in double_scalars
demod_bits[i * self.num_bits_symbol + self.num_bits_symbol - 1 - bit_index] = log(llr_num / llr_den)

This is maybe related with the problem with MCS 6 at SNRs over 25, I can easily replicate this issue.

Also, I can add code if it helps. :-)

Failed to check for complex number in `channels.awgn`

Hi @veeresht ,

Thank you for creating this library, I am currently using cp.channels.awgn. However, it failed to check for complex ndarray.

Code to reproduce bug:

import commpy as cp
import numpy as np

np.random.seed(2019)
qpsk_modem = cp.modulation.QAMModem(m=4)

# Random message_bits length 5
msg_bits = np.random.randint(0, 2, size=5)

# A modulated signal, type = complex np.complex128
moded_signal = qpsk_modem.modulate(msg_bits)

# Corrupt signal
moded_singal = cp.channels.awgn(moded_signal, snr_dB=10)

• The awgn failed to add noise to imaginary part as you see ( it is still [.. -1j, ... +1j, ... -1j]

print(msg_bits)      # [0 0 1 1 0]
print(moded_singal)   # [-0.76902489-1.j  2.36032112+1.j -0.20336901-1.j]

Proposed Solution:

• Update type checking in channels.py

if type(input_signal[0]) == complex:

In order to fix it, I have to change it to:

if isinstace(input_singal[0], complex):

Output after the change

I have created a pull request. Feel free to check it out.

Thanks

Modulating with PSKModem throws a TypeError

I have been testing modulation with the commpy package. Apparently, modulating with the QAMModem works, whereas PSKModem (irrelevant of the M value) throws a TypeError. Find the sample code here

Coverage is not computed since commit 2b102b7

I noticed that the coverage reports are no more updated since commit 2b102b7. PR #50 shows the last coverage output. The section Recent builds on https://coveralls.io/github/veeresht/CommPy confirms that there are no more builds.

@veeresht any clues or fixes?

Puncturing convolutional code at rate 1/2 to get rate 3/4

Hi,
I am trying to execute the readme.md example for convolutional coding at code rate 3/4. The generator matrix [5, 7] for R=1/2 is used and punctured with [[1,0,1],[1,1,0]] to get R=3/4. But the BER values of viterbi decoded outputs are higher than that of uncoded output. Also the size of the coded bits is as per 1/2 rate and not 3/4. Is it an issue with the 'conv_encode' function or am I missing any steps in between?

Also is there any method to get the coding rate 3/4 using generator matrix?

import numpy as np
import commpy.channelcoding.convcode as cc
import commpy.modulation as modulation

def BER_calc(a, b):
  num_ber = np.sum(np.abs(a - b))
  ber = np.mean(np.abs(a - b))
  return int(num_ber), ber

N = 100 #number of symbols per the frame
message_bits = np.random.randint(0, 2, N) # message

M = 2 # modulation order (BPSK)
k = np.log2(M) #number of bit per modulation symbol
modem = modulation.PSKModem(M) # M-PSK modem initialization

generator_matrix = np.array([[5, 7]]) # generator branches
trellis = cc.Trellis(np.array([2]), generator_matrix) # Trellis structure
punctureMatrix=np.array([[1,0,1],[1,1,0]])

rate = 3/4 # code rate
L = 7 # constraint length
m = np.array([L-1]) # number of delay elements

tb_depth = 5*(m.sum() + 1) # traceback depth

EbNo = 5 # energy per bit to noise power spectral density ratio (in dB)
snrdB = EbNo + 10*np.log10(k*rate) # Signal-to-Noise ratio (in dB)
noiseVar = 10**(-snrdB/10) # noise variance (power)

N_c = 10 # number of trials

BER_soft = np.zeros(N_c)
BER_hard = np.zeros(N_c)
BER_uncoded = np.zeros(N_c)

for cntr in range(N_c):
  
  message_bits = np.random.randint(0, 2, N) # message
  coded_bits = cc.conv_encode(message_bits, trellis,puncture_matrix=punctureMatrix) # encoding
  modulated = modem.modulate(coded_bits) # modulation
  modulated_uncoded = modem.modulate(message_bits) # modulation (uncoded case)

  Es = np.mean(np.abs(modulated)**2) # symbol energy
  No = Es/((10**(EbNo/10))*np.log2(M)) # noise spectrum density

  noisy = modulated + np.sqrt(No/2)*\
      (np.random.randn(modulated.shape[0])+\
       1j*np.random.randn(modulated.shape[0])) # AWGN
  
  noisy_uncoded = modulated_uncoded + np.sqrt(No/2)*\
      (np.random.randn(modulated_uncoded.shape[0])+\
       1j*np.random.randn(modulated_uncoded.shape[0])) # AWGN (uncoded case)

  demodulated_soft = modem.demodulate(noisy, demod_type='soft', noise_var=noiseVar) # demodulation (soft output)
  demodulated_hard = modem.demodulate(noisy, demod_type='hard') # demodulation (hard output)
  demodulated_uncoded = modem.demodulate(noisy_uncoded, demod_type='hard') # demodulation (uncoded case)

  decoded_soft = cc.viterbi_decode(demodulated_soft, trellis, tb_depth, decoding_type='unquantized') # decoding (soft decision)
  decoded_hard = cc.viterbi_decode(demodulated_hard, trellis, tb_depth, decoding_type='hard') # decoding (hard decision)

  NumErr, BER_soft[cntr] = BER_calc(message_bits, decoded_soft[:message_bits.size]) # bit-error ratio (soft decision)
  NumErr, BER_hard[cntr] = BER_calc(message_bits, decoded_hard[:message_bits.size]) # bit-error ratio (hard decision)
  NumErr, BER_uncoded[cntr] = BER_calc(message_bits, demodulated_uncoded[:message_bits.size]) # bit-error ratio (uncoded case)

mean_BER_soft = BER_soft.mean() # averaged bit-error ratio (soft decision)
mean_BER_hard = BER_hard.mean() # averaged bit-error ratio (hard decision)
mean_BER_uncoded = BER_uncoded.mean() # averaged bit-error ratio (uncoded case)

print("Soft decision:\n{}\n".format(mean_BER_soft))
print("Hard decision:\n{}\n".format(mean_BER_hard))
print("Uncoded message:\n{}\n".format(mean_BER_uncoded))

PS: The package and libraries are from the Github cloned version.

[Potential Bug] Normalize real `input` in `awgn`.

Hi @veeresht again,

According to Wikipedia on SNR and dsp.stackexchange, it seems that we might need to normalize real inputs to be zero-mean as: inputs = 2.0 * inputs - 1.0, when adding White Gaussian Noise to signals.

I wonder if this is a potential bug. I ran two tests (decoding convolutional codes over AWGN Channel using Viterbi) to validate my speculation:

• First test: Original implementation of awgn
• Second test : Modified version as

 def awgn(input_signal, snr_dB, rate=1.0):
    avg_energy = sum(abs(input_signal) * abs(input_signal))/len(input_signal)
    snr_linear = 10**(snr_dB/10.0)
    noise_variance = avg_energy/(2*rate*snr_linear)
    if isinstance(input_signal[0], complex):
        noise = (sqrt(noise_variance) * randn(len(input_signal))) + (sqrt(noise_variance) * randn(len(input_signal))*1j)
    else:
        noise = sqrt(2*noise_variance) * randn(len(input_signal))
        input_signal = 2.0 * input_signal - 1.0  # Zero-mean normalization HERE

    output_signal = input_signal + noise
    return output_signal

Test file:

import multiprocessing as mp
import numpy as np
import commpy as cp
import matplotlib.pyplot as plt
import commpy.channelcoding.convcode as cc

from commpy.channelcoding import Trellis

# For reproducability
np.random.seed(2018)

NUM_EXAMPLES = 50
SNRs = np.arange(-1.0, 5.0, 1)
BLOCK_LEN = 100

M = np.array([2])
G =  np.array([[0o7, 0o5]])
TB_DEPTH = 15
trellis = Trellis(M, G, feedback=7)

def simulation_func(snr):
    # Generate random message bits to be encoded
    message_bits = np.random.randint(0, 2, BLOCK_LEN)

    # Encode message bits
    coded_bits = cc.conv_encode(message_bits, trellis)

    # Add Additive White Gaussian noise
    noisy_signal = cp.channels.awgn(coded_bits, snr, rate=1/2)
    
    # Estimate original message bit using Viterbit
    decoded_bits = cc.viterbi_decode(noisy_signal, trellis, TB_DEPTH, decoding_type='unquantized')     
    
    num_bit_errors = cp.utilities.hamming_dist(message_bits, decoded_bits[:-int(M)])
    return num_bit_errors


if __name__ == '__main__':
    BERs, BLERs =  [], []
    for snr in SNRs:
        # Test 100 examples
        with mp.Pool(2) as pool:
            hamming_dists = pool.map(simulation_func, [(snr) for _ in range(NUM_EXAMPLES)])
        # Save result
        BERs.append(np.sum(hamming_dists)/ (BLOCK_LEN * NUM_EXAMPLES))
        BLERs.append(np.count_nonzero(hamming_dists)/NUM_EXAMPLES)
        print('SNR = %d BER: %.2f BLER: %.2f'% (snr, BERs[-1], BLERs[-1]))

    # Plot the result
    fig, (left, right) = plt.subplots(1, 2, figsize=(10, 5))
    left.plot(SNRs, BERs)
    left.set_ylabel('Bit Error Rate (BER)')
    left.set_xlabel('Signal-to-Noise Ratio (SNR in dB)')
    left.semilogy()
    left.grid(True, which='both')
    right.plot(SNRs, BLERs)
    right.set_ylabel('Bit Block Error Rate (BLER)')
    right.set_xlabel('Signal-to-Noise Ratio (SNR in dB)')
    right.semilogy()
    right.grid(True, which='both')
    fig.tight_layout()
    plt.show()

BEFORE (First Test)

SNR = -1 BER: 0.36 BLER: 1.00
SNR = 0 BER: 0.34 BLER: 1.00
SNR = 1 BER: 0.33 BLER: 1.00
SNR = 2 BER: 0.31 BLER: 1.00
SNR = 3 BER: 0.30 BLER: 1.00
SNR = 4 BER: 0.28 BLER: 1.00

AFTER (Second Test)

SNR = -1 BER: 0.13 BLER: 1.00
SNR = 0 BER: 0.07 BLER: 0.96
SNR = 1 BER: 0.04 BLER: 0.66
SNR = 2 BER: 0.02 BLER: 0.44
SNR = 3 BER: 0.00 BLER: 0.18
SNR = 4 BER: 0.00 BLER: 0.12

I have not submitted my pull request yet because I am not sure. What you you think?

ImportError: cannot import name 'gcd' from 'fractions' (/usr/lib64/python3.9/fractions.py)

In Python 3.9 , "gcd" is deprecated : https://docs.python.org/3.5/library/fractions.html#fractions.gcd

Fix :

$ diff commpy/channelcoding/gfields.py*
6c6
< #from fractions import gcd
---
> from fractions import gcd
112c112
<             orders[idx] = (2**self.m - 1)/(math.gcd(i, 2**self.m-1))
---
>             orders[idx] = (2**self.m - 1)/(gcd(i, 2**self.m-1))

Requirement on pip

It would be usefull to include the dependencies in pip so that pip install scikit-commpy also installs the required packages.

Using the trellis

Hello,
This is not a issue with the code, more with trying to understand and using it.
I'm trying to simulate 802.11, and having difficulties with which parameters are the right to use to build the Trellis instance.
I have found a paper indicating the generator polynomials g0= 1011011 and g1= 1111001 [1], but to use it in this framework?
Like this g_matrix = np.array(((1, 0, 1, 1, 0, 1, 1), (1, 1, 1, 1, 0, 0, 1)), ndmin=2)?

Also, how to build the different coding rates of 802.11 (1/2, 3/4, 5/6)?

Thanks for the help.

[1] - (page 19) https://pdfs.semanticscholar.org/c63b/71e43dc23b17ca57267f3b769224c64d5e33.pdf

v0.5 release

Hi @veeresht,

I drafted a new release as the new content seems enough to justify it. Does it look good to you?
If so, I'll make the release and try to update pip. If I remember well, you haded me as a contributor on pip.

This version abandons Python 2 compatibility and all tests are now validated on Python 3.
The Conda version is no longer updated but pip should be updated accordingly as soon as possible.

• The convolutional code module now supports puncturing, all code ratios and provides a method for tracing FSM.
• LDPC encoder, design saving utilities, WiMAX designs and BP decoder enhancements.
• A tool to estimate the BER performance of a link model using a Monte Carlo estimator.
• K-best MIMO detector.
• Updates and corrections in the doc.
• Several bug fixes.

Add all contributors

@all-contributors please add @veeresht @BastienTr @kirlf @akou97 @eSoares @rtucker @prsudheer @NinjaComics @matchius @mborgerding @datlife

Adding AWGN to a signal

Adding Gaussian noise to a signal of a specific SNR. The signal may or may not have noise already

Import error prevent Travis to work

As pointed out in the recent PRs, there is a new import error in Travis. It deals with the conda environment setup: the line conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION atlas numpy scipy matplotlib nose used to work 2 months ago and now returns a PackagesNotFoundError for atlas.

confused about turbo.py

I use parameters in LTE 36.212 to realize a turbo code. Then I find in line 47 (turbo.py):

stream = conv_encode(msg_bits, trellis1, 'rsc')
sys_stream = stream[::2]
non_sys_stream_1 = stream[1::2]

interlv_msg_bits = interleaver.interlv(sys_stream)
puncture_matrix = array([[0, 1]])
non_sys_stream_2 = conv_encode(interlv_msg_bits, trellis2, 'rsc', puncture_matrix)

Here 'rsc' is not the right parameter in conv_encode, since the third parameter is termination(selected from {'term','cont'}). And 'rsc' should be set in Trellis like this:

memory = np.array([3])
g_matrix = np.array([[1,13]])
feedback = 11
trellis = cc.Trellis(memory, g_matrix, feedback, code_type='rsc')

So am I right? (Does there exist some wrong in turbo.py?)