Strange results in wifi 802.11

eSoares commented 4 years ago

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. :-)

BastienTr commented 4 years ago

Hi Eduardo,

I agree that there does appear to be a glitch somewhere. I don't feel competent enough on WiFi to spot a bug on that area. I can try to check if it is coming from the remainder (simulation, coding, modulation, etc).

Could you share the file that produces the faulty graphs? Besides, do you have any reference that gives the expected curves so that we can tell if it works as expected or not?

eSoares commented 4 years ago

Hi BastienTr, Sorry for the time to reply, I was reproducing the problem in a easier way to show and debug. The issue is the strange behaviour with the parameters: modem 64 QAM, and performing puncturing 3/4 ([1, 1, 1, 0, 0, 1]). Different modems (16 QAM) or different puncturing levels (i.e. 5/6 [1, 1, 1, 0, 0, 1, 1, 0, 0, 1]) don't manifest the issue.

To illustrate, I show where the behaviour of 16QAM with different coding rates (1/2 was no puncturing, the other coding rates are achieved with first coding at 1/2 and then puncturing):

The expected behaviour is starting with a BER of about 0.5 and with the increase of SNR, it will decrease until zero. The higher coding rate, the higher the BER for a same SNR.

The same test (adjusted for higher SNRs) with 64QAM the behaviour is the following:

With coding rates of 1/2 and 5/6 the behaviour is the expected. Coding 3/4 manifests the strange behaviour.

The code to test the results is bellow (used multiprocessing just to speed up, but behaviour also occurs if no multiprocess is used):

# Authors: CommPy contributors
# License: BSD 3-Clause

import math
import time
from datetime import timedelta
from itertools import product

import matplotlib.pyplot as plt
import numpy as np
from p_tqdm import p_map

import commpy.channelcoding.convcode as cc
import commpy.channels as chan
import commpy.links as lk
import commpy.modulation as mod
import sys

# =============================================================================
# Convolutional Code 1: G(D) = [1+D^2, 1+D+D^2]
# Standard code with rate 1/2
# =============================================================================

# Number of delay elements in the convolutional encoder
memory = np.array(7, ndmin=1)

# Generator matrix
g_matrix = np.array((133, 171), ndmin=2)

# Create trellis data structure
trellis1 = cc.Trellis(memory, g_matrix)

is16 = False
if len(sys.argv) > 1 and sys.argv[1] == "16":
    is16 = True

modem = mod.QAMModem(64)
if is16:
    modem = mod.QAMModem(16)
# modem = mod.PSKModem(2)

# AWGN channel
channels = chan.SISOFlatChannel(None, (1 + 0j, 0j))

puncture_matrix34 = [1, 1, 1, 0, 0, 1]
puncture_matrix56 = [1, 1, 1, 0, 0, 1, 1, 0, 0, 1]
code_rate34 = 3 / 4
code_rate56 = 5 / 6

# Modulation function
def _modulate(bits):
    bits2 = cc.conv_encode(bits, trellis1, 'cont')
    return modem.modulate(bits2)

# Modulation function
def modulate_punturing34(bits):
    # bits2 = cc.conv_encode(bits, trellis1, 'cont', puncture_matrix=puncture_matrix)
    bits2 = cc.conv_encode(bits, trellis1, 'cont')
    bits2 = cc.puncturing(bits2, puncture_matrix34)
    return modem.modulate(bits2)

# Modulation function
def modulate_punturing56(bits):
    # bits2 = cc.conv_encode(bits, trellis1, 'cont', puncture_matrix=puncture_matrix)
    bits2 = cc.conv_encode(bits, trellis1, 'cont')
    bits2 = cc.puncturing(bits2, puncture_matrix56)
    return modem.modulate(bits2)

# Receiver function (no process required as there are no fading)
def receiver(y, h, constellation, noise_var):
    return modem.demodulate(y, 'soft', noise_var)

# Decoder function
def _decoder(msg):
    return cc.viterbi_decode(msg, trellis1, decoding_type='soft')

# Decoder function
def decoder_puncturing34(msg):
    msg2 = cc.depuncturing(msg, puncture_matrix34, math.ceil(len(msg) * code_rate34 * 2))
    return cc.viterbi_decode(msg2, trellis1, decoding_type='soft')

# Decoder function
def decoder_puncturing56(msg):
    msg2 = cc.depuncturing(msg, puncture_matrix56, math.ceil(len(msg) * code_rate56 * 2))
    return cc.viterbi_decode(msg2, trellis1, decoding_type='soft')

def perf(model_):
    if model_[1] == 0:
        code_rate = code_rate34
        modulate = modulate_punturing34
        decoder = decoder_puncturing34
    elif model_[1] == 1:
        code_rate = code_rate56
        modulate = modulate_punturing56
        decoder = decoder_puncturing56
    else:
        code_rate = 1 / 2
        modulate = _modulate
        decoder = _decoder

    model = lk.LinkModel(modulate, channels, receiver,
                         modem.num_bits_symbol, modem.constellation, modem.Es,
                         decoder, code_rate)

    SNRs = np.arange(model_[0], model_[0] + 1) + 10 * math.log10(modem.num_bits_symbol)
    return {"code_rate": model_[1], "SNR": SNRs[0],
            "res": model.link_performance_full_metrics(SNRs, 10, 0, 600, code_rate, stop_on_surpass_error=False)}

models = list(product(np.arange(5, 20), [0, 1, 2]))
if is16:
    models = list(product(np.arange(0, 15), [0, 1, 2]))

start = time.time()
data = p_map(perf, models)
print(str(timedelta(seconds=(time.time() - start))))

print("16 QAM" if is16 else "64QAM")
print("SNRs", list(map(lambda x: x[0] + 10 * math.log10(modem.num_bits_symbol), filter(lambda x: x[1] == 0, models))))
rate12 = list(map(lambda x: x["res"][0][0], sorted(filter(lambda x: x["code_rate"] == 2, data), key=lambda x: x["SNR"])))
rate34 = list(map(lambda x: x["res"][0][0], sorted(filter(lambda x: x["code_rate"] == 0, data), key=lambda x: x["SNR"])))
rate56 = list(map(lambda x: x["res"][0][0], sorted(filter(lambda x: x["code_rate"] == 1, data), key=lambda x: x["SNR"])))
print("1/2", rate12)
print("3/4", rate34)
print("5/6", rate56)

# plt.semilogy(SNRs, BERs_no_punturing[0], 'o-', SNRs, BERs_puncturing[0], 'o-')
# plt.grid()
# plt.xlabel('Signal to Noise Ration (dB)')
# plt.ylabel('Bit Error Rate')
# plt.legend(('W/o Punturing', 'W/punturing'))
# plt.show()

BastienTr commented 4 years ago

Thanks for the minimal example Eduardo. I can't have a look this week. I'm overwhelmed by some PhD staff... I'll try to look closer in the next weeks.

kirlf commented 4 years ago

Hi, guys!

Frankly speaking, I am not so confident, but it seems like issue of convolutional codes implementation. I saw something similar when use MatLab: https://www.mathworks.com/matlabcentral/fileexchange/72392-exact-vs-approximate-llr-calculation?s_tid=prof_contriblnk

Hopefully, I have not increased discussion entropy now! :)

eSoares commented 4 years ago

Hi kirlf, every guess is a good as any! It was in some way related with the convolutional codes, but not an issue with the implementation (as far as I know).

I found the issue and I'm making a PR as soon as travis passes the tests (they pass in my machine, so it should be soon).

The issue was with the send chunk after coding and puncturing not being divisible by the number of bits of the modulation! There was code that tried to select the correct send chunk in order to avoid the issue, but was only taking in to account the modulation, without looking at the code rate! (Creating issues even if the send chunk chosen was perfect for after coding rate being divisible by the modulation!)

The fix is slightly adjusting the code to select a good send chunk. Making sure the chunk is both divisible by the modulation after applying the code rate, and divisible by the code rate after demodulating at the receiver side.

I have written the following test to check it actually fixed it (click to expand):

```python # Authors: CommPy contributors # License: BSD 3-Clause import math import time from datetime import timedelta from fractions import Fraction from itertools import product import matplotlib.pyplot as plt import numpy as np from p_tqdm import p_map import commpy.channelcoding.convcode as cc import commpy.channels as chan import commpy.links as lk import commpy.modulation as mod import sys # ============================================================================= # Convolutional Code 1: G(D) = [1+D^2, 1+D+D^2] # Standard code with rate 1/2 # ============================================================================= # Number of delay elements in the convolutional encoder memory = np.array(7, ndmin=1) # Generator matrix g_matrix = np.array((133, 171), ndmin=2) # Create trellis data structure trellis1 = cc.Trellis(memory, g_matrix) is16 = False if len(sys.argv) > 1 and sys.argv[1] == "16": is16 = True modem = mod.QAMModem(64) if is16: modem = mod.QAMModem(16) # modem = mod.PSKModem(2) # AWGN channel channels = chan.SISOFlatChannel(None, (1 + 0j, 0j)) puncture_matrix34 = [1, 1, 1, 0, 0, 1] puncture_matrix56 = [1, 1, 1, 0, 0, 1, 1, 0, 0, 1] code_rate34 = Fraction(3, 4) code_rate56 = Fraction(5, 6) # Modulation function def _modulate(bits): bits2 = cc.conv_encode(bits, trellis1, 'cont') return modem.modulate(bits2) # Modulation function def modulate_punturing34(bits): # bits2 = cc.conv_encode(bits, trellis1, 'cont', puncture_matrix=puncture_matrix) bits2 = cc.conv_encode(bits, trellis1, 'cont') bits2 = cc.puncturing(bits2, puncture_matrix34) return modem.modulate(bits2) # Modulation function def modulate_punturing56(bits): # bits2 = cc.conv_encode(bits, trellis1, 'cont', puncture_matrix=puncture_matrix) bits2 = cc.conv_encode(bits, trellis1, 'cont') bits2 = cc.puncturing(bits2, puncture_matrix56) return modem.modulate(bits2) # Receiver function (no process required as there are no fading) def receiver(y, h, constellation, noise_var): return modem.demodulate(y, 'soft', noise_var) # Decoder function def _decoder(msg): return cc.viterbi_decode(msg, trellis1, decoding_type='soft') # Decoder function def decoder_puncturing34(msg): msg2 = cc.depuncturing(msg, puncture_matrix34, math.ceil(len(msg) * code_rate34 * 2)) return cc.viterbi_decode(msg2, trellis1, decoding_type='soft') # Decoder function def decoder_puncturing56(msg): msg2 = cc.depuncturing(msg, puncture_matrix56, math.ceil(len(msg) * code_rate56 * 2)) return cc.viterbi_decode(msg2, trellis1, decoding_type='soft') def perf(model_): if model_[1] == 0: code_rate = code_rate34 modulate = modulate_punturing34 decoder = decoder_puncturing34 elif model_[1] == 1: code_rate = code_rate56 modulate = modulate_punturing56 decoder = decoder_puncturing56 else: code_rate = Fraction(1, 2) modulate = _modulate decoder = _decoder model = lk.LinkModel(modulate, channels, receiver, modem.num_bits_symbol, modem.constellation, modem.Es, decoder, code_rate) SNRs = np.arange(model_[0], model_[0] + 1) + 10 * math.log10(modem.num_bits_symbol) nbits = model_[2] return {"code_rate": model_[1], "SNR": SNRs[0], "size": nbits, "res": model.link_performance_full_metrics(SNRs, 10, 0, nbits, code_rate, stop_on_surpass_error=False)} SNR = 30 if is16: SNR = 25 psizes = list(range(10, 1500, 10)) # psizes = [10] models = list(product([SNR], [0, 1, 2], psizes)) start = time.time() data = p_map(perf, models) # data = [perf([20, 0, 10])] print(str(timedelta(seconds=(time.time() - start)))) for psize in psizes: rate12 = list(map(lambda x: x["res"][0][0], sorted(filter(lambda x: x["code_rate"] == 2 and x["size"] == psize, data), key=lambda x: x["SNR"]))) if rate12[0] > 0: print("FAILED, rate 1/2, size %d, BER %f" % (psize, rate12[0])) rate34 = list(map(lambda x: x["res"][0][0], sorted(filter(lambda x: x["code_rate"] == 0 and x["size"] == psize, data), key=lambda x: x["SNR"]))) if rate34[0] > 0: print("FAILED, rate 3/4, size %d, BER %f" % (psize, rate34[0])) rate56 = list(map(lambda x: x["res"][0][0], sorted(filter(lambda x: x["code_rate"] == 1 and x["size"] == psize, data), key=lambda x: x["SNR"]))) if rate56[0] > 0: print("FAILED, rate 5/6, size %d, BER %f" % (psize, rate56[0])) ```

I didn't add to the suite of tests, because it takes a huge amount of time.

BastienTr commented 4 years ago

it takes a huge amount of time.

I'm running some tests and... well you were right :sweat_smile:. I keep in touch when computations and tests end.

eSoares commented 4 years ago

it takes a huge amount of time.

I'm running some tests and... well you were right sweat_smile. I keep in touch when computations and tests end.

If it helps, in my 48 core running machine was something like 20 minutes! (all cores used)

BastienTr commented 4 years ago

Then I must trust you, this will be too much for my 4-cores laptop.

eSoares commented 4 years ago

You can run a smaller set, instead of calculating for every send chunk between 10, 20, ...1500, just select a few!

veeresht / CommPy

Strange results in wifi 802.11 #85