# HW 3 Q 6 smoothspec.R JLCY, 19 Dec 2013

############################################################ 

smoothspec = function(y, deltaT) {

    N = length(y)
    freqs = 2 * pi * (1:round((N/2)))/N
    plotfreqs = freqs/(2 * pi * (deltaT))
    nfreq = length(freqs)
    a = 1:nfreq
    b = 1:nfreq

    nlag = N - 1
    # nlag = 150
    c = acf(y, type = "covariance", lag.max = nlag, plot = F)$acf

    # lag window width
    M = round(N/10)  # M should be even
    if (M%%2 > 0) 
        M = M + 1
    # try different Window widths and different window types....

    # Parzen Window..
    lambdas = 1:M
    lambdas[1:(M/2)] = 1 - (6 * ((lambdas[1:(M/2)]/M)^2)) + (6 * ((lambdas[1:(M/2)]/M)^3))
    lambdas[((M/2) + 1):M] = 2 * ((1 - (lambdas[((M/2) + 1):M]/M))^3)
    lambda0 = 1
    degf = 3.709 * N/M

    # Tukey window.. lambdas=1:M lambdas[1:M]=0.5*(1+cos(pi*(1:M)/M)) lambda0=1
    # degf=8*N/(3*M)

    # number of frequencies possible..
    freqs = 2 * pi * (1:round((N/2)))/N
    plotfreqs = freqs/(2 * pi * (deltaT))
    nfreq = length(freqs)

    tordf = 1:round((N/2))
    tord = 1:N

    # compute the smoothed periodgram/spectrum.. - Equation 12.3.15 Wei

    pgrams = 1:nfreq

    for (i in 1:nfreq) {
        pgrams[i] = ((lambda0 * c[1]) + (2 * sum(lambdas * c[2:(M + 1)] * cos(freqs[i] * 
            (1:M)))))/pi
    }

    pgrams
}