Geek asks and answers

Convolution for processing digital signals in R

I have a simple digital system that has the input x (n) = u (n) - u (n-4). System 1

I am trying to find the output y (n) with the conv () function from the signal package or convolve () from the "stats" package and plot 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 recycle of the vectors when I do the convolution, and the output of the conv () function does not seem to match the original time index. I just can't figure out how to fix my logic here. I understand that the function conv (n, m) returns a vector of length (m + n) -1, is there a good way to easily match this vector with a time index vector?

, R, , - R . .

+5
1

. conv() . :

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)

, , . - , , .:)

UPDATE

conv() . , , , .

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)

2

FFT. , sapply(), . . .

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))]
}

, , , .

+4

More articles:


All Articles