Benjamini and Hochberg False Discovery Rate

//*********************************************************************
// statsBHFDR
//
// Implements Benjamini and Hochberg False Discovery Rate. Tolerates
// more false positives but less false negatives. Least stringent
// correction.
//
// Parameters:
// w - input wave of p-values
//
// Optional Parameters:
// ncol - column number of p-values if input wave w is a 2D matrix
//*********************************************************************
function statsBHFDR(w,[ncol])

    //-----------------------------------------------------------------
    // parameters
    //-----------------------------------------------------------------
    wave w;
    variable ncol;
   
    //-----------------------------------------------------------------
    // variables
    //-----------------------------------------------------------------
    variable npnts,ndim;
   
    //-----------------------------------------------------------------
    // init variables
    //-----------------------------------------------------------------
   
    if(paramIsDefault(ncol))
        ncol = 0;
    endif
   
    npnts = dimsize(w,0);
    ndim = dimsize(w,1);

    //-----------------------------------------------------------------
    // main function
    //-----------------------------------------------------------------
   
    // copy p-values out to new wave and make index wave
    make/free/n=(npnts) wpv,wpi;
   
    if(ndim)
        wpv[] = w[p][ncol];
    else
        wpv[] = w[p];
    endif
   
    // populate index wave
    wpi[] = p;
   
    // sort from highest to lowest
    sort/r wpv,wpv,wpi;
   
    // perform correction
    wpv[] = wpv[p]*(npnts/(npnts-p));
   
    // restore original sorting
    sort wpi,wpv,wpi;
   
    // copy adjusted p-values to original wave
    if(ndim)
        w[][ncol] = wpv[p];
    else
        w[] = wpv[p];
    endif
   
    // set p-values greater than 1 to 1
    w[] = w[p] > 1 ? 1 : w[p];

end