R中数字信号处理的卷积

Question

我有一个简单的数字系统,它有一个输入x(n)= u(n) - u(n-4).

我试图用'signal'包中的conv()函数或'stats'包中的convolve()函数找到输出y(n),并绘制y(n)对n为-10≤n≤ 10.

到目前为止,我有以下代码:

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)  

我检查了这个情节,这是非常错误的.我认为在进行卷积时有一些向量的循环,而conv()函数的输出似乎与原始时间索引不一致.我似乎无法弄清楚如何修复我的逻辑.我意识到conv(n,m)函数返回一个长度为(m + n)-1的向量,有没有一种很好的方法可以很容易地将这个向量与时间索引向量相匹配？

这需要一些数字信号处理的知识以及R中的编码,如果有人有为此目的使用R的经验并且可以提供一些指示,那将是很好的.提前致谢.

Answer 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)

我将制定一个通用表格来完成此任务，因为我将定期执行此操作，并且不想每次都手动执行此操作。如果有人打败了我，请分享。:)

更新

创建了 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))]
}

接下来，我认为我需要努力使其变得多维。