## 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)