y2<-20*cos(2*pi*(1/22)*xvals)+cos(2*pi*0.15*xvals)
I2<-abs(fft(y2)/sqrt(n))^2
P2<-I2*4/n

par(mfrow = c(2, 1))
par(mar = c(5, 5, 1, 1))

plot(P2[1:51], type = "l", xlab = "Frequency", ylab = "", lwd = 3,
     xaxt = "n", yaxt = "n", cex.lab = 1.6)
axis(1, at = seq(0, 50, by = 10),
     labels = seq(0, 0.5, by = 0.1), cex.axis = 1.4)
axis(2, at = seq(0, 180, by = 60),
     labels = seq(0, 180, by = 60), cex.axis = 1.4)



plot(log10(P2[1:51]), type = "l", xlab = "Frequency", ylab = "",
     lwd = 3, xaxt = "n", yaxt = "n", cex.lab = 1.6)
axis(1, at = seq(0, 50, by = 10),
     labels = seq(0, 0.5, by = 0.1), cex.axis = 1.4)
axis(2, at = seq(-1, 2, by = 1),
     labels = seq(-1, 2, by = 1), cex.axis = 1.4)

dev.print(device = postscript, "18.14.eps", horizontal = TRUE)