phaplot1 <- function(x, ci = 0.95)
{
        alfa <- 1 - (1 - ci)/2
        plot(x$freq, (x$phase+pi)%%pi-pi, type = "l", ylim = c( - pi, pi), xlab = "", ylab
                 = "", main = "Phase", sub = paste("Phase spectra with ", ci, 
                "% pointwise confidence interval"), sep = "")
        gg <- 2/x$df
        se <- sqrt(gg/2)
        cl <- asin(pmin(0.9999, qt(ci, 2/gg - 2) * sqrt((gg * (x$coh^{-2} - 1))/(2 * (1 - gg)))))
        lines(x$freq, (x$phase+pi)%%pi-pi + cl, lty = 3)
        lines(x$freq, (x$phase+pi)%%pi-pi - cl, lty = 3)
}