Does anyone know how to figure out a workaround to avoid computing overflow/u​nderflow/N​aN/inf in this algorithm?

        M14 = Signal.^14;
        M12 = Signal.^12    ; M10 = Signal.^10;
        M8  = Signal.^8     ; M6  = Signal.^6;
        M4  = Signal.^4     ; M2  = Signal.^2;
        S14 = Sigma.^14;
        S12 = Sigma.^12     ; S10 = Sigma.^10;
        S8  = Sigma.^8      ; S6  = Sigma.^6;
        S4  = Sigma.^4      ; S2  = Sigma.^2;

        nPiD2               = pi/2;
        sqrtNpiD2           = sqrt(nPiD2);
        n1D2                = 1/2;
        n1D4                = 1/4;

        n1DM10Sig           = 1./(M10.*Sigma);
        n1DM12Sig           = 1./(M12.*Sigma);

        alpha               = M2./S2;
        nAlphaD4            = n1D4*alpha;

        FirstTerm  = n1DM10Sig.*(M12 + 9*M10.*S2 - 15*M8.*S4 + 90*M6.*S6 - 495*M4.*S8 + 2160*M2.*S10 - 5760*S12).*besseli(0,nAlphaD4);
        SecondTerm = n1DM12Sig.*(M14 + 7*M12.*S2 - 27*M10.*S4 + 150*M8.*S6 - 855*M6.*S8 + 4320*M4.*S10 - 17280*M2.*S12 + 46080*S14).*besseli(1,nAlphaD4);

        biasedSignal = n1D2*sqrtNpiD2*exp(-nAlphaD4).*(FirstTerm + SecondTerm);

As you can imagine, because of the powers of these numbers being rather high, I am running into issues with computing inf/NaN where I don't actually want it. Is there a way to avoid computing these values?

Answer 1

[Jan]

You can calculate the logarithm of all equations to keep the ranges of the values inside the limits. Replace besseli by its taylor series to build its log.

[Eric Diaz]

Great thinking! I'll try that now and see if it works! Thanks!

[Eric Diaz]

It turns out that most of the overflow problem was occuring in the besseli function.

I actually was having some problems with overflow of the besseli function two nights ago and I found a solution which works really well without having to use a taylor approximation.

Instead of explaining it, I will give you the link of where I found the solution. As you may know, Cleve Moler, who is the person that is providing the solution, is the person that founded MATLAB.

https://www.mathworks.com/matlabcentral/newsreader/view_thread/101943

I found that solution #2, works really well!