function y = rice(X, sigma, alpha)
    // This function computes the Rice distribution, which gives
    // the statistical distribution of peaks of band-filtered 
    // random noise.
    // Usage:
    //       y = rice(X, sigma, alpha)
    // where
    //       X:     Vector containing values of the random variable 
    //       alpha: Rice parameter, equal to (z/2/m)^2, where z is 
    //              the number of zero crossings per unit time and 
    //              m is the number of maxima per unit time. alpha 
    //              is 0 for wide band noise and 1 for very narrow
    //              band noise.    
    //       sigma: RMS value of signal
    //       y:     vector containing the probability density of the
    //              Rice distribution correspnding to the variable 
    //              values
    //
    // Reference:
    // Rice, Stephen O. (1944). Mathematical Analysis of Random
    //       Noise. Bell System Tech. Journal, 23(3), 282–332.
    //       doi:10.1002/j.1538-7305.1944.tb00874.x        
    //
    // ------------------------------
    // Author: Federico Miyara
    // Date:   2020-01-13
    
    K1 = sqrt(1-alpha)/sqrt(2*%pi)/sigma;
    K2 = sqrt(alpha)/2/sigma^2;
    y1 = exp(-X.^2/2/sigma^2/(1-alpha));
    y2 = X.*(1 + erf(X/sigma*sqrt(alpha/2/(1-alpha)))).* ...
         exp(-X.^2/2/sigma^2);
    y = K1*y1 + K2*y2;
    
    // Trap nan's caused by alpha being 1 and X being 0, which causes
    // erf to be nan and y2 to be 0*nan
    q = 1 - isnan(y);
    y(find(q==0)) = 0;
    
endfunction
