source("moving_average_functions.R")

# Frequency Response for filter with a=1 and coefficients b given (e.g. Moving Average)
frequencyResponseExp <- function(b, ws = seq(-pi, pi, by=0.01))
{
  H <- complex(length(ws))
  k <- 1
  for (w in ws)
  {
    n <- 0
    for (c in b)
    {
      H[k] <- H[k] + c * exp(0 - 1i * w * n)
      n <- n+1
    }
    k <- k + 1
  }
  return(abs(H))
}

# Frequency Response for filter with a=1 and coefficients b given (e.g. Moving Average)
frequencyResponseCosSin <- function(b, ws = seq(-pi, pi, by=0.01))
{
  C <- numeric(length(ws))
  i <- 1
  for (w in ws)
  {
    cossum <- 0
    sinsum <- 0
    k <- 0
    for (bk in b)
    {
      cossum <- cossum + bk * cos(k*w)
      sinsum <- sinsum + bk * sin(k*w)
      k <- k + 1
    }
    C[i] <- sqrt(cossum^2+sinsum^2)
    i <- i + 1
  }
  return(C)
}

# Frequency Response (absolute value) by cosine response
cosInOut <- function(b, a = 1, ws = seq(-pi, pi, by=0.01)) 
{
  tic <- Sys.time()
  H <- numeric(length(ws))
  k <- 1
  for (w in ws)
  {
    inVal <- cos(w*seq(0, 200*pi))
    outVal <- filter(b, a, inVal)
    #outVal <- filterMA(inVal, length(b))
    #outVal <- filterCMA(inVal, length(b))
    #outVal <- filterWMA(inVal, b)
    #outVal <- filterCWMA(inVal, b)
    # H[k] <- max(outVal[(length(b)+1):(length(outVal)-length(b)+1)]) / max(inVal)
    H[k] <- max(outVal) / max(inVal)
    k <- k + 1
  }
  print(Sys.time() - tic)
  return(H)
}

# Frequency Response with fft for FIR
frequencyResponseFFT_FIR <- function(b, ws = seq(-pi, pi, by=0.01))
{
  Nfft = 2^ceiling(log2(length(ws)))
  H <- waved::fftshift(fft(c(b, numeric(Nfft-length(b)))))
  wfft <- seq(-Nfft/2, Nfft/2-1)*(2*pi)/(Nfft-1)
  Habs <- approx(wfft, abs(H), ws, rule=2)
  return(Habs[[2]])
}

# Frequency Response with fft for IIR
frequencyResponseFFT_IIR <- function(b, a, ws = seq(-pi, pi, by=0.01))
{
  Nfft = 2^ceiling(log2(length(ws)))
  H <- waved::fftshift(fft(c(b, numeric(Nfft-length(b))))/fft(c(a, numeric(Nfft-length(a)))))
  wfft <- seq(-Nfft/2, Nfft/2-1)*(2*pi)/(Nfft-1)
  Habs <- approx(wfft, abs(H), ws, rule=2)
  return(Habs[[2]])
}