# -*- coding: utf-8 -*-
"""
Created on Sat May 18 23:52:28 2024

@author: shara
"""

import numpy as np
from scipy.stats import norm
from scipy.integrate import quad

def integrand( y, sigma ):
    p = 0.5 * ( norm.pdf( y + 1, loc=0, scale = sigma ) + norm.pdf( y - 1, loc=0, scale = sigma ) )
    if p < 1e-20:
        return 0
    return p * np.log2(p)

def biawgn_capacity( sigma ):
    integral_eval = quad(integrand, -np.inf, np.inf, args=(sigma))
    marginal_entropy = -1 * integral_eval[0]
    error = integral_eval[1]
    return marginal_entropy - 0.5 * np.log2( 2 * np.pi * np.exp(1) * sigma**2 )