## SPDX-FileCopyrightText: Copyright (c) 2021-2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved.# SPDX-License-Identifier: Apache-2.0#"""Utility functions and blocks for the FEC package."""importtensorflowastfimportnumpyasnpimportmatplotlib.pyplotaspltimportwarningsfromimportlib_resourcesimportfiles,as_filefromsionna.phyimportBlockfromsionna.phy.fec.ldpcimportcodesfromsionna.phy.utilsimportlog2fromsionna.phyimportconfig
[docs]classGaussianPriorSource(Block):r"""Generates synthetic Log-Likelihood Ratios (LLRs) for Gaussian channels. Generates synthetic Log-Likelihood Ratios (LLRs) as if an all-zero codeword was transmitted over a Binary Additive White Gaussian Noise (Bi-AWGN) channel. The LLRs are generated based on either the noise variance ``no`` or mutual information. If mutual information is used, it represents the information associated with a binary random variable observed through an AWGN channel. .. image:: ../figures/GaussianPriorSource.png The generated LLRs follow a Gaussian distribution with parameters: .. math:: \sigma_{\text{llr}}^2 = \frac{4}{\sigma_\text{ch}^2} .. math:: \mu_{\text{llr}} = \frac{\sigma_\text{llr}^2}{2} where :math:`\sigma_\text{ch}^2` is the noise variance specified by ``no``. If the mutual information is provided as input, the J-function as described in [Brannstrom]_ is used to relate the mutual information to the corresponding LLR distribution. Parameters ---------- precision : `None` (default) | 'single' | 'double' Precision used for internal calculations and outputs. If set to `None`, :py:attr:`~sionna.phy.config.precision` is used. Input ----- output_shape : tf.int Integer tensor or Python list defining the shape of the generated LLR tensor. no : None (default) | tf.float Scalar defining the noise variance for the synthetic AWGN channel. mi : None (default) | tf.float Scalar defining the mutual information for the synthetic AWGN channel. Only used of ``no`` is None. Output ------ tf.Tensor of dtype ``dtype`` (defaults to `tf.float32`) Tensor with shape defined by ``output_shape``. """def__init__(self,*,precision=None,**kwargs):super().__init__(precision=precision,**kwargs)defcall(self,output_shape,no=None,mi=None):"""Generate Gaussian distributed fake LLRs as if the all-zero codeword was transmitted over an Bi-AWGN channel. Args: output_shape : tf.int Integer tensor or Python list defining the shape of the generated LLR tensor. no : None (default) | tf.float Scalar defining the noise variance for the synthetic AWGN channel. mi : None (default) | tf.float Scalar defining the mutual information for the synthetic AWGN channel. Only used of ``no`` is None. Returns: 1+D Tensor (``dtype``): Shape as defined by ``output_shape``. """ifnoisNone:ifmiisNone:raiseValueError("Either no or mi must be provided.")# clip Ia to range (0,1)mi=tf.cast(mi,self.rdtype)mi=tf.maximum(mi,1e-7)mi=tf.minimum(mi,1.)mu_llr=j_fun_inv(mi)sigma_llr=tf.math.sqrt(2*mu_llr)else:# noise_var must be positiveno=tf.cast(no,self.rdtype)no=tf.maximum(no,1e-7)sigma_llr=tf.math.sqrt(4/no)mu_llr=sigma_llr**2/2# generate LLRs with Gaussian approximation (BPSK, all-zero cw)# Use negative mean as we generate logits with definition p(b=1)/p(b=0)llr=config.tf_rng.normal(output_shape,mean=-1.*mu_llr,stddev=sigma_llr,dtype=self.rdtype)returnllr
[docs]defllr2mi(llr,s=None,reduce_dims=True):# pylint: disable=line-too-longr"""Approximates the mutual information based on Log-Likelihood Ratios (LLRs). This function approximates the mutual information for a given set of ``llr`` values, assuming an `all-zero codeword` transmission as derived in [Hagenauer]_: .. math:: I \approx 1 - \sum \operatorname{log_2} \left( 1 + \operatorname{e}^{-\text{llr}} \right) The approximation relies on the `symmetry condition`: .. math:: p(\text{llr}|x=0) = p(\text{llr}|x=1) \cdot \operatorname{exp}(\text{llr}) For cases where the transmitted codeword is not all-zero, this method requires knowledge of the original bit sequence ``s`` to adjust the LLR signs accordingly, simulating an all-zero codeword transmission. Note that the LLRs are defined as :math:`\frac{p(x=1)}{p(x=0)}`, which reverses the sign compared to the solution in [Hagenauer]_. Parameters ---------- llr : tf.float Tensor of arbitrary shape containing LLR values. s : None | tf.float Tensor of the same shape as ``llr`` representing the signs of the transmitted sequence (assuming BPSK), with values of +/-1. reduce_dims : `bool`, (default `True`) If `True`, reduces all dimensions and returns a scalar. If False, averages only over the last dimension. Returns ------- mi : tf.float If ``reduce_dims`` is `True`, returns a scalar tensor. Otherwise, returns a tensor with the same shape as ``llr`` except for the last dimension, which is removed. Contains the approximated mutual information. """ifllr.dtypenotin(tf.float16,tf.bfloat16,tf.float32,tf.float64):raiseTypeError("Dtype of llr must be a real-valued float.")ifsisnotNone:# ensure compatible typess=tf.cast(s,llr.dtype)# scramble sign as if all-zero cw was transmittedllr_zero=tf.multiply(s,llr)else:llr_zero=llr# clip for numerical stabilityllr_zero=tf.clip_by_value(llr_zero,-100.,100.)x=log2(1.+tf.exp(1.*llr_zero))ifreduce_dims:x=1.-tf.reduce_mean(x)else:x=1.-tf.reduce_mean(x,axis=-1)returnx
[docs]defj_fun(mu):# pylint: disable=line-too-longr"""Computes the `J-function` The `J-function` relates mutual information to the mean of Gaussian-distributed Log-Likelihood Ratios (LLRs) using the Gaussian approximation. This function implements the approximation proposed in [Brannstrom]_: .. math:: J(\mu) \approx \left( 1 - 2^{H_\text{1}(2\mu)^{H_\text{2}}}\right)^{H_\text{3}} where :math:`\mu` represents the mean of the LLR distribution, and the constants are defined as :math:`H_\text{1}=0.3073`, :math:`H_\text{2}=0.8935`, and :math:`H_\text{3}=1.1064`. Input values are clipped to [1e-10, 1000] for numerical stability. The output is clipped to a maximum LLR of 20. Parameters ---------- mu : tf.float32 Tensor of arbitrary shape, representing the mean of the LLR values. Returns ------- tf.float32 Tensor of the same shape and dtype as ``mu``, containing the calculated mutual information values. """# input must be positive for numerical stabilitymu=tf.maximum(mu,1e-10)mu=tf.minimum(mu,1000)h1=0.3073h2=0.8935h3=1.1064mi=(1-2**(-h1*(2*mu)**h2))**h3returnmi
[docs]defj_fun_inv(mi):# pylint: disable=line-too-longr"""Computes the inverse of the `J-function` The `J-function` relates mutual information to the mean of Gaussian-distributed Log-Likelihood Ratios (LLRs) using the Gaussian approximation. This function computes the inverse `J-function` based on the approximation proposed in [Brannstrom]_: .. math:: J(\mu) \approx \left( 1 - 2^{H_\text{1}(2\mu)^{H_\text{2}}}\right)^{H_\text{3}} where :math:`\mu` is the mean of the LLR distribution, and constants are defined as :math:`H_\text{1}=0.3073`, :math:`H_\text{2}=0.8935`, and :math:`H_\text{3}=1.1064`. Input values are clipped to [1e-10, 1] for numerical stability. The output is clipped to a maximum LLR of 20. Parameters ---------- mi : tf.float32 Tensor of arbitrary shape, representing mutual information values. Returns ------- tf.float32 Tensor of the same shape and dtype as ``mi``, containing the computed mean values of the LLR distribution. """# input must be positive for numerical stabilitymi=tf.maximum(mi,1e-10)# ensure that I>0mi=tf.minimum(mi,1.)# ensure that I=<1h1=0.3073h2=0.8935h3=1.1064mu=0.5*((-1/h1)*log2((1-mi**(1/h3))))**(1/(h2))returntf.minimum(mu,20)# clip the output to mu_max=20
[docs]defplot_trajectory(plot,mi_v,mi_c,ebno=None):"""Plots the trajectory of an EXIT-chart. This utility function plots the trajectory of mutual information values in an EXIT-chart, based on variable and check node mutual information values. Parameters ---------- plot : matplotlib.figure.Figure A handle to a matplotlib figure where the trajectory will be plotted. mi_v : numpy.ndarray Array of floats representing the variable node mutual information values. mi_c : numpy.ndarray Array of floats representing the check node mutual information values. ebno : float The Eb/No value in dB, used for the legend entry. """assert(len(mi_v)==len(mi_c)),"mi_v and mi_c must have same length."# number of decoding iterations to plotiters=np.shape(mi_v)[0]-1x=np.zeros([2*iters])y=np.zeros([2*iters])# iterate between VN and CN MI valuey[1]=mi_v[0]foriinrange(1,iters):x[2*i]=mi_c[i-1]y[2*i]=mi_v[i-1]x[2*i+1]=mi_c[i-1]y[2*i+1]=mi_v[i]label_str="Actual trajectory"ifebnoisnotNone:label_str+=f" @ {ebno} dB"# plot trajectoryplot.plot(x,y,"-",linewidth=3,color="g",label=label_str)# and show the legendplot.legend(fontsize=18)
[docs]defplot_exit_chart(mi_a=None,mi_ev=None,mi_ec=None,title="EXIT-Chart"):"""Plots an EXIT-chart based on mutual information curves [tenBrinkEXIT]_. This utility function generates an EXIT-chart plot. If all inputs are `None`, an empty EXIT chart is created; otherwise, mutual information curves are plotted. Parameters ---------- mi_a : numpy.ndarray, optional Array of floats representing the a priori mutual information. mi_v : numpy.ndarray, optional Array of floats representing the variable node mutual information. mi_c : numpy.ndarray, optional Array of floats representing the check node mutual information. title : str Title of the EXIT chart. Returns ------- matplotlib.figure.Figure A handle to the generated matplotlib figure. """assertisinstance(title,str),"title must be a string."ifnot(mi_evisNoneandmi_ecisNone):ifmi_aisNone:raiseValueError("mi_a cannot be None if mi_e is provided.")ifmi_evisnotNone:assert(len(mi_a)==len(mi_ev)),"mi_a and mi_ev must have same length."ifmi_ecisnotNone:assert(len(mi_a)==len(mi_ec)),"mi_a and mi_ec must have same length."plt.figure(figsize=(10,10))plt.title(title,fontsize=25)plt.xlabel("$I_{a}^v$, $I_{e}^c$",fontsize=25)plt.ylabel("$I_{e}^v$, $I_{a}^c$",fontsize=25)plt.grid(visible=True,which='major')# for MI, the x,y limits are always (0,1)plt.xlim(0,1)plt.ylim(0,1)plt.xticks(fontsize=18)plt.yticks(fontsize=18)# and plot EXIT curvesifmi_ecisnotNone:plt.plot(mi_ec,mi_a,"r",linewidth=3,label="Check node")plt.legend()ifmi_evisnotNone:plt.plot(mi_a,mi_ev,"b",linewidth=3,label="Variable node")plt.legend()returnplt
[docs]defget_exit_analytic(pcm,ebno_db):"""Calculates analytic EXIT curves for a given parity-check matrix. This function extracts the degree profile from the provided parity-check matrix ``pcm`` and calculates the EXIT (Extrinsic Information Transfer) curves for variable nodes (VN) and check nodes (CN) decoders. Note that this approach relies on asymptotic analysis, which requires a sufficiently large codeword length for accurate predictions. It assumes transmission over an AWGN channel with BPSK modulation at an SNR specified by ``ebno_db``. For more details on the equations, see [tenBrink]_ and [tenBrinkEXIT]_. Parameters ---------- pcm : numpy.ndarray The parity-check matrix. ebno_db : float Channel SNR in dB. Returns ------- mi_a : numpy.ndarray Array of floats containing the a priori mutual information. mi_ev : numpy.ndarray Array of floats containing the extrinsic mutual information of the variable node decoder for each ``mi_a`` value. mi_ec : numpy.ndarray Array of floats containing the extrinsic mutual information of the check node decoder for each ``mi_a`` value. Notes ----- This function assumes random, unstructured parity-check matrices. Thus, applying it to parity-check matrices with specific structures or constraints may result in inaccurate EXIT predictions. Additionally, this function is based on asymptotic properties and performs best with large parity-check matrices. For more information, refer to [tenBrink]_. """# calc coderaten=pcm.shape[1]k=n-pcm.shape[0]coderate=k/n# calc mean and noise_var of Gaussian distributed LLRs for given channel SNRebno=10**(ebno_db/10)snr=ebno*coderatenoise_var=1/(2*snr)# For BiAWGN channels the LLRs follow a Gaussian distr. as given below [1]sigma_llr=np.sqrt(4/noise_var)mu_llr=sigma_llr**2/2# calculate max node degree# "+1" as the array indices later directly denote the node degrees and we# have to account the array start at position 0 (i.e., we need one more# element)c_max=int(np.max(np.sum(pcm,axis=1))+1)v_max=int(np.max(np.sum(pcm,axis=0))+1)# calculate degree profile (node perspective)c=np.histogram(np.sum(pcm,axis=1),bins=c_max,range=(0,c_max),density=False)[0]v=np.histogram(np.sum(pcm,axis=0),bins=v_max,range=(0,v_max),density=False)[0]# calculate degrees from edge perspectiver=np.zeros([c_max])foriinrange(1,c_max):r[i]=(i-1)*c[i]r=r/np.sum(r)l=np.zeros([v_max])foriinrange(1,v_max):l[i]=(i-1)*v[i]l=l/np.sum(l)mi_a=np.arange(0.002,0.998,0.001)# quantize Ia with 0.01 resolution# Exit function of check node updatemi_ec=np.zeros_like(mi_a)foriinrange(1,c_max):mi_ec+=r[i]*j_fun((i-1.)*j_fun_inv(1-mi_a).numpy()).numpy()mi_ec=1-mi_ec# Exit function of variable node updatemi_ev=np.zeros_like(mi_a)foriinrange(1,v_max):mi_ev+=l[i]*j_fun(mu_llr+(i-1.)*j_fun_inv(mi_a)).numpy()returnmi_a,mi_ev,mi_ec
[docs]defload_parity_check_examples(pcm_id,verbose=False):# pylint: disable=line-too-long"""Loads parity-check matrices of built-in example codes. This utility function loads predefined example codes, including Hamming, BCH, and LDPC codes. The following codes are available: - ``pcm_id`` =0 : `(7,4)` Hamming code with `k=4` information bits and `n=7` codeword length. - ``pcm_id`` =1 : `(63,45)` BCH code with `k=45` information bits and `n=63` codeword length. - ``pcm_id`` =2 : `(127,106)` BCH code with `k=106` information bits and `n=127` codeword length. - ``pcm_id`` =3 : Random LDPC code with variable node degree 3 and check node degree 6, with `k=50` information bits and `n=100` codeword length. - ``pcm_id`` =4 : 802.11n LDPC code with `k=324` information bits and `n=648` codeword length. Parameters ---------- pcm_id : int An integer identifying the code matrix to load. verbose : `bool`, (default `False`) If `True`, prints the code parameters. Returns ------- pcm : numpy.ndarray Array containing the parity-check matrix (values are `0` and `1`). k : int Number of information bits. n : int Number of codeword bits. coderate : float Code rate, assuming full rank of the parity-check matrix. """source=files(codes).joinpath("example_codes.npy")withas_file(source)ascode:pcms=np.load(code,allow_pickle=True)pcm=np.array(pcms[pcm_id])# load parity-check matrixn=int(pcm.shape[1])# number of codeword bits (codeword length)k=int(n-pcm.shape[0])# number of information bits k per codewordcoderate=k/nifverbose:print(f"\nn: {n}, k: {k}, coderate: {coderate:.3f}")returnpcm,k,n,coderate
[docs]defbin2int(arr):"""Converts a binary array to its integer representation. This function converts an iterable binary array to its equivalent integer. For example, `[1, 0, 1]` is converted to `5`. Parameters ---------- arr : iterable of int or float An iterable that contains binary values (0's and 1's). Returns ------- int The integer representation of the binary array. """iflen(arr)==0:returnNonereturnint(''.join([str(x)forxinarr]),2)
[docs]defbin2int_tf(arr):""" Converts a binary tensor to an integer tensor. Interprets the binary representation across the last dimension of ``arr``, from most significant to least significant bit. For example, an input of `[0, 1, 1]` is converted to `3`. Parameters ---------- arr : tf.Tensor Tensor of integers or floats containing binary values (0's and 1's) along the last dimension. Returns ------- tf.Tensor Tensor with the integer representation of ``arr``. """len_=tf.shape(arr)[-1]shifts=tf.range(len_-1,-1,-1)# (2**len_-1)*arr[0] +... 2*arr[len_-2] + 1*arr[len_-1]op=tf.reduce_sum(tf.bitwise.left_shift(arr,shifts),axis=-1)returnop
[docs]defint2bin(num,length):# pylint: disable=line-too-long""" Converts an integer to a binary list of specified length. This function converts an integer ``num`` to a list of 0's and 1's, with the binary representation padded to a length of ``length``. Both ``num`` and ``length`` must be non-negative. For example: - If ``num = 5`` and ``length = 4``, the output is `[0, 1, 0, 1]`. - If ``num = 12`` and ``length = 3``, the output is `[1, 0, 0]` (truncated to length). Parameters ---------- num : int The integer to be converted into binary representation. length : int The desired length of the binary output list. Returns ------- list of int A list of 0's and 1's representing the binary form of ``num``, padded or truncated to a length of ``length``. """assertnum>=0,"Input integer should be non-negative"assertlength>=0,"length should be non-negative"bin_=format(num,f'0{length}b')binary_vals=[int(x)forxinbin_[-length:]]iflengthelse[]returnbinary_vals
[docs]defint2bin_tf(ints,length):""" Converts an integer tensor to a binary tensor with specified bit length. This function converts each integer in the input tensor ``ints`` to a binary representation, with an additional dimension of size ``length`` added at the end to represent the binary bits. The ``length`` parameter must be non-negative. Parameters ---------- ints : tf.Tensor Tensor of arbitrary shape `[..., k]` containing integers to be converted into binary representation. length : int An integer specifying the bit length of the binary representation for each integer. Returns ------- tf.Tensor A tensor of the same shape as ``ints`` with an additional dimension of size ``length`` at the end, i.e., shape `[..., k, length]`. This tensor contains the binary representation of each integer in ``ints``. """assertlength>=0shifts=tf.range(length-1,-1,delta=-1)bits=tf.math.floormod(tf.bitwise.right_shift(tf.expand_dims(ints,-1),shifts),2)returnbits
[docs]defalist2mat(alist,verbose=True):# pylint: disable=line-too-longr"""Converts an `alist` [MacKay]_ code definition to a NumPy parity-check matrix. This function converts an `alist` format representation of a code's parity-check matrix to a NumPy array. Many example codes in `alist` format can be found in [UniKL]_. About the `alist` format (see [MacKay]_ for details): - Row 1: Defines the parity-check matrix dimensions `m x n`. - Row 2: Contains two integers, `max_CN_degree` and `max_VN_degree`. - Row 3: Lists the degrees of all `n` variable nodes (columns). - Row 4: Lists the degrees of all `m` check nodes (rows). - Next `n` rows: Non-zero entries of each column, zero-padded as needed. - Following `m` rows: Non-zero entries of each row, zero-padded as needed. Parameters ---------- alist : list Nested list in `alist` format [MacKay]_ representing the parity-check matrix. verbose : `bool`, (default `True`) If `True`, prints the code parameters. Returns ------- pcm : numpy.ndarray NumPy array of shape `[n - k, n]` representing the parity-check matrix. k : int Number of information bits. n : int Number of codeword bits. coderate : float Code rate of the code. Notes ----- Use :class:`~sionna.phy.fec.utils.load_alist` to import an `alist` from a text file. Example ------- The following code snippet imports an `alist` from a file called ``filename``: .. code-block:: python al = load_alist(path=filename) pcm, k, n, coderate = alist2mat(al) """assertlen(alist)>4,"Invalid alist format."n=alist[0][0]m=alist[0][1]v_max=alist[1][0]c_max=alist[1][1]k=n-mcoderate=k/nvn_profile=alist[2]cn_profile=alist[3]# plausibility checksassertnp.sum(vn_profile)==np.sum(cn_profile),"Invalid alist format."assertnp.max(vn_profile)==v_max,"Invalid alist format."assertnp.max(cn_profile)==c_max,"Invalid alist format."iflen(alist)==len(vn_profile)+4:print("Note: .alist does not contain (redundant) CN perspective.")print("Recovering parity-check matrix from VN only.")print("Please verify the correctness of the results manually.")vn_only=Trueelse:assertlen(alist)==len(vn_profile)+len(cn_profile)+4, \
"Invalid alist format."vn_only=Falsepcm=np.zeros((m,n))num_edges=0# count number of edgesforidx_vinrange(n):foridx_iinrange(vn_profile[idx_v]):# first 4 rows of alist contain meta informationidx_c=alist[4+idx_v][idx_i]-1# "-1" as this is pythonpcm[idx_c,idx_v]=1num_edges+=1# count number of edges (=each non-zero entry)# validate results from CN perspectiveifnotvn_only:foridx_cinrange(m):foridx_iinrange(cn_profile[idx_c]):# first 4 rows of alist contain meta information# following n rows contained VN perspectiveidx_v=alist[4+n+idx_c][idx_i]-1# "-1" as this is pythonassertpcm[idx_c,idx_v]==1# entry must already existifverbose:print("Number of variable nodes (columns): ",n)print("Number of check nodes (rows): ",m)print("Number of information bits per cw: ",k)print("Number edges: ",num_edges)print("Max. VN degree: ",v_max)print("Max. CN degree: ",c_max)print("VN degree: ",vn_profile)print("CN degree: ",cn_profile)returnpcm,k,n,coderate
[docs]defload_alist(path):"""Reads an `alist` file and returns a nested list describing a code's parity-check matrix. This function reads a file in `alist` format [MacKay]_ and returns a nested list representing the parity-check matrix. Numerous example codes in `alist` format are available in [UniKL]_. Parameters ---------- path : str Path to the `alist` file to be loaded. Returns ------- list A nested list containing the imported `alist` data representing the parity-check matrix. """alist=[]withopen(path,"r")asreader:# pylint: disable=unspecified-encoding# read list line by line (different length)forlineinreader:l=[]# append all entriesforwordinline.split():l.append(int(word))ifl:# ignore empty linesalist.append(l)returnalist
[docs]defmake_systematic(mat,is_pcm=False):r"""Converts a binary matrix to its systematic form. This function transforms a binary matrix into systematic form, where the first `k` columns (or last `k` columns if ``is_pcm`` is True) form an identity matrix. Parameters ---------- mat : numpy.ndarray Binary matrix of shape `[k, n]` to be transformed to systematic form. is_pcm : `bool`, (default `False`) If `True`, ``mat`` is treated as a parity-check matrix, and the identity part will be placed in the last `k` columns. Returns ------- mat_sys : numpy.ndarray Binary matrix in systematic form, where the first `k` columns (or last `k` columns if ``is_pcm`` is True) form the identity matrix. column_swaps : list of tuple of int A list of integer tuples representing the column swaps performed to achieve the systematic form, in order of execution. Notes ----- This function may swap columns of the input matrix to achieve systematic form. As a result, the output matrix represents a permuted version of the code, defined by the ``column_swaps`` list. To revert to the original column order, apply the inverse permutation in reverse order of the swaps. If ``is_pcm`` is `True`, indicating a parity-check matrix, the identity matrix portion will be arranged in the last `k` columns. """m=mat.shape[0]n=mat.shape[1]assertm<=n,"Invalid matrix dimensions."# check for all-zero columns (=unchecked nodes)ifis_pcm:c_node_deg=np.sum(mat,axis=0)ifnp.any(c_node_deg==0):warnings.warn("All-zero column in parity-check matrix detected. " \
"It seems as if the code contains unprotected nodes.")mat=np.copy(mat)column_swaps=[]# store all column swaps# convert to bool for faster arithmeticsmat=mat.astype(bool)# bring in upper triangular formforidx_cinrange(m):success=mat[idx_c,idx_c]ifnotsuccess:# skip if leading "1" already occurred# step 1: find next leading "1"foridx_rinrange(idx_c+1,m):# skip if entry is "0"ifmat[idx_r,idx_c]:mat[[idx_c,idx_r]]=mat[[idx_r,idx_c]]# swap rowssuccess=Truebreak# Could not find "1"-entry for column idx_c# => swap with columns from non-sys part# The task is to find a column with index idx_cc that has a "1" at# row idx_cifnotsuccess:foridx_ccinrange(idx_c+1,n):ifmat[idx_c,idx_cc]:# swap columnsmat[:,[idx_c,idx_cc]]=mat[:,[idx_cc,idx_c]]column_swaps.append([idx_c,idx_cc])success=Truebreakifnotsuccess:raiseValueError("Could not succeed; mat is not full rank?")# we can now assume a leading "1" at row idx_cforidx_rinrange(idx_c+1,m):ifmat[idx_r,idx_c]:mat[idx_r,:]^=mat[idx_c,:]# bin. add of row idx_c to idx_r# remove upper triangle part in inverse orderforidx_cinrange(m-1,-1,-1):foridx_rinrange(idx_c-1,-1,-1):ifmat[idx_r,idx_c]:mat[idx_r,:]^=mat[idx_c,:]# bin. add of row idx_c to idx_r# verify resultsassertnp.array_equal(mat[:,:m],np.eye(m)), \
"Internal error, could not find systematic matrix."# bring identity part to end of matrix if parity-check matrix is providedifis_pcm:# individual column swaps instead of copying entire block# this simplifies the tracking of column swaps.foriinrange(n-1,(n-1)-m,-1):j=i-(n-m)mat[:,[i,j]]=mat[:,[j,i]]column_swaps.append([i,j])# return integer arraymat=mat.astype(int)returnmat,column_swaps
[docs]defgm2pcm(gm,verify_results=True):r"""Generates the parity-check matrix for a given generator matrix. This function converts the generator matrix ``gm`` (denoted as :math:`\mathbf{G}`) to systematic form and uses the following relationship to compute the parity-check matrix :math:`\mathbf{H}` over GF(2): .. math:: \mathbf{G} = [\mathbf{I} | \mathbf{M}] \Rightarrow \mathbf{H} = [\mathbf{M}^T | \mathbf{I}]. \tag{1} This is derived from the requirement for an all-zero syndrome, such that: .. math:: \mathbf{H} \mathbf{c}^T = \mathbf{H} * (\mathbf{u} * \mathbf{G})^T = \mathbf{H} * \mathbf{G}^T * \mathbf{u}^T = \mathbf{0}, where :math:`\mathbf{c}` represents an arbitrary codeword and :math:`\mathbf{u}` the corresponding information bits. This leads to: .. math:: \mathbf{G} * \mathbf{H}^T = \mathbf{0}. \tag{2} It can be seen that (1) satisfies (2), as in GF(2): .. math:: [\mathbf{I} | \mathbf{M}] * [\mathbf{M}^T | \mathbf{I}]^T = \mathbf{M} + \mathbf{M} = \mathbf{0}. Parameters ---------- gm : numpy.ndarray Binary generator matrix of shape `[k, n]`. verify_results : `bool`, (default `True`) If `True`, verifies that the generated parity-check matrix is orthogonal to the generator matrix in GF(2). Returns ------- numpy.ndarray Binary parity-check matrix of shape `[n - k, n]`. Notes ----- This function requires ``gm`` to have full rank. An error is raised if ``gm`` does not meet this requirement. """k=gm.shape[0]n=gm.shape[1]assertk<n,"Invalid matrix dimensions."# bring gm in systematic formgm_sys,c_swaps=make_systematic(gm,is_pcm=False)m_mat=np.transpose(np.copy(gm_sys[:,-(n-k):]))i_mat=np.eye(n-k)pcm=np.concatenate((m_mat,i_mat),axis=1)# undo column swapsforlinc_swaps[::-1]:# reverse ordering when going through listpcm[:,[l[0],l[1]]]=pcm[:,[l[1],l[0]]]# swap columnsifverify_results:assertverify_gm_pcm(gm=gm,pcm=pcm), \
"Resulting parity-check matrix does not match to generator matrix."returnpcm
[docs]defpcm2gm(pcm,verify_results=True):r"""Generates the generator matrix for a given parity-check matrix. This function converts the parity-check matrix ``pcm`` (denoted as :math:`\mathbf{H}`) to systematic form and uses the following relationship to compute the generator matrix :math:`\mathbf{G}` over GF(2): .. math:: \mathbf{G} = [\mathbf{I} | \mathbf{M}] \Rightarrow \mathbf{H} = [\mathbf{M}^T | \mathbf{I}]. \tag{1} This derivation is based on the requirement for an all-zero syndrome: .. math:: \mathbf{H} \mathbf{c}^T = \mathbf{H} * (\mathbf{u} * \mathbf{G})^T = \mathbf{H} * \mathbf{G}^T * \mathbf{u}^T = \mathbf{0}, where :math:`\mathbf{c}` represents an arbitrary codeword and :math:`\mathbf{u}` the corresponding information bits. This leads to: .. math:: \mathbf{G} * \mathbf{H}^T = \mathbf{0}. \tag{2} It can be shown that (1) satisfies (2), as in GF(2): .. math:: [\mathbf{I} | \mathbf{M}] * [\mathbf{M}^T | \mathbf{I}]^T = \mathbf{M} + \mathbf{M} = \mathbf{0}. Parameters ---------- pcm : numpy.ndarray Binary parity-check matrix of shape `[n - k, n]`. verify_results : `bool`, (default `True`) If `True`, verifies that the generated generator matrix is orthogonal to the parity-check matrix in GF(2). Returns ------- numpy.ndarray Binary generator matrix of shape `[k, n]`. Notes ----- This function requires ``pcm`` to have full rank. An error is raised if ``pcm`` does not meet this requirement. """n=pcm.shape[1]k=n-pcm.shape[0]assertk<n,"Invalid matrix dimensions."# bring pcm in systematic formpcm_sys,c_swaps=make_systematic(pcm,is_pcm=True)m_mat=np.transpose(np.copy(pcm_sys[:,:k]))i_mat=np.eye(k)gm=np.concatenate((i_mat,m_mat),axis=1)# undo column swapsforlinc_swaps[::-1]:# reverse ordering when going through listgm[:,[l[0],l[1]]]=gm[:,[l[1],l[0]]]# swap columnsifverify_results:assertverify_gm_pcm(gm=gm,pcm=pcm), \
"Resulting parity-check matrix does not match to generator matrix."returngm
[docs]defverify_gm_pcm(gm,pcm):r"""Verifies that the generator matrix :math:`\mathbf{G}` (``gm``) and parity-check matrix :math:`\mathbf{H}` (``pcm``) are orthogonal in GF(2). For a valid code with an all-zero syndrome, the following condition must hold: .. math:: \mathbf{H} \mathbf{c}^T = \mathbf{H} * (\mathbf{u} * \mathbf{G})^T = \mathbf{H} * \mathbf{G}^T * \mathbf{u}^T = \mathbf{0}, where :math:`\mathbf{c}` represents an arbitrary codeword and :math:`\mathbf{u}` the corresponding information bits. Since :math:`\mathbf{u}` can be arbitrary, this leads to the condition: .. math:: \mathbf{H} * \mathbf{G}^T = \mathbf{0}. Parameters ---------- gm : numpy.ndarray Binary generator matrix of shape `[k, n]`. pcm : numpy.ndarray Binary parity-check matrix of shape `[n - k, n]`. Returns ------- bool `True` if ``gm`` and ``pcm`` define a valid pair of orthogonal parity-check and generator matrices in GF(2). """# check for valid dimensionsk=gm.shape[0]n=gm.shape[1]n_pcm=pcm.shape[1]k_pcm=n_pcm-pcm.shape[0]assertk==k_pcm,"Inconsistent shape of gm and pcm."assertn==n_pcm,"Inconsistent shape of gm and pcm."# check that both matrices are binaryassert((gm==0)|(gm==1)).all(),"gm is not binary."assert((pcm==0)|(pcm==1)).all(),"pcm is not binary."# check for zero syndromes=np.mod(np.matmul(pcm,np.transpose(gm)),2)# mod2 to account for GF(2)returnnp.sum(s)==0# Check for Non-zero syndrome of H*G'
[docs]defgenerate_reg_ldpc(v,c,n,allow_flex_len=True,verbose=True):r"""Generates a random regular (v, c) LDPC code. This function generates a random Low-Density Parity-Check (LDPC) parity-check matrix of length ``n`` where each variable node (VN) has degree ``v`` and each check node (CN) has degree ``c``. Note that the generated LDPC code is not optimized to avoid short cycles, which may result in a non-negligible error floor. For encoding, the :class:`~sionna. fec.utils.LinearEncoder` block can be used, but the construction does not guarantee that the parity-check matrix (``pcm``) has full rank. Parameters ---------- v : int Desired degree of each variable node (VN). c : int Desired degree of each check node (CN). n : int Desired codeword length. allow_flex_len : `bool`, (default `True`) If `True`, the resulting codeword length may be slightly increased to meet the degree requirements. verbose : `bool`, (default `True`) If `True`, prints code parameters. Returns ------- pcm : numpy.ndarray Parity-check matrix of shape `[n - k, n]`. k : int Number of information bits per codeword. n : int Number of codeword bits. coderate : float Code rate of the LDPC code. Notes ----- This algorithm is designed only for regular node degrees. To achieve state-of-the-art bit-error-rate performance, optimizing irregular degree profiles is usually necessary (see [tenBrink]_). """# check input values for consistencyassertisinstance(allow_flex_len,bool), \
'allow_flex_len must be bool.'# allow slight change in n to keep num edges# from CN and VN perspective an integerifallow_flex_len:forn_modinrange(n,n+2*c):ifnp.mod((v/c)*n_mod,1.)==0:n=n_modifverbose:print("Setting n to: ",n)break# calculate number of nodescoderate=1-(v/c)n_v=nn_c=int((v/c)*n)k=n_v-n_c# generate socketsv_socks=np.tile(np.arange(n_v),v)c_socks=np.tile(np.arange(n_c),c)ifverbose:print("Number of edges (VN perspective): ",len(v_socks))print("Number of edges (CN perspective): ",len(c_socks))assertlen(v_socks)==len(c_socks),"Number of edges from VN and CN " \
"perspective does not match. Consider to (slightly) change n."# apply random permutationsconfig.np_rng.shuffle(v_socks)config.np_rng.shuffle(c_socks)# and generate matrixpcm=np.zeros([n_c,n_v])idx=0shuffle_max=200# stop if no successshuffle_counter=0cont=Truewhilecont:# if edge is available, take itifpcm[c_socks[idx],v_socks[idx]]==0:pcm[c_socks[idx],v_socks[idx]]=1idx+=1# and go to next socketshuffle_counter=0# reset counterifidx==len(v_socks):cont=Falseelse:# shuffle socketsshuffle_counter+=1ifshuffle_counter<shuffle_max:config.np_rng.shuffle(v_socks[idx:])config.np_rng.shuffle(c_socks[idx:])else:print("Stopping - no solution found!")cont=Falsev_deg=np.sum(pcm,axis=0)c_deg=np.sum(pcm,axis=1)assert((v_deg==v).all()),"VN degree not always v."assert((c_deg==c).all()),"CN degree not always c."ifverbose:print(f"Generated regular ({v},{c}) LDPC code of length n={n}")print(f"Code rate is r={coderate:.3f}.")plt.spy(pcm)returnpcm,k,n,coderate
[docs]defint_mod_2(x):r"""Modulo 2 operation and implicit rounding for floating point inputs Performs more efficient modulo-2 operation for integer inputs. Uses `tf.math.floormod` for floating inputs and applies implicit rounding for floating point inputs. Parameters ---------- x : tf.Tensor Tensor to which the modulo 2 operation is applied. Returns ------- x_mod: tf.Tensor Binary Tensor containing the result of the modulo 2 operation, with the same shape as ``x``. """ifx.dtypein(tf.int8,tf.int32,tf.int64):x_mod=tf.bitwise.bitwise_and(x,tf.constant(1,x.dtype))else:# round to next integerx_=tf.math.abs(tf.math.round(x))# tf.math.mod seems deprecatedx_mod=tf.math.floormod(x_,2)returnx_mod