I'm trying to compute using matlab the mutual information for an $ \infty $-PAM input (the limit of a very dense PAM constellation) for a range of snr and I got stuck.
I'm working with a real-valued scalar Gaussian channel with input power constrained to unity and standard Gaussian noise $ N \sim (0, 1) $. The output of the channel is given by, $$ y = \sqrt{ \gamma } x + n $$
I know that under the power input constraint the input is uniformly distributed on $ \left[ -\sqrt{ 3 }, \sqrt{ 3 } \right] $.
I found that in HIGH-SNR regime, there exists a gap between the Shannon capacity and the capacity achieved with $ \infty $-PAM around 0.254 bits. I know this is called the shaping loss and represents the loss of using an uniform input rather than a Gaussian input.
I know how to compute that. But, I got stuck trying to compute the capacity for a range of snr or $ \gamma $. I tried using the formula, $$ I(x;y) = h(y) - h(n) $$ with $ f_{y}(y) = \int f_{x}(x)f_{y|x}(y|x) dx $ and I got $$ f_{y} (y) = \frac{ 1 }{ 2 \sqrt{ 3 \gamma } } \frac{ 1 }{ 2 } \int_{ \frac{ y - \sqrt{3} }{ \sqrt{2} } }^{ \frac{ y + \sqrt{3} }{ \sqrt{2} } } \frac{ 2 }{ \sqrt{ \pi } } e^{ -t^{2} } dt \\ = \frac{ 1 }{ 2 \sqrt{ 3 \gamma } } \frac{ 1 }{ 2 } \left[ \text{erf} \left( \frac{ y + \sqrt{ 3 \gamma } }{ \sqrt{ 2 } } \right) - \text{erf} \left( \frac{ y - \sqrt{ 3 \gamma } }{ \sqrt{ 2 } } \right) \right] $$
I tried to compute $ h(y) = \int_{-\infty}^{\infty} f_{y} (y) \log_{2} f_{y} (y) $ and the integral does not converge. I tried with int$(\cdot)$ and integral$(\cdot)$ methods in matlab.
I also tried approximating the $ \text{erf}(\cdot) $ function with Taylor series expansion $$ erf(x) = \frac{ 2 }{ \sqrt{ \pi } } \sum_{ n = 0 }^{ \infty } \frac{ (-1)^{n} x^{ 2n + 1 } }{ n!(2n + 1) } $$ but it didn't work.
How can I obtain a more manageable integral? Is my approach correct?
I appreciate any help!! Thanks.