require (signal)

dt <- 0.01
t <- seq (-10, 50, by=dt)

s <- ifelse (t>=0, 1, 0)

CR.RC4 <- function (s, tau) {
	my.HP <- Arma (b = c(tau / dt, -tau / dt),
			a = c(1+tau/dt, -tau/dt))
	dummy <- filter (my.HP, s)
	my.LP <- Arma (b = 1,
			a = c(1+tau/dt, -tau/dt))
	dummy <- filter (my.LP, dummy)
	dummy <- filter (my.LP, dummy)
	dummy <- filter (my.LP, dummy)
	filter (my.LP, dummy)
}

out <- CR.RC4(s, 1)
plot (t, out, type="l", lwd=3)
grid()
lines (t, 1/factorial(4)*(t)^4*exp(-t), col="red")

plot (c(-10,50), c(0,0.2), type="n", xlab="t", ylab="simple CR-RC4 response")
grid()

tau <- 1
lines (t, CR.RC4(s,tau), type="l", lwd=5)
lines (t, ifelse(t>=0, 1/factorial(4)*(t/tau)^4*exp(-t/tau), 0), col="white")

tau <- 2
lines (t, CR.RC4(s,tau), type="l", lwd=5)
lines (t, ifelse(t>=0, 1/factorial(4)*(t/tau)^4*exp(-t/tau), 0), col="white")

tau <- 4
lines (t, CR.RC4(s,tau), type="l", lwd=5)
lines (t, ifelse(t>=0, 1/factorial(4)*(t/tau)^4*exp(-t/tau), 0), col="white")

##
## true CSP response
##
par (mar=c(4,4,1,1))

plot (t, s, type="l", lwd=3)
grid ()

tau <- 1
plot (t, CR.RC4(s,tau), type="l", lwd=5)
grid ()

s <- ifelse (t>=0, exp(-t/100), 0)

plot (t, s, type="l", lwd=3)
grid ()

tau <- 1
plot (t, CR.RC4(s,tau), type="l", lwd=5)
grid ()

s <- ifelse (t>=0, exp(-t/100), 0) + ifelse (t>=15, exp(-(t-15)/100), 0)
 
plot (t, s, type="l", lwd=3)
grid ()

tau <- 1
plot (t, CR.RC4(s,tau), type="l", lwd=5)
grid ()

##
## with PZajd
##

CR.RC4 <- function (s, tau, pz.adj) {
	my.PZ <- Arma (b = c(1+pz.adj / dt, -pz.adj / dt),
			a = 1)
	dummy <- filter (my.PZ, s) / pz.adj
	my.LP <- Arma (b = 1,
			a = c(1+tau/dt, -tau/dt))
	dummy <- filter (my.LP, dummy)
	dummy <- filter (my.LP, dummy)
	dummy <- filter (my.LP, dummy)
	dummy <- filter (my.LP, dummy)
	filter (my.LP, dummy)
}

s <- ifelse (t>=0, exp(-t/100), 0) + ifelse (t>=15, exp(-(t-15)/100), 0)
 
plot (t, s, type="l", lwd=3)
grid ()

tau <- 1
plot (t, CR.RC4(s,tau,100), type="l", lwd=5)
grid ()