# Plot the spectrum of a harmonic series
# a: cosine coefficients 
# b: sine coefficients
# f: frequencies

harsp <- function(a=c(2,1,6,7),
              b=c(3,3,1,4),f=c(1/12,1/36,1/6,1/3))
{
 y <- (a^2+b^2)/2

 plot(f,y,type="h",xlim=c(0,.5),xlab="Frequency",
                  ylab="Power spectrum")
}