###Define the tukey window and constants
postscript("test.ps")
###first compute the contants
auxwindow <- function(x)(1-x^3)^3
TT <- 5001
N <- (TT-1)/2
waux <-  apply(as.array(c(N:0,1:N)/N),1,auxwindow)
aux <- c(1:TT)/TT - .5/TT
W0 <- mean(waux)
W1 <-  mean(aux*waux)
W2 <-  mean(aux^2*waux)
U0 <-  mean(waux^2)
U1 <-  mean(aux*waux^2)
U2 <-  mean(aux^2*waux^2)
a0 <- (W0*W2-W1^2)^{-2}
a1 <- (U0*U2 - U1^2)
a2 <- W0^(-2)*(W0*U1 - W1*U0)^2
b0 <- (U2*W0^2 - 2*W1*U1*W0+U0*W1^2)
b1 <-  (U2*W1^2 - 2*W2*U1*W1+U0*W2^2)
 
C0 <- a0*b0
C1 <- U0/W0^2
C2 <- a0*b1
C3 <- a0*W1/(W0^2)*(W0^2*W1*U2 - W1^3*U0 - 2*W0^2*W2*U1 + 2*W0*W1*W2*U0)
C4 <- a0*(W0*W1*U2 + W1*W2*U0 - U1*(W1^2  + W0*W2))
rm(a0,a1,a2,b0,b1,auxwindow)

Window <- function(n,h){
  aux <- seq(-1,1,len=round(n*h))
  Naux1 <- round(n- length(aux))/2
  Naux2 <- n- length(aux) - Naux1
  c(rep(0,Naux1),(1-abs(aux)^3)^3,rep(0,Naux2))
}

##now we run the program that computes all the variances

cycles <- c(1:30)
par(mfrow=c(1,2))
plot(cycles,rep(0,length(cycles)),ylab="percent difference",main="Comparison of asymptotic vars\nof fund freq",ylim=c(0,100),type="n",las=1)
plot(cycles,rep(0,length(cycles)),ylab="percent difference",main="Comparison of asymptotic vars\nof total amp",ylim=c(0,100),type="n",las=1)
Ns <- c(50,100,500,1000)
A <- 2
B<-1
sigma2 <- .01
for(j in 1:4){
N <- Ns[j]
tt <- 1:N
w <- Window(N,1)
vars <- sapply(cycles,function(i){
  l<- i*pi/N
  ltt <- l*tt
  J1<-matrix(c(A^2*sum(w*tt^2*sin(ltt)^2)+B^2*sum(w*tt^2*cos(ltt)^2)
		-2*A*B*sum(w*tt^2*cos(ltt)*sin(ltt)),
	       -A*sum(w*tt*cos(ltt)*sin(ltt))+B*sum(w*tt*cos(ltt)^2),
	       -A*sum(w*tt*sin(ltt)^2)+B*sum(w*tt*cos(ltt)*sin(ltt)),
	       -A*sum(w*tt*cos(ltt)*sin(ltt))+B*sum(w*tt*cos(ltt)^2),
	       sum(w*cos(ltt)^2),
	       sum(w*cos(ltt)*sin(ltt)),
	       -A*sum(w*tt*sin(ltt)^2)+B*sum(w*tt*cos(ltt)*sin(ltt)),
	       sum(w*cos(ltt)*sin(ltt)),
	       sum(w*sin(ltt)^2)),3,3)
  J1 <- solve(J1)
  J2<-matrix(c(A^2*sum(w^2*tt^2*sin(ltt)^2)+B^2*sum(w^2*tt^2*cos(ltt)^2)
		-2*A*B*sum(w^2*tt^2*cos(ltt)*sin(ltt)),
	        -A*sum(w^2*tt*cos(ltt)*sin(ltt))+B*sum(w^2*tt*cos(ltt)^2),
	       -A*sum(w^2*tt*sin(ltt)^2)+B*sum(w^2*tt*cos(ltt)*sin(ltt)),
	       -A*sum(w^2*tt*cos(ltt)*sin(ltt))+B*sum(w^2*tt*cos(ltt)^2),
	       sum(w^2*cos(ltt)^2),
	       sum(w^2*cos(ltt)*sin(ltt)),
	       -A*sum(w^2*tt*sin(ltt)^2)+B*sum(w^2*tt*cos(ltt)*sin(ltt)),
	       sum(w^2*cos(ltt)*sin(ltt)),
	       sum(w^2*sin(ltt)^2)),3,3)
  aux <- (J1%*%J2%*%J1)
  c(aux[1,1],2/N^3*C0/(A^2+B^2),1/(A^2+B^2)*(A^2*aux[2,2]+B^2*aux[3,3]+2*A*B*aux[2,3]),2/N*(U0/W0^2))
})
par(mfg=c(1,1,1,2))
lines(cycles,100*abs(vars[1,]-vars[2,])/vars[1,],lty=j)
par(mfg=c(1,2,1,2))
lines(cycles,100*abs(vars[3,]-vars[4,])/vars[3,],lty=j)
}

par(mfg=c(1,1,1,2))
legend(5,90,sapply(as.character(Ns),function(x) paste("N=",x,sep="")),lty=1:length(Ns))
par(mfg=c(1,2,1,2))
legend(5,90,sapply(as.character(Ns),function(x) paste("N=",x,sep="")),lty=1:length(Ns))
dev.off()