diff --git a/.travis.yml b/.travis.yml
index 4afe9c3..a989699 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -23,7 +23,7 @@ install:
     - conda info -a
 
     # Replace dep1 dep2 ... with your dependencies
-    - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib nose
+    - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib nose sympy
     - source activate test-environment
     #- conda install --yes -c dan_blanchard python-coveralls nose-cov
     - pip install coveralls
diff --git a/README.md b/README.md
index 1c34b76..e3b47a3 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
 
 
 [![Build Status](https://secure.travis-ci.org/veeresht/CommPy.svg?branch=master)](https://secure.travis-ci.org/veeresht/CommPy)
-[![Coverage](https://coveralls.io/repos/veeresht/CommPy/badge.svg)](https://coveralls.io/r/veeresht/CommPy)
+[![Coverage](https://coveralls.io/repos/veeresht/CommPy/badge.svg?branch=master)](https://coveralls.io/r/veeresht/CommPy)
 [![PyPi](https://badge.fury.io/py/scikit-commpy.svg)](https://badge.fury.io/py/scikit-commpy)
 [![Docs](https://readthedocs.org/projects/commpy/badge/?version=latest)](http://commpy.readthedocs.io/en/latest/?badge=latest)
 
@@ -37,6 +37,9 @@ Available Features
 - Binary Symmetric Channel (BSC)
 - Binary AWGN Channel (BAWGNC)
 
+[Wifi 802.11 Simulation Class](https://github.com/veeresht/CommPy/blob/master/commpy/wifi80211.py)
+- A class to simulate the transmissions and receiving parameters of physical layer 802.11 (currently till VHT (ac)).
+
 [Filters](https://github.com/veeresht/CommPy/blob/master/commpy/filters.py)
 -------
 - Rectangular
@@ -99,6 +102,7 @@ Requirements/Dependencies
 - scipy 0.15 or above
 - matplotlib 1.4 or above
 - nose 1.3 or above
+- sympy 1.7 or above
 
 Installation
 ------------
@@ -122,9 +126,9 @@ 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
+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: veeresht@gmail.com
 
-For more details on CommPy, please visit http://veeresht.github.com/CommPy
+For more details on CommPy, please visit https://veeresht.info/CommPy/
diff --git a/THANKS.txt b/THANKS.txt
index c9ec33d..68c5492 100644
--- a/THANKS.txt
+++ b/THANKS.txt
@@ -2,6 +2,7 @@ Please add names as needed so that we can keep up with all the contributors.
 
 Veeresh Taranalli for initial creation and contribution of CommPy.
 Bastien Trotobas for adding some features, bugfixes, maintenance.
+Eduardo Soares (@eSoares) for several enhancements including WiFi 802.11 class and code speedups.
 Vladimir Fadeev for bugfixes, addition of convolutional code puncturing and readme explanation of codings.
 Youness Akourim for adding features and fixing some bugs.
 Rey Tucker for python3 compatibility fixes.
diff --git a/commpy/channelcoding/README.md b/commpy/channelcoding/README.md
index 57eb338..14fd9f4 100644
--- a/commpy/channelcoding/README.md
+++ b/commpy/channelcoding/README.md
@@ -78,11 +78,12 @@ The [following](https://en.wikipedia.org/wiki/File:Conv_code_177_133.png) convol
 Convolutional encoder parameters:
 
 ```python
+generator_matrix = np.array([[5, 7]]) # generator branches
+trellis = cc.Trellis(np.array([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
 ```
 
 Viterbi decoder parameters:
@@ -106,9 +107,9 @@ 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,))
+BER_soft = np.zeros(N_c)
+BER_hard = np.zeros(N_c)
+BER_uncoded = np.zeros(N_c)
 
 for cntr in range(N_c):
     
@@ -136,31 +137,30 @@ for cntr in range(N_c):
     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)
 
-    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)
+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"+str(mean_BER_soft)+"\n")
-print("Hard decision:\n"+str(mean_BER_hard)+"\n")
-print("Uncoded message:\n"+str(mean_BER_uncoded)+"\n")
+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))
 ```
 
 Outputs:
 
 ```python
->>> Soft decision:
->>> 0.0
->>>
->>> Hard decision:
->>> 3.0000000000000004e-05
->>>
->>> Uncoded message:
->>> 0.008782
+Soft decision:
+2.8000000000000003e-05
+
+Hard decision:
+0.0007809999999999999
+
+Uncoded message:
+0.009064
 ```
 
 ### Reference
diff --git a/commpy/channelcoding/convcode.py b/commpy/channelcoding/convcode.py
index dc530f7..d1dcf16 100644
--- a/commpy/channelcoding/convcode.py
+++ b/commpy/channelcoding/convcode.py
@@ -695,12 +695,13 @@ def viterbi_decode(coded_bits, trellis, tb_depth=None, decoding_type='hard'):
     rate = k/n
     total_memory = trellis.total_memory
 
-    if tb_depth is None:
-        tb_depth = 5*total_memory
-
     # Number of message bits after decoding
     L = int(len(coded_bits)*rate)
 
+    if tb_depth is None:
+        tb_depth = min(5 * total_memory, L)
+
+
     path_metrics = np.full((trellis.number_states, 2), np.inf)
     path_metrics[0][0] = 0
     paths = np.empty((trellis.number_states, tb_depth), 'int')
@@ -747,7 +748,8 @@ def viterbi_decode(coded_bits, trellis, tb_depth=None, decoding_type='hard'):
 
     return decoded_bits[:L]
 
-def puncturing(message, punct_vec):
+
+def puncturing(message: np.ndarray, punct_vec: np.ndarray) -> np.ndarray:
     """
     Applying of the punctured procedure.
     Parameters
@@ -771,7 +773,8 @@ def puncturing(message, punct_vec):
             shift = shift + 1
     return np.array(punctured)
 
-def depuncturing(punctured, punct_vec, shouldbe):
+
+def depuncturing(punctured: np.ndarray, punct_vec: np.ndarray, shouldbe: int) -> np.ndarray:
     """
     Applying of the inserting zeros procedure.
     Parameters
@@ -797,5 +800,5 @@ def depuncturing(punctured, punct_vec, shouldbe):
         else:
             shift2 = shift2 + 1
         if idx%N == 0:
-            shift = shift + 1;
+            shift = shift + 1
     return depunctured
diff --git a/commpy/channelcoding/gfields.py b/commpy/channelcoding/gfields.py
index 49cc680..6bd0e1d 100644
--- a/commpy/channelcoding/gfields.py
+++ b/commpy/channelcoding/gfields.py
@@ -3,7 +3,7 @@
 
 """ Galois Fields """
 
-from fractions import gcd
+from math import gcd
 
 from numpy import array, zeros, arange, convolve, ndarray, concatenate
 
@@ -52,7 +52,7 @@ def __init__(self, x, m):
         if type(x) is int and x >= 0 and x < pow(2, m):
             self.elements = array([x])
         elif type(x) is ndarray and len(x) >= 1:
-            self.elements = x
+            self.elements = x.astype(int)
 
     # Overloading addition operator for Galois Field
     def __add__(self, x):
diff --git a/commpy/channelcoding/ldpc.py b/commpy/channelcoding/ldpc.py
index f843f68..b8fe73c 100644
--- a/commpy/channelcoding/ldpc.py
+++ b/commpy/channelcoding/ldpc.py
@@ -241,11 +241,11 @@ def ldpc_bp_decode(llr_vec, ldpc_code_params, decoder_algorithm, n_iters):
 
             # Variable Node Update
             msg_sum = np.array(message_matrix.sum(0)).squeeze()
-            message_matrix *= -1
-            message_matrix += parity_check_matrix.multiply(msg_sum + llr_vec[i_start:i_stop])
+            message_matrix.data *= -1
+            message_matrix.data += parity_check_matrix.multiply(msg_sum + llr_vec[i_start:i_stop]).data
 
-            out_llrs = msg_sum + llr_vec[i_start:i_stop]
-            np.signbit(out_llrs, out=dec_word[i_start:i_stop])
+            out_llrs[i_start:i_stop] = msg_sum + llr_vec[i_start:i_stop]
+            np.signbit(out_llrs[i_start:i_stop], out=dec_word[i_start:i_stop])
 
     # Reformat outputs
     n_blocks = llr_vec.size // ldpc_code_params['n_vnodes']
@@ -304,7 +304,7 @@ def triang_ldpc_systematic_encode(message_bits, ldpc_code_params, pad=True):
     Encode bits using the LDPC code specified. If the  generator matrix is not computed, this function will build it
     and add it to the dictionary. It will also add the parity check matrix.
 
-    This function work only for LDPC specified by a triangular parity check matrix.
+    This function work only for LDPC specified by a approximate triangular parity check matrix.
 
     Parameters
     ----------
@@ -312,7 +312,8 @@ def triang_ldpc_systematic_encode(message_bits, ldpc_code_params, pad=True):
         Message bit to encode.
 
     ldpc_code_params : dictionary that at least contains one of these options:
-        Option 1: generator matrix is available.
+        Option 1: generator matrix and parity-check matrix are available.
+                parity_check_matrix (CSC sparse matrix of int8) - parity check matrix.
                 generator_matrix (2D-array or sparse matrix) - generator matrix of the code.
         Option 2: generator and parity check matrices will be added as sparse matrices.
                 n_vnodes (int) - number of variable nodes.
@@ -337,7 +338,7 @@ def triang_ldpc_systematic_encode(message_bits, ldpc_code_params, pad=True):
             If the message length is not a multiple of block length and pad is False.
     """
 
-    if ldpc_code_params.get('generator_matrix') is None:
+    if ldpc_code_params.get('generator_matrix') is None or ldpc_code_params.get('parity_check_matrix') is None:
         build_matrix(ldpc_code_params)
 
     block_length = ldpc_code_params['generator_matrix'].shape[1]
diff --git a/commpy/channels.py b/commpy/channels.py
index 37866af..d374755 100644
--- a/commpy/channels.py
+++ b/commpy/channels.py
@@ -19,7 +19,7 @@
 
 from __future__ import division, print_function  # Python 2 compatibility
 
-from numpy import complex, abs, sqrt, sum, zeros, identity, hstack, einsum, trace, kron, absolute, fromiter, array, exp, \
+from numpy import abs, sqrt, sum, zeros, identity, hstack, einsum, trace, kron, absolute, fromiter, array, exp, \
     pi, cos
 from numpy.random import randn, random, standard_normal
 from scipy.linalg import sqrtm
@@ -54,7 +54,7 @@ def generate_noises(self, dims):
         else:
             self.noises = standard_normal(dims) * self.noise_std
 
-    def set_SNR_dB(self, SNR_dB, code_rate=1, Es=1):
+    def set_SNR_dB(self, SNR_dB, code_rate: float = 1., Es=1):
 
         """
         Sets the the noise standard deviation based on SNR expressed in dB.
@@ -449,7 +449,7 @@ def specular_compo(self, thetat, dt, thetar, dr):
         H = zeros((self.nb_rx, self.nb_tx), dtype=complex)
         for n in range(self.nb_rx):
             for m in range(self.nb_tx):
-                H[n, m] = exp(1j * 2 * pi * (n * dr * cos(thetar) - m * dt * cos(thetat)))
+                H[n, m] = exp(1j * 2 * pi * (n * dr * cos(thetar) + m * dt * cos(thetat)))
         return H
 
     @property
diff --git a/commpy/examples/wifi80211_conv_encode_decode.py b/commpy/examples/wifi80211_conv_encode_decode.py
new file mode 100644
index 0000000..5806cc9
--- /dev/null
+++ b/commpy/examples/wifi80211_conv_encode_decode.py
@@ -0,0 +1,34 @@
+# Authors: CommPy contributors
+# License: BSD 3-Clause
+
+import math
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import commpy.channels as chan
+# ==================================================================================================
+# Complete example using Commpy Wifi 802.11 physical parameters
+# ==================================================================================================
+from commpy.wifi80211 import Wifi80211
+
+# AWGN channel
+channels = chan.SISOFlatChannel(None, (1 + 0j, 0j))
+
+w2 = Wifi80211(mcs=2)
+w3 = Wifi80211(mcs=3)
+
+# SNR range to test
+SNRs2 = np.arange(0, 6) + 10 * math.log10(w2.get_modem().num_bits_symbol)
+SNRs3 = np.arange(0, 6) + 10 * math.log10(w3.get_modem().num_bits_symbol)
+
+BERs_mcs2 = w2.link_performance(channels, SNRs2, 10, 10, 600, stop_on_surpass_error=False)
+BERs_mcs3 = w3.link_performance(channels, SNRs3, 10, 10, 600, stop_on_surpass_error=False)
+
+# Test
+plt.semilogy(SNRs2, BERs_mcs2, 'o-', SNRs3, BERs_mcs3, 'o-')
+plt.grid()
+plt.xlabel('Signal to Noise Ration (dB)')
+plt.ylabel('Bit Error Rate')
+plt.legend(('MCS 2', 'MCS 3'))
+plt.show()
diff --git a/commpy/links.py b/commpy/links.py
index a124956..356f32d 100644
--- a/commpy/links.py
+++ b/commpy/links.py
@@ -16,6 +16,7 @@
 from __future__ import division  # Python 2 compatibility
 
 import math
+from fractions import Fraction
 from inspect import getfullargspec
 
 import numpy as np
@@ -49,7 +50,7 @@ def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_
                   so it should be large enough regarding the code type.
                   *Default*: send_chunck = err_min
 
-    code_rate : float in (0,1]
+    code_rate : float or Fraction in (0,1]
                 Rate of the used code.
                 *Default*: 1 i.e. no code.
 
@@ -58,7 +59,8 @@ def link_performance(link_model, SNRs, send_max, err_min, send_chunk=None, code_
     BERs : 1d ndarray
            Estimated Bit Error Ratio corresponding to each SNRs
     """
-
+    if not send_chunk:
+        send_chunk = err_min
     return link_model.link_performance(SNRs, send_max, err_min, send_chunk, code_rate)
 
 
@@ -93,12 +95,12 @@ class LinkModel:
          *Default* Es=1.
 
     decoder : function with prototype decoder(array) or decoder(y, H, constellation, noise_var, array) that return a
-                binary array.
+                binary ndarray.
               *Default* is no process.
 
-    rate : float
-        Code rate.
-        *Default* is 1.
+    rate : float or Fraction in (0,1]
+           Rate of the used code.
+           *Default*: 1 i.e. no code.
 
     Attributes
     ----------
@@ -125,7 +127,7 @@ class LinkModel:
     Es : float
          Average energy per symbols.
 
-    decoder : function with prototype decoder(binary array) that return a binary array.
+    decoder : function with prototype decoder(binary array) that return a binary ndarray.
               *Default* is no process.
 
     rate : float
@@ -133,21 +135,25 @@ class LinkModel:
         *Default* is 1.
     """
 
-    def __init__(self, modulate, channel, receive, num_bits_symbol, constellation, Es=1, decoder=None, rate=1):
+    def __init__(self, modulate, channel, receive, num_bits_symbol, constellation, Es=1, decoder=None, rate=Fraction(1, 1)):
         self.modulate = modulate
         self.channel = channel
         self.receive = receive
         self.num_bits_symbol = num_bits_symbol
         self.constellation = constellation
         self.Es = Es
+        if type(rate) is float:
+            rate = Fraction(rate).limit_denominator(100)
         self.rate = rate
 
         if decoder is None:
             self.decoder = lambda msg: msg
         else:
             self.decoder = decoder
+        self.full_simulation_results = None
 
-    def link_performance(self, SNRs, send_max, err_min, send_chunk=None, code_rate=1):
+    def link_performance_full_metrics(self, SNRs, tx_max, err_min, send_chunk=None, code_rate: Fraction = Fraction(1, 1),
+                                      number_chunks_per_send=1, stop_on_surpass_error=True):
         """
         Estimate the BER performance of a link model with Monte Carlo simulation.
 
@@ -157,8 +163,8 @@ def link_performance(self, SNRs, send_max, err_min, send_chunk=None, code_rate=1
                Signal to Noise ratio in dB defined as :math:`SNR_{dB} = (E_b/N_0)_{dB} + 10 \log_{10}(R_cM_c)`
                where :math:`Rc` is the code rate and :math:`Mc` the modulation rate.
 
-        send_max : int
-                   Maximum number of bits send for each SNR.
+        tx_max : int
+                 Maximum number of transmissions for each SNR.
 
         err_min : int
                   link_performance send bits until it reach err_min errors (see also send_max).
@@ -168,10 +174,117 @@ def link_performance(self, SNRs, send_max, err_min, send_chunk=None, code_rate=1
                       so it should be large enough regarding the code type.
                       *Default*: send_chunck = err_min
 
-        code_rate : float in (0,1]
+        code_rate : Fraction in (0,1]
                     Rate of the used code.
                     *Default*: 1 i.e. no code.
 
+        number_chunks_per_send : int
+                                 Number of chunks per transmission
+
+        stop_on_surpass_error : bool
+                                Controls if during simulation of a SNR it should break and move to the next SNR when
+                                the bit error is above the err_min parameter
+
+        Returns
+        -------
+        List[BERs, BEs, CEs, NCs]
+           BERs : 1d ndarray
+                  Estimated Bit Error Ratio corresponding to each SNRs
+           BEs : 2d ndarray
+                 Number of Estimated Bits with Error per transmission corresponding to each SNRs
+           CEs : 2d ndarray
+                 Number of Estimated Chunks with Errors per transmission corresponding to each SNRs
+           NCs : 2d ndarray
+                 Number of Chunks transmitted per transmission corresponding to each SNRs
+        """
+
+        # Initialization
+        BERs = np.zeros_like(SNRs, dtype=float)
+        BEs = np.zeros((len(SNRs), tx_max), dtype=int)  # Bit errors per tx
+        CEs = np.zeros((len(SNRs), tx_max), dtype=int)  # Chunk Errors per tx
+        NCs = np.zeros((len(SNRs), tx_max), dtype=int)  # Number of Chunks per tx
+        # Set chunk size and round it to be a multiple of num_bits_symbol* nb_tx to avoid padding taking in to account the coding rate
+        if send_chunk is None:
+            send_chunk = err_min
+        if type(code_rate) is float:
+            code_rate = Fraction(code_rate).limit_denominator(100)
+        self.rate = code_rate
+        divider = (Fraction(1, self.num_bits_symbol * self.channel.nb_tx) * 1 / code_rate).denominator
+        send_chunk = max(divider, send_chunk // divider * divider)
+
+        receive_size = self.channel.nb_tx * self.num_bits_symbol
+        full_args_decoder = len(getfullargspec(self.decoder).args) > 1
+
+        # Computations
+        for id_SNR in range(len(SNRs)):
+            self.channel.set_SNR_dB(SNRs[id_SNR], float(code_rate), self.Es)
+            total_tx_send = 0
+            bit_err = np.zeros(tx_max, dtype=int)
+            chunk_loss = np.zeros(tx_max, dtype=int)
+            chunk_count = np.zeros(tx_max, dtype=int)
+            for id_tx in range(tx_max):
+                if stop_on_surpass_error and bit_err.sum() > err_min:
+                    break
+                # Propagate some bits
+                msg = np.random.choice((0, 1), send_chunk * number_chunks_per_send)
+                symbs = self.modulate(msg)
+                channel_output = self.channel.propagate(symbs)
+
+                # Deals with MIMO channel
+                if isinstance(self.channel, MIMOFlatChannel):
+                    nb_symb_vector = len(channel_output)
+                    received_msg = np.empty(int(math.ceil(len(msg) / float(self.rate))))
+                    for i in range(nb_symb_vector):
+                        received_msg[receive_size * i:receive_size * (i + 1)] = \
+                            self.receive(channel_output[i], self.channel.channel_gains[i],
+                                         self.constellation, self.channel.noise_std ** 2)
+                else:
+                    received_msg = self.receive(channel_output, self.channel.channel_gains,
+                                                self.constellation, self.channel.noise_std ** 2)
+                # Count errors
+                if full_args_decoder:
+                    decoded_bits = self.decoder(channel_output, self.channel.channel_gains,
+                                                self.constellation, self.channel.noise_std ** 2,
+                                                received_msg, self.channel.nb_tx * self.num_bits_symbol)
+                else:
+                    decoded_bits = self.decoder(received_msg)
+                # calculate number of error frames
+                for i in range(number_chunks_per_send):
+                    errors = np.bitwise_xor(msg[send_chunk * i:send_chunk * (i + 1)],
+                                            decoded_bits[send_chunk * i:send_chunk * (i + 1)].astype(int)).sum()
+                    bit_err[id_tx] += errors
+                    chunk_loss[id_tx] += 1 if errors > 0 else 0
+
+                chunk_count[id_tx] += number_chunks_per_send
+                total_tx_send += 1
+            BERs[id_SNR] = bit_err.sum() / (total_tx_send * send_chunk)
+            BEs[id_SNR] = bit_err
+            CEs[id_SNR] = np.where(bit_err > 0, 1, 0)
+            NCs[id_SNR] = chunk_count
+            if BEs[id_SNR].sum() < err_min:
+                break
+        self.full_simulation_results = BERs, BEs, CEs, NCs
+        return BERs, BEs, CEs, NCs
+
+    def link_performance(self, SNRs, send_max, err_min, send_chunk=None, code_rate=1):
+        """
+        Estimate the BER performance of a link model with Monte Carlo simulation.
+        Parameters
+        ----------
+        SNRs : 1D arraylike
+               Signal to Noise ratio in dB defined as :math:`SNR_{dB} = (E_b/N_0)_{dB} + 10 \log_{10}(R_cM_c)`
+               where :math:`Rc` is the code rate and :math:`Mc` the modulation rate.
+        send_max : int
+                   Maximum number of bits send for each SNR.
+        err_min : int
+                  link_performance send bits until it reach err_min errors (see also send_max).
+        send_chunk : int
+                      Number of bits to be send at each iteration. This is also the frame length of the decoder if available
+                      so it should be large enough regarding the code type.
+                      *Default*: send_chunck = err_min
+        code_rate : float or Fraction in (0,1]
+                    Rate of the used code.
+                    *Default*: 1 i.e. no code.
         Returns
         -------
         BERs : 1d ndarray
@@ -183,7 +296,10 @@ def link_performance(self, SNRs, send_max, err_min, send_chunk=None, code_rate=1
         # Set chunk size and round it to be a multiple of num_bits_symbol*nb_tx to avoid padding
         if send_chunk is None:
             send_chunk = err_min
-        divider = self.num_bits_symbol * self.channel.nb_tx
+        if type(code_rate) is float:
+            code_rate = Fraction(code_rate).limit_denominator(100)
+        self.rate = code_rate
+        divider = (Fraction(1, self.num_bits_symbol * self.channel.nb_tx) * 1 / code_rate).denominator
         send_chunk = max(divider, send_chunk // divider * divider)
 
         receive_size = self.channel.nb_tx * self.num_bits_symbol
@@ -191,7 +307,7 @@ def link_performance(self, SNRs, send_max, err_min, send_chunk=None, code_rate=1
 
         # Computations
         for id_SNR in range(len(SNRs)):
-            self.channel.set_SNR_dB(SNRs[id_SNR], code_rate, self.Es)
+            self.channel.set_SNR_dB(SNRs[id_SNR], float(code_rate), self.Es)
             bit_send = 0
             bit_err = 0
             while bit_send < send_max and bit_err < err_min:
@@ -203,7 +319,7 @@ def link_performance(self, SNRs, send_max, err_min, send_chunk=None, code_rate=1
                 # Deals with MIMO channel
                 if isinstance(self.channel, MIMOFlatChannel):
                     nb_symb_vector = len(channel_output)
-                    received_msg = np.empty(int(math.ceil(len(msg) / self.rate)), dtype=np.int8)
+                    received_msg = np.empty(int(math.ceil(len(msg) / float(self.rate))))
                     for i in range(nb_symb_vector):
                         received_msg[receive_size * i:receive_size * (i + 1)] = \
                             self.receive(channel_output[i], self.channel.channel_gains[i],
@@ -216,9 +332,9 @@ def link_performance(self, SNRs, send_max, err_min, send_chunk=None, code_rate=1
                     decoded_bits = self.decoder(channel_output, self.channel.channel_gains,
                                                 self.constellation, self.channel.noise_std ** 2,
                                                 received_msg, self.channel.nb_tx * self.num_bits_symbol)
-                    bit_err += np.bitwise_xor(msg, decoded_bits[:len(msg)]).sum()
+                    bit_err += np.bitwise_xor(msg, decoded_bits[:len(msg)].astype(int)).sum()
                 else:
-                    bit_err += np.bitwise_xor(msg, self.decoder(received_msg)[:len(msg)]).sum()
+                    bit_err += np.bitwise_xor(msg, self.decoder(received_msg)[:len(msg)].astype(int)).sum()
                 bit_send += send_chunk
             BERs[id_SNR] = bit_err / bit_send
             if bit_err < err_min:
diff --git a/commpy/modulation.py b/commpy/modulation.py
index 5389daf..4f90a52 100644
--- a/commpy/modulation.py
+++ b/commpy/modulation.py
@@ -21,14 +21,14 @@
 
 """
 from bisect import insort
-from itertools import product
 
 import matplotlib.pyplot as plt
-from numpy import arange, array, zeros, pi, cos, sin, sqrt, log2, argmin, \
+from numpy import arange, array, zeros, pi, sqrt, log2, argmin, \
     hstack, repeat, tile, dot, shape, concatenate, exp, \
-    log, vectorize, empty, eye, kron, inf, full, abs, newaxis, minimum, clip
+    log, vectorize, empty, eye, kron, inf, full, abs, newaxis, minimum, clip, fromiter
 from numpy.fft import fft, ifft
 from numpy.linalg import qr, norm
+from sympy.combinatorics.graycode import GrayCode
 
 from commpy.utilities import bitarray2dec, dec2bitarray, signal_power
 
@@ -65,10 +65,16 @@ class Modem:
                         If the constellation is changed to an array-like with length that is not a power of 2.
         """
 
-    def __init__(self, constellation):
+    def __init__(self, constellation, reorder_as_gray=True):
         """ Creates a custom Modem object. """
 
-        self.constellation = constellation
+        if reorder_as_gray:
+            m = log2(len(constellation))
+            gray_code_sequence = GrayCode(m).generate_gray()
+            gray_code_sequence_array = fromiter((int(g, 2) for g in gray_code_sequence), int, len(constellation))
+            self.constellation = array(constellation)[gray_code_sequence_array.argsort()]
+        else:
+            self.constellation = constellation
 
     def modulate(self, input_bits):
         """ Modulate (map) an array of bits to constellation symbols.
@@ -197,10 +203,11 @@ class PSKModem(Modem):
     def __init__(self, m):
         """ Creates a Phase Shift Keying (PSK) Modem object. """
 
-        def _constellation_symbol(i):
-            return cos(2 * pi * (i - 1) / m) + sin(2 * pi * (i - 1) / m) * (0 + 1j)
+        num_bits_symbol = log2(m)
+        if num_bits_symbol != int(num_bits_symbol):
+            raise ValueError('Constellation length must be a power of 2.')
 
-        self.constellation = list(map(_constellation_symbol, arange(m)))
+        super().__init__(exp(1j * arange(0, 2 * pi, 2 * pi / m)))
 
 
 class QAMModem(Modem):
@@ -229,6 +236,7 @@ class QAMModem(Modem):
         ------
         ValueError
                         If the constellation is changed to an array-like with length that is not a power of 2.
+                        If the parameter m would lead to an non-square QAM during initialization.
     """
 
     def __init__(self, m):
@@ -237,16 +245,21 @@ def __init__(self, m):
         Parameters
         ----------
         m : int
-            Size of the QAM constellation.
+            Size of the QAM constellation. Must lead to a square QAM (ie sqrt(m) is an integer).
 
+        Raises
+        ------
+        ValueError
+                        If m would lead to an non-square QAM.
         """
 
-        def _constellation_symbol(i):
-            return (2 * i[0] - 1) + (2 * i[1] - 1) * (1j)
+        num_symb_pam = sqrt(m)
+        if num_symb_pam != int(num_symb_pam):
+            raise ValueError('m must lead to a square QAM.')
 
-        mapping_array = arange(1, sqrt(m) + 1) - (sqrt(m) / 2)
-        self.constellation = list(map(_constellation_symbol,
-                                      list(product(mapping_array, repeat=2))))
+        pam = arange(-num_symb_pam + 1, num_symb_pam, 2)
+        constellation = tile(hstack((pam, pam[::-1])), int(num_symb_pam) // 2) * 1j + pam.repeat(num_symb_pam)
+        super().__init__(constellation)
 
 
 def ofdm_tx(x, nfft, nsc, cp_length):
diff --git a/commpy/sequences.py b/commpy/sequences.py
index ba1a97b..2862290 100644
--- a/commpy/sequences.py
+++ b/commpy/sequences.py
@@ -15,6 +15,7 @@
 """
 __all__ = ['pnsequence', 'zcsequence']
 
+import numpy as np
 from numpy import empty, exp, pi, arange, int8, fromiter, sum
 
 def pnsequence(pn_order, pn_seed, pn_mask, seq_length):
@@ -72,23 +73,40 @@ def pnsequence(pn_order, pn_seed, pn_mask, seq_length):
 
     return pnseq
 
-def zcsequence(u, seq_length):
+def zcsequence(u, seq_length, q=0):
     """
     Generate a Zadoff-Chu (ZC) sequence.
 
     Parameters
     ----------
     u : int
-        Root index of the the ZC sequence.
+        Root index of the the ZC sequence: u>0.
 
     seq_length : int
-        Length of the sequence to be generated. Usually a prime number.
+        Length of the sequence to be generated. Usually a prime number:
+        u<seq_length, greatest-common-denominator(u,seq_length)=1.
+
+    q : int
+        Cyclic shift of the sequence (default 0).
 
     Returns
     -------
     zcseq : 1D ndarray of complex floats
         ZC sequence generated.
     """
-    zcseq = exp((-1j * pi * u * arange(seq_length) * (arange(seq_length)+1)) / seq_length)
+
+    for el in [u,seq_length,q]:
+        if not float(el).is_integer():
+            raise ValueError('{} is not an integer'.format(el))
+    if u<=0:
+        raise ValueError('u is not stricly positive')
+    if u>=seq_length:
+        raise ValueError('u is not stricly smaller than seq_length')
+    if np.gcd(u,seq_length)!=1:
+        raise ValueError('the greatest common denominator of u and seq_length is not 1')
+
+    cf = seq_length%2
+    n = np.arange(seq_length)
+    zcseq = np.exp( -1j * np.pi * u * n * (n+cf+2.*q) / seq_length)
 
     return zcseq
diff --git a/commpy/tests/test_links.py b/commpy/tests/test_links.py
index 9a4ae6b..5a8c548 100644
--- a/commpy/tests/test_links.py
+++ b/commpy/tests/test_links.py
@@ -8,44 +8,94 @@
 from numpy.testing import run_module_suite, assert_allclose, dec
 from scipy.special import erfc
 
+from commpy.channelcoding.ldpc import get_ldpc_code_params, triang_ldpc_systematic_encode, ldpc_bp_decode
 from commpy.channels import MIMOFlatChannel, SISOFlatChannel
 from commpy.links import link_performance, LinkModel
-from commpy.modulation import QAMModem, kbest
+from commpy.modulation import QAMModem, kbest, best_first_detector
 
 
 @dec.slow
 def test_link_performance():
     # Set seed
-    seed(17121996)
+    seed(8071996)
+    ######################################
+    # Build models & desired solutions
+    ######################################
+    models = []
+    desired_bers = []
+    snr_range = []
+    labels = []
+    rtols = []
+    code_rates = []
 
-    # Apply link_performance to SISO QPSK and AWGN channel
+    # SISO QPSK and AWGN channel
     QPSK = QAMModem(4)
 
     def receiver(y, h, constellation, noise_var):
         return QPSK.demodulate(y, 'hard')
-    model = LinkModel(QPSK.modulate, SISOFlatChannel(fading_param=(1 + 0j, 0)), receiver,
-                      QPSK.num_bits_symbol, QPSK.constellation, QPSK.Es)
 
-    BERs = link_performance(model, range(0, 9, 2), 600e4, 600)
-    desired = erfc(sqrt(10**(arange(0, 9, 2) / 10) / 2)) / 2
-    assert_allclose(BERs, desired, rtol=0.25,
-                    err_msg='Wrong performance for SISO QPSK and AWGN channel')
+    models.append(LinkModel(QPSK.modulate, SISOFlatChannel(fading_param=(1 + 0j, 0)), receiver,
+                            QPSK.num_bits_symbol, QPSK.constellation, QPSK.Es))
+    snr_range.append(arange(0, 9, 2))
+    desired_bers.append(erfc(sqrt(10 ** (snr_range[-1] / 10) / 2)) / 2)
+    labels.append('SISO QPSK and AWGN channel')
+    rtols.append(.25)
+    code_rates.append(1)
 
-    # Apply link_performance to MIMO 16QAM and 4x4 Rayleigh channel
+    # MIMO 16QAM, 4x4 Rayleigh channel and hard-output K-Best
     QAM16 = QAMModem(16)
     RayleighChannel = MIMOFlatChannel(4, 4)
     RayleighChannel.uncorr_rayleigh_fading(complex)
 
     def receiver(y, h, constellation, noise_var):
         return QAM16.demodulate(kbest(y, h, constellation, 16), 'hard')
-    model = LinkModel(QAM16.modulate, RayleighChannel, receiver,
-                      QAM16.num_bits_symbol, QAM16.constellation, QAM16.Es)
-    SNRs = arange(0, 21, 5) + 10 * log10(QAM16.num_bits_symbol)
-
-    BERs = link_performance(model, SNRs, 600e4, 600)
-    desired = (2e-1, 1e-1, 3e-2, 2e-3, 4e-5)  # From reference
-    assert_allclose(BERs, desired, rtol=1.25,
-                    err_msg='Wrong performance for MIMO 16QAM and 4x4 Rayleigh channel')
+
+    models.append(LinkModel(QAM16.modulate, RayleighChannel, receiver,
+                            QAM16.num_bits_symbol, QAM16.constellation, QAM16.Es))
+    snr_range.append(arange(0, 21, 5) + 10 * log10(QAM16.num_bits_symbol))
+    desired_bers.append((2e-1, 1e-1, 3e-2, 2e-3, 4e-5))  # From reference
+    labels.append('MIMO 16QAM, 4x4 Rayleigh channel and hard-output K-Best')
+    rtols.append(1.25)
+    code_rates.append(1)
+
+    # MIMO 16QAM, 4x4 Rayleigh channel and soft-output best-first
+    QAM16 = QAMModem(16)
+    RayleighChannel = MIMOFlatChannel(4, 4)
+    RayleighChannel.uncorr_rayleigh_fading(complex)
+    ldpc_params = get_ldpc_code_params('commpy/channelcoding/designs/ldpc/wimax/1440.720.txt', True)
+
+    def modulate(bits):
+        return QAM16.modulate(triang_ldpc_systematic_encode(bits, ldpc_params, False).reshape(-1, order='F'))
+
+    def decoder(llrs):
+        return ldpc_bp_decode(llrs, ldpc_params, 'MSA', 15)[0][:720].reshape(-1, order='F')
+
+    def demode(symbs):
+        return QAM16.demodulate(symbs, 'hard')
+
+    def receiver(y, h, constellation, noise_var):
+        return best_first_detector(y, h, constellation, (1, 3, 5), noise_var, demode, 500)
+
+    models.append(LinkModel(modulate, RayleighChannel, receiver,
+                            QAM16.num_bits_symbol, QAM16.constellation, QAM16.Es,
+                            decoder, 0.5))
+    snr_range.append(arange(17, 20, 1))
+    desired_bers.append((1.7e-1, 1e-1, 2.5e-3))  # From reference
+    labels.append('MIMO 16QAM, 4x4 Rayleigh channel and soft-output best-first')
+    rtols.append(2)
+    code_rates.append(.5)
+
+    ######################################
+    # Make tests
+    ######################################
+
+    for test in range(len(models)):
+        BERs = link_performance(models[test], snr_range[test], 5e5, 200, 720, models[test].rate)
+        assert_allclose(BERs, desired_bers[test], rtol=rtols[test],
+                        err_msg='Wrong performance for ' + labels[test])
+        full_metrics = models[test].link_performance_full_metrics(snr_range[test], 2500, 200, 720, models[test].rate)
+        assert_allclose(full_metrics[0], desired_bers[test], rtol=rtols[test],
+                        err_msg='Wrong performance for ' + labels[test])
 
 
 if __name__ == "__main__":
diff --git a/commpy/tests/test_modulation.py b/commpy/tests/test_modulation.py
index 3dfc663..52d5780 100644
--- a/commpy/tests/test_modulation.py
+++ b/commpy/tests/test_modulation.py
@@ -3,16 +3,21 @@
 
 from itertools import product
 
-from numpy import zeros, identity, arange, concatenate, log2, array, inf
+from numpy import zeros, identity, arange, concatenate, log2, log10, array, inf, sqrt, sin, pi
 from numpy.random import seed
 from numpy.testing import run_module_suite, assert_allclose, dec, assert_raises, assert_array_equal
+from scipy.special import erf
 
-from commpy.channels import MIMOFlatChannel
+from commpy.channels import MIMOFlatChannel, SISOFlatChannel
 from commpy.links import *
 from commpy.modulation import QAMModem, mimo_ml, bit_lvl_repr, max_log_approx, PSKModem, Modem
 from commpy.utilities import signal_power
 
 
+def Qfunc(x):
+    return 0.5 - 0.5 * erf(x / sqrt(2))
+
+
 @dec.slow
 def test_bit_lvl_repr():
     # Set seed
@@ -124,8 +129,33 @@ def do_custom(self, modem):
         pass
 
 
+@dec.slow
 class TestModulateHardDemodulate(ModemTestcase):
 
+    @staticmethod
+    def check_BER(modem, EbN0dB, BERs_expected):
+        seed(8071996)
+        model = LinkModel(modem.modulate,
+                          SISOFlatChannel(fading_param=(1 + 0j, 0)),
+                          lambda y, _, __, ___: modem.demodulate(y, 'hard'),
+                          modem.num_bits_symbol, modem.constellation, modem.Es)
+        BERs = model.link_performance(EbN0dB + 10 * log10(log2(modem.m)), 5e5, 400, 720)
+        assert_allclose(BERs, BERs_expected, atol=1e-4, rtol=.1,
+                        err_msg='Wrong BER for a standard modulation with {} symbols'.format(modem.m))
+
+    def do_qam(self, modem):
+        EbN0dB = arange(8, 25, 4)
+        nb_symb_pam = sqrt(modem.m)
+        BERs_expected = 2 * (1 - 1 / nb_symb_pam) / log2(nb_symb_pam) * \
+                        Qfunc(sqrt(3 * log2(nb_symb_pam) / (nb_symb_pam ** 2 - 1) * (2 * 10 ** (EbN0dB / 10))))
+        self.check_BER(modem, EbN0dB, BERs_expected)
+
+    def do_psk(self, modem):
+        EbN0dB = arange(15, 25, 4)
+        SERs_expected = 2 * Qfunc(sqrt(2 * modem.num_bits_symbol * 10 ** (EbN0dB / 10)) * sin(pi / modem.m))
+        BERs_expected = SERs_expected / modem.num_bits_symbol
+        self.check_BER(modem, EbN0dB, BERs_expected)
+
     def do(self, modem):
         for bits in product(*((0, 1),) * modem.num_bits_symbol):
             assert_array_equal(bits, modem.demodulate(modem.modulate(bits), 'hard'),
diff --git a/commpy/tests/test_sequences.py b/commpy/tests/test_sequences.py
index f48f332..51a5bcb 100644
--- a/commpy/tests/test_sequences.py
+++ b/commpy/tests/test_sequences.py
@@ -1,16 +1,17 @@
 # Authors: CommPy contributors
 # License: BSD 3-Clause
 
+import numpy as np
 from numpy import array
-from numpy.testing import run_module_suite, assert_raises, assert_equal
-
-from commpy.sequences import pnsequence
+from numpy.testing import run_module_suite, assert_raises, assert_equal, assert_almost_equal
 
+from commpy.sequences import pnsequence, zcsequence
 
 def test_pnsequence():
     # Test the raises of errors
     with assert_raises(ValueError):
         pnsequence(4, '001', '1101', 2**4 - 1)
+    with assert_raises(ValueError):
         pnsequence(4, '0011', '110', 2 ** 4 - 1)
 
     # Test output with
@@ -19,6 +20,34 @@ def test_pnsequence():
     assert_equal(pnsequence(4, (0, 0, 1, 1), array((1, 1, 0, 1)), 7), array((1, 1, 0, 0, 1, 0, 1), int),
                  err_msg='Pseudo-noise sequence is not the one expected.')
 
+def test_zcsequence():
+    # Test the raises of errors
+    with assert_raises(ValueError):
+        zcsequence(u=-1, seq_length=20, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=20, seq_length=0, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=3, seq_length=18, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=3.1, seq_length=11, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=3, seq_length=11.1, q=0)
+    with assert_raises(ValueError):
+        zcsequence(u=3, seq_length=11, q=0.1)
+
+    # Test output with
+    assert_almost_equal(zcsequence(u=1, seq_length=2, q=0), array([1.000000e+00+0.j, 6.123234e-17-1.j]),
+        err_msg='CAZAC sequence is not the expected one.')
+
+    # Test if output cross-correlation is valid
+    seqCAZAC = zcsequence(u=3, seq_length=20, q=0)
+    x = np.fft.fft(seqCAZAC) / np.sqrt(seqCAZAC.size)
+    h = (np.fft.ifft(np.conj(x) * x)*np.sqrt(seqCAZAC.size)).T
+    corr = np.absolute(h)**2/h.size
+    assert_almost_equal(corr[0], 1.,
+        err_msg='CAZAC sequence auto-correlation is not valid, first term is not 1')
+    assert_almost_equal(corr[1:], np.zeros(corr.size-1),
+        err_msg='CAZAC sequence auto-correlation is not valid, all terms except first are not 0')
 
 if __name__ == "__main__":
     run_module_suite()
diff --git a/commpy/wifi80211.py b/commpy/wifi80211.py
new file mode 100644
index 0000000..58ee435
--- /dev/null
+++ b/commpy/wifi80211.py
@@ -0,0 +1,216 @@
+# Authors: CommPy contributors
+# License: BSD 3-Clause
+
+"""
+============================================
+Wifi 802.11 simulation (:mod:`commpy.wifi80211`)
+============================================
+
+.. autosummary::
+   :toctree: generated/
+
+   Wifi80211    -- Class to simulate the transmissions and receiving parameters of physical layer 802.11
+"""
+
+import math
+from typing import List
+
+import numpy as np
+
+import commpy.channelcoding.convcode as cc
+import commpy.links as lk
+import commpy.modulation as mod
+from commpy.channels import _FlatChannel
+
+
+# =============================================================================
+# Convolutional Code
+# =============================================================================
+
+
+class Wifi80211:
+    """
+    This class aims to simulate the transmissions and receiving parameters of physical layer 802.11 (currently till VHT (ac))
+
+    First the chunk is coded according to the generator matrix from the standard, having a rate of 1/2.
+    Then, depending on the Modulation Coding Scheme (MCS) used, puncturing is applied to achieve other coding rates.
+    For more details of which MCS map to which modulation and each coding the standard is *the* recommended place,
+    but for a lighter and faster source to check  https://mcsindex.com is a good place.
+    Finally the bits are then mapped to the modulation scheme in conformity to the MCS (BPSK, QPSK, 16-QAM, 64-QAM, 256-QAM).
+
+    On the receiving the inverse operations are perform, with depuncture when MCS needs it.
+
+    """
+    # Build memory and generator matrix
+    # Number of delay elements in the convolutional encoder
+    # "The encoder uses a 6-stage shift register."
+    # (https://pdfs.semanticscholar.org/c63b/71e43dc23b17ca57267f3b769224c64d5e33.pdf p.19)
+    memory = np.array(6, ndmin=1)
+    generator_matrix = np.array((133, 171), ndmin=2)  # from 802.11 standard, page 2295
+
+    def get_modem(self) -> mod.Modem:
+        """
+        Gets the modem that is going to be used for this particular WiFi simulation according to the MCS
+        """
+        bits_per_symbol = [
+            2,
+            4,
+            4,
+            16,
+            16,
+            64,
+            64,
+            64,
+            256,
+            256
+        ]
+        if self.mcs <= 2:
+            # BPSK for mcs 0
+            # QPSK for mcs 1 a 2
+            return mod.PSKModem(bits_per_symbol[self.mcs])
+        else:
+            # Modem : QAMModem
+            return mod.QAMModem(bits_per_symbol[self.mcs])
+
+    @staticmethod
+    def _get_puncture_matrix(numerator: int, denominator: int) -> List:
+        if numerator == 1 and denominator == 2:
+            return None
+        # from the standard 802.11 2016
+        if numerator == 2 and denominator == 3:
+            # page 2297
+            return [1, 1, 1, 0]
+        if numerator == 3 and denominator == 4:
+            # page 2297
+            return [1, 1, 1, 0, 0, 1]
+        if numerator == 5 and denominator == 6:
+            # page 2378
+            return [1, 1, 1, 0, 0, 1, 1, 0, 0, 1]
+        return None
+
+    def _get_coding(self):
+        coding = [
+            (1, 2),
+            (1, 2),
+            (3, 4),
+            (1, 2),
+            (3, 4),
+            (2, 3),
+            (3, 4),
+            (5, 6),
+            (3, 4),
+            (5, 6),
+        ]
+        return coding[self.mcs]
+
+    @staticmethod
+    def _get_trellis():
+        return cc.Trellis(Wifi80211.memory, Wifi80211.generator_matrix)
+
+    def __init__(self, mcs: int):
+        """
+        Build WiFi 802.11 simulation class
+        Parameters
+        ----------
+        mcs : int
+              The Modulation Coding Scheme (MCS) to simulate.
+              A list of MCS and which coding and modulations they correspond to is bellow:
+                 - 0 : BPSK 1/2
+                 - 1 : QPSK 1/2
+                 - 2 : QPSK 3/4
+                 - 3 : 16-QAM 1/2
+                 - 4 : 16-QAM 3/4
+                 - 5 : 64-QAM 2/3
+                 - 6 : 64-QAM 3/4
+                 - 7 : 64-QAM 5/6
+                 - 8 : 256-QAM 3/4
+                 - 9 : 256-QAM 5/6
+        """
+        self.mcs = mcs
+        self.modem = None
+
+    def link_performance(self, channel: _FlatChannel, SNRs, tx_max, err_min, send_chunk=None,
+                         frame_aggregation=1, receiver=None, stop_on_surpass_error=True):
+        """
+        Estimate the BER performance of a link model with Monte Carlo simulation as in commpy.links.link_performance
+
+        Parameters
+        ----------
+        channel : _FlatChannel
+                  The channel to be used for the simulation
+
+        SNRs : 1D arraylike
+               Signal to Noise ratio in dB defined as :math:`SNR_{dB} = (E_b/N_0)_{dB} + 10 \log_{10}(R_cM_c)`
+               where :math:`Rc` is the code rate and :math:`Mc` the modulation rate.
+
+        tx_max : int
+                 Maximum number of transmissions for each SNR.
+
+        err_min : int
+                  link_performance send bits until it reach err_min errors (see also send_max).
+
+        send_chunk : int
+                     Number of bits to be send at each frame. This is also the frame length of the decoder if available
+                     so it should be large enough regarding the code type.
+                     *Default*: send_chunck = err_min
+
+        frame_aggregation : int
+                            Number of frames aggregated per transmission (each frame with size send_chunk)
+
+        receiver  : function
+                    Specify a custom receiver function to be used in the simulation.
+                    This is particular useful for MIMO simulations.
+
+        stop_on_surpass_error : bool
+                                Controls if during simulation of a SNR it should break and move to the next SNR when
+                                the bit error is above the err_min parameter
+
+        Returns
+        -------
+        BERs : 1d ndarray
+               Estimated Bit Error Ratio corresponding to each SNRs
+        """
+
+        trellis1 = Wifi80211._get_trellis()
+        coding = self._get_coding()
+        modem = self.get_modem()
+
+        def modulate(bits):
+            res = cc.conv_encode(bits, trellis1, 'cont')
+            puncture_matrix = Wifi80211._get_puncture_matrix(coding[0], coding[1])
+            res_p = res
+            if puncture_matrix:
+                res_p = cc.puncturing(res, puncture_matrix)
+
+            return modem.modulate(res_p)
+
+        # Receiver function (no process required as there are no fading)
+        def _receiver(y, h, constellation, noise_var):
+            return modem.demodulate(y, 'soft', noise_var)
+
+        if not receiver:
+            receiver = _receiver
+
+        # Decoder function
+        def decoder_soft(msg):
+            msg_d = msg
+            puncture_matrix = Wifi80211._get_puncture_matrix(coding[0], coding[1])
+            if puncture_matrix:
+                try:
+                    msg_d = cc.depuncturing(msg, puncture_matrix, math.ceil(len(msg) * coding[0] / coding[1] * 2))
+                except IndexError as e:
+                    print(e)
+                    print("Decoded message size %d" % (math.ceil(len(msg) * coding[0] / coding[1] * 2)))
+                    print("Encoded message size %d" % len(msg))
+                    print("Coding %d/%d" % (coding[0], coding[1]))
+            return cc.viterbi_decode(msg_d, trellis1, decoding_type='soft')
+
+        self.model = lk.LinkModel(modulate, channel, receiver,
+                                  modem.num_bits_symbol, modem.constellation, modem.Es,
+                                  decoder_soft, coding[0] / coding[1])
+        return self.model.link_performance_full_metrics(SNRs, tx_max,
+                                                        err_min=err_min, send_chunk=send_chunk,
+                                                        code_rate=coding[0] / coding[1],
+                                                        number_chunks_per_send=frame_aggregation,
+                                                        stop_on_surpass_error=stop_on_surpass_error
+                                                        )
diff --git a/requirements.txt b/requirements.txt
index 4991962..731c11c 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,3 +2,4 @@ numpy>=1.9.2
 scipy>=0.15.0
 matplotlib>=1.4.3
 nose>=1.3.4
+sympy>=1.7.1
diff --git a/setup.py b/setup.py
index 080e8c0..7d3d2ce 100644
--- a/setup.py
+++ b/setup.py
@@ -1,22 +1,22 @@
 # Authors: CommPy contributors
 # License: BSD 3-Clause
 
-from distutils.core import setup
+from setuptools import setup
 
 # Taken from scikit-learn setup.py
 DISTNAME = 'scikit-commpy'
 DESCRIPTION = 'Digital Communication Algorithms with Python'
-LONG_DESCRIPTION = open('README.md').read()
+LONG_DESCRIPTION = open('README.md', encoding="utf8").read()
 MAINTAINER = 'Veeresh Taranalli & Bastien Trotobas'
-MAINTAINER_EMAIL = 'veeresht@gmail.com'
+MAINTAINER_EMAIL = 'bastien.trotobas@gmail.com'
 URL = 'http://veeresht.github.com/CommPy'
 LICENSE = 'BSD 3-Clause'
-VERSION = '0.5.0'
+VERSION = '0.8.0'
 
 #This is a list of files to install, and where
 #(relative to the 'root' dir, where setup.py is)
 #You could be more specific.
-files = ["channelcoding/*, channelcoding/tests/*, tests/*"]
+files = ["channelcoding/*", "channelcoding/tests/*", "tests/*", "channelcoding/designs/ldpc/gallager/*", "channelcoding/designs/ldpc/wimax/*"]
 
 setup(
     name=DISTNAME,
@@ -35,6 +35,7 @@
           'numpy',
           'scipy',
           'matplotlib',
+          'sympy'
     ],
     #'package' package must contain files (see list above)
     #This dict maps the package name =to=> directories
@@ -52,6 +53,7 @@
         'Intended Audience :: Science/Research',
         'Intended Audience :: Telecommunications Industry',
         'Operating System :: Unix',
+        'Operating System :: Microsoft :: Windows',
         'Programming Language :: Python',
         'Topic :: Scientific/Engineering',
         'Topic :: Software Development',