## contributed by: ## Professor Douglas P. Wiens ## Department of Mathematical and Statistical Sciences ## University of Alberta ## http://www.stat.ualberta.ca/%7Ewiens/wiens.html ## ###################################### ###################################### ## These lines set the series, smoothing parameter (L), ## number of estimates (M) and desired shape of the real, symmetric ## frequency response function A(nu). After that the script should ## work for any series to be filtered in this manner. series = ts(scan("/mydata/soi.dat")) L = 5 # Smoothing M = 80 # Number of estimates # Compute the spectrum series = ts(series, frequency = 1) # This script assumes frequency = 1 (so 0 < nu < .5) spectrum = spec.pgram(series, kernel("daniell", (L-1)/2)) mat = cbind(spectrum\$freq, spectrum\$spec) mat # plot.ts(mat[1:50,1], mat[1:50,2]) # From the plot and numerical display we decide on the # desired (symmetric)shape of A(nu) for 0 < nu <.5. # In this case I want it to = 1 if .01 < nu < .05 and to = 0 otherwise. A <- function(nu) { qwe = ifelse((nu > .01 && nu < .05), 1, 0) # See help(if) and help(ifelse) qwe # Sets A(nu) to be qwe } ###################################### ## Compute the filter coefficients N = 2*length(spectrum\$freq) # This will be T'. # First look only at frequencies 1/M, 2/M, ... .5*M/M # Currently the frequencies used are are 1/N, 2/N, ... .5*N/N sampled.indices = (N/M)*(1:(M/2)) # These are the indices of the frequencies we want fr.N = spectrum\$freq fr.M = fr.N[sampled.indices] # These will be the frequencies we want spec.N = spectrum\$spec spec.M = spec.N[sampled.indices] # Power at these frequencies # Evaluate the desired frequency reponse A(nu) at these frequencies: A.desired = vector(length = length(fr.M)) for(k in 1:length(fr.M)) A.desired[k] = A(fr.M[k]) # Invert A.desired, by discretizing the defining integral, to get the coefficients a: delta = 1/M Omega = seq(from = 1/M, to = .5, length = M/2) aa = function(s) 2*delta*sum(exp(2i*pi*Omega*s)*A.desired) S = ((-M/2+1):(M/2-1)) a = vector(length = length(S)) for(k in 1:length(S)) a[k] = aa(S[k]) a = Re(a) # The filter coefficients # Apply a cosine taper h = .5*(1+cos(2*pi*S/length(S))) a = a*h # Comment out this line, to see the effect of NOT tapering cat("The filter coefficients are", "\n") qwe = cbind(S, a) colnames(qwe) = c("s", "a(s)") print(qwe[((length(S)+1)/2):1,]) cat("for s >=0; and a(-s) = a(s).", "\n") # Compute the realized frequency response function, and the filtered series: A.M = function(nu) Re(sum(exp(-2i*pi*nu*S)*a)) A.attained = vector(length = length(fr.N)) A.theoretical = vector(length = length(fr.N)) for(k in 1:length(fr.N)) { A.attained[k] = A.M(fr.N[k]) # The attained freq. resp. A.theoretical[k] = A(fr.N[k]) } series.filt = filter(series, a, sides = 2) # The filtered series par(mfrow=c(2,1)) plot.ts(series, main = "Original series") plot.ts(series.filt, main = "Filtered series") x11() par(mfrow=c(2,1)) spectrum = spec.pgram(series, kernel("daniell", (L-1)/2), main = "Spectrum of original series") spectrum.filtered = spec.pgram(na.omit(series.filt), kernel("daniell", (L-1)/2), main = "Spectrum of filtered series") x11() par(mfrow=c(2,1)) plot(S, a, xlab = "s", ylab = "a(s)", main = "Filter coefficients") plot(fr.N, A.theoretical, type = "l", lty = 2, xlab = "freq", ylab = "freq. response", main = "Desired and attained frequency response functions") lines(fr.N, A.attained, lty = 1, col = 2)