Convolution for Digital Signal Processing in R

Question

I have a simple digital system which has an input x(n) = u(n) - u(n-4).

I am trying to find the output y(n) with the conv() function from the 'signal' package or the convolve() function from the 'stats' package and plot the y(n) versus n for -10 ≤ n ≤ 10.

So far I have the following code:

library(signal)

n <- c(-10:10)                           # Time index
x <- c(rep(0, 10), rep(1, 4), rep(0, 7)) # Input Signal
h1 <- c(rep(0, 11), 0.5, rep(0, 9))      # Filter 1
h2 <- 0.8^n                              # Filter 2
h2[0:11] <- 0                            #

system <- data.frame(n, x, h1, h2)


y <- conv(x + conv(x, h1), h2)           # Output Signal

system <- transform(system, y=y[1:21]) 

plot(system$n, system$y)  

I checked this plot and it is very wrong. I think there is some recycling of the vectors when I do the convolution and the output of the conv() function doesn't seem to line up with the original time index. I just can't seem to figure out how to fix my logic here. I realize the conv(n, m) function returns a vector of length (m+n)-1, is there a good way to easily match this vector to a time index vector?

This would require some knowledge of Digital Signal Processing as well as coding in R, and it would be great if someone had experience in using R for this purpose and could give a few pointers. Thanks in advance.

Answer 1

I figured it out.. The center of the output of the conv() function lines up with the center of the time index vector. As such:

library(signal)

n <- c(-10:10)                           # Time index
x <- c(rep(0, 10), rep(1, 4), rep(0, 7)) # Input Signal, square pulse
h1 <- c(rep(0, 11), 0.5, rep(0, 9))      # Filter 1
h2 <- 0.8^n                              # Filter 2
h2[1:10] <- 0                            #

system <- data.frame(n, x, h1, h2)

y <- conv(x + conv(x, h1)[11:31], h2)    # Output Signal

system <- transform(system, y=y[11:31]) 

plot(system$n, system$y)

I'll work on a general form to accomplish this, as I will be doing this regularly and wouldn't want to do this manually every time. If someone beats me to it, please share. :)

UPDATE

Created a general form of the conv() function to automatically line up indices of the input and output vectors. This comes at the cost of not getting the full convolution, so you would have to set up your input as depicting the full area of interest first.

library(signal) # Should this be inside the func. with attach(), detach()?

conv2 <- function(x, y){
    conv(x, y)[ceiling(length(x)/2):(length(x)+floor(length(x)/2))]
}

# so 
y <- conv2(x + conv2(x, h1), h2)

UPDATE 2

I wanted a function to compare to FFT. I'm not exactly happy with this version, i wanted to use sapply(), but it works. For now, it'll do.. I'll work on improvements.

conv3 <- function(x, h){
m <- length(x)
n <- length(h)
X <- c(x, rep(floor(n/2), 0, floor(n/2)))   
H <- c(h, rep(floor(m/2), 0, floor(m/2)))
   Y <- vector()

for(i in 1:n+m-1){
    Y[i] <- 0 
    for(j in 1:m){
        Y[i] <- ifelse(i-j+1>0, Y[i] + X[j]*H[i-j+1], 0)
    }
}
Y[is.na(Y)] <- 0
Y[ceiling(m/2):(m+floor(m/2))]
}

Next, I think I need to work on making it multidimensional.