# Script with example of impulse response and step response estimation

# First time install package for reading audio files
#install.packages("tuneR")
# Load the package
library(tuneR)

# Two plots
par(mfrow=c(2,1))
iplot <- 1:300
pst <- paste0


# ----------------------------------------------------------------
# Plot'em
# Impulse reponse
x <- readWave("wavfiles/input_impulse.wav")@left[iplot]
plot(x, type="h", xlab="", main="Impulse input")
for(i in 1:4){
  y <- readWave(pst("wavfiles/impulseresponse",i,".wav"))@left[iplot]
  if(i == 1){ plot(y, type="n", xlab="", main="Impulse Response") }
  lines(y, col=i)
}

# Step response
x <- readWave("wavfiles/input_step.wav")@left[iplot]
plot(x, type="h", xlab="", main="Step input", ylim=c(0,max(x)))
for(i in 1:4){
  y <- readWave(pst("wavfiles/stepresponse",i,".wav"))@left[iplot]
  if(i == 1){ plot(y, type="n", xlab="", main="Step response") }
  lines(y, col=i)
}

# Noise response
x <- readWave("wavfiles/input_noise.wav")@left[iplot]
plot(x, type="h", xlab="", main="Noise input")
for(i in 1:4){
  y <- readWave(pst("wavfiles/noiseresponse",i,".wav"))@left[iplot]
  if(i == 1){ plot(y, type="n", xlab="", main="Noise response") }
  lines(y, col=i)
}
# ----------------------------------------------------------------


# ----------------------------------------------------------------
yir <- readWave("wavfiles/impulseresponse1.wav")@left[iplot]
ystep <- readWave("wavfiles/stepresponse1.wav")@left[iplot]

# Can we get the step response from the impulse response?
plot(ystep, type="n")
lines(cumsum(yir), col="blue")
lines(ystep, col="red")
# Back to slides, return for estimation
# ----------------------------------------------------------------


# ----------------------------------------------------------------
# Can we estimate the impulse response from the noise!?
i <- 1000:2000
ynoise <- readWave("wavfiles/noiseresponse1.wav")@left
y <- ynoise[i]
plot(y)

# Load the WAV file
wav <- readWave("wavfiles/input_noise.wav")

# Access the audio data
x <- wav@left[i]  # For mono or left channel

# Use LS
library(onlineforecast)
p <- 100
D <- as.data.frame(y=y, lagdf(x, 0:p))
(frml <- paste0("y ~ 0+",paste0("k",0:p,collapse="+")))
fit <- lm(frml, D)
summary(fit)

# See the result for the LS
par(mfrow=c(2,1))
plot(0:p, fit$coef, type="h", ylim=c(-1,1), xlab="lag", ylab="Impulse response Te to Ph", col="red", main="LS method")
lines(yir*0.00005)

# See the result and for the CCF
ccf(y,x, xlim=c(0,p), lag.max=p, col="red", ylim=c(-0.5,0.5), main="CCF method (works because of white noise input, i.e. pre-whitening not needed)")
lines(1:length(yir), yir*0.00002)
# ----------------------------------------------------------------